给出T个询问,询问(sum_{i=1}^Nsum_{j=1}^M(gcd(i,j)in prime)),T = 10000,N, M <= 10000000。
解
显然质数是需要枚举的,设N<M,于是
[ans=sum_{pin prime}sum_{i=1}^Nsum_{j=1}^M(gcd(i,j)==p)
]
于是设
[f(p)=sum_{i=1}^Nsum_{j=1}^M(gcd(i,j)==p)
]
[F(p)=sum_{i=1}^Nsum_{j=1}^M(p|gcd(i,j))=[N/p][M/p]
]
由Mobius反演定理,我们有
[f(p)=sum_{p|d}F(d)mu(d/p)
]
于是
[ans=sum_{pin prime}sum_{p|d}F(d)mu(d/p)=sum_{d=1}^NF(d)sum_{p|d,pin prime}mu(d/p)
]
显然后式是可以(O(nlon(n)))维护的,对F(d)进行整除分块即可。
参考代码
#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define ll long long
#define limit 10000000
#define swap(x,y) x^=y^=x^=y
using namespace std;
bool check[limit+1];
int mb[limit+1],prime[700000],pt,
opt[limit+1];ll ans;
il void read(int&);
il int min(int,int);
void pen(ll),prepare();
int main(){
int lsy,a,b,i,j;read(lsy);
prepare();while(lsy--){
read(a),read(b),ans&=0;
if(a>b)swap(a,b);
for(i=1;i<=a;i=j+1){
j=min(a/(a/i),b/(b/i));
ans+=(ll)(a/i)*(b/i)*(opt[j]-opt[i-1]);
}pen(ans),putchar('
');
}
return 0;
}
void prepare(){
check[1]=mb[1]=1;
for(ri int i(2),j;i<=limit;++i){
if(!check[i])prime[++pt]=i,mb[i]=-1;
for(j=1;j<=pt&&prime[j]<=limit/i;++j){
check[i*prime[j]]|=true;
if(!(i%prime[j]))break;
mb[i*prime[j]]=-mb[i];
}
}for(ri int i(1),j;i<=pt;++i)
for(j=1;j*prime[i]<=limit;++j)
opt[j*prime[i]]+=mb[j];
for(ri int i(1);i<=limit;++i)opt[i]+=opt[i-1];
}
il int min(int a,int b){
return a<b?a:b;
}
void pen(ll x){
if(x>9)pen(x/10);putchar(x%10+48);
}
il void read(int &x){
x&=0;ri char c;while(c=getchar(),c<'0'||c>'9');
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
}