# BZOJ4407 : 于神之怒加强版

$\begin{eqnarray*} &&\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k\\ &=&\sum_d d^k\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]\\ &=&\sum_d d^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1]\\ &=&\sum_d d^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{u|i,u|j}\mu(u)\\ &=&\sum_d d^k\sum_{u}\mu(u)\lfloor\frac{n}{du}\rfloor\lfloor\frac{m}{du}\rfloor\\ &=&\sum_D \lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor\sum_{d|D}d^k\mu(\frac{D}{d})\\ &=&\sum_D \lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor F(D) \end{eqnarray*}$

#include<cstdio>
const int N=5000001,P=1000000007;
int T,n,m,i,j,k,tot,p[N],f[N],g[N],F[N],ans;bool v[N];
inline int pow(int a,int b){int t=1;for(;b;b>>=1,a=1LL*a*a%P)if(b&1)t=1LL*t*a%P;return t;}
inline int min(int a,int b){return a<b?a:b;}
int main(){
scanf("%d%d",&T,&k);
for(F[1]=1,i=2;i<N;i++){
if(!v[i])f[i]=pow(i,k),g[i]=i,F[i]=f[i]-1,p[tot++]=i;
for(j=0;j<tot&&i*p[j]<N;j++){
v[i*p[j]]=1;
if(i%p[j]){
g[i*p[j]]=p[j];
F[i*p[j]]=1LL*F[i]*F[p[j]]%P;
}else{
g[i*p[j]]=g[i]*p[j];
F[i*p[j]]=g[i]!=i?1LL*F[i/g[i]]*F[g[i]*p[j]]%P:1LL*F[i]*f[p[j]]%P;
break;
}
}
}
for(i=2;i<N;i++)F[i]=(F[i-1]+F[i])%P;
while(T--){
scanf("%d%d",&n,&m);
for(ans=0,i=1;i<=n&&i<=m;i=j+1){
j=min(n/(n/i),m/(m/i));
ans=(1LL*(F[j]-F[i-1]+P)*(n/i)%P*(m/i)+ans)%P;
}
printf("%d\n",ans);
}
return 0;
}


posted @ 2016-02-15 23:06 Claris 阅读(...) 评论(...) 编辑 收藏