刷题笔记：莫比乌斯反演

入门紫题：yy的GCD

  一句话题意：
  求：

\[\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j) \in Prime] \]

  首先强制令\(n \le m\)。
  考虑枚举\(gcd\)：

\[\sum_{k=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=k]*[k \in Prime] \]

\[= \]

\[\sum_{k=1}^{n} \sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{k} \rfloor} [gcd(i,j)=1]*[k \in Prime] \]

  可以知道莫比乌斯函数有一个性质：

\[\sum_{d|n} \mu(d)=[n=1] \]

  所以原式等价于：

\[\sum_{k=1}^{n} \sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{k} \rfloor}\sum_{d|gcd(i,j)}\mu(d) *[k \in Prime] \]

  \(d\)|\(gcd(i,j)\),那么，\(i,j\)都是\(d\)的倍数，考虑枚举\(d\)，则原式等价于：

\[\sum_{k=1}^{n}\sum_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d)*{\lfloor \frac{n}{kd} \rfloor}*{\lfloor \frac{m}{kd} \rfloor}*[k \in Prime] \]

  这里是考虑\(d\)会给范围内所有的倍数贡献，所以计算出倍数的个数即可将循环省去。
  这里的复杂度会炸掉。
  有一个套路，就是考虑变换枚举的对象。
  另\(t=kd\)，枚举\(t\)，易知\(t \in [1,n]\)。
  那么原式等价于：

\[\sum_{t=1}^{n} {\lfloor \frac{n}{t} \rfloor}*{\lfloor \frac{m}{t} \rfloor} \sum_{x|t,x \in Prime} \mu(\frac{t}{x}) \]

  后者可以处理前缀和，前者可以除法分块。
  总复杂度：\(O(q\sqrt{n})\).

Code

#include<bits/stdc++.h>
using namespace std;
namespace STD
{
	#define rr register
	using ll=long long ;
	const int N=1e7+4;
	int t,n,m;
	inline void swap(int & x,int & y){x^=y,y^=x,x^=y;}
	int read()
	{
		rr int x_read=0,y_read=1;
		rr char c_read=getchar();
		while(c_read<'0'||c_read>'9')
		{
			if(c_read=='-') y_read=-1;
			c_read=getchar();
		}
		while(c_read<='9'&&c_read>='0')
		{
			x_read=(x_read<<3)+(x_read<<1)+(c_read^48);
			c_read=getchar();
		}
		return x_read*y_read;
	}
	int cnt;
	int mu[N],f[N];
	ll sum[N];
	int prime[N/10];
	bool is_not_prime[N];
	void pre()
	{
		mu[1]=1;
		for(int i=2;i<=N-4;i++)
		{
			if(!is_not_prime[i]) prime[++cnt]=i,mu[i]=-1;
			for(int j=1;j<=cnt&&i*prime[j]<=N-4;j++)
			{
				is_not_prime[i*prime[j]]=1;
				if(i%prime[j]) mu[i*prime[j]]=-mu[i];
				else
				{
					mu[i*prime[j]]=0;
					break;
				}
			}
		}
		for(int i=1;i<=cnt;i++)
			for(int j=1;j*prime[i]<=N-4;j++)
				f[j*prime[i]]+=mu[j];
		for(int i=1;i<=N-4;i++)
			sum[i]=sum[i-1]+f[i];
	}
};
using namespace STD;
int main()
{
	pre();
	t=read();
	while(t--)
	{
		n=read(),m=read();
		if(n>m) swap(n,m);
		ll ans=0;
		int l=1,r=0;
		while(r<n)
		{
			r=min(n/(n/l),m/(m/l));
			ans+=((ll)n/l)*((ll)m/l)*(ll)(sum[r]-sum[l-1]);
			l=r+1;
		}
		printf("%lld\n",ans);
	}
}

SDOI2014： 数表

  离线询问+莫比乌斯反演+\(BIT\);
  一句话题意：
  求：

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

  强制令\(n \le m\)
  首先记\(f(x)\)表示\(x\)所有的因数的和。
  那么原式就是：

