莫比乌斯反演基础

tonyfang posted @ 2017年1月16日 10:22 in 随笔 with tags c++ OI , 314 阅读
本来想写一篇数论基础的,后来发现太长了省选之前填不完就弃了
所以现在是莫比乌斯反演基础。

1 定义

$F(n) = \sum_{d|n}f(d)$可以反演至$f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d)$或$f(n)=\sum_{d|n}\mu(d)F(\frac{n}{d})$。
可以理解为,前后等价。
$\mu(i)$为莫比乌斯函数。
1. 如果$i=1$,则$\mu(i)=1$
2. 如果$i=p_1p_2p_3...p_k$, $\{p_k\}$各不相同的素数,那么$\mu(i)=(-1)^k$
3. 否则$\mu(i)=0$
 
P.S 以下预备知识需要你熟练掌握:线性筛,数学推导。
 

2 例题

2.1 BZOJ 2301 Problem B

$n$次询问,每次询问多少个$(x,y)$满足,$a\leq x\leq b, c\leq y\leq d$,且$(x,y)=k$。
$1 \leq n, a, b, c, d \leq 5\times 10^4$
【题解】
首先询问拆分。
也就是现在要求求出多少个$(x,y)$满足$1\leq x\leq n, 1\leq y \leq m$,且$(x,y)=k$。
接下来就是莫比乌斯反演的事情啦!
那我们令$f(i)$为$1 \leq x \leq n, 1 \leq y \leq m$且$(x,y)=i$的数对$(x,y)$的个数。
令$F(i)$为为$1 \leq x \leq n, 1 \leq y \leq m$且$i|(x,y)$的数对$(x,y)$的个数。
那么就可以反演啦!
显然可以得到$F(i)=\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{i}\rfloor$
反演可得
$f(i) = \sum_{i|d}\mu(\frac{d}{i})F(d)$
$f(i) = \sum_{i|d}\mu(\frac{d}{i})\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$
显然原题要求$f(k)$,那么枚举$k$的每一个倍数$d$即可。
但是这样的复杂度是$O(n)$的,考虑如何优化。
观察发现,$\lfloor \frac{n}{d}\rfloor$最多有$2\sqrt{n}$个取值。
为什么呢?
分析$\lfloor \frac{n}{d} \rfloor (1 \leq d \leq n)$的取值:
1. 当$1\leq d\leq \lfloor \sqrt{n} \rfloor$时,$\lfloor \frac{n}{d} \rfloor$最多有$\sqrt{n}$种不同取值。
2. 当$\lfloor \sqrt{n} \rfloor +1 \leq d \leq n$时,由于$\lfloor \frac{n}{d} \rfloor \leq \lfloor \sqrt{n} \rfloor$,所以也最多有$\sqrt{n}$种不同取值。
综上,最多共有$2\sqrt{n}$个取值。
同样的,$\lfloor \frac{m}{d}\rfloor$最多有$2\sqrt{m}$个取值。
那么问题来啦!
$\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$有多少个取值呢?
啥?你说有$4\sqrt{nm}$个?Naive!
所以我们枚举这$2(\sqrt{n}+\sqrt{m})$个取值,对莫比乌斯函数维护前缀和,就能在$O(\sqrt{n})$内解答了。
那么,如何枚举除法呢?
if(n>m) swap(n, m);
for (int i=1, last; i<=n; i=last+1) {
    last = min(n/(n/i), m/(m/i));
    ret += (n/i)*(m/i)*(sum[last]-sum[i-1]);
}
return ret;
 
是不是很简单啦?
那么代码实现也就水到渠成了。
# include <stdio.h>
# include <algorithm>

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;

const int M = 500010;

int q, a, b, c, d, k;
int p[M], pn=0;
int mu[M], s[M];
bool wait[M];

inline ll solve(int n, int m) {
	ll res = 0;
	if(n>m) swap(n, m);
	for (int i=1, last=0; i<=n; i=last+1) {
		last = min(n/(n/i), m/(m/i));
		res += (ll)(s[last] - s[i-1]) * (n/i) * (m/i);
	}
	return res;
}

int main() {
	scanf("%d", &q);
	mu[1] = 1;
	for (int i=2; i<=500000; ++i) {
		if(!wait[i]) {
			p[++pn] = i;
			mu[i] = -1;
		}
		for (int j=1; j<=pn&&i*p[j]<=500000; ++j) {
			wait[i*p[j]] = 1;
			if(i%p[j]==0) {
				mu[i*p[j]] = 0;
				break;
			} else mu[i*p[j]] = -mu[i];
		}
	}
	for (int i=1; i<=500000; ++i)
		s[i] = s[i-1] + mu[i];
	while(q--) {
		scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
		long long ans = solve(b/k, d/k) - solve((a-1)/k, d/k) - solve(b/k, (c-1)/k) + solve((a-1)/k, (c-1)/k);
		printf("%lld\n", ans);
	}
	return 0;
}

