【洛谷2257】YY的GCD（莫比乌斯反演）

点此看题面

大致题意： 求${\sum }_{x=1}^N{\sum }_{y=1}^{M}IsPrime(gcd(x,y))$。

---

莫比乌斯反演

听说此题是莫比乌斯反演入门题？

$Link$

莫比乌斯反演 详见博客 初学莫比乌斯反演

---

一些定义

首先，我们可以定义$f(d)$和$F(d)$如下：

$$f(d)=\sum _{i=1}^N\sum _{j=1}^M[gcd(i,j)==d]$$

$$F(d)=\sum _{i=1}^N\sum _{j=1}^M[d|gcd(i,j)]$$

通过定义，不难发现：

$$F(n)=\sum _{n|d}f(d)=\lfloor \frac{N}{n}\rfloor \lfloor \frac{M}{n}\rfloor$$

由于莫比乌斯反演的某些性质，我们又可以得到：

$$f(n)=\sum _{n|d}\mu (\lfloor \frac{d}{n}\rfloor )F(d)$$

---

公式化简

首先，我们应该不难想到：

$$answer=\sum _{IsPrime(p)}\sum _{i=1}^N\sum _{j=1}^M[gcd(i,j)==p]$$

然后就是一波化简。

应该挺容易看出，由于$f(p)$的定义，上面的式子其实就相当于下面这个式子：

$answer={\sum }_{IsPrime(p)}f(p)$

然后是莫比乌斯反演：

$$answer=\sum _{IsPrime(p)}\sum _{p|d}\mu (\lfloor \frac{d}{p}\rfloor )F(d)$$

但是，这样有点难以处理。

于是，我们改成枚举$\lfloor \frac{d}{p}\rfloor$，于是原式就变成了这样：

$$answer=\sum _{IsPrime(p)}\sum _{d=1}^{min(\lfloor \frac{N}{p}\rfloor ,\lfloor \frac{M}{p}\rfloor )}\mu (d)F(d\cdot p)$$

将$F(n)=\lfloor \frac{N}{n}\rfloor \lfloor \frac{M}{n}\rfloor$代入进一步化简，可以得到：

$$answer=\sum _{IsPrime(p)}\sum _{d=1}^{min(\lfloor \frac{N}{p}\rfloor ,\lfloor \frac{M}{p}\rfloor )}\mu (d)\lfloor \frac{N}{d\cdot p}\rfloor \lfloor \frac{M}{d\cdot p}\rfloor$$

如果我们用$G$来表示$d\cdot p$，则$d=\frac{G}{p}$，原式就变成了这个样子：

$$answer=\sum _{G=1}^{min(N,M)}\sum _{IsPrime(p),p|G}\mu (\frac{G}{p})\lfloor \frac{N}{G}\rfloor \lfloor \frac{M}{G}\rfloor$$

通过乘法交换律和乘法结合律，我们可以再一步转化，得：

$$answer=\sum _{G=1}^{min(N,M)}\lfloor \frac{N}{G}\rfloor \lfloor \frac{M}{G}\rfloor (\sum _{IsPrime(p),p|G}\mu (\frac{G}{p}))$$

然后就可以$O(n)$求解了。

$But$，多组数据，$T\le 10000$… …

---

求解答案

首先，我们用定义一个$g(n)$，它的定义如下：

$$g(n)=\sum _{IsPrime(p),p|n}\mu (\frac{n}{p})$$

于是，我们便能将上面的式子进一步转化：

$$answer=\sum _{G=1}^{min(N,M)}\lfloor \frac{N}{G}\rfloor \lfloor \frac{M}{G}\rfloor g(n)$$

然后，我们可以用除法分块。

