bzoj 4407: 于神之怒加强版【莫比乌斯反演+线性筛】

看着就像反演，所以先推式子（默认n<m）：

\[\sum_{d=1}^{n}d^k\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d] \]

\[=\sum_{d=1}^{n}d^k\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor }\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor }[gcd(i,j)==d] \]

\[=\sum_{d=1}^{n}d^k\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor }\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor }\sum_{g|i,g|j}\mu(g) \]

\[=\sum_{d=1}^{n}d^k\sum_{g=1}^{\lfloor \frac{n}{d} \rfloor }\mu(g)\lfloor\frac{n}{dg}\rfloor\lfloor\frac{m}{dg}\rfloor \]

一般这样就行了，但是这里T很大，所以看看有没有能预处理的东西，枚举p=dg

\[=\sum_{d=1}^{n}d^k\sum_{d|p}\mu(\frac{p}{d})\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor \]

\[=\sum_{d=1}^{n}\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor\sum_{d|p}\mu(\frac{p}{d})d^k \]

前面那段和nm有关，分块来做；考虑怎么预处理后面的
显然是个积性的，所以考虑线性筛出来前缀和即可

#include<iostream>
#include<cstdio>
using namespace std;
const int N=5000005,mod=1e9+7;
int T,k,n,m,p[N],tot,s[N],f[N],sm[N],ans;
bool v[N];
int read()
{
	int r=0,f=1;
	char p=getchar();
	while(p>'9'||p<'0')
	{
		if(p=='-')
			f=-1;
		p=getchar();
	}
	while(p>='0'&&p<='9')
	{
		r=r*10+p-48;
		p=getchar();
	}
	return r*f;
}
int ksm(int a,int b)
{
	int r=1;
	while(b)
	{
		if(b&1)
			r=1ll*r*a%mod;
		a=1ll*a*a%mod;
		b>>=1;
	}
	return r;
}
int main()
{
	T=read(),k=read();
	f[1]=1;
	for(int i=2;i<=5000000;i++)
	{
		if(!v[i])
		{
			p[++tot]=i;
			s[i]=ksm(i,k);
			f[i]=s[i]-1;
		}
		for(int j=1;j<=tot&&i*p[j]<=5000000;j++)
		{
			v[i*p[j]]=1;
			if(i%p[j]==0)
			{
				f[i*p[j]]=1ll*f[i]*s[p[j]]%mod;
				break;
			}
			f[i*p[j]]=1ll*f[i]*f[p[j]]%mod;
		}
	}
	for(int i=1;i<=5000000;i++)
		sm[i]=(sm[i-1]+f[i])%mod;
	while(T--)
	{
		n=read(),m=read(),ans=0;
		if(n>m)
			swap(n,m);
		for(int i=1,la;i<=n;i=la+1)
		{
			la=min(n/(n/i),m/(m/i));
			ans=(ans+1ll*(n/i)*(m/i)%mod*(sm[la]-sm[i-1])%mod)%mod;
		}
		printf("%d\n",(ans+mod)%mod);
	}
	return 0;
}