2.2 BZOJ2820 YY的GCD

求有多少个数对$(x,y)$满足$1\leq x\leq n, 1\leq y\leq m$且$(x,y)$是质数。
$n \leq 10^7, T\leq 10^4$
【题解】
设$(x,y)=p$,那么问题转化为对于每个**质数**$p$,求出有多少个数对$(x,y)$满足$1 \leq x \leq \lfloor \frac{n}{p} \rfloor, 1 \leq y \leq \lfloor \frac{m}{p} \rfloor$,且$(x,y)=1$。
那么答案就是
$ans=\sum_p^{min(n,m)} \sum_{d=1}^{min(n,m)} \mu(d) \lfloor \frac{n}{pd} \rfloor \lfloor \frac{m}{pd} \rfloor$
令$pd=T$,则
$ans = \sum_T^{min(n,m)} \sum_{p|T} \mu(\frac{T}{p}) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor$
定值可以提取,则
$ans = \sum_T^{min(n,m)} \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T}\rfloor \sum_{p|T} \mu(\frac{T}{p})$
所以暴力枚举每个质数,更新倍数,处理前缀和即可。复杂度$O(n)$。
总复杂度$O(n)$。
需要较强的公式推导和代换能力。
# include <stdio.h>
# include <algorithm>
// # include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
const int M = 10000010;
const int P = 4000010;
bool wait[M];
int p[P], pn=0;
int mu[M], n, m;
ll s[M];

inline ll solve(int n, int m) {
	ll res=0;
	if(n>m) swap(n, m);
	for (int i=1, last=0; i<=n; i=last+1) {
		last=min(n/(n/i), m/(m/i));
		res=res+(ll)(n/i)*(m/i)*(s[last]-s[i-1]);
	}
	return res;
}

int main() {
	mu[1] = 1;
	for (int i=2; i<=10000000; ++i) {
		if(!wait[i]) {
			p[++pn] = i;
			mu[i] = -1;
		}
		for (int j=1; j<=pn&&i*p[j]<=10000000; ++j) {
			wait[i*p[j]] = 1;
			if(i%p[j]==0) {
				mu[i*p[j]] = 0;
				break;
			} else mu[i*p[j]] = -mu[i];
		}
	}
	for (int i=1; i<=pn; ++i) 
		for (int j=1; j*p[i]<=10000000; j++) s[j*p[i]] += mu[j];
	for (int i=1; i<=10000000; ++i) s[i] += s[i-1];
	int T;
	scanf("%d", &T);
	while(T--) {
		scanf("%d%d", &n, &m);
		printf("%lld\n", solve(n, m));
	}
	return 0;
}

2.3 BZOJ 3529 数表

令$F(i)$为$i$的约数和,$q$次给定$n,m,a$,求$\sum_{1\leq i\leq n, 1\leq j\leq m, F(gcd(i,j)) \leq a} F(gcd(i,j))~mod~ 2^{31}$
$n, m \leq 10^5, q \leq 20000, a \leq 10^9$
【题解】
假定没有$a$的限制,只是求$1 \leq x \leq n, 1 \leq y \leq m, (x,y)=i$的对数,那么:
$f(i) = \sum_{i|d}\mu(\frac{d}{i})\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$
$F(i)$可以利用线性筛处理出来。
考虑答案怎么表示:
$ans = \sum_{i=1}^{min(n,m)} F(i) \sum_{i|d} \mu(\frac{d}{i})\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$
把$\sum$换位,
$ans = \sum_{d=1}^{min(n,m)} \lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor \sum_{i|d}F(i)\mu(\frac{d}{i})$
对比上一题,同理暴力枚举每个质数,更新倍数,处理前缀和即可。
那么考虑了$a$的限制了呢?
离线,对$a$排序,按照$F(i)$将$i$排序,发现从小到大每次只要维护前缀和的增加即可,用树状数组。
$O(nlog^2n + q\sqrt{n}logn)$
# include <stdio.h>
# include <string.h>
# include <algorithm>

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;

const int M = 500010, Ms = 200000;

