YY的GCD

时间:2022-02-14 07:58:43

YY的GCD

给出T个询问,询问\(\sum_{i=1}^N\sum_{j=1}^M(gcd(i,j)\in prime)\),T = 10000,N, M <= 10000000。

显然质数是需要枚举的,设N<M,于是

\[ans=\sum_{p\in prime}\sum_{i=1}^N\sum_{j=1}^M(gcd(i,j)==p)
\]

于是设

\[f(p)=\sum_{i=1}^N\sum_{j=1}^M(gcd(i,j)==p)
\]

\[F(p)=\sum_{i=1}^N\sum_{j=1}^M(p|gcd(i,j))=[N/p][M/p]
\]

由Mobius反演定理,我们有

\[f(p)=\sum_{p|d}F(d)\mu(d/p)
\]

于是

\[ans=\sum_{p\in prime}\sum_{p|d}F(d)\mu(d/p)=\sum_{d=1}^NF(d)\sum_{p|d,p\in 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('\n');
}
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();
}