洛谷P2257 YY的GCD（莫比乌斯反演）

传送门

 

原来……莫比乌斯反演是这么用的啊……（虽然仍然不是很明白）

首先，题目所求如下$$\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=prim]$$

我们设$f(d)$表示$gcd(i,j)=d$的$(i,j)$的对数，$g(d)$表示存在公因数为$d$的$(i,j)$的对数

那么就有$$f(d)=\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=d]$$

$$g(d)=\sum_{d|k}f(k)=\lfloor\frac{N}{d}\rfloor\lfloor\frac{M}{d}\rfloor$$

那么根据莫比乌斯反演定理，则有$$f(n)=\sum_{n\mid d}\mu(\frac{d}{n})g(d)$$

然后就可以化简了$$ans=\sum_{p\in prim}\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=p]$$

将$f(p)$代入，得$$ans=\sum_{p\in prim}f(p)$$

$$ans=\sum_{p\in prim}\sum_{p\mid d}\mu(\frac{d}{p})g(d)$$

我们考虑换一个枚举项，不枚举$p$，枚举$\frac{d}{p}$

$$ans=\sum_{p\in prim}\sum_{d=1}^{min(\lfloor\frac{n}{p}\rfloor,\lfloor\frac{m}{p}\rfloor)}\mu(d)g(dp)\\=\sum_{p\in prim}\sum_{d=1}^{min(\lfloor\frac{n}{p}\rfloor,\lfloor\frac{m}{p}\rfloor)}\mu(d)\lfloor\frac{n}{dp}\rfloor\lfloor\frac{m}{dp}\rfloor$$

然后我们把$dp$给换成$T$

$$ans=\sum_{T=1}^{min(n,m)}\sum_{t|T,t\in prim}\mu(\lfloor\frac{T}{t}\rfloor)\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor$$

$$ans=\sum_{T=1}^{min(n,m)}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor(\sum_{t|T,t\in prim}\mu(\lfloor\frac{T}{t}\rfloor))$$

然后对$(\sum_{t\mid T,t\in prim}\mu(\frac{T}{t}))$求一个前缀和，就好了

这数学公式真的是打的我心力憔悴……

//minamoto
#include<cstdio>
#define ll long long
//#define min(a,b) ((a)<(b)?(a):(b))
inline int min(int a,int b){return a<b?a:b;}
const int N=1e7+5;
int vis[N],mu[N],p[N],g[N],m;ll sum[N],ans;
void init(int n){
    mu[1]=1;
    for(int i=2;i<=n;++i){
        if(!vis[i]) p[++m]=i,mu[i]=-1;
        for(int j=1;j<=m&&p[j]*i<=n;++j){
            vis[i*p[j]]=1;
            if(i%p[j]==0) break;
            mu[i*p[j]]=-mu[i];
        }
    }
    for(int j=1;j<=m;++j)
    for(int i=1;i*p[j]<=n;++i)
    g[i*p[j]]+=mu[i];
    for(int i=1;i<=n;++i)
    sum[i]=sum[i-1]+g[i];
}
int main(){
//    freopen("testdata.in","r",stdin);
    init(1e7);
    int T;scanf("%d",&T);
    while(T--){
        int n,m;
        scanf("%d%d",&n,&m);
        if(n>m) n^=m^=n^=m;
        ans=0;
        for(int l=1,r;l<=n;l=r+1){
            r=min(n/(n/l),m/(m/l));
            ans+=1ll*(n/l)*(m/l)*(sum[r]-sum[l-1]);
        }
        printf("%lld\n",ans);
    }
    return 0;
}