bzoj 2820 luogu 2257 yy的gcd (莫比乌斯反演)

Posted 2021-01-20 guapisolo

题目大意：求$gcd(i,j)==k,iin[1,n],jin[1,m] ,kin prime,n,m<=10^{7}$的有序数对个数，不超过10^{4}次询问

莫比乌斯反演入门题

为方便表述，由于n和m等价，以下内容均默认n<=m

题目让我们求：$sum_{k=1}^{n}sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)==k]$

容易变形为：$sum_{k=1}^{n}sum_{i=1}^{left lfloor frac{n}{k}  ight floor}sum_{j=1}^{left lfloor frac{m}{k}  ight floor}[gcd(i,j)==1]$

很显然$gcd(i,j)==1$这个东西不太好处理，因为它是最大公约数，直接统计gcd的数量很困难

所以，根据莫比乌斯反演的常规套路：

把$gcd(i,j)==1$变形为$sum _{d|gcd(i,j)}mu(d)$

可得$sum_{k=1}^{n}sum_{i=1}^{left lfloor frac{n}{k}  ight floor}sum_{j=1}^{left lfloor frac{m}{k}  ight floor}sum _{d|gcd(i,j)}mu(d)$

此时，统计$d$的数量就很容易了，即$left lfloor frac{n}{k}  ight floorcdot left lfloor frac{m}{k}  ight floor$，数学含义是"[1,n]中能整除d的数的数量*[1,m]中能整除d的数的数量"

继续变形$sum_{k=1}^{n}sum_{d=1}^{left lfloor frac{n}{k} ight floor}left lfloor frac{n}{kd} ight floor cdot left lfloor frac{m}{kd} ight floorcdot mu(d)$

我们看kd这个东西很不爽，把它换掉，令Q=kd

$sum_{Q}^{n}left lfloor frac{n}{Q} ight floor cdot left lfloor frac{m}{Q} ight floorcdot sum_{k|Q,kin prime}mu(frac{Q}{k})$

$sum_{k|Q,kin prime}mu(frac{Q}{k})$可以用线性筛筛出，将其命名为$g(i)$,再维护一个前缀和

在线性筛内分三种情况讨论，对于数i，我们遍历到了一个质数$p$，$i$中$p$的幂次为$k$

$k=0$，$i$中不含$p$，根据莫比乌斯函数的性质,$g(icdot p)=-g(i)+mu(i)$，即$g/i$的部分全都反过来，再加上$p$带来的贡献

$k=1$，$i$中含一个$p$，那么在$p$加入后，除了$p$的其它质因子的贡献全部为0，仅计算$p$的贡献，即$g(icdot p)=g(i)$

$k>=2$，$i$中含两个$p$，加入$p$后，所有质因子的贡献均为0

而$sum_{Q}^{n}left lfloor frac{n}{Q} ight floor cdot left lfloor frac{m}{Q} ight floor$可以使用整除分块（套路的味道）在$O(sqrt n)$的时间内算出来，再用前缀和统计即可

代码

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <algorithm>
 4 #define N 10010000
 5 #define maxn 10000000
 6 #define ll long long 
 7 using namespace std;
 8 
 9 int T,n,m,cnt;
10 int mu[N],pr[N],use[N],qm[N];
11 int pm[N];
12 
13 void Pre()
14 {
15     mu[1]=1,qm[1]=1;
16     for(int i=2;i<=maxn;i++)
17     {
18         if(!use[i]) pr[++cnt]=i,mu[i]=-1,qm[i]=1;
19         for(int j=1;j<=cnt&&i*pr[j]<=maxn;j++){
20             use[i*pr[j]]=1;
21             if(i%pr[j]==0){
22                 mu[i*pr[j]]=0;
23                 if((i/pr[j])%pr[j]==0) qm[i*pr[j]]=0;
24                 else qm[i*pr[j]]=mu[i];
25                 break;
26             }else{
27                 mu[i*pr[j]]=-mu[i];
28                 qm[i*pr[j]]=-qm[i]+mu[i];
29             }
30         }
31     }
32     for(int i=1;i<=maxn;i++)
33         pm[i]=pm[i-1]+qm[i];
34 }
35 ll solve(int n,int m)
36 {
37     ll ans=0;//int nd,md;
38     for(int i=2,la,mi=min(n,m);i<=mi;i=la+1)
39     {
40         la=min(n/(n/i),m/(m/i));
41         ans+=1ll*(n/i)*(m/i)*(pm[la]-pm[i-1]);
42     }
43     return ans;
44 }
45 int main()
46 {
47     scanf("%d",&T);
48     Pre();
49     for(int t=1;t<=T;t++)
50     {
51         scanf("%d%d",&n,&m);
52         if(n>m) swap(n,m);
53         printf("%lld
",solve(n,m));
54     }
55     return 0;
56 }