bool wait[M];
int p[M], pn=0, mu[M];
struct Fs {
	int s, id;
	inline void init(int x) {
		s = 0; id = x;
	}
	friend bool operator <(Fs a, Fs b) {
		return a.s<b.s || (a.s==b.s && a.id<b.id);
	}
}F[M];

int Q;
struct ques {
	int n, m, a, id;
	friend bool operator <(ques a, ques b) {
		return a.a<b.a;
	}
}q[M];

int ans[M];

struct BIT {
	int c[M], n;
	inline void init(int x) {
		n = x;
		memset(c, 0, sizeof(c));
	}
	inline int lb(int x) {return x&(-x);}
	inline void change(int x, int del) {
		for (; x<=n; x+=lb(x)) c[x]+=del;
	}
	inline int ask(int x) {
		int ret=0;
		for (; x; x-=lb(x)) ret+=c[x];
		return ret;
	}
}T;

inline int solve(int n, int m) {
	if(n>m) swap(n, m);
	int ret=0;
	for (int i=1, last=0; i<=n; i=last+1) {
		last = min(n/(n/i), m/(m/i));
		ret += (n/i)*(m/i)*(T.ask(last)-T.ask(i-1));
	}
	return ret;
}

int main() {
	
	mu[1] = 1;
	for (int i=2; i<=Ms; ++i) {
		if(!wait[i]) {
			p[++pn] = i;
			mu[i] = -1;
		}
		for (int j=1; j<=pn&&i*p[j]<=Ms; ++j) {
			wait[i*p[j]] = 1;
			if(i%p[j]==0) {
				mu[i*p[j]]=0;
				break;
			} 
			mu[i*p[j]] = -mu[i];
		}
	}
	
	for (int i=1; i<=Ms; ++i)
		F[i].init(i);
	
	for (int i=1; i<=Ms; ++i)
		for (int j=i; j<=Ms; j+=i)
			F[j].s += i;
	
	sort(F+1, F+Ms+1);
	
	scanf("%d", &Q);
	for (int i=1; i<=Q; ++i) {
		scanf("%d%d%d", &q[i].n, &q[i].m, &q[i].a);
		q[i].id=i;
	}
	
	sort(q+1, q+Q+1);
	
	T.init(Ms);
	
	int cur = 0;
	for (int i=1; i<=Q; ++i) {
		while(cur+1 <= Ms && F[cur+1].s<=q[i].a) {
			++cur;
			for (int j=1; j*F[cur].id<=Ms; ++j)
				T.change(j*F[cur].id, F[cur].s*mu[j]);
		}
		ans[q[i].id] = solve(q[i].n, q[i].m);
	}
	
	for (int i=1; i<=Q; ++i) printf("%d\n", ans[i]&0x7fffffff);

	return 0;
}

2.4 BZOJ2154 Crash的数字表格 (BZOJ2693 jzptab)

给$n,m$,求$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)$。
$1 \leq n, m \leq 10^7$
【题解】
$ans = \sum_{i=1}^{n}\sum_{j=1}^{m} \frac{ij}{(i,j)}$
枚举$(i,j)=d$。
令$F(x,y) = \sum_{1\leq i \leq x, 1\leq j \leq y}^{(i,j)=1}ij$
那么$ans = \sum_{d=1}^{min(n,m)}\frac{d^2F(\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor)}{d} = \sum_{d=1}^{min(n,m)}dF(\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor)$
$S(x,y) = \sum_{i=1}^{x} \sum_{j=1}^{y} ij = \frac{xy(x+1)(y+1)}{4}$
$F(x,y) = \sum_{i=1}^{min(x,y)}\mu(i)\times i^2\times S(\lfloor \frac{x}{i} \rfloor, \lfloor \frac{y}{i}\rfloor)$
对两个式子进行计算,复杂度分别为$O(\sqrt{n})$和$O(\sqrt{n})$。乘起来就是$O(n)$。
进一步,如果有多组询问呢?
$ans = \sum_{d=1}^{min(n,m)}d \sum_{i=1}^{min(\lfloor \frac{n}{d} \rfloor,\lfloor \frac{m}{d} \rfloor)}\mu(i)\times i^2\times S(\lfloor \frac{x}{di} \rfloor, \lfloor \frac{y}{di}\rfloor)$
令$T = di$
那么$ans = \sum_{d=1}^{min(n,m)}d \sum_{i=1}^{min(\lfloor \frac{n}{d} \rfloor,\lfloor \frac{m}{d} \rfloor)}\mu(i)\times i^2\times S(\lfloor \frac{x}{di} \rfloor, \lfloor \frac{y}{di}\rfloor)$
$ans = \sum_{D=1}^{min(n,m)}\sum_{i|D}\frac{D}{i}\times\mu(i)\times i^2 \times S(\lfloor\frac{x}{D}\rfloor\lfloor\frac{y}{D}\rfloor)$
发现后面那一坨跟$i$无关,提前。
$ans = \sum_{D=1}^{min(n,m)}S(\lfloor\frac{x}{D}\rfloor\lfloor\frac{y}{D}\rfloor)\sum_{i|D}\frac{D}{i}\times\mu(i)\times i^2$
$ans = \sum_{D=1}^{min(n,m)}S(\lfloor\frac{x}{D}\rfloor\lfloor\frac{y}{D}\rfloor)\sum_{i|D}D\times\mu(i)\times i$。
令$g(D)=\sum_{i|D}D\times\mu(i)\times i$。
则$g(D)=D\sum_{i|D}i\times \mu_i$。
对于后面维护前缀和即可。
如果$D$是质数,那么$g(D)=D\times (1\times 1 + D \times (-1)) = -D^2+D$。
如果我们已经知道了$g(D')$,$g(D=D' \times p)$怎么求呢?
首先很好证明$g(D)$是积性函数,那么如果$(D', p)=1$,则$g(D)=g(D')\times g(p)$。
如果$(D', p) = 0$,那么根据线性筛的性质,$g(D)=g(D')\times p$。
那么就能线性筛啦!
预处理$O(n)$,单次询问$O(\sqrt{n})$。
# include <stdio.h>
# include <algorithm>

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;

