题面
题目分析
[egin{split}
sumlimits_{i=1}^nsumlimits_{j=1}^mgcd(i,j)^k&=sumlimits_{d=1}^nd^ksumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)==d]\
end{split}
]
设(f(x))表示(gcd(i,j)=x)，(g(x))表示(gcd(i,j)==kx,kin Z)。
[egin{split}
g(x)&=sumlimits_{x|d}^nf(d)\
&=sumlimits_{i=1}^nsumlimits_{j=1}^m[x|gcd(i,j)]\
&=sumlimits_{i=1}^{lfloorfrac n x
floor}sumlimits_{j=1}^{lfloorfrac m x
floor}lfloorfrac n x
floorlfloorfrac m x
floor\
f(x)&=sumlimits_{x|d}^nmu(frac dx)g(d)=sumlimits_{x|d}^nmu(frac dx)lfloorfrac n d
floorlfloorfrac m d
floor
end{split}
]
[egin{split}
ans&=sumlimits_{d=1}^nd^kcdot f(d)\
&=sumlimits_{d=1}^nd^ksumlimits_{d|T}^nmu(frac Td)lfloorfrac n T
floorlfloorfrac m T
floor\
&=sumlimits_{T=1}^nlfloorfrac n T
floorlfloorfrac m T
floorsumlimits_{d|T}mu(frac Td)d^k
end{split}
]
由于(mu)和(d^k)均为积性函数，所以(sumlimits_{d|T}mu(frac Td)d^k)也为积性函数，可以在线性筛中(O(nlog n))预处理。
前面部分用整除分块加速。
代码实现
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<iomanip>
#include<cstdlib>
#define MAXN 0x7fffffff
typedef long long LL;
const int N=5000005,mod=1e9+7;
using namespace std;
inline int Getint(){register int x=0,f=1;register char ch=getchar();while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}return x*f;}
int g[N],mu[N],prime[N];
bool vis[N];
LL ksm(LL x,LL k){
LL ret=1;
while(k){
if(k&1)ret=ret*x%mod;
x=x*x%mod;
k>>=1;
}
return ret;
}
int low[N];
int main(){
int T=Getint(),K=Getint();
mu[1]=g[1]=1;
for(int i=2;i<=5e6;i++){
if(!vis[i]){
prime[++prime[0]]=i,mu[i]=-1;
low[i]=i,g[i]=ksm(i,K)-1;
}
for(int j=1;j<=prime[0]&&1ll*prime[j]*i<=5e6;j++){
vis[i*prime[j]]=1;
if(i%prime[j]==0){
low[i*prime[j]]=low[i]*prime[j];
if(low[i*prime[j]]==i*prime[j])
g[i*prime[j]]=g[i]*ksm(prime[j],K)%mod;
else
g[i*prime[j]]=(1ll*g[low[i*prime[j]]]*g[i*prime[j]/low[i*prime[j]]])%mod;
break;
}
low[i*prime[j]]=prime[j];
g[i*prime[j]]=(1ll*g[i]*g[prime[j]])%mod;
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=5e6;i++)g[i]=(g[i]+g[i-1])%mod;
while(T--){
int n=Getint(),m=Getint();
if(n>m)swap(n,m);
int ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+1ll*(n/l)*(m/l)%mod*(g[r]-g[l-1])%mod+mod)%mod;
}
cout<<ans<<'
';
}
return 0;
}