\[\sum_{i=1}^{n}\sum_{j=1}^{m} f(gcd(i,j)) \]

\[= \]

\[\sum_{k=1}^{n}\sum_{i=1}^{\lfloor \frac{n}{k} \rfloor }\sum_{j=1}^{\lfloor \frac{m}{k} \rfloor } [gcd(i,j)=1]*f(k)*[f(k) \le a] \]

\[= \]

\[\sum_{k=1}^{n}\sum_{i=1}^{\lfloor \frac{n}{k} \rfloor }\sum_{j=1}^{\lfloor \frac{m}{k} \rfloor } \sum_{d|gcd(i,j)} \mu(d) * f(k) * [f(k) \le a] \]

  和上提一样的套路可以化为：

\[\sum_{t=1}^{n} {\lfloor \frac{n}{t} \rfloor}*{\lfloor \frac{m}{t} \rfloor} \sum_{d|t} {\frac{t}{d}} f(d) * [f(d) \le a] \]

  可以考虑一件事就是更小的\(a\)的值是可以给更大的\(a\)用的，将询问离线，然后按\(a\)升序排序。
  每次将新加入的值加给它所有比\(n_{max}\)小的倍数，用树状数组维护前缀和，然后除法分块即可。
  由于每个数只会被加一次，每次枚举倍数，所以这一步复杂度为\(N*ln N*log N\)
  每次回答的复杂度为\(\sqrt{N}*logN\)。
  所以总复杂度为\(O(N*ln N*log N)\)

Code

#include<bits/stdc++.h>
using namespace std;
namespace STD
{
	#define rr register
	using ll=long long;
	const ll mod=1ll<<31;
	const int N=1e5;
	const int Q=2e4+4;
	int q,last;
	ll ans[Q];
	void my_swap(int & x,int & y){x^=y,y^=x,x^=y;}
	struct req
	{
		int n,m,a,id;
		req():n(0),m(0),a(0),id(0) {}
		bool operator<(const req x) const {return a<x.a;}
	}a[Q];
	template<typename type>
	inline type cmin(type x,type y){return x<y?x:y;}
	int read()
	{
		rr int x_read=0,y_read=1;
		rr char c_read=getchar();
		while(c_read<'0'||c_read>'9')
		{
			if(c_read=='-') y_read=-1;
			c_read=getchar();
		}
		while(c_read<='9'&&c_read>='0')
		{
			x_read=(x_read<<3)+(x_read<<1)+(c_read^48);
			c_read=getchar();
		}
		return x_read*y_read;
	}
	class BIT
	{
		private:
			ll t[N+4];
			inline int lowbit(int x){return x&(-x);}
		public:
			void insert(int,int);
			ll query(int,int);
	}t;
	void BIT::insert(int pos,int val)
	{
		while(pos<=N)
		{
			t[pos]=t[pos]+val;
			if(t[pos]>=mod) t[pos]-=mod;
			pos+=lowbit(pos);
		}
	}
	ll BIT::query(int l,int r)
	{
		ll sum=0;
		while(r)
		{
			sum=sum+t[r];
			if(sum>=mod) sum-=mod;
			r-=lowbit(r);
		}
		while(l)
		{
			sum=sum-t[l];
			if(sum<0) sum+=mod;
			l-=lowbit(l);
		}
		return sum;
	}
	struct pack
	{
		ll f;
		int num;
		pack():f(0ll),num(0) {}
		bool operator<(const pack x) const {return f<x.f;}
	}b[N+4];
	int cnt;
	int prime[(N+4)/10],mu[N+4];
	bool tag[N+4];
	void pre()
	{
		mu[1]=1;
		for(int i=2;i<=N;i++)
		{
			if(!tag[i]) prime[++cnt]=i,mu[i]=-1;
			for(int j=1;j<=cnt&&prime[j]*i<=N;j++)
			{
				tag[i*prime[j]]=1;
				if(i%prime[j]) mu[i*prime[j]]=-mu[i];
				else
				{
					mu[i*prime[j]]=0;
					break;
				}
			}
		}
		for(int i=1;i<=N;i++)
		{
			b[i].num=i;
			for(int j=1;(ll)(j*i)<=(ll)N;j++)
			{
				b[i*j].f+=i;
				if(b[i*j].f>=mod)
					b[i*j].f-=mod;
			}
		}
		sort(b+1,b+1+N);
	}
	int find(ll x)
	{
		int l=1,r=N;
		while(l<r)
		{
			int mid=(l+r+1)>>1;
			if(b[mid].f<=x) l=mid;
			else r=mid-1;
		}
		return l;
	}
};
using namespace STD;
int main()
{
	q=read();
	pre();
	for(int i=1;i<=q;i++)
	{
		a[i].n=read(),a[i].m=read();
		a[i].a=read();
		a[i].id=i;
		if(a[i].n>a[i].m) my_swap(a[i].n,a[i].m);
	}
	sort(a+1,a+1+q);
	for(int i=1;i<=q;i++)
	{
		if(a[i].a<=0) continue;
		int pos=find(a[i].a);
		for(int j=last+1;j<=pos;j++)
			for(int k=1;(ll)(k*b[j].num)<=(ll)N;k++)
				t.insert(k*b[j].num,mu[k]*b[j].f);		
		int l=1,r=0;
		while(r<a[i].n)
		{
			r=cmin(a[i].n/(a[i].n/l),a[i].m/(a[i].m/l));
			ans[a[i].id]=(ans[a[i].id]+((a[i].n/l)*(a[i].m/l)%mod)*t.query(l-1,r)%mod)%mod;
			l=r+1;
		}	
		last=pos;
	}
	for(int i=1;i<=q;i++) printf("%lld\n",ans[i]);
}

