莫比乌斯反演例题集 ^_^

【luogu 2257】YY的GCD

Problem Here

预备知识

除法分块、莫比乌斯反演

最终公式：

[ans= sum_{T=1}^{n} ( frac{n}{T}) ( frac{m}{T}) sum_{p|t,isprime[p]=1} mu ( frac{T}{p} ) ]

#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return x*f;
}
#define MN 10000000
int T,n,m;
int prime[MN+5],cnt,mu[MN+5];
bool mark[MN+5];
ll sum[MN+5],ans;
inline void init(){
	mark[1]=1;mu[1]=1;register int i,j;
	for(i=2;i<=MN;i++){
		if(!mark[i]) {mu[i]=-1;prime[++cnt]=i;}
		for(j=1;j<=cnt&&prime[j]*i<=MN;j++){
			mark[i*prime[j]]=1;
			if(i%prime[j]==0) break;
			else mu[i*prime[j]]=-mu[i];
		}
	}
	for(i=1;i<=cnt;i++)
	for(j=1;j*prime[i]<=MN;j++) sum[j*prime[i]]+=mu[j];
	for(i=2;i<=MN;i++) sum[i]+=sum[i-1];
}
int main(){
	T=read();init();
	register int i,r;
	while(T--){
		n=read();m=read();int MIN=min(n,m);
		ans=0ll;
		for(i=1;i<=MIN;i=r+1){
			r=min(n/(n/i),m/(m/i));
			ans+=(1ll)*(n/i)*(m/i)*(sum[r]-sum[i-1]);
		}
		printf("%lld
",ans);
	}
	return 0;
}

【HAOI2011】Problem b

Problem Here

预备知识

容斥

#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return x*f;
}
#define MN 50005
int prime[MN+5],cnt,mu[MN+5];
bool mark[MN+5];
inline void init()
{
	mark[1]=1;mu[1]=1;register int i,j;
	for(i=2;i<=MN;i++)
	{
		if(!mark[i]) {mu[i]=-1;prime[++cnt]=i;}
		for(j=1;j<=cnt&&prime[j]*i<=MN;j++)
		{
			mark[i*prime[j]]=1;
			if(i%prime[j]==0) break;
			else mu[i*prime[j]]=-mu[i];
		}
	}
	for(i=2;i<=MN;i++) mu[i]+=mu[i-1];
}
ll q(int n,int m,int k)
{
	ll ans=0ll;
	register int i,r;
	for(i=1;i*k<=n&&i*k<=m;i=r+1){
		r=min(n/(n/(i*k))/k,m/(m/(i*k))/k);
		ans+=1ll*(mu[r]-mu[i-1])*(n/(i*k))*(m/(i*k));
	}
	return ans;
}
int main()
{
	int n=read(),a,b,c,d,k;
	init();
	while(n--)
	{
		a=read(),b=read();
		c=read(),d=read();
		k=read();
		printf("%lld
",q(b,d,k)-q(a-1,d,k)-q(b,c-1,k)+q(a-1,c-1,k));
	}
	return 0;
}

【NOI2010】能量采集

Problem Here

预备知识

求gcd的新姿势

[sum_{d|n}^{} phi(d) =n ]

prove: 考虑小于等于n的数中，与n的gcd为(frac{n}{d})的数的个数是(phi(d)),而所有的(phi(d))相加正好是n

所有满足 i|x,i|y的数，均满足i|gcd(x,y)

//一个正整数，等于它所有因数的欧拉函数之和。
#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return x*f;
}
#define MN 100000
int prime[MN+5],cnt;
ll phi[MN+5],ans;
bool mark[MN+5];
inline void init()
{
	mark[1]=1;phi[1]=1;register int i,j;
	for(i=2;i<=MN;i++)
	{
		if(!mark[i]) {phi[i]=i-1;prime[++cnt]=i;}
		for(j=1;j<=cnt&&prime[j]*i<=MN;j++)
		{
			mark[i*prime[j]]=1;
			if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;}
            else phi[i*prime[j]]=(prime[j]-1)*phi[i];
		}
	}
	for(i=2;i<=MN;i++) phi[i]+=phi[i-1];
}
int main(){
	register int n,m;
	n=read(),m=read();init();
	register int i,r;
	for(i=1;i<=n&&i<=m;i=r+1){
		r=min(n/(n/i),m/(m/i));
		ans+=1ll*(n/i)*(m/i)*(phi[r]-phi[i-1]);
	}
	(ans<<=1)-=1ll*n*m;
	printf("%lld
",ans);
	return 0;
}

【luogu 1829】Crash的数字表格

Problem Here

预备知识

有一些莫比乌斯反演的的常见套路：

• 对于求和gcd的，把gcd提到最前面枚举
• 对于判断一个数是否为1，用$$sum_{d|n}^{} mu(d)=[n==1]$$来实现，之后在想办法把(mu(d))给提出来
• 把能够进行除法分块的部分提出来

最终公式:

[ans=sum_{T=1}^{n} Sum( frac{n}{T}) Sum( frac{m}{T}) T sum_{d|T} d mu(d) ]

而函数 $$g(x)= sum_{d|x}^{} d mu(d) $$是个积性函数，可以用欧拉筛来求。

//把能够除法分块的部分拖到最前面
//积性函数可以用欧拉筛来求 
#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return x*f;
}
#define MN 10000007
#define mod 20101009
int prime[MN],cnt;
ll f[MN],ans;
bool mark[MN];
inline void init()
{
	mark[1]=1;f[1]=1;register int i,j;
	for(i=2;i<=MN-7;i++)
	{
		if(!mark[i]==1){prime[++cnt]=i;f[i]=1-i+mod;}
		for(j=1;j<=cnt&&i*prime[j]<=MN-7;j++)
		{
			mark[i*prime[j]]=1;
			if(i%prime[j]==0){f[i*prime[j]]=f[i];break;}
			else f[i*prime[j]]=f[i]*f[prime[j]]%mod;
		}
	}
	for(i=1;i<=MN-7;i++) f[i]=f[i]*i%mod;
	for(i=2;i<=MN-7;i++) (f[i]+=f[i-1])%=mod; 
}
inline ll sum(int x){return (1ll*x*(x+1)/2)%mod;}
int main()
{
	register int i,r,n,m,M;
	n=read(),m=read();M=min(n,m);
	init();
	for(i=1;i<=M;i=r+1){
		r=min(n/(n/i),m/(m/i));
		(ans+=(sum(n/i)*sum(m/i)%mod*(mod+f[r]-f[i-1])%mod)%mod)%=mod;
	}
	printf("%lld
",ans);
	return 0;
}

【luogu 3327】[SDOI2005] 约数个数和

Problem Here

预备知识

• [d(nm)=sum_{i|n}^{} sum_{j|m} [gcd(i,j)==1] ]

  prove:

  考虑分别理解左式和右式：

  1. 取(nm)的约数个数和相当于各个质因子的（次数(+1)）的累乘
  2. 试想该如何分配(i)，(j)中包含质因子(p)的次数？显然，只有（次数(+1)）种做法，累乘即为答案

• [sum _{i=1}^{n} frac{n}{i}= sum_{j=1}^{n} d(j)$$ ，从右往左理解要方便的多 具体可以参见 [约束研究](https://www.luogu.org/problemnew/show/P1403) ]

最终公式:

令 $$g(n)=sum _{i=1}^{n}frac{n}{i}$$

[ans=sum_{d=1}^{min(n,m)} mu(d) g( frac{n}{d}) g ( frac{m}{d}) ]

#include<bits/stdc++.h>
#define ll long long
using namespace std;
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return x*f;
}
#define MN 50005
int mu[MN],cnt,prime[MN];
ll d[MN],p[MN];
bool mark[MN];
inline void init()
{
	mark[1]=true;
	p[1]=d[1]=mu[1]=1;
	register int i,j;
	for(i=2;i<=MN-5;++i)
	{
		if(!mark[i]){prime[++cnt]=i;mu[i]=-1;d[i]=2;p[i]=2;}
		for(j=1;j<=cnt&&prime[j]*i<=MN-5;++j)
		{
			mark[i*prime[j]]=true;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;d[i*prime[j]]=d[i]/p[i]*(p[i]+1);
				p[i*prime[j]]=p[i]+1;break; 
			}
			else
			{
				mu[i*prime[j]]=-mu[i];d[i*prime[j]]=d[i]*2;
				p[i*prime[j]]=2;
			}
		}
	}
	for(i=2;i<=MN-5;i++) mu[i]+=mu[i-1],d[i]+=d[i-1];
} 
int main()
{
	register int i,r,n,m,T;
	T=read(),init();
	while(T--)
	{
		register ll ans=0ll;
		n=read(),m=read();
		for(i=1;i<=n&&i<=m;i=r+1)
		{
			r=min(n/(n/i),m/(m/i));
			ans+=1ll*d[n/i]*d[m/i]*(mu[r]-mu[i-1]);
		}
		printf("%lld
",ans);
	}
	return 0;
}

【SDOI2014】数表

Problem Here

预备知识

莫比乌斯反演+树状数组！

• “同时整除i和j的所有自然数之和”其实可以理解为(sum_{d|gcd(i,j)}^{} d),也就是(gcd(i,j))的约数和

  而每个数的约数和是可以在O(n log n) 的时间内算出的

• 记F[i]为i的约数和

• 根据莫比乌斯反演的基本操作，我们要求的是$$sum_{d=1}^{n} left lfloor frac{n}{d} ight floorleft lfloor frac{m}{d} ight floorsum_{x|d}^{} F(x)mu (frac{d}{x})$$

• 因为加上了a的限制，先对询问按a排序，采用树状数组维护前缀和

#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
using namespace std;
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return x*f;
}
#define MN 100005
#define MQ 20005
struct ques
{
	int n,m,a,id;
	bool operator < (const ques&o) const{return a<o.a;}
}q[MQ];
struct pac
{
	int a,id;
	bool operator < (const pac&o) const{return a<o.a;}
}F[MN];
//F[x]:x的约数和 
//mu[x]：莫比乌斯函数 
int mu[MN],prime[MN],cnt;
bool mark[MN];
inline void init()
{
	mark[1]=1;mu[1]=1;
	register int i,j;
	for(i=1;i<=MN-5;++i) for(j=1;j*i<=MN-5;++j) F[i*j].a+=i;
	for(i=1;i<=MN-5;++i) F[i].id=i; 
	std::sort(F+1,F+MN-4);
	for(i=2;i<=MN-5;++i)
	{
		if(!mark[i]){prime[++cnt]=i;mu[i]=-1;}
		for(j=1;j<=cnt&&prime[j]*i<=MN-5;++j)
		{
			#define k i*prime[j] 
			mark[k]=1;
			if(i%prime[j]==0){mu[k]=0;break;}
			else mu[k]=-mu[i];
			#undef k
		}
	}
}
#define lowbit(x) (x&(-x))
int t[MN];
inline void C(int x,int val){for(;x<=MN-5;x+=lowbit(x)) t[x]+=val;}
inline int G(int x){int res=0;for(;x;x-=lowbit(x)) res+=t[x];return res;}
int ans[MQ];
int main()
{
	register int i,r,Q,now=1;
	Q=read();init();
	for(i=1;i<=Q;++i) q[i].n=read(),q[i].m=read(),q[i].a=read(),q[i].id=i;
	std::sort(q+1,q+Q+1);
	for(int t=1;t<=Q;++t)
	{
		#define N q[t].n
		#define M q[t].m
		#define A q[t].a
		for(;F[now].a<=A&&now<=MN-5;++now)
			for(i=1;i*F[now].id<=MN-5;++i)
				C(i*F[now].id,F[now].a*mu[i]);
		for(i=1;i<=N&&i<=M;i=r+1)
		{
			r=min(N/(N/i),M/(M/i));
			ans[q[t].id]+=(M/i)*(N/i)*(G(r)-G(i-1));
		}
		#undef N
		#undef M
		#undef A 
	}
	for(i=1;i<=Q;++i) printf("%d
",ans[i]&0x7fffffff);
	return 0; 
}

【SDOI2017】数字表格

Problem Here

预备知识

其实没有预备知识der~，都是套路

最终公式:

[ans=prod_{T=1}^{n} left ( prod_{d|T}^{ } f[d]^{mu (frac{n}{T})} ight )^{left [ frac{n}{T} ight ]left [ frac{m}{T} ight ]} ]

还有就是要小心别T了......

#include<bits/stdc++.h>
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return x*f;
}
#define MN 1000005
#define mod 1000000007
#define swap(x,y) (x^=y^=x^=y)
int mu[MN],prime[MN],cnt/*,fib[MN]*/,f[MN],invf[MN];
bool mark[MN];
inline int fpow(int x,int m)
{
	int res=1;
	while(m)
	{
		if(m&1) res=1ll*res*x%mod;
		x=1ll*x*x%mod;m>>=1;
	}
	return res;
}
inline void init()
{
	register int i,j,k;
	//fib[1]=fib[2]=1;
	//for(i=3;i<=MN-5;++i) fib[i]=fib[i-1]+fib[i-2],fib[i]>mod?fib[i]-=mod:0;
	mark[1]=1;mu[1]=1;
	for(i=2;i<=MN-5;++i)
	{
		if(!mark[i]){prime[++cnt]=i;mu[i]=-1;}
		for(j=1;j<=cnt&&prime[j]*i<=MN-5;++j)
		{
			#define K i*prime[j] 
			mark[K]=1;
			if(i%prime[j]==0){mu[K]=0;break;}
			else mu[K]=-mu[i];
			#undef K
		}
	}
	for(i=1;i<=MN-5;++i) f[i]=invf[i]=1;
    int A=1,B=0;
    for(i=1;i<=MN-5;++i){
        B=(A+B)%mod;A=(B-A+mod)%mod;
        int G[3]={fpow(B,mod-2),1,B};
        for(j=i,k=1;j<=MN-5;j+=i,++k){
            f[j]=(ll)f[j]*G[mu[k]+1]%mod,
            invf[j]=(ll)invf[j]*G[1-mu[k]]%mod;
        }
    }
    f[0]=invf[0]=1;
    for(i=1;i<=MN-5;++i)
        f[i]=(ll)f[i-1]*f[i]%mod,
        invf[i]=(ll)invf[i-1]*invf[i]%mod; 
}
int main()
{
	register int t=read(),i,r,n,m;init();
	while(t--)
	{
		int ans=1;
		n=read(),m=read();if(n>m) swap(n,m);
		for(i=1;i<=n;i=r+1)
		{
			r=min(n/(n/i),m/(m/i));
			#define Fall 1ll*f[r]*invf[i-1]%mod 
            //注意!幂应该对mod-1取模 
			ans=1ll*ans*fpow(Fall,1ll*(n/i)*(m/i)%(mod-1))%mod;
			#undef Fall
		}
		printf("%d
",ans);
	}
	return 0;
}

---

Blog来自PaperCloud，未经允许，请勿转载，TKS！