莫比乌斯反演

经典柿子

[[gcd(i,j)=1]=sum_{wmid gcd(i,j)}mu(w)=sum_{wmid i,wmid j}mu(w) ]

POJ3904 Sky Code

description

给定(n) 个不超过(10000) 的正整数(a_1,a_2,cdots,a_n) ，求

[sum_{1le i<j<k<lle n}[gcd(a_i,a_j,a_k,a_l)=1] ]

(nle 10000,Tle10000)

solution

大力推柿子

[egin{aligned} ans&=sum_{1le i<j<k<lle n}[gcd(a_i,a_j,a_k,a_l)=1]\ &=sum_{1le i<j<k<lle n}sum_{wmid a_i,wmid a_j,wmid a_k,wmid a_l}mu(w)\ &=sum_{w=1}^{10000}mu(w)sum_{1le i<j<k<lle n}[wmid a_i][wmid a_j][wmid a_k][wmid a_l] end{aligned} ]

我们不妨假设有(c(w)) 个数是(w) 的倍数，那么后面那些东西就等于(dbinom{c(w)}{4}) 。

复杂度(mathcal O(nT)) 。

code

#include<cstdio>
#include<vector>
using namespace std;
const int N=1e4+5;
typedef long long ll;ll C[N];
bool flag[N];int mu[N];
vector<int>pr;
inline void pre()
{
	int n=N-5;
	flag[1]=1;mu[1]=1;
	for(int i=2;i<=n;++i)
	{
		if(!flag[i]){pr.push_back(i);mu[i]=-1;}
		for(int j=0;j<pr.size();++j)
		{
			int num=i*pr[j];if(num>n)break;
			flag[num]=1;
			if(i%pr[j])mu[num]=-mu[i];
			else{mu[num]=0;break;}
		}
	}
	for(int i=4;i<=n;++i)
		C[i]=1ll*i*(i-1)*(i-2)*(i-3)/24;
}
int n,cnt[N];
int main()
{
	pre();
	while(scanf("%d",&n)==1)
	{
		int mx=0;
		for(int i=1,d;i<=n;++i)
		{
			scanf("%d",&d);
			for(int j=1;j*j<=d;++j)
				if(d%j==0)
				{
					int w1=d/j,w2=j;
					++cnt[w1];
					if(w1^w2)++cnt[w2];
				}
			mx=max(mx,d);
		}ll ans=0;
		for(int i=1;i<=mx;++i)
		{
			if(!mu[i])continue;
			ans+=mu[i]*C[cnt[i]];
		}
		printf("%lld
",ans);
		for(int i=1;i<=mx;++i)cnt[i]=0;
	}
	return 0;
}

可见格点

description

三维直角坐标系中三个坐标均非负且不超过(n) 的整点中有多少点可以从原点看到。

(nle 10^6)

solution

原题相当于是求

[sum_{xge 0,yge0,zge0}[gcd(x,y,z)=1] ]

这是经典的莫反柿子。不过需要考虑(x,y,z) 中部分为零的情况，于是就有

[egin{aligned} &=3+3sum_{x>0,y>0}[gcd(x,y)=1]+sum_{x>0,y>0,z>0}[gcd(x,y,z)=1]\ &=3+3sum_{x>0,y>0}sum_{wmid x,wmid y}mu(w)+sum_{x>0,y>0,z>0}sum_{wmid x,wmid y,wmid z}mu(w)\ &=3+3sum_{w=1}^nmu(w)lfloorfrac nw floor^2+sum_{w=1}^nmu(w)lfloorfrac nw floor^3\ &=3+sum_{w=1}^nmu(w)lfloorfrac nw floor^2(lfloorfrac nw floor+3) end{aligned} ]

复杂度(mathcal O(nT)) 。

code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1E+6+5,V=1e+6;
bool flag[N];vector<int>pr;
int mu[N];
inline void pre()
{
	mu[1]=1,flag[1]=1;
	for(int i=2;i<=V;++i)
	{
		if(!flag[i])pr.push_back(i),mu[i]=-1;
		for(int c:pr)
		{
			int num=c*i;if(num>V)break;
			flag[num]=1;
			if(i%c)mu[num]=-mu[i];
			else{mu[num]=0;break;}
		}
	}
}
int main()
{
	pre();
	int T;scanf("%d",&T);
	while(T-->0)
	{
		int n;scanf("%d",&n);
		ll ans=3;
		for(int i=1;i<=n;++i)
		{
			if(!mu[i])continue;
			int k=n/i;
			ans+=1ll*mu[i]*k*k*(k+3);
		}
		printf("%lld
",ans);
	}
	return 0;
}

数一数a x b

description

令

[f(n)=sum_{i=1}^nsum_{j=1}^n[n mid ij] ]

求

[g(n)=sum_{dmid n}f(d) ]

(nle 10^9,Tle 20000)

solution

不整除显然是难以计算的，我们考虑简单容斥，即

[f(n)=n^2-h(n)=n^2-sum_{i=1}^nsum_{j=1}^n[nmid ij] ]

因此就有

[g(n)=sum_{dmid n}d^2-sum_{dmid n}sum_{i=1}^dsum_{j=1}^d[dmid ij] ]

分开考虑。我们假设(n) 的质因数分解为(n=prod_limits{i=1}^kp_i^{a_i}) ，那么就有

[sum_{dmid n}d^2=prod_{i=1}^ksum_{j=0}^{a_i}(p_i^j)^2 ]

因为次数之和大概是(mathcal O(log n)) 级别的，因此这部分复杂度为(mathcal O(log n)) 。

对于后面部分，我们枚举(w=gcd(d,i)) ，则有

[egin{aligned} &=sum_{dmid n}sum_{i=1}^dsum_{j=1}^d[dmid ij]\ &=sum_{w=1}^nsum_{dmid n,wmid d}sum_{wmid i}^dsum_{j=1}^d[frac dwmid frac iwcdot j][gcd(frac dw,frac iw)=1]\ &=sum_{wmid n}sum_{dmid frac nw}sum_{j=1}^{wd}[dmid j]sum_{i=1}^{d}[gcd(d,i)=1]\ &=sum_{wmid n}sum_{dmid frac nw}sum_{j=1}^{wd}[dmid j]varphi(d)\ &=sum_{wmid n}sum_{dmid frac nw}sum_{j=1}^{w}varphi(d)\ &=sum_{wmid n}wsum_{dmid frac nw}varphi(d)\ &=sum_{wmid n}wfrac nw\ &=ncdot sigma_0(n) end{aligned} ]

可以通过预处理质数来优化后续复杂度。至此，本题完结。

总复杂度为(mathcal O(sqrt n+Tlog n)) 。

code

#include<bits/stdc++.h>
using namespace std;
const int N=1e5+5;
typedef unsigned long long ll;
bool flag[N];vector<int>pr;
inline void pre(int n)
{
	flag[1]=1;
	for(int i=2;i<=n;++i)
	{
		if(!flag[i])pr.push_back(i);
		for(int j=0;j<pr.size();++j)
		{
			int num=i*pr[j];if(num>n)break;
			flag[num]=1;
			if(i%pr[j]==0)break;
		}
	}
}
int main()
{
	pre(4e4);int T;scanf("%d",&T);
	while(T-->0)
	{
		int n;scanf("%d",&n);ll g=1,h=n;
		for(int i=0;i<pr.size();++i)
		{
			int p=pr[i];ll ret=0;
			if(n<p)break;if(n%p)continue;
			int k=0;while(n%p==0)n/=p,++k;
			for(int j=0,cur=1;j<=k;++j,cur*=p)ret+=1ull*cur*cur;
			g*=ret;h*=k+1;
		}
		if(n>1)h<<=1,g*=(1+1ull*n*n);
		printf("%llu
",g-h);
	}
	return 0;
}

墨菲斯

description

给定(n,m,p) ，求

[sum_{i=1}^nsum_{j=1}^m[h(gcd(i,j))le p] ]

其中(h(x)) 表示(x) 的唯一分解中质数个数。

solution

首先套路地推柿子，枚举(gcd)

[egin{aligned} ans&=sum_{d=1}^n[h(d)le p]sum_{dmid i}^nsum_{dmid j}^m[gcd(frac id,frac jd)=1]\ &=sum_{d=1}^n[h(d)le p]sum_isum_jsum_{wmid i,wmid j}mu(w)\ &=sum_{d=1}^n[h(d)le p]sum_{w=1}^{lfloor frac nd floor}mu(w)lfloor frac n{dw} floorlfloor frac m{dw} floor\ &=sum_{T=1}^nlfloor frac nT floorlfloor frac mT floorsum_{dmid T}[h(d)le p]mu(frac Td) end{aligned} ]

似乎遇到了障碍。不过我们注意到(2^{19}>5 imes 10^5) ，因此当(pge 18) 时答案就是(nm) 。我们只用考虑(p<18) 的情况。

令

[f(p,d)=sum_{dmid T}[h(d)le p]mu(frac Td) ]

我们可以在(mathcal O(18nlog n)) 的时间内求出所有的(f(p,d)) ，然后对其求前缀和以便在整除分块时可以快速计算区间和。

复杂度(mathcal O(18nlog n+qsqrt n)) 。

code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e5+5,MX=18;
bool flag[N];int mu[N],h[N];ll f[MX+5][N];
vector<int>pr;
inline void pre()
{
	int n=N-5;
	flag[1]=1;mu[1]=1;h[1]=0;
	for(int i=2;i<=n;++i)
	{
		if(!flag[i])pr.push_back(i),mu[i]=-1,h[i]=1;
		for(int j=0;j<pr.size();++j)
		{
			int num=i*pr[j];if(num>n)break;
			flag[num]=1;h[num]=h[i]+1;
			if(i%pr[j])mu[num]=-mu[i];
			else{mu[num]=0;break;}
		}
	}
	for(int p=0;p<=MX;++p)
	{
		for(int d=1;d<=n;++d)
		{
			if(h[d]>p)continue;
			for(int t=d,c=1;t<=n;t+=d,++c)f[p][t]+=mu[c];
		}
		for(int i=1;i<=n;++i)f[p][i]+=f[p][i-1];
	}
}
int main()
{
	pre();int T;scanf("%d",&T);
	while(T-->0)
	{
		int n,m,p;scanf("%d%d%d",&n,&m,&p);
		if(p>MX){printf("%lld
",1ll*n*m);continue;}
		if(n>m)swap(n,m);ll ans=0;
		for(int l=1,r;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			ans+=1ll*(n/l)*(m/l)*(f[p][r]-f[p][l-1]);
		}
		printf("%lld
",ans);
	}
	return 0;
}