DZY Loves Math

  线性筛的拓展使用+莫比乌斯反演。
  首先考虑如何快速地求出题中给出的\(f\)。
  传统做法比较无脑，就是直接\(O(n \sqrt n)\)枚举所有约数，边除边加加
  有一个\(O(n)\)做法，就是利用线性筛，记录两个数组\(f\)与\(f'\)。\(f\)就是要求的数组，\(f'\)就是自己最小的质因子的幂数。
  转移就有：

f[iprime[j]]=max(f'[i]+1,f[i]),f'[iprime[j]]=f'[i]+1;

  其中\(prime[j]\)是线性筛枚举的素数。
  最后的公式就是：

\[\sum_{t=1}^{n}\sum_{d|t} \mu(d)*f[\frac{t}{d}]*{\lfloor \frac{n}{t} \rfloor}{\lfloor \frac{m}{t} \rfloor } \]

  \(m,n\)老规矩。
  推导略，因为太水了。

Code

#include<bits/stdc++.h>
using namespace std;
namespace STD
{
	#define rr register
	using ll=long long;
	const int N=1e7;
	int t;
	int f[N+4],f_[N+4];
	ll g[N+4];
	template<typename type>
	inline type cmax(type x,type y){return x>y?x:y;}
	template<typename type>
	inline type cmin(type x,type y){return x<y?x:y;}
	inline void my_swap(int & x,int & y){x^=y,y^=x,x^=y;}
	inline int read()
	{
		rr int x_read=0,y_read=1;
		rr char c_read=getchar();
		while(c_read<'0'||c_read>'9')
		{
			if(c_read=='-') y_read=-1;
			c_read=getchar();
		}
		while(c_read<='9'&&c_read>='0')
		{
			x_read=(x_read<<3)+(x_read<<1)+(c_read^48);
			c_read=getchar();
		}
		return x_read*y_read;
	}
	int cnt;
	int prime[N+4],mu[N+4];
	bool tag[N+4];
	void pre()
	{
		mu[1]=1;
		int i,j,en;
		for(i=2;i<=N;i++)
		{
			if(!tag[i]) prime[++cnt]=i,mu[i]=-1,f[i]=f_[i]=1;
			for(j=1;j<=cnt&&prime[j]*i<=N;j++)
			{
				tag[i*prime[j]]=1;
				if(i%prime[j]) mu[i*prime[j]]=-mu[i],f[i*prime[j]]=f[i],f_[i*prime[j]]=1;
				else
				{
					f[i*prime[j]]=cmax(f[i],f_[i]+1);
					f_[i*prime[j]]=f_[i]+1;
					break;
				}
			}
		}	
		for(i=2;i<=N;i++)
		{
			en=cmin(N/i,i);
			for(j=1;j<=en;j++)
				g[i*j]+=mu[j]*f[i]+mu[i]*f[j];
			if(i==en&&i*i<=N)
				g[i*i]-=mu[i]*f[i];
			g[i]+=g[i-1];
		}
	}
};
using namespace STD;
int main()
{
	t=read();
	pre();
	int a,b;
	while(t--)
	{
		a=read(),b=read();
		if(a>b) my_swap(a,b);
		int l=1,r=0;
		ll ans=0ll;
		while(r<a)
		{
			r=cmin(a/(a/l),b/(b/l));
			ans+=1ll*(a/l)*(b/l)*(g[r]-g[l-1]);
			l=r+1;
		}
		printf("%lld\n",ans);
	}
}

约数个数和

  首先要知道\(d(ij)\)的性质：

\[d(ij)=\sum_{x|i} \sum_{y|j} [gcd(x,y)=1] \]

  证明见这里。

Code

#include<bits/stdc++.h>
using namespace std;
namespace STD
{
	#define rr register
	typedef long long ll;
	const int N=50000;
	int t;
	int read()
	{
		rr int x_read=0,y_read=1;
		rr char c_read=getchar();
		while(c_read<'0'||c_read>'9')
		{
			if(c_read=='-') y_read=-1;
			c_read=getchar();
		}
		while(c_read<='9'&&c_read>='0')
		{
			x_read=(x_read<<3)+(x_read<<1)+(c_read^48);
			c_read=getchar();
		}
		return x_read*y_read;
	}
	int cnt;
	int mu[N+4],prime[N+4];
	bool tag[N+1];ll f[N+1];
	void pre()
	{
		mu[1]=1;
		for(int i=2;i<=N;i++)
		{
			if(!tag[i]) prime[++cnt]=i,mu[i]=-1;
			for(int j=1;j<=cnt&&prime[j]*i<=N;j++)
			{
				tag[i*prime[j]]=1;
				if(i%prime[j]) mu[prime[j]*i]=-mu[i];
				else
				{
					mu[prime[j]*i]=0;
					break;
				}
			}
			mu[i]+=mu[i-1];
		}	
		for(int i=1;i<=N;i++)
		{
			int l=1,r=0;
			while(r<i)
			{
				r=i/(i/l);
				f[i]+=(r-l+1)*(i/l);
				l=r+1;
			}
		}
	}
	inline void my_swap(int & x,int & y){x^=y,y^=x,x^=y;}
	template<typename type>
	inline type cmin(type x,type y){return x<y?x:y;}
};
using namespace STD;
int main()
{
	pre();
	t=read();
	while(t--)
	{
		int n=read(),m=read();
		ll ans=0;
		if(n>m) my_swap(m,n);
		int l=1,r=0;
		while(r<n)
		{
			r=cmin(m/(m/l),n/(n/l));
			ans+=(mu[r]-mu[l-1])*f[m/l]*f[n/l];
			l=r+1;
		}
		printf("%lld\n",ans);
	}
}

数字表格

  人生第一道黑题，虽然说是黑里发紫。
  首先\(O(n)\)搞出\(fib\)。
  然后颓式子：

\[\prod_{i=1}^{n} \prod_{j=1}^{m} f[gcd(i,j) \]

\[= \]

\[\prod_{k=1}^{n} \prod_{i=1}^{n} \prod_{j=1}^{m} f[k]^{[gcd(i,j)=k]} \]

\[= \]

\[\prod_{k=1}^{n} \prod_{i=1}^{\lfloor \frac{n}{k} \rfloor} \prod_{j=1}^{\lfloor \frac{m}{k} \rfloor} f[k]^{[gcd(i,j)=1]} \]

\[= \]

\[\prod_{k=1}^{n} \prod_{i=1}^{\lfloor \frac{n}{k} \rfloor} \prod_{j=1}^{\lfloor \frac{m}{k} \rfloor} f[k]^{\sum_{d|gcd(i,j)} \mu(d)} \]

\[= \]

\[\prod_{k=1}^{n} \prod_{i=1}^{\lfloor \frac{n}{k} \rfloor} \prod_{j=1}^{\lfloor \frac{m}{k} \rfloor} \prod_{d|gcd(i,j)} f[k]^{\mu(d)} \]

\[= \]

\[\prod_{k=1}^{n} \prod_{d=1}^{\lfloor \frac{n}{k} \rfloor} (f[k]^{\mu(d)})^{{\lfloor \frac{n}{kd} \rfloor}{\lfloor \frac{m}{kd} \rfloor}} \]

\[= \]

\[\prod_{t=1}^{n} \prod_{d|t} (f[d]^{\mu(\frac{t}{d})})^{{\lfloor \frac{n}{t} \rfloor }{\lfloor \frac{m}{t} \rfloor}} \]

  处理一下前缀积即可。

Code

#include<bits/stdc++.h>
using namespace std;
namespace STD
{
	#define rr register
	using ll=long long;
	const int mod=1e9+7;
	const int N=1e6;
	int t;
	int read()
	{
		rr int x_read=0,y_read=1;
		rr char c_read=getchar();
		while(c_read<'0'||c_read>'9')
		{
			if(c_read=='-') y_read=-1;
			c_read=getchar();
		}
		while(c_read<='9'&&c_read>='0')
		{
			x_read=(x_read<<3)+(x_read<<1)+(c_read^48);
			c_read=getchar();
		}
		return x_read*y_read;
	}
	ll qpow(ll base,ll exp)
	{
		if(exp<0) return qpow(base,mod-2);
		ll ret=1ll;
		while(exp)
		{
			if(exp&1) ret=ret*base%mod;
			base=base*base%mod;
			exp>>=1;
		}
		return ret;
	}
	int cnt;
	int prime[N+1],mu[N+1];
	ll f[N+1],g[N+1];
	bool tag[N+1];
	void pre()
	{
		mu[1]=1;
		for(int i=2;i<=N;i++)
		{
			if(!tag[i]) prime[++cnt]=i,mu[i]=-1;
			for(int j=1;j<=cnt&&prime[j]*i<=N;j++)
			{
				tag[i*prime[j]]=1;
				if(!(i%prime[j])) break;
				mu[i*prime[j]]=-mu[i];
			}
		}
		f[1]=f[2]=1;
		for(int i=3;i<=N;i++)
			f[i]=(f[i-1]+f[i-2])%mod;
		for(int i=0;i<=N;i++) g[i]=1ll;
		for(int i=1;i<=N;i++)
		{
			for(int j=1;j*i<=N;j++)
				g[i*j]=g[i*j]*qpow(f[j],mu[i])%mod;
			g[i]=g[i]*g[i-1]%mod;
		}
	}
	inline void my_swap(ll & x,ll & y){x^=y,y^=x,x^=y;}
	template<typename type>
	inline type cmin(type x,type y){return x<y?x:y;}
};
using namespace STD;
int main()
{
	pre();
	t=read();
	while(t--)
	{
		ll n=read(),m=read();
		if(n>m) my_swap(m,n);
		int l=1,r=0;
		ll ans=1ll,temp;
		while(r<n)
		{
			r=cmin(n/(n/l),m/(m/l));
			temp=g[r]*qpow(g[l-1],mod-2)%mod;
			ans=(ans*qpow(temp,(n/l)*(m/l)))%mod;
			l=r+1;
		}
		printf("%lld\n",ans);
	}
}