傳送門

原來……莫比烏斯反演是這麽用的啊……（雖然仍然不是很明白）

首先，題目所求如下$$\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}))$求一個前綴和，就好了

這數學公式真的是打的我心力憔悴……

 1 //minamoto
 2 #include<cstdio>
 3 #define ll long long
 4 //#define min(a,b) ((a)<(b)?(a):(b))
 5 inline int min(int a,int b){return a<b?a:b;}
 6 const int N=1e7+5;
 7 int vis[N],mu[N],p[N],g[N],m;ll sum[N],ans;
 8 void init(int n){
 9     mu[1]=1;
10     for(int i=2;i<=n;++i){
11         if(!vis[i]) p[++m]=i,mu[i]=-1;
12         for(int j=1;j<=m&&p[j]*i<=n;++j){
13             vis[i*p[j]]=1;
14             if(i%p[j]==0) break;
15             mu[i*p[j]]=-mu[i];
16         }
17     }
18     for(int j=1;j<=m;++j)
19     for(int i=1;i*p[j]<=n;++i)
20     g[i*p[j]]+=mu[i];
21     for(int i=1;i<=n;++i)
22     sum[i]=sum[i-1]+g[i];
23 }
24 int main(){
25 //    freopen("testdata.in","r",stdin);
26     init(1e7);
27     int T;scanf("%d",&T);
28     while(T--){
29         int n,m;
30         scanf("%d%d",&n,&m);
31         if(n>m) n^=m^=n^=m;
32         ans=0;
33         for(int l=1,r;l<=n;l=r+1){
34             r=min(n/(n/l),m/(m/l));
35             ans+=1ll*(n/l)*(m/l)*(sum[r]-sum[l-1]);
36         }
37         printf("%lld\n",ans);
38     }
39     return 0;
40 }

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