[学习笔记]莫比乌斯反演

1、(Crash)的数字表格

(assume n<m)

(sum_{i=1}^{n}sum_{j=1}^{m}lcm(i,j))

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

(sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]ij)

(sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}sum_{x|gcd(i,j)}mu(x)ij)

(sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}sum_{x|i,x|j}mu(x)ij)

(sum_{x=1}^{n}sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}mu(x)ij[x|i,x|j])

(assume i=ix,j=jx)

(sum_{x=1}^{n}mu(x)x^{2}sum_{d=1}^{n}dsum_{i=1}^{frac{n}{dx}}sum_{j=1}^{frac{m}{dx}}ij)

(g(x)=sum_{i=1}^{frac{n}{dx}}sum_{j=1}^{frac{n}{dx}}ij=[frac{frac{n}{dx}(frac{n}{dx}+1)}{2}]^{2})

(sum_{d=1}^{n}dsum_{x=1}^{n/d}mu(x)x^{2}g(x))

用数论分块优化(O(nsqrt{n}))

记得开(O_{2})

(Code Below:)

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=10000000+10;
const int p=20101009;
ll n,m,prim[maxn],vis[maxn],mu[maxn],sum[maxn],cnt,ans;
ll a,b,d;

void getmu(ll n){
	mu[1]=1;
	register ll i,j;
	for(i=2;i<=n;i++){
		if(!vis[i]){prim[++cnt]=i;mu[i]=-1;}
		for(j=1;i*prim[j]<=n&&j<=cnt;j++){
			vis[i*prim[j]]=1;
			if(i%prim[j]==0) break;
			mu[i*prim[j]]=-mu[i];
		}
	}
	for(i=1;i<=n;i++) sum[i]=(sum[i-1]+mu[i]*i*i%p)%p;
}

ll g(ll x){
	return ((a/x)*(a/x+1)/2%p)*((b/x)*(b/x+1)/2%p)%p;
}

int main()
{
	scanf("%lld%lld",&n,&m);
	if(n>m) swap(n,m);
	getmu(n);
	register ll l,r,num;
	for(d=1;d<=n;d++){
		num=0;a=n/d;b=m/d;
		for(l=1;l<=a;l=r+1){
			r=min(a/(a/l),b/(b/l));
			num=(num+(sum[r]-sum[l-1])*g(l)%p)%p;
		}
		ans=(ans+d*num%p)%p;
	}
	printf("%lld
",(ans+p)%p);
    return 0;
}

2、简单的数学题

(sum_{i=1}^{n}sum_{j=1}^{n}ijgcd(i,j))

(sum_{d=1}^{n}dsum_{i=1}^{n}sum_{j=1}^{n}[gcd(i,j)==d]ij)

(sum_{d=1}^{n}d^{3}sum_{i=1}^{n/d}sum_{j=1}^{n/d}[gcd(i,j)==1]ij)

(sum_{d=1}^{n}d^{3}sum_{i=1}^{n/d}sum_{j=1}^{n/d}sum_{x|gcd(i,j)}mu(x)ij)

(sum_{d=1}^{n}d^{3}sum_{i=1}^{n/d}sum_{j=1}^{n/d}sum_{x=1}^{n/d}mu(x)[x|i,x|j]ij)

(assume i=ix,j=jx)

(sum_{d=1}^{n}d^{3}sum_{x=1}^{n/d}mu(x)x^{2}sum_{i=1}^{frac{n}{dx}}sum_{j=1}^{frac{n}{dx}}ij)

(g(x)=sum_{i=1}^{frac{n}{dx}}sum_{j=1}^{frac{n}{dx}}ij=[frac{frac{n}{dx}(frac{n}{dx}+1)}{2}]^{2})

(sum_{d=1}^{n}d^{3}sum_{x=1}^{n/d}mu(x)x^{2}g(x))

突然发现时间复杂度(O(nsqrt{n}))，解法是伪的

所以(mobius)函数没用了，要用(phi)函数

(sum_{i=1}^{n}sum_{j=1}^{n}ijgcd(i,j))

(sum_{i=1}^{n}sum_{j=1}^{n}ijsum_{k|i,k|j}varphi(k))

(sum_{k=1}^{n}varphi(k)sum_{k|i}sum_{k|j}ij)

(sum_{k=1}^{n}varphi(k)sum_{i=1}^{n}sum_{j=1}^{n}[k|i,k|j]ij)

(sum_{k=1}^{n}varphi(k)k^{2}sum_{i=1}^{n/k}sum_{j=1}^{n/k}ij)

(g(x)=(sum_{i=1}^{x}i)^{2}=[frac{x(x+1)}{2}]^{2})

(sum_{k=1}^{n}varphi(k)k^{2}g(n/k)^{2})

#include <bits/stdc++.h>
#define ll long long
#define res register ll
using namespace std;
const int maxn=7500000+10;
ll p,n,prim[maxn],vis[maxn],phi[maxn],cnt,ans,inv2,inv6;
map<ll,ll> Phi;

ll fast_pow(ll a,ll b){
	ll ret=1;
	for(;b;b>>=1,a=a*a%p)
		ret=ret*(b&1?a:1)%p;
	return ret;
}

void getphi(ll n){
	phi[1]=1;
	res i,j;
	for(i=2;i<=n;i++){
		if(!vis[i]){prim[++cnt]=i;phi[i]=i-1;}
		for(j=1;i*prim[j]<=n&&j<=cnt;j++){
			vis[i*prim[j]]=1;
			if(i%prim[j]==0){
				phi[i*prim[j]]=phi[i]*prim[j]%p;
				break;
			}
			phi[i*prim[j]]=phi[i]*(prim[j]-1)%p;
		}
	}
	for(i=1;i<=n;i++) phi[i]=(i*i%p*phi[i]%p+phi[i-1])%p;
}

ll g(ll x){
	x%=p;
	return (x*(x+1)/2%p)*(x*(x+1)/2%p)%p;
}

ll f(ll x){
	x%=p;
	return x*(x+1)%p*(2*x+1)%p*inv6%p;
}


ll calc_phi(res n){
	if(n<=7500000) return phi[n];
	if(Phi[n]) return Phi[n];
	res ans=g(n);
	for(res l=2,r;l<=n;l=r+1){
		r=n/(n/l);
		ans=(ans-(f(r)-f(l-1)+p)%p*calc_phi(n/l)%p)%p;
	}
	return Phi[n]=(ans+p)%p;
}

int main()
{
	scanf("%lld%lld",&p,&n);
	getphi(7500000);
	inv2=fast_pow(2,p-2);
	inv6=fast_pow(6,p-2);
	for(res l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		ans=(ans+(calc_phi(r)-calc_phi(l-1)+p)%p*g(n/l)%p)%p;
	}
	printf("%lld
",(ans+p)%p);
	return 0;
}

3、于神之怒加强版

(assume n<m)

(sum_{i=1}^{n}sum_{j=1}^{m}gcd(i,j)^{k})

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

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

(sum_{d=1}^{n}d^{k}sum_{i=1}^{n/d}sum_{j=1}^{m/d}sum_{x|gcd(i,j)}mu(x))

(sum_{d=1}^{n}d^{k}sum_{x=1}^{n}mu(x)sum_{i=1}^{n/d}sum_{j=1}^{m/d}[x|i,x|j])

(sum_{d=1}^{n}d^{k}sum_{x=1}^{n/d}mu(x)lfloorfrac{n}{dx} floorlfloorfrac{m}{dx} floor)

(sum_{t=1}^{n}lfloorfrac{n}{t} floorlfloorfrac{m}{t} floorsum_{d|t}d^{k}mu(frac{t}{d}))

(f(x)=sum_{d|t}d^{k}mu(frac{t}{d}))

(sum_{t=1}^{n}f(t)lfloorfrac{n}{t} floorlfloorfrac{m}{t} floor)

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=5000000+10;
const int p=1e9+7;
ll n,m,k,prim[maxn],vis[maxn],f[maxn],cnt,ans;

ll fast_pow(ll a,ll b){
	ll ret=1;
	for(;b;b>>=1,a=a*a%p)
		ret=ret*(b&1?a:1)%p;
	return ret;
}

void sieve(ll n){
	f[1]=1;
	ll i,j;
	for(i=2;i<=n;i++){
		if(!vis[i]){prim[++cnt]=i;f[i]=fast_pow(i,k)-1;}
		for(j=1;i*prim[j]<=n&&j<=cnt;j++){
			vis[i*prim[j]]=1;
			if(i%prim[j]==0){
				f[i*prim[j]]=f[i]*(f[prim[j]]+1)%p;
				break;
			}
			f[i*prim[j]]=f[i]*f[prim[j]]%p;
		}
	}
	for(i=1;i<=n;i++) f[i]=(f[i]+f[i-1])%p;
}

int main()
{
	ll T,l,r;
	scanf("%lld%lld",&T,&k);
	sieve(5000000);
	while(T--){
		scanf("%lld%lld",&n,&m);
		if(n>m) swap(n,m);
		ans=0;
		for(l=1,r;l<=n;l=r+1){
			r=min(n/(n/l),m/(m/l));
			ans=(ans+(n/l)*(m/l)%p*(f[r]-f[l-1]+p)%p)%p;
		}
		printf("%lld
",ans);
	}
	return 0;
}