【bzoj4407】【于神之怒加強版】【莫比烏斯反演】

Description

給下N,M,K.求

∑i=1n∑j=1mgcd(i,j)kmod(109+7)

Input

輸入有多組資料，輸入資料的第一行兩個正整數T,K，代表有T組資料，K的意義如上所示，下面第二行到第T+1行，每行為兩個正整數N,M，其意義如上式所示。

Output

如題

Sample Input

1 2

3 3

Sample Output

20

HINT

1<=N,M,K<=5000000,1<=T<=2000

題解：

首先可以化成

∑d=1ndk∑i=1⌊nd⌋⌊nid⌋⌊mid⌋u(i)

設T=i*d

則可以化成

∑T=1n⌊nT⌋⌊mT⌋∑d|Tdku(Td)

設f(T)= ∑d|Tdku(Td)

隻要求出f數組,處理每組查詢就是o( n√) 了。

考慮 dk 是積性函數,是以f(T)也是積性函數。

是以我們就可以線上篩中處理f數組。

然後搞個字首和即可。

#include<iostream>
#include<cstdio>
#include<cstring>
#define N 5000010
#define P 1000000007 
#define LL long long
using namespace std;
int T,k,p[N],f[N],pos;
LL n,m,g[N],u[N],ans;
LL power(LL a,int b){
 LL ans;a%=P;
 for (ans=;b;a=a*a%P,b>>=) if (b&) ans=ans*a%P;
 return ans;
}
void pre(){
  u[]=;g[]=;
  for (int i=;i<=n;i++){
   if (!f[i]){p[++p[]]=i;u[i]=-;g[i]=(power(i,k)-+P)%P;}
   for (int j=;j<=p[]&&i*p[j]<=n;j++){
    f[i*p[j]]=; 
    if (i%p[j]==){g[i*p[j]]=g[i]*power(p[j],k)%P;break;}
    u[i*p[j]]=-u[i];g[i*p[j]]=(g[i]*g[p[j]])%P;
   } 
  }
  for (int i=;i<=n;i++) g[i]=(g[i]+g[i-])%P;
}
int main(){
 scanf("%d%d",&T,&k);n=N-;pre();
 while (T--){
  scanf("%d%d",&n,&m);
  if (n>m) swap(n,m);ans=;pos=; 
  for (int i=;i<=n;i=pos+){
    pos=min(n/(n/i),m / (m / i));
    ans+=(LL)(n / i)*(LL)(m / i)%P*(LL)(g[pos]-g[i-]+P)%P;
    ans%=P;
  }
  printf("%lld\n",ans);  
 }  
}