【系列】莫比乌斯反演

几个有用的结论：

记狄利克雷卷积： $(f*g)(n)=sumlimits_{d|n} f(d) * g(frac{n}{d})$

则有$mu * 1 = [n=1]$与$phi *1 = id$

若$F(n)=sumlimits_{d|n} f(d)$

则$f(n)=sumlimits_{d|n} mu(d) * F(frac{n}{d})$

bzoj 2190 仪仗队

题目大意：

求$n*n$的矩形中能从左下角被直接看到的点的个数

思路：

设左下角$(0,0)$ 则相当于求$sumlimits_{i=1}^{n-1} sumlimits_{j=1}^{n-1} [gcd(i,j)==1]+2$（$(1,0),(0,1)$

将式子转化为$1+ 2* sumlimits_{i=1}^{n-1} sumlimits_{j=1}^{i} [gcd(i,j)==1] $（由于有$(1,0),(0,1),(1,1)$三个特殊点存在

然后发现满足欧拉函数的定义，则所求为$1+2* sumlimits_{i=1}^{n-1} phi(i) $ 筛出欧拉函数即可

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstring>
 4 #include<cstdlib>
 5 #include<cmath>
 6 #include<algorithm>
 7 #include<queue>
 8 #include<vector>
 9 #include<map>
10 #include<set>
11 #define ll long long
12 #define db double
13 #define inf 2139062143
14 #define MAXN 50010
15 #define rep(i,s,t) for(register int i=(s),i##__end=(t);i<=i##__end;++i)
16 #define dwn(i,s,t) for(register int i=(s),i##__end=(t);i>=i##__end;--i)
17 #define ren for(register int i=fst[x];i;i=nxt[i])
18 #define pb(i,x) vec[i].push_back(x)
19 #define pls(a,b) (a+b)%MOD
20 #define mns(a,b) (a-b+MOD)%MOD
21 #define mul(a,b) (1LL*(a)*(b))%MOD
22 using namespace std;
23 inline int read()
24 {
25     int x=0,f=1;char ch=getchar();
26     while(!isdigit(ch)) {if(ch=='-') f=-1;ch=getchar();}
27     while(isdigit(ch)) {x=x*10+ch-'0';ch=getchar();}
28     return x*f;
29 }
30 int n,tot,p[MAXN+100],ntp[MAXN+100],phi[MAXN+100];
31 void mem()
32 {
33     phi[1]=1;
34     rep(i,2,MAXN)
35     {
36         if(!ntp[i]) p[++tot]=i,phi[i]=i-1;
37         rep(j,1,tot) if(i*p[j]>MAXN) break;
38             else {ntp[i*p[j]]=1;if(i%p[j]) phi[i*p[j]]=phi[i]*phi[p[j]];else {phi[i*p[j]]=phi[i]*p[j];break;}}
39     }
40     rep(i,2,n-1) phi[i]+=phi[i-1];
41 }
42 int main()
43 {
44     n=read();mem();printf("%d
",n>1?phi[n-1]<<1|1:0);
45 }

View Code

bzoj 2818 Gcd

题目大意：

求有序数对$(i,j),i,j leq n$满足$gcd(i,j)$为素数的数对个数

思路：

法1：和上一题很相似，设$p$为素数则$sumlimits_{i=1}^n sumlimits_{j=1}^n [gcd(i,j)==p] =sumlimits_{i=1}^{n/p} sumlimits_{j=1}^{n/p} [gcd(i,j)==1]$

这样我们枚举每一个素数就可以解决了

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstring>
 4 #include<cstdlib>
 5 #include<cmath>
 6 #include<algorithm>
 7 #include<queue>
 8 #include<vector>
 9 #include<map>
10 #include<set>
11 #define ll long long
12 #define db double
13 #define inf 2139062143
14 #define MAXN 10001000
15 #define rep(i,s,t) for(register int i=(s),i##__end=(t);i<=i##__end;++i)
16 #define dwn(i,s,t) for(register int i=(s),i##__end=(t);i>=i##__end;--i)
17 #define ren for(register int i=fst[x];i;i=nxt[i])
18 #define pb(i,x) vec[i].push_back(x)
19 #define pls(a,b) (a+b)%MOD
20 #define mns(a,b) (a-b+MOD)%MOD
21 #define mul(a,b) (1LL*(a)*(b))%MOD
22 using namespace std;
23 inline int read()
24 {
25     int x=0,f=1;char ch=getchar();
26     while(!isdigit(ch)) {if(ch=='-') f=-1;ch=getchar();}
27     while(isdigit(ch)) {x=x*10+ch-'0';ch=getchar();}
28     return x*f;
29 }
30 int n,tot,p[MAXN],ntp[MAXN+5],phi[MAXN];
31 ll ans,sum[MAXN];
32 void mem()
33 {
34     phi[1]=1;
35     rep(i,2,n)
36     {
37         if(!ntp[i]) p[++tot]=i,phi[i]=i-1;
38         rep(j,1,tot) if(i*p[j]>n) break;
39             else {ntp[i*p[j]]=1;if(i%p[j]) phi[i*p[j]]=phi[i]*phi[p[j]];else {phi[i*p[j]]=phi[i]*p[j];break;}}
40     }
41     rep(i,1,n) sum[i]=sum[i-1]+phi[i];
42 }
43 int main()
44 {
45     n=read();mem();
46     rep(i,1,tot) ans+=sum[n/p[i]]*2-1;printf("%lld
",ans);
47 }

View Code

法2：因为$gcd(i,j)==1$的形式可以转化为$sumlimits_{d|gcd(i,j)} mu(d)$

由于对于每一个$d$，相对应的$i,j$有$n/i,n/j$个。所以前两个枚举是没有必要的，设$m$为$n/p$，式子为$sumlimits_{d=1}^m mu(d)*frac{m}{d}*frac{m}{d}$

枚举素数后数论分块

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstring>
 4 #include<cstdlib>
 5 #include<cmath>
 6 #include<algorithm>
 7 #include<queue>
 8 #include<vector>
 9 #include<map>
10 #include<set>
11 #define ll long long
12 #define db double
13 #define inf 2139062143
14 #define MAXN 10001000
15 #define rep(i,s,t) for(register int i=(s),i##__end=(t);i<=i##__end;++i)
16 #define dwn(i,s,t) for(register int i=(s),i##__end=(t);i>=i##__end;--i)
17 #define ren for(register int i=fst[x];i;i=nxt[i])
18 #define pb(i,x) vec[i].push_back(x)
19 #define pls(a,b) (a+b)%MOD
20 #define mns(a,b) (a-b+MOD)%MOD
21 #define mul(a,b) (1LL*(a)*(b))%MOD
22 using namespace std;
23 inline int read()
24 {
25     int x=0,f=1;char ch=getchar();
26     while(!isdigit(ch)) {if(ch=='-') f=-1;ch=getchar();}
27     while(isdigit(ch)) {x=x*10+ch-'0';ch=getchar();}
28     return x*f;
29 }
30 int n,tot,p[MAXN],ntp[MAXN+5],mu[MAXN];
31 ll ans,sum[MAXN];
32 void mem()
33 {
34     mu[1]=1;
35     rep(i,2,n)
36     {
37         if(!ntp[i]) p[++tot]=i,mu[i]=-1;
38         rep(j,1,tot) if(i*p[j]>n) break;
39             else {ntp[i*p[j]]=1;if(i%p[j]) mu[i*p[j]]=-mu[i];else {mu[i*p[j]]=0;break;}}
40     }
41     rep(i,1,n) sum[i]=sum[i-1]+mu[i];
42 }
43 ll solve(int x,ll res=0)
44 {
45     int pos,lmt=n/x;rep(i,1,lmt) {pos=lmt/(lmt/i);res+=1LL*(sum[pos]-sum[i-1])*(lmt/i)*(lmt/i);i=pos;}
46     return res;
47 }
48 int main()
49 {
50     n=read();mem();    
51     rep(i,1,tot) ans+=solve(p[i]);printf("%lld
",ans);
52 }

View Code

bzoj 1101 Zap

题目大意：

求有序数对$(i,j),ileq a,j leq b$满足$gcd(i,j)==d$的数对个数

思路：

和上一题基本一样

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstring>
 4 #include<cstdlib>
 5 #include<cmath>
 6 #include<algorithm>
 7 #include<queue>
 8 #include<vector>
 9 #include<map>
10 #include<set>
11 #define ll long long
12 #define db double
13 #define inf 2139062143
14 #define MAXN 60100
15 #define rep(i,s,t) for(register int i=(s),i##__end=(t);i<=i##__end;++i)
16 #define dwn(i,s,t) for(register int i=(s),i##__end=(t);i>=i##__end;--i)
17 #define ren for(register int i=fst[x];i;i=nxt[i])
18 #define pb(i,x) vec[i].push_back(x)
19 #define pls(a,b) (a+b)%MOD
20 #define mns(a,b) (a-b+MOD)%MOD
21 #define mul(a,b) (1LL*(a)*(b))%MOD
22 using namespace std;
23 inline int read()
24 {
25     int x=0,f=1;char ch=getchar();
26     while(!isdigit(ch)) {if(ch=='-') f=-1;ch=getchar();}
27     while(isdigit(ch)) {x=x*10+ch-'0';ch=getchar();}
28     return x*f;
29 }
30 int n,tot,p[MAXN],ntp[MAXN+5],mu[MAXN];
31 ll ans,sum[MAXN];
32 void mem(int n)
33 {
34     mu[1]=1;
35     rep(i,2,n)
36     {
37         if(!ntp[i]) p[++tot]=i,mu[i]=-1;
38         rep(j,1,tot) if(i*p[j]>n) break;
39             else {ntp[i*p[j]]=1;if(i%p[j]) mu[i*p[j]]=-mu[i];else {mu[i*p[j]]=0;break;}}
40     }
41     rep(i,1,n) sum[i]=sum[i-1]+mu[i];
42 }
43 ll solve(int x,int y,ll res=0)
44 {
45     int pos,lmt=min(x,y);rep(i,1,lmt)
46         {pos=min(x/(x/i),(y/(y/i)));res+=1LL*(sum[pos]-sum[i-1])*(x/i)*(y/i);i=pos;}
47     return res;
48 }
49 int main()
50 {
51     int T=read(),a,b,d;mem(50010);while(T--)
52         {a=read(),b=read(),d=read();printf("%lld
",solve(a/d,b/d));}
53 }

View Code

bzoj 4804 欧拉心算

题目大意：

求$sumlimits_i^n sumlimits_j^n phi(gcd(i,j))$

思路：

枚举$gcd$ 得到$sumlimits_{k=1}^n phi(k) sumlimits_{i=1}^n sumlimits_{j=1}^n [gcd(i,j)==k]$

后面的式子转化为第一题的样子即原式为$sumlimits_{k=1}^n phi(k) (-1+2*sumlimits_{i=1}^{n/k}phi(i))$

设$f(k)=-1+2*sumlimits_{i=1}^k phi(i)$ 则原式为$sumlimits_{k=1}^n phi(k)*f(frac{n}{k})$

使用数论分块即可

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstring>
 4 #include<cstdlib>
 5 #include<cmath>
 6 #include<algorithm>
 7 #include<queue>
 8 #include<vector>
 9 #include<map>
10 #include<set>
11 #define ll long long
12 #define db double
13 #define inf 2139062143
14 #define MAXN 10000001
15 #define rep(i,s,t) for(register int i=(s),i##__end=(t);i<=i##__end;++i)
16 #define dwn(i,s,t) for(register int i=(s),i##__end=(t);i>=i##__end;--i)
17 #define ren for(register int i=fst[x];i;i=nxt[i])
18 #define pb(i,x) vec[i].push_back(x)
19 #define pls(a,b) (a+b)%MOD
20 #define mns(a,b) (a-b+MOD)%MOD
21 #define mul(a,b) (1LL*(a)*(b))%MOD
22 using namespace std;
23 inline int read()
24 {
25     int x=0,f=1;char ch=getchar();
26     while(!isdigit(ch)) {if(ch=='-') f=-1;ch=getchar();}
27     while(isdigit(ch)) {x=x*10+ch-'0';ch=getchar();}
28     return x*f;
29 }
30 int n,tot,p[MAXN],ntp[MAXN];
31 ll ans,phi[MAXN],f[MAXN];
32 void mem(int n)
33 {
34     phi[1]=1;
35     rep(i,2,n)
36     {
37         if(!ntp[i]) p[++tot]=i,phi[i]=i-1;
38         rep(j,1,tot) if(i*p[j]>n) break;
39             else {ntp[i*p[j]]=1;if(i%p[j]) phi[i*p[j]]=phi[i]*phi[p[j]];else {phi[i*p[j]]=phi[i]*p[j];break;}}
40     }
41     rep(i,1,n) phi[i]+=phi[i-1],f[i]=-1+(phi[i]<<1);
42 }
43 ll solve(int n,ll res=0)
44 {
45     int pos;rep(i,1,n) {pos=n/(n/i);res+=1LL*(phi[pos]-phi[i-1])*f[n/i];i=pos;}return res;
46 }
47 int main()
48 {
49     int T=read(),a,b,d;mem(10000000);while(T--)
50         {a=read();printf("%lld
",solve(a));}
51 }

View Code

（卡我一个数组的空间常数怕不是有毒

bzoj 3529 数表

题目大意：

有一张 n×m 的数表，其第 i 行第 j 列（1 <= i <= n, 1 <= j <= m）的数值为能同时整除 i 和 j 的所有自然数之和

q次询问每次给出n m a 求数表中不大于 a 的数之和

思路：

首先设$f(i)=sumlimits_{d|i} d$

则所求为$sumlimits_{i=1}^n sumlimits_{j=1}^m f(gcd(i,j)) *[f(gcd(i,j))<=a]$

先不管最后那一个限制，继续推式子，我们枚举$gcd$：$sumlimits_{k=1}^n f(k) sumlimits_{i=1}^{n/k} sumlimits_{j=1}^{m/k} [gcd(i,j)==1]$

后面的式子好熟啊，变成莫比乌斯函数形式：$sumlimits_{k=1}^n f(k) sumlimits_d^{min(n/k,m/k)} mu(d) frac{n/k}{d} frac{m/k}{d}$

如果后面把$k$写在分数线下面，发现我们可以枚举$k*d$设为$t$ 得到：$sumlimits_{t=1}^{min(n,m)} frac{n}{t} frac{m}{t} sumlimits_{i|t} f(i)*mu(frac{t}{i})$

把后面的拿出来设$h(t)=sumlimits_{i|t} f(i)*mu(frac{t}{i})$ 原式为$sumlimits_{t=1}^{min(n,m)} h(t) frac{n}{t} frac{m}{t}$我们可以数论分块

对于$h(i)$，我们可以把询问的$a$限制从小到大排序，将所有$i$按照$f(i)$排序后一次加入

每一个$i$，枚举其倍数$t=k *i $即对$h(t)$造成的影响，其复杂度为调和级数，暴力使用树状数组维护对每个倍数的位置$pos(i*j)+=f(i)*mu(j)$

数论分块的时候前缀和用树状数组查询即可

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstring>
 4 #include<cstdlib>
 5 #include<cmath>
 6 #include<algorithm>
 7 #include<queue>
 8 #include<vector>
 9 #include<map>
10 #include<set>
11 #define ll long long
12 #define db double
13 #define inf 2139062143
14 #define MAXN 100100
15 #define MOD 2147483647
16 #define rep(i,s,t) for(register int i=(s),i##__end=(t);i<=i##__end;++i)
17 #define dwn(i,s,t) for(register int i=(s),i##__end=(t);i>=i##__end;--i)
18 #define ren for(register int i=fst[x];i;i=nxt[i])
19 #define pb(i,x) vec[i].push_back(x)
20 #define pls(a,b) (a+b)%MOD
21 #define mns(a,b) (a-b+MOD)%MOD
22 #define mul(a,b) (1LL*(a)*(b))%MOD
23 using namespace std;
24 inline int read()
25 {
26     int x=0,f=1;char ch=getchar();
27     while(!isdigit(ch)) {if(ch=='-') f=-1;ch=getchar();}
28     while(isdigit(ch)) {x=x*10+ch-'0';ch=getchar();}
29     return x*f;
30 }
31 int n,tot,p[MAXN],ntp[MAXN],mu[MAXN],lmt;
32 ll phi[MAXN],f[MAXN],g[MAXN],ans[MAXN],af[MAXN];
33 struct Ask{int id,n,m,a;}q[MAXN];
34 bool cmpa(Ask x,Ask y) {return x.a<y.a;}
35 bool cmp(int a,int b) {return f[a]<f[b];}
36 ll c[MAXN];
37 int lowbit(int x) {return x&(-x);}
38 void mdf(int x,int val) {for(;x<=lmt;x+=lowbit(x)) (c[x]+=val)&=MOD;}
39 ll query(int x,ll res=0) {for(;x;x-=lowbit(x)) (res+=c[x])&=MOD;return res;}
40 void mem(int n)
41 {
42     mu[1]=1;rep(i,2,n)
43     {
44         if(!ntp[i]) p[++tot]=i,mu[i]=-1;
45         rep(j,1,tot) if(i*p[j]>n) break;
46             else {ntp[i*p[j]]=1;if(i%p[j]) mu[i*p[j]]=-mu[i];else {mu[i*p[j]]=0;break;}}
47     }
48     rep(i,1,n) for(int j=af[i]=i;j<=n;j+=i) f[j]=(f[j]+i)&MOD;
49 }
50 ll solve(int n,int m,ll res=0)
51 {
52     int pos;rep(i,1,min(n,m)) 
53         {pos=min(n/(n/i),m/(m/i));(res+=1LL*(query(pos)-query(i-1))*(((n/i)*(m/i))&MOD))&=MOD;i=pos;}
54     return res;
55 }
56 int main()
57 {
58     int num=read(),a,b,d;mem(lmt=100000);
59     rep(i,1,num) q[i].n=read(),q[i].m=read(),q[i].a=read(),q[i].id=i;
60     sort(q+1,q+num+1,cmpa);sort(af+1,af+lmt+1,cmp);int pos=1;
61     rep(i,1,num)
62     {
63         while(pos<lmt&&f[af[pos]]<=q[i].a)
64             {rep(j,1,lmt/af[pos]) mdf(af[pos]*j,f[af[pos]]*mu[j]);pos++;}
65         ans[q[i].id]=solve(q[i].n,q[i].m);
66     }
67     rep(i,1,num) printf("%d
",ans[i]);
68 }

View Code

bzoj 2154 Crash的数字表格

题目大意：

求$sumlimits_{i=1}^n sumlimits_{j=1}^m lcm(i,j)$

思路：

先将$lcm$转化为$gcd$：$sumlimits_{i=1}^n sumlimits_{j=1}^m frac{i*j}{gcd(i,j)}$

令$n$为较小的那个数，然后枚举一波gcd：$sumlimits_{d=1}^n sumlimits_{i=1}^{n/d} sumlimits_{j=1}^{m/d} frac{id*jd}{d} *[gcd(i,j)==1]$

把$d$提到前面去接化简后面那个很熟的式子：$sumlimits_{d=1}^n d* sumlimits_{i=1}^{n/d} sumlimits_{j=1}^{m/d} sumlimits_{t|gcd(i,j)} i*j*mu(t)$

对于第二个$gcd$继续枚举：$sumlimits_{d=1}^n d * sumlimits_{t=1}^{n/d} sumlimits_{i=1}^{n/(dt)} sumlimits_{j=1}^{m/(dt)} it *jt*mu(t)$

提出$t$： $sumlimits_{d=1}^n d *( sumlimits_{t=1}^{n/d} t^2 * mu(t)) sumlimits_{i=1}^{n/(dt)} i *sumlimits_{j=1}^{m/(dt)} j$

设$s(t)=sumlimits_{i=1}^t i$，则：$sumlimits_{d=1}^n d * sumlimits_{t=1}^{n/d} t^2 * mu(t) * s(frac{n}{dt}) *s(frac{m}{dt})$

与前一题类似的，我们设$g=d*t$，得到：$sumlimits_{g=1}^n sumlimits_{d|g} d * mu(frac{g}{d}) * (frac{g}{d})^2 * s(frac{n}{g}) * s(frac{m}{g})$

把$g$提出来：$sumlimits_{g=1}^n g * s(frac{n}{g}) * s(frac{m}{g})sumlimits_{d|g} frac{g}{d}* mu(frac{g}{d}) $

把后面的部分单独拿出来，设$f(g)=sumlimits_{d|g} d * mu(d)$，发现$f$是$d*mu(d)$与$1$的狄利克雷卷积结果

由于那两个都是积性函数，因此$f$也是积性函数，可以通过筛法来求

对于素数$p$，两个因数分别$1,p$，因此$f(p)=-p+1$

对于线筛中$i%p[j]==0$的情况，考虑新增添的因数一定含有平方因子，因此其莫比乌斯函数一定为0，不需要考虑

最终式子为$sumlimits_{g=1}^n g * s(frac{n}{g}) * s(frac{m}{g}) *f(g) $然后愉快的数论分块即可

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstring>
 4 #include<cstdlib>
 5 #include<cmath>
 6 #include<algorithm>
 7 #include<queue>
 8 #include<vector>
 9 #include<map>
10 #include<set>
11 #define ll long long
12 #define db double
13 #define inf 2139062143
14 #define MAXN 10001000
15 #define MOD 20101009
16 #define rep(i,s,t) for(register int i=(s),i##__end=(t);i<=i##__end;++i)
17 #define dwn(i,s,t) for(register int i=(s),i##__end=(t);i>=i##__end;--i)
18 #define ren for(register int i=fst[x];i;i=nxt[i])
19 #define pb(i,x) vec[i].push_back(x)
20 #define pls(a,b) (a+b+MOD)%MOD
21 #define mns(a,b) (a-b+MOD)%MOD
22 #define mul(a,b) (1LL*(a)*(b))%MOD
23 using namespace std;
24 inline int read()
25 {
26     int x=0,f=1;char ch=getchar();
27     while(!isdigit(ch)) {if(ch=='-') f=-1;ch=getchar();}
28     while(isdigit(ch)) {x=x*10+ch-'0';ch=getchar();}
29     return x*f;
30 }
31 int n,m,tot,p[MAXN],ntp[MAXN],lmt;
32 ll f[MAXN],s[MAXN];
33 void mem(int n)
34 {
35     f[1]=1;rep(i,2,n)
36     {
37         if(!ntp[i]) p[++tot]=i,f[i]=-i+1;
38         rep(j,1,tot) if(i*p[j]>n) break;
39             else {ntp[i*p[j]]=1;if(i%p[j]) f[i*p[j]]=mul(f[i],f[p[j]]);else {f[i*p[j]]=f[i];break;}}
40     }
41     rep(i,1,n) s[i]=pls(s[i-1],i),f[i]=pls(f[i-1],mul(f[i],i));
42 }
43 ll solve(int n,int m,ll res=0)
44 {
45     int pos;rep(i,1,min(n,m)) 
46         {pos=min(n/(n/i),m/(m/i));res=pls(res,mul(f[pos]-f[i-1],mul(s[n/i],s[m/i])));i=pos;}
47     return res;
48 }
49 int main()
50 {
51     n=read(),m=read();mem(max(n,m));printf("%lld
",solve(n,m));
52 }

View Code

【系列】 莫比乌斯反演

【系列】莫比乌斯反演