不难发现，在一定范围内$\lfloor \frac{N}{i}\rfloor$的值是保持不变的（$\lfloor \frac{M}{i}\rfloor$同理），则$\lfloor \frac{N}{G}\rfloor \lfloor \frac{M}{G}\rfloor$其实最多只有$\sqrt{N}+\sqrt{M}$，而对于$\lfloor \frac{N}{G}\rfloor \lfloor \frac{M}{G}\rfloor$相同的值，我们可以一起求解，于是就能想到用$su{m}_i$来表示${\sum }_{i=1}^{i}g(i)$，这样就能快速求解了。

---

代码

#include<bits/stdc++.h> #define max(x,y) ((x)>(y)?(x):(y)) #define min(x,y) ((x)<(y)?(x):(y)) #define uint unsigned int #define LL long long #define ull unsigned long long #define swap(x,y) (x^=y,y^=x,x^=y) #define abs(x) ((x)<0?-(x):(x)) #define INF 1e9 #define Inc(x,y) ((x+=(y))>=MOD&&(x-=MOD)) #define ten(x) (((x)<<3)+((x)<<1)) #define N 10000000 using namespace std; int n,m; class FIO { private: #define Fsize 100000 #define tc() (FinNow==FinEnd&&(FinEnd=(FinNow=Fin)+fread(Fin,1,Fsize,stdin),FinNow==FinEnd)?EOF:*FinNow++) #define pc(ch) (FoutSize<Fsize?Fout[FoutSize++]=ch:(fwrite(Fout,1,FoutSize,stdout),Fout[(FoutSize=0)++]=ch)) int f,FoutSize,OutputTop;char ch,Fin[Fsize],*FinNow,*FinEnd,Fout[Fsize],OutputStack[Fsize]; public: FIO() {FinNow=FinEnd=Fin;} inline void read(int &x) {x=0,f=1;while(!isdigit(ch=tc())) f=ch^'-'?1:-1;while(x=ten(x)+(ch&15),isdigit(ch=tc()));x*=f;} inline void read_char(char &x) {while(isspace(x=tc()));} inline void read_string(string &x) {x="";while(isspace(ch=tc()));while(x+=ch,!isspace(ch=tc())) if(!~ch) return;} inline void write(LL x) {if(!x) return (void)pc('0');if(x<0) pc('-'),x=-x;while(x) OutputStack[++OutputTop]=x%10+48,x/=10;while(OutputTop) pc(OutputStack[OutputTop]),--OutputTop;} inline void write_char(char x) {pc(x);} inline void write_string(string x) {register int i,len=x.length();for(i=0;i<len;++i) pc(x[i]);} inline void end() {fwrite(Fout,1,FoutSize,stdout);} }F; class Class_Mobius//莫比乌斯反演 { private: int Prime_cnt,mu[N+5],Prime[N+5];bool IsNotPrime[N+5]; public: LL sum[N+5]; Class_Mobius()//预处理 { register int i,j; for(mu[1]=1,i=2;i<=N;++i)//求出莫比乌斯函数 { if(!IsNotPrime[i]) Prime[++Prime_cnt]=i,mu[i]=-1; for(j=1;j<=Prime_cnt&&i*Prime[j]<=N;++j) if(IsNotPrime[i*Prime[j]]=true,i%Prime[j]) mu[i*Prime[j]]=-mu[i];else break; } for(j=1;j<=Prime_cnt;++j) for(i=Prime[j];i<=N;i+=Prime[j]) sum[i]+=mu[i/Prime[j]];//计算g(i) for(i=1;i<=N;++i) sum[i]+=sum[i-1];//求前缀和 } }Mobius; int main() { register int i,nxt,T;register LL ans;F.read(T); while(T--) { if(F.read(n),F.read(m),n>m) swap(n,m); for(ans=0,i=1;i<=n;i=nxt+1) nxt=min(n/(n/i),m/(m/i)),ans+=1LL*(n/i)*(m/i)*(Mobius.sum[nxt]-Mobius.sum[i-1]);//除法分块 F.write(ans),F.write_char('\n');//输出答案 } return F.end(),0; }