const int M = 10000000, mod = 20101009, P = 4000010;

bool wait[M+10];
int p[P], pn=0, n, m;
int g[M+10];

inline int S(int x, int y) {
	return 1ll * (1ll*x*(x+1)/2%mod) * (1ll*y*(y+1)/2%mod) % mod;
}

int main() {
	g[0] = 0, g[1] = 1;
	for (int i=2; i<=M; ++i) {
		if(!wait[i]) {
			p[++pn] = i;
			g[i] = 1ll * i * (1-i+mod) % mod;
		}
		for (int j=1; j<=pn && i*p[j] <= M; ++j) {
			wait[i*p[j]] = 1;
			if(i%p[j] == 0) {
				g[i*p[j]] = 1ll * g[i] * p[j] % mod;
				break;
			} else g[i*p[j]] = 1ll * g[i] * g[p[j]] % mod;
		}
	}
	for (int i=1; i<=M; ++i) g[i] = (g[i-1] + g[i]) % mod;
	
	scanf("%d%d", &n, &m);
	int ans = 0;
	if(n>m) swap(n, m);
	for (int i=1, last=0; i<=n; i=last+1) {
		last = min(n/(n/i), m/(m/i));
		ans = ans + 1ll * S(n/i, m/i) * (g[last] - g[i-1] + mod) % mod;
		ans %= mod;
	}
	printf("%d\n", ans);
	return 0;
}
# include <stdio.h>
# include <algorithm>

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;

const int M = 10000000, mod = 100000009, P = 4000010;

bool wait[M+10];
int p[P], pn=0, n, m;
int g[M+10];

inline int S(int x, int y) {
	return 1ll * (1ll*x*(x+1)/2%mod) * (1ll*y*(y+1)/2%mod) % mod;
}

int main() {
	g[0] = 0, g[1] = 1;
	for (int i=2; i<=M; ++i) {
		if(!wait[i]) {
			p[++pn] = i;
			g[i] = 1ll * i * (1-i+mod) % mod;
		}
		for (int j=1; j<=pn && i*p[j] <= M; ++j) {
			wait[i*p[j]] = 1;
			if(i%p[j] == 0) {
				g[i*p[j]] = 1ll * g[i] * p[j] % mod;
				break;
			} else g[i*p[j]] = 1ll * g[i] * g[p[j]] % mod;
		}
	}
	for (int i=1; i<=M; ++i) g[i] = (g[i-1] + g[i]) % mod;
	int T; scanf("%d", &T);
	while(T--) {
		scanf("%d%d", &n, &m);
		int ans = 0;
		if(n>m) swap(n, m);
		for (int i=1, last=0; i<=n; i=last+1) {
			last = min(n/(n/i), m/(m/i));
			ans = ans + 1ll * S(n/i, m/i) * (g[last] - g[i-1] + mod) % mod;
			ans %= mod;
		}
		printf("%d\n", ans);
	}
	return 0;
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter