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;
}