传送门
这个题真的也是有点难度啊（应该是因为我太菜了）
[ans=sum_{i=1}^{n}sum_{j=1}^{m}gcd(i,j)^k\
]
可以设
[f(d)=sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)==d]d^k\
ans=sum_{d=1}^{min(n,m)}f(d)=sum_{d=1}^{min(n,m)}sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)==d]d^k\
ans=sum_{d=1}^{min(n,m)}d^ksum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)==d]\
]
设
[F(d)=sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)==d]\
g(x)=sum_{x|d}F(d)=sum_{i=1}^{n}sum_{j=1}^{m}[x|gcd(i,j)]=lfloor frac{n}{x}
floorlfloor frac{m}{x}
floor\
F(d)=sum_{d|x}g(x)mu(frac{x}{d})=sum_{d|x}lfloor frac{n}{x}
floorlfloor frac{m}{x}
floormu(frac{x}{d})\
]
设(T=frac{x}{d})
[F(d)=sum_{T=1}^{min(n,m)/d}lfloor frac{n}{Td}
floorlfloor frac{m}{Td}
floormu(T)\
]
那么答案就是
[ans=sum_{d=1}^{min(n,m)}d^{k}F(d)\
ans=sum_{d=1}^{min(n,m)}d^{k}sum_{T=1}^{min(n,m)/d}lfloor frac{n}{Td}
floorlfloor frac{m}{Td}
floormu(T)\
]
然后这个式子已经可以预处理(d^k)，然后数论分块来做了
不过会(TLE)，此时时间复杂度为(O(nlogn+Tn))
所以我们还需要优化
考虑设(k=Td)
[ans=sum_{k=1}^{min(n,m)}lfloor frac{n}{k}
floorlfloor frac{m}{k}
floorsum_{d|k}mu(frac{k}{d})d^k
]
考虑设(s(d)=sum_{d|k}mu(frac{k}{d})d^k)
可以发现(s)是积性函数，因为(mu)和(d^k)都是积性函数，积性函数和积性函数的狄利克雷卷积也是积性函数
所以(s)可以通过线筛求出
考虑如何线筛，对于(a)与(b)不互质的情况下，直接用(s(ab)=s(a)*s(b))
当(x)为质数时，除了(mu(1))以及(mu(x))不为(0)，其余都为(0)，那么(s(x)=mu(1)*x^k+mu(x)*1^k)
考虑(a)有(x)个质因子(prime_i)，那么对于(prime_i)的贡献(s(prime_{i}^{x})=mu(1)*prime_i^x+mu(prime_i)*prime_i^{x-1})
那么对于(s(a*prime_i))，将相当于整体乘上一个(prime_i^k)
然后就可以筛了，我讲的有点不清楚，可以看代码理解一下
#include<cstdio>
#include<iostream>
#include<algorithm>
using namespace std;
void read(int &x) {
char ch; bool ok;
for(ok=0,ch=getchar(); !isdigit(ch); ch=getchar()) if(ch=='-') ok=1;
for(x=0; isdigit(ch); x=x*10+ch-'0',ch=getchar()); if(ok) x=-x;
}
#define rg register
const int maxn=5e6+10,mod=1e9+7;bool vis[maxn];
int T,k,n,m,pri[maxn],tot,ans,f[maxn],g[maxn];
int mi(int a,int b)
{
int ans=1;
while(b)
{
if(b&1)ans=1ll*ans*a%mod;
b>>=1,a=1ll*a*a%mod;
}
return ans;
}
void prepare()
{
f[1]=1;
for(rg int i=2;i<=5e6;i++)
{
if(!vis[i])pri[++tot]=i,g[i]=mi(i,k),f[i]=(g[i]-1+mod)%mod;
for(rg int j=1;pri[j]*i<=5e6&&j<=tot;j++)
{
vis[pri[j]*i]=1;
if(!(i%pri[j])){f[i*pri[j]]=1ll*f[i]*g[pri[j]]%mod;break;}
else f[i*pri[j]]=1ll*f[i]*f[pri[j]]%mod;
}
}
for(rg int i=1;i<=5e6;i++)f[i]=(f[i]+f[i-1])%mod;
}
int main()
{
read(T),read(k),prepare();
while(T--)
{
read(n),read(m);if(n>m)swap(n,m);ans=0;
for(rg int i=1,j;i<=n;i=j+1)
{
j=min(n/(n/i),m/(m/i));
ans=(ans+1ll*(f[j]-f[i-1]+mod)%mod*(n/i)%mod*(m/i)%mod)%mod;
}
printf("%d
",ans);
}
}