傳送門

公式太長了……我就直接抄一下這位大佬好了……實在懶得打了

首先據說$d(ij)$有個性質$$d(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]$$

我們所求的答案為$$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}d(ij)$$

$$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]$$

考慮一下$gcd(x,y)=1$，我們可以考慮莫比烏斯函數的性質，那麽即$\sum_{d\mid n}\mu(d)$與$[n=1]$的結果相同

則有$$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x|i}\sum_{y|j}\sum_{d|gcd(x,y)}\mu(d)$$

然後我們由枚舉$gcd(x,y)$的約數改為直接枚舉$d$

$$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x|i}\sum_{y|j}\sum_{d=1}^{min(n,m)}\mu(d)*[d|gcd(x,y)]$$

然後把$\mu(d)$提取出來

$$ans=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x|i}\sum_{y|j}[d|gcd(x,y)]$$

然後，我們把枚舉$i,j$和約數改為直接枚舉約數，然後每個約數都會對他所有的倍數產生貢獻

$$ans=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^{n}\sum_{y=1}^{m}[d|gcd(x,y)]\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor$$

然後我們把枚舉$x,y$改為枚舉$dx,dy$，那麽就可以把$[d|gcd(x,y)]$這個條件給消掉

$$ans=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{{\lfloor\frac{m}{d}\rfloor}}\lfloor\frac{n}{dx}\rfloor\lfloor\frac{m}{dy}\rfloor$$

$$ans=\sum_{d=1}^{min(n,m)}\mu(d)(\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{dx}\rfloor)(\sum_{y=1}^{{\lfloor\frac{m}{d}\rfloor}}\lfloor\frac{m}{dy}\rfloor)$$

然後$\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{dx}\rfloor$和$\sum_{y=1}^{{\lfloor\frac{m}{d}\rfloor}}\lfloor\frac{m}{dy}\rfloor$的前綴和都可以預處理，直接上整除分塊就可以了

 1 //minamoto
 2 #include<iostream>
 3 #include<cstdio>
 4 #define ll long long
 5 using namespace std;
 6 #define getc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
 7 char buf[1<<21],*p1=buf,*p2=buf;
 8 inline int read(){
 9     #define num ch-‘0‘
10     char ch;bool flag=0;int res;
11     while(!isdigit(ch=getc()))
12     (ch==‘-‘)&&(flag=true);
13     for(res=num;isdigit(ch=getc());res=res*10+num);
14     (flag)&&(res=-res);
15     #undef num
16     return res;
17 }
18 const int N=50005;
19 int vis[N],p[N],mu[N],sum[N],m;
20 ll g[N],ans;
21 void init(int n){
22     mu[1]=1;
23     for(int i=2;i<=n;++i){
24         if(!vis[i]) p[++m]=i,mu[i]=-1;
25         for(int j=1;j<=m&&p[j]*i<=n;++j){
26             vis[i*p[j]]=1;
27             if(i%p[j]==0) break;
28             mu[i*p[j]]=-mu[i];
29         }
30     }
31     for(int i=1;i<=n;++i) sum[i]=sum[i-1]+mu[i];
32     for(int i=1;i<=n;++i){
33         ans=0;
34         for(int l=1,r;l<=i;l=r+1){
35             r=(i/(i/l));
36             ans+=1ll*(r-l+1)*(i/l);
37         }
38         g[i]=ans;
39     }
40 }
41 int main(){
42 //    freopen("testdata.in","r",stdin);
43     init(50000);
44     int n,m,T,lim;scanf("%d",&T);
45     while(T--){
46         scanf("%d%d",&n,&m);
47         lim=min(n,m),ans=0;
48         for(int l=1,r;l<=lim;l=r+1){
49             r=min(n/(n/l),m/(m/l));
50             ans+=(sum[r]-sum[l-1])*g[n/l]*g[m/l];
51         }
52         printf("%lld\n",ans);
53     }
54     return 0;
55 }

洛谷P3327 [SDOI2015]約數個數和（莫比烏斯反演）