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*}\]
设$f(n)=n^k$,则$F$为$f$和$\mu$的狄利克雷卷积。
对于$F$的计算,质数时直接快速幂,质数的幂递推计算,其它数可以通过线性筛得到。
对于每次询问,分块求和即可。
时间复杂度$O(n+T\sqrt{n})$。
#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; }