【于神之怒加强版】

求

[sum_{i=1}^Nsum_{j=1}^M(i,j)^k ]

多组询问，但是每次的(k)都是不变的

先是套路

[f(n)=sum_{i=1}^Nsum_{j=1}^M[(i,j)=n] ]

[F(n)=sum_{n|d}f(d)=left lfloor frac{N}{n} ight floor imesleft lfloor frac{M}{n} ight floor ]

我们要求的就是

[Ans=sum_{i=1}^Ni^kf(i) ]

[=sum_{i=1}^Ni^ksum_{i|d}mu(frac{d}{i}) imes left lfloor frac{N}{d} ight floor imesleft lfloor frac{M}{d} ight floor ]

交换一下枚举顺序

[Ans=sum_{d=1}left lfloor frac{N}{d} ight floor imesleft lfloor frac{M}{d} ight floorsum_{i|d}mu(frac{d}{i})i^k ]

我们设(h(d)=sum_{i|d}mu(frac{d}{i})i^k)

这显然是一个积性函数

考虑一下如何线筛这个积性函数

如果(d)为质数，显然(h(d)=d^k-1)

发现这个形式很像欧拉函数，可以尝试搞一下(h(d^t))和(h(d^{t-1}))之间的关系

发现从(d^t)到(d^{t-1})只不过是多了一个约数，但是显然这个约数本身并没有什么贡献

可以直接不用考虑这个约数

于是

[h(d^{t-1})=sum_{i|d^{t-1}}mu(frac{d^{t-1}}{i})i^k ]

因为新多出来的那一个约数并不需要考虑，因为其(mu)值等于0

所以(h(d^t))和(h(d^{t-1}))的有用的约数并没有什么变化

所以

[h(d^t)=sum_{i|d^{t-1}}mu(frac{d^{t-1}}{i})(i imes d)^k ]

[=d^ksum_{i|d^{k-1}}mu(frac{d^{k-1}}{i})i^k=d^k imes h(d^{t-1}) ]

发现这个跟欧拉函数好像啊，我们完全可以按照线性筛欧拉函数的方式来线筛(h)

所以线筛之后求一个前缀和之后分块就好了

#include<iostream>
#include<cstring>
#include<cstdio>
#define re register
#define maxn 5000005
#define LL long long
#define min(a,b) ((a)<(b)?(a):(b))
#define max(a,b) ((a)>(b)?(a):(b))
const LL mod=1000000007;
inline int read()
{
	char c=getchar();
	int x=0;
	while(c<'0'||c>'9') c=getchar();
	while(c>='0'&&c<='9')
		x=(x<<3)+(x<<1)+c-48,c=getchar();
	return x;
}
int p[maxn>>1],f[maxn],T,n,m;
LL pre[maxn];
int N[2005],M[2005];
inline LL quick(LL a,LL b)
{
	LL S=1;
	while(b){if(b&1) S=S*a%mod;b>>=1;a=a*a%mod;}
	return S;
}
LL K;
int main()
{
	T=read(),K=read();
	for(re int i=1;i<=T;i++) N[i]=read(),M[i]=read();
	for(re int i=1;i<=T;i++) if(N[i]>M[i]) std::swap(N[i],M[i]);
	for(re int i=1;i<=T;i++) n=max(n,M[i]);
	f[1]=1,pre[1]=1;
	for(re int i=2;i<=n;i++)
	{
		if(!f[i]) p[++p[0]]=i,pre[i]=quick(i,K)-1;
		for(re int j=1;j<=p[0]&&p[j]*i<=n;j++)
		{
			f[p[j]*i]=1;
			if(i%p[j]==0)
			{
				pre[p[j]*i]=pre[i]*(pre[p[j]]+1)%mod;
				break;
			}
			pre[p[j]*i]=pre[p[j]]*pre[i]%mod;
		}
	}
	for(re int i=1;i<=n;i++) pre[i]=(pre[i]+pre[i-1])%mod;
	for(re int t=1;t<=T;t++)
	{
		n=N[t],m=M[t];
		LL ans=0;
		for(re LL l=1,r;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			ans=(ans+(n/l)*(m/l)%mod*(pre[r]-pre[l-1]+mod)%mod)%mod;
		}
		printf("%lld
",ans);
	}
	return 0;
}