莫比乌斯反演练习

线性筛板子

const int N = 1e6+10;
int phi[N], mu[N], p[N], cnt, vis[N];
void init() {
	phi[1] = mu[1] = 1;
	for (int i=2; i<N; ++i) {
		if (!vis[i]) p[++cnt]=i,mu[i]=-1,phi[i]=i-1;
		for (int j=1,t; j<=cnt&&i*p[j]<N; ++j) {
			vis[t=i*p[j]] = 1;
			if (i%p[j]==0) {mu[t]=0,phi[t]=phi[i]*p[j];break;}
			mu[t]=-mu[i],phi[t]=phi[i]*phi[p[j]];
		}
	}
}

$O(nloglogn)$求$g=mu * f$ 板子

void calc(int n, int *f, int *g) {
	REP(i,1,n) g[i] = f[i];
	REP(i,1,cnt) PER(j,1,n/p[i]) g[j*p[i]] -= g[j];
}

整除分块板子, 优化了一下除法次数

int mx = sqrt(n);
REP(i,1,mx) ans += (ll)(s[i]-s[i-1])*f(n/i);
for (int i=mx+1,j,k=n/i; i<=n; i=j+1,--k) {
    j = n/k;
    ans += (ll)(s[j]-s[i-1])*f(k);
}

练习1. luogu P2257 YY的GCD

大意: 求$sumlimits_{i=1}^n sumlimits_{j=1}^m[gcd(i,j)为素数]$

这种反演题目一般都是先转化为某个$g$函数的卷积, 再交换求和次序.

设$f(x)=[x为素数]$,  $g=mu *f=sumlimits_{p|n}mu(frac{n}{p})$, 则有

$egin{align*} sumlimits_{i=1}^n sumlimits_{j=1}^mf(gcd(i,j)) &= sumlimits_{i=1}^n sumlimits_{j=1}^m sumlimits_{d|gcd(i,j)} g(d) \ &= sumlimits_{d=1}^{min(n,m)}g(d)sumlimits_{i=1}^nsumlimits_{j=1}^m [d|gcd(i,j)] \ &= sumlimits_{d=1}^{min(n,m)}g(d) lfloor frac{n}{d} floor lfloor frac{m}{d} floor end{align*}$

所以求出$g$以后整除分块即可, 考虑如何求$g$, $O(nloglogn)$当然很容易求出

这里思考一下$O(n)$的求法

对于$g(x)$, 假设$y$为$x$的最小素因子, 即$x=iy$

若$x$为素数, 显然有$g(x)=mu(1)=1$

若$i%y==0$, $g(x)=mu(i)$

若$i%y!=0$, $g(x)=-g(i)+mu(i)$

欧拉筛一次即可$O(n)-O(qsqrt{n})$

const int N = 1e7+10;
int mu[N], p[N], cnt, vis[N], g[N];
void init() {
	mu[1] = 1;
	for (int i=2; i<N; ++i) {
		if (!vis[i]) p[++cnt]=i,mu[i]=-1,g[i]=1;
		for (int j=1,t; j<=cnt&&i*p[j]<N; ++j) {
			vis[t=i*p[j]] = 1;
			if (i%p[j]==0) {mu[t]=0,g[t]=mu[i];break;}
			mu[t]=-mu[i],g[t]=-g[i]+mu[i];
		}
		g[i]+=g[i-1];
	}
}

int main() {
	init();
	int t;
	scanf("%d", &t);
	REP(i,1,t) {
		int n, m;
		scanf("%d%d", &n, &m);
		ll ans = 0;
		if (n>m) swap(n,m);
		int lim = min(int(sqrt(n)*2),n);
		REP(i,1,lim) ans+=(ll)(g[i]-g[i-1])*(n/i)*(m/i);
		for (int i=lim+1,j; i<=n; i=j+1) {
			j = min(n/(n/i),m/(m/i));
			ans+=(ll)(g[j]-g[i-1])*(n/i)*(m/i);
		}
		printf("%lld
", ans);
	}
}

练习2. Luogu P4213 杜教筛

大意: 求$sumlimits_{i=1}^{n}varphi(i),sumlimits_{i=1}^{n}mu(i)$
杜教筛是用来求一类积性函数前缀和的算法, 核心思想是利用公式
$$g(1)S(n)=sumlimits_{i=1}^n(f*g)(i)-sumlimits_{i=2}^n g(i)S(lfloor frac{n}{i} floor)$$
其中$S(n)$是我们要求的$f(n)$的前缀和, $g(n)$是我们需要构造的一个函数, 需要满足$f*g$和$g$的前缀和可以$O(1)$求出.
可以证明上面这个递归式调用的不同的$S(n)$不超过$2sqrt{n}$个, 即
$$S(1),S(2),...,S(sqrt{n}),S(frac{n}{1}),S(frac{n}{2}),...,S(frac{n}{sqrt{n}})$$
对于$ile n^{frac{2}{3}}$的$S(i)$, 先暴力求出.
而对于$i>n^{frac{2}{3}}$的$S(i)$, 即$S(frac{n}{1}),S(frac{n}{2}),...,S(frac{n}{n^{frac{1}{3}}})$, 只有$n^{frac{1}{3}}$项, 可以每一项$sqrt{n}$求出
所以总复杂度为$O(n^{frac{2}{3}})$

对于本题来说, 对于$sumlimits_{i=1}^nmu(i)$取$g=oldsymbol 1$, 就有$f*g=varepsilon$, 对于$sumlimits_{i=1}^{n}varphi(i)$同理

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <math.h>
#include <set>
#include <map>
#include <queue>
#include <string>
#include <string.h>
#include <bitset>
#define REP(i,a,n) for(int i=a;i<=n;++i)
#define PER(i,a,n) for(int i=n;i>=a;--i)
#define hr putchar(10)
#define pb push_back
#define lc (o<<1)
#define rc (lc|1)
#define mid ((l+r)>>1)
#define ls lc,l,mid
#define rs rc,mid+1,r
#define x first
#define y second
#define io std::ios::sync_with_stdio(false)
#define endl '
'
using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
const int P = 1e9+7, INF = 0x3f3f3f3f;
ll gcd(ll a,ll b) {return b?gcd(b,a%b):a;}
ll qpow(ll a,ll n) {ll r=1%P;for (a%=P;n;a=a*a%P,n>>=1)if(n&1)r=r*a%P;return r;}
ll inv(ll x){return x<=1?1:inv(P%x)*(P-P/x)%P;}
//head


const int N1 = ~0u>>1, N2 = 2e6, N3 = 1300;

int N;
ll phi[N2], phi2[N3];
int mu[N2], mu2[N3];
int vis[N3], v[N2];
int p[N2], cnt;

inline void init(int n) {
    phi[1] = mu[1] = 1;
    REP(i,2,n) {
        if (!v[i]) p[++cnt]=i,mu[i]=-1,phi[i]=i-1;
        for (int j=1,t; j<=cnt&&i*p[j]<=n; ++j) {
            v[t=i*p[j]] = 1;
            if (i%p[j]==0) {mu[t]=0,phi[t]=phi[i]*p[j];break;}
            mu[t]=-mu[i],phi[t]=phi[i]*phi[p[j]];
        }
    }
    REP(i,2,N2-1) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}

inline ll sum_phi(int n) {
    if (n<N2) return phi[n];
    int x = N/n;
    if (vis[x]) return phi2[x];
    vis[x] = 1;
    ll &ans = phi2[x];
    ans = (ll)n*(n+1)/2;
    int mx = sqrt(n);
    REP(i,2,mx) ans -= sum_phi(n/i);
    for (int i=mx+1,j,k=n/i; i<=n; i=j+1,--k) {
        j = n/k;
        ans -= (j-i+1)*sum_phi(k);
    }
    return ans;
}
inline int sum_mu(int n) {
    if (n<N2) return mu[n];
    int x = N/n;
    if (vis[x]) return mu2[x];
    vis[x] = 1;
    int &ans = mu2[x];
    ans = 1;
    int mx = sqrt(n);
    REP(i,2,mx) ans -= sum_mu(n/i);
    for (int i=mx+1,j,k=n/i; i<=n; i=j+1,--k) {
        j = n/k;
        ans -= (j-i+1)*sum_mu(k);
    }
    return ans;
}

int main() {
    init(N2-1);
    int t;
    scanf("%d", &t);
    while (t--) {
        scanf("%d", &N);
        if (!N) {puts("0
0");continue;}
        memset(vis,0,sizeof vis);
        ll ans1 = sum_phi(N);
        memset(vis,0,sizeof vis);
        int ans2 = sum_mu(N);
        printf("%lld %d
", ans1, ans2);
    }
}

练习3. Luogu P3768 简单数学题

大意: 求$sumlimits_{i=1}^nsumlimits_{j=1}^nijgcd(i,j)$

又是个杜教板子题, 简单推一下.

$egin{align*} sumlimits_{i=1}^nsumlimits_{j=1}^nijgcd(i,j) &= sumlimits_{i=1}^nsumlimits_{j=1}^nijsumlimits_{d|gcd(i,j)} varphi(d) \&= sumlimits_{d=1}^nvarphi(d)sumlimits_{i=1}^nsumlimits_{j=1}^nij[d|gcd(i,j)] \ &= sumlimits_{d=1}^nd^2varphi(d)sumlimits_{i=1}^{lfloor frac{n}{d} floor}isumlimits_{j=1}^{lfloor frac{n}{d} floor}jend{align*}$

对这个式子用整除分块, 转化为求$sumlimits_{i=1}^ni^2varphi(i)$

按照杜教的套路, 这里可以凑一个函数$g=Id^2$, 这样$f*g=Id^3$, 总的复杂度仍然是$O(n^{frac{2}{3}})$

const int N2 = 5e6+10, N3 = 2500;
int P, inv6, inv2, cnt;
int phi[N2], vis[N2], p[N2];
int S1[N2], S2[N3], v[N3];
ll N;

ll inv(ll x){return x<=1?1:inv(P%x)*(P-P/x)%P;}
void init(int n) {
	phi[1] = 1;
	REP(i,2,n-1) {
		if (!vis[i]) p[++cnt]=i,phi[i]=i-1;
		for (int j=1,t;j<=cnt&&i*p[j]<=n; ++j) {
			vis[t=i*p[j]]=1;
			if (i%p[j]==0) {phi[t]=phi[i]*p[j];break;}
			phi[t]=phi[i]*phi[p[j]];
		}
	}
	REP(i,1,n) S1[i]=(S1[i-1]+(ll)i*i%P*phi[i]%P)%P;
}
ll sqr(ll n) {n%=P;return n*n%P;}
ll s_1(ll n) {n%=P;return n*(n+1)%P*inv2%P;}
ll s_2(ll n) {n%=P;return n*(n+1)%P*(2*n+1)%P*inv6%P;}
ll s_3(ll n) {n%=P;return sqr(s_1(n));}
int sum(ll n) {
	if (n<N2) return S1[n]; int x = N/n;
	if (v[x]) return S2[x]; v[x] = 1;
	int &ans = S2[x] = s_3(n), mx = sqrt(n);
	REP(i,2,mx) (ans-=sqr(i)*sum(n/i)%P)%=P;
	for (ll i=mx+1,j,k=n/i; i<=n; i=j+1,--k) {
		(ans-=(s_2(j=n/k)-s_2(i-1))*sum(k)%P)%=P;
	}
	return ans;
}

int main() {
	scanf("%d%lld", &P, &N);
	inv2 = inv(2), inv6 = inv(6);
	init(N2-1);
	int mx = sqrt(N), ans = 0;
	REP(i,1,mx) (ans+=(S1[i]-S1[i-1])*s_3(N/i)%P)%=P;
	for (ll i=mx+1,j,k=N/i; i<=N; i=j+1,--k) {
		(ans+=(sum(j=N/k)-sum(i-1))*s_3(k)%P)%=P;	
	}
	if (ans<0) ans+=P;
	printf("%d
", ans);
}

练习4. luogu P3172 [CQOI2015]选数

大意: 求选择$n$个范围[L,R]的数, 且$gcd=k$的方案数

相当于求$sumlimits_{Lle a_1 le R}cdotssumlimits_{Lle a_n le R}[gcd=k]$

等价于$sumlimits_{lle a_1 le r}cdotssumlimits_{lle a_n le r}[gcd=1]$, 其中$l=lfloor frac{L-1}{K} floor+1, r=frac{R}{K}$

再反演一下就没了.....$sumlimits_{d=1}^{r}mu(d)(lfloor frac{r}{d} floor-lfloor frac{l-1}{d} floor)^n$

这题因为有$lfloor frac{r}{d} floor$和$lfloor frac{l-1}{d} floor$两部分, 所以为了方便直接用map记忆化了

这题刚开始想错了, 反演成$varphi$了, 结果改了半天................

const int N2 = 2e6+10, N3 = 2000, P = 1e9+7;
int mu[N2], vis[N2], p[N2];
int n, k, l, r, cnt;
unordered_map<int,int> S2;
ll qpow(ll a,ll n) {ll r=1%P;for (a%=P;n;a=a*a%P,n>>=1)if(n&1)r=r*a%P;return r;} 
void init(int n) {
    mu[1] = 1;
    REP(i,2,n-1) {
        if (!vis[i]) p[++cnt]=i,mu[i]=-1;
        for (int j=1,t;j<=cnt&&i*p[j]<=n; ++j) {
            vis[t=i*p[j]]=1;
            if (i%p[j]==0) {mu[t]=0;break;}
            mu[t] = -mu[i];
        }
        mu[i] += mu[i-1];
    }
}
int sum(int n) {
    if (n<N2) return mu[n];
    if (S2.count(n)) return S2[n];
    int ans = 1, mx = sqrt(n);
    REP(i,2,mx) (ans-=sum(n/i))%=P;
    for (int i=mx+1,j,k=n/i; i<=n; i=j+1,--k) {
        (ans-=((j=n/k)-i+1)*sum(k)%P)%=P;
    }
    return S2[n]=ans;
}

int main() {
    scanf("%d%d%d%d", &n, &k, &l, &r);
    l = (l-1)/k, r /= k;
    init(N2-1);
    int ans = 0;
    for (int i=1,j; i<=r; i=j+1) {
        j = r/(r/i);
        if (l/i) j = min(j,l/(l/i));
        (ans+=(sum(j)-sum(i-1))*qpow(r/i-l/i,n)%P)%=P;
    }
    if (ans<0) ans+=P;
    printf("%d
", ans);
}

练习5. Luogu P3327 [SDOI2015]约数个数和

大意: 求$sumlimits_{i=1}^nsumlimits_{j=1}^md(ij)$, $d(n)$为除数函数

对于除数函数可以这样处理$d(ij)=sumlimits_{x|i}sumlimits_{y|i}[xperp y]=sumlimits_{x|i}sumlimits_{y|i}sumlimits_{d|gcd(x,y)}mu(d)$

所以就可以得到$sumlimits_{i=1}^nsumlimits_{j=1}^md(ij)=sumlimits_{d=1}^{min(n,m)}varphi(d)sumlimits_{i=1}^{lfloor frac{n}{d} floor}lfloor frac{n}{id} floor sumlimits_{j=1}^{lfloor frac{m}{d} floor}lfloor frac{m}{jd} floor$

然后就没了....主要这题的数据范围太小, 不需要杜教筛也行