📚【笔记】莫比乌斯反演
前置
-
积性函数
-
积性函数：满足\(f_1 = 1\)且\(\forall x,y\in \mathbb{N}^*,\gcd(x,y) = 1\)都有\(f(x\cdot y)=f(x)\cdot f(y)\)的数论函数\(f_n\)。
-
完全积性函数：满足\(f_1 = 1\)且\(\forall x,y \in \mathbb{N}^*\)都有\(f(x\cdot y) = f(x)\cdot f(y)\)的数论函数\(f_n\)。
-
-
常见数论函数
-
单位函数：\(\epsilon(n) = (n = 1)\)。（完全积性函数）
-
恒等函数：\(\operatorname{id}_k(n) = n^k\)。常用\(\operatorname{id}_1(n)\)，可记为\(\operatorname{id}(n)\)。（完全积性函数）
-
常数函数：\(\operatorname{1}(n) = 1\)。（完全积性函数）
-
因数函数：\(\sigma_k(n)=\sum\limits_{d\mid n} d^k\)。常用\(\sigma_0(n)\)，记作\(\tau(n)\)（因数个数）。以及\(\sigma_1(n)\)，记作\(\sigma(n)\)（因数和）。
-
欧拉函数：\(\varphi(n) = \sum\limits_{i = 1}^{n}\left[\gcd(i,n) = 1\right]\)
-
-
数论卷积（狄利克雷卷积）
-
\(\left(f*g\right)(n) = \sum\limits_{d\mid n}f(d)\cdot g(\frac{n}{d})\)。
-
有：交换律、结合律、分配律。
-
单位元：狄利克雷卷积的单位元是单位函数\(\epsilon(n)\)。
-
逆：称满足\(f^{-1}*f = \epsilon\)的数论函数\(f^{-1}\)为数论函数\(f\)的狄利克雷逆。
-
有\(\operatorname{id}_k * 1 = \sigma_k\) ， $ \varphi * \operatorname{1} = \operatorname{id} $ ， \(\operatorname{1}*\operatorname{1}=\operatorname{id}\)
-
-
数论分块
用于处理类似\(\sum\limits_{i = 1}^{n} f(i)\cdot \left\lfloor\frac{n}{i}\right\rfloor\)的式子。详见这里。
莫比乌斯反演
- 莫比乌斯函数
-
定义莫比乌斯函数\(\mu(n)=\begin{cases} 1 & n = 1.\\ 0 & \exists \ d\in \mathbb{N}^*,d^2\mid n.\\ (-1)^k & \operatorname{otherwise}.\end{cases}\)，其中\(k\)为\(n\)本质不同质因子个数。
-
性质1：\(\mu * \operatorname{1} = \epsilon\)，或者说\(\sum\limits_{d\mid n}\mu(d) = (n=1)\)，通过二项式定理可证。（其实实践中，形式\(\left[\gcd(i,j) = 1\right] = \sum\limits_{d\mid \gcd(i,j)}\mu(d)\)极为常见）
-
性质2：\(\sum\limits_{d\mid n}\frac{\mu(d)}{d} = \frac{\varphi(n)}{n}\)。
- 莫比乌斯反演
- 莫比乌斯反演：\(f = g*\operatorname{1}\implies g = f*\mu\)。或者说\(f(n) = \sum\limits_{d\mid n}g(d)\implies g(n) = \sum\limits_{d\mid n}\mu(d)\cdot f(\frac{n}{d})\).
- 莫比乌斯反演乘积形式：\(f(n)=\prod\limits_{d\mid n}g(d)\implies g(n)=\prod\limits_{d\mid n}{f(d)}^{\mu(\tfrac{n}{d})}\)
keys
- \(\large\sum\limits_{i = 1}^{n}i^2 = \frac{n\cdot(n+1)\cdot(2\cdot n+1)}{6}\)
我导的很奇怪：
\(\text{设}F_n = \sum\limits_{i = 1}^{n}i^2\)
\(\begin{aligned}\text{则}F_n &= \sum\limits_{i = 1}^{n}\sum\limits_{j = i}^{n} j \\ &= \sum\limits_{i = 1}^{n}\frac{n \cdot (n+1)}{2}-\frac{i \cdot (i-1)}{2} \\ &= \frac{n^2 \cdot (n+1)}{2}-\sum\limits_{i = 1}^{n}\frac{i \cdot (i-1)}{2} \\ &= \frac{n^2 \cdot (n+1)}{2}-\frac{\sum\limits_{i = 1}^{n}i^2}{2}+\frac{n \cdot (n+1)}{4}\end{aligned}\)
\(\text{于是}F_n = \frac{n \cdot (n+1)}{4}+\frac{n^2 \cdot (n+1)}{2}-\frac{F_n}{2}\)。
□
易得。
-
\(\large\tau(n\times m) =\sum\limits_{i \mid n}\sum\limits_{j \mid m}\left[\gcd(i,j) = 1\right]\)
-
\(\sum\limits_{a_1 = 1}^{n}\sum\limits_{a_2 = 1}^{n}\ldots\sum\limits_{a_m = 1}^{n}\left[\gcd\limits_{i = 1}^{m}\{a_i\} = 1\right] = \sum\limits_{d = 1}^{n}\mu(d) \times \left\lfloor\frac{n}{d}\right\rfloor^m\)
例题
YY的GCD
给你两个整数 \(n,m\) ，试求 \(\sum\limits_{i = 1}^{n}\sum\limits_{j = 1}^{m}\left[\gcd(i,j) \in prime\right]\) 。多组询问。
題解
一种方法：
然后莫反一下：
再化：
此时的时间复杂度在 \(O(T\times \pi(\min(n,m))\times\sqrt n)\) ，仍不可接受。
我们设 \(T = k\times d\) ，于是：
可得：
只要筛出来 \(g(n) = \sum\limits_{\substack{k\mid n\\k \in prime}}\mu\left(\left\lfloor\frac{n}{k}\right\rfloor\right)\) 即可。
#include <iostream>
const int N = 1e7+10;
const int P = 1e6+10;
int prime[P];
int minprime[N];
int mu[N], g[N];
long long pre[N];
void line(int max_size) {
register int upd;
register int m_h;
mu[1] = 1;
for(register int i = 2;i <= max_size;++i) {
if(!minprime[i]) {
prime[++prime[0]] = i;
minprime[i] = i;
mu[i] = -1;
}
if((i<<1) <= max_size) {
upd = std :: min(max_size/i,minprime[i]);
for(register int* j = prime+1;*j&&*j <= upd;++j) {
minprime[m_h = i*(*j)] = *j;
if(*j != minprime[i])
mu[m_h] = -mu[i];
}
}
}
for(register int i = 1;i <= prime[0];++i)
for(register int j = 1;j*prime[i] <= max_size;++j)
g[j*prime[i]] += mu[j];
for(register int i = 1;i <= max_size;++i)
pre[i] = pre[i-1]+g[i];
}
int T, n, m;
long long ans;
int main() {
line(10000000);
scanf("%d",&T);
for(register int outs = 1;outs <= T;++outs) {
scanf("%d %d",&n,&m);
if(n > m)
std :: swap(n,m);
ans = 0;
for(register int l = 1, r;l <= n;l = r+1) {
r = std :: min(n/(n/l),m/(m/l));
ans += ((long long)(n/l)*(m/l))*(pre[r]-pre[l-1]);
}
printf("%lld\n",ans);
}
return 0;
}
数字表格
给你两个整数 \(n,m\;,\;(n,m \in \left[1,10^7\right])\)
试求：\(\sum\limits_{i = 1}^{n}\sum\limits_{j = 1}^{m}\operatorname{lcm}(i,j)\)
題解
也是好题。
首先 \(\operatorname{lcm}\) 不是很好化，先改成 \(\gcd\)。
老套路，枚举 \(\gcd\)：
再利用经典公式 \(\left[\gcd(i,j) = 1\right] = \sum\limits_{d \mid \gcd(i,j)} \mu(d)\)：
首先对于前面的 \(k\) 这一部分数论分块一次，然后再拿 \(\left\lfloor\frac{n}{k}\right\rfloor\) 和 \(\left\lfloor\frac{m}{k}\right\rfloor\) 再做一次数论分块。预处理出来 \(\mu(d) \cdot d^2\) 的前缀和即可。
时间复杂度 \(O(n+m)\)。
#include <iostream>
#define llt long long int
const int N = 1e7+10;
const int P = 4e6+10;
const llt moyn = 20101009LL;
int mu[N];
int minprime[N];
int prime[P];
void init(int init_size) {
const int half_size = init_size/2;
register int i, *j, upd, m_h;
register int *p = prime+1;
mu[1] = 1;
for(i = 2;i <= half_size;++i) {
if(!minprime[i]) {
minprime[i] = i;
*p++ = i;
mu[i] = -1;
}
if((upd = init_size/i) >= minprime[i]) {
for(j = prime+1;*j < minprime[i];++j) {
minprime[m_h = i*(*j)] = *j;
mu[m_h] = -mu[i];
}
minprime[m_h = i*(*j)] = *j;
mu[m_h] = 0;
} else {
for(j = prime+1;*j <= upd;++j) {
minprime[m_h = i*(*j)] = *j;
mu[m_h] = -mu[i];
}
}
mu[i] = mu[i]*(llt)i%moyn*i%moyn;
mu[i] = (mu[i]+mu[i-1]+moyn)%moyn;
}
for(;i <= init_size;++i) {
if(!minprime[i])
mu[i] = -1;
mu[i] = mu[i]*(llt)i%moyn*i%moyn;
mu[i] = (mu[i]+mu[i-1]+moyn)%moyn;
}
}
int n, m;
int sum(int x,int y) {
return (x*1LL*(x+1LL)/2%moyn)*(y*1LL*(y+1LL)/2%moyn)%moyn;
}
int f(int x,int y) {
int res = 0;
for(register int l = 1, r;l <= x;l = r+1) {
r = std :: min(x/(x/l),y/(y/l));
res = (res+(mu[r]-mu[l-1]+moyn)%moyn*(llt)sum(x/l,y/l))%moyn;
}
return res;
}
int work(int x,int y) {
int res= 0;
for(register int l = 1, r;l <= n;l = r+1) {
r = std :: min(n/(n/l),m/(m/l));
res = (res+(((r-l+1LL)*(llt)(l+r))/2)%moyn*f(n/l,m/l))%moyn;
}
return res;
}
int main() {
scanf("%d %d",&n,&m);
if(n > m)
std :: swap(n,m);
init(n);
printf("%d\n",work(n,m));
return 0;
}
DZY Loves Math
題解
好題++
首先原式不難化：
然後考慮如何線性遞推 \(h(n) = \sum\limits_{d \mid n}\operatorname{f}(d) \cdot \mu\left(\frac{n}{d}\right)\)。
先擺個結論吧（笑）：
證明：
首先由於對於一個有平方因子的數 \(n\)，\(\mu(n) = 0\)。也就是這個 \(n\) 是沒有貢獻的。所以當我們選擇 \(d = \prod\limits_{i = 1}^{k}p_i^{\beta_i}\) 時，一定要保證 \(\forall i \in \left[1,k\right],0 \le \alpha_i-\beta_i \le 1\)。
-
對於 \(\exists\ i,j \in \left[1,k\right],\alpha_i \lt \alpha_j\) 的情況，我們考慮對於 \(\beta_i\) 而言，無論選擇 \(\alpha_i\) 還是 \(\alpha_{i}-1\) 都沒有影響 \(\operatorname{f}(d)\)，但是 \(\mu\left(\frac{n}{d}\right)\) 會分別取 \(1\) 和 \(-1\)，所以兩者相消，不會產生任何貢獻。最後 \(\operatorname{h}(n) = 0\)。
-
否則 \(\operatorname{h}(n)\) 中除了 \(\beta_i\) 都選擇相應的 \(\alpha_i - 1\) 這一種情況外，其它都有 \(\operatorname{f}(d)\) 相等。最後兩兩相消，只剩下一個 \(\operatorname{f}\left(\frac{n}{\prod_{i = 1}^{k}p_i} \times p_j\right) \cdot \mu\left(\frac{\prod_{i = 1}^{k}p_i}{p_j}\right) - \operatorname{f}\left(\frac{n}{\prod_{i = 1}^{k}p_i}\right) \cdot \mu\left(\prod_{i = 1}^{k}p_i\right)\)，其實就是 \(k \cdot (-1)^{k-1} - (k-1) \cdot (-1)^{k} = (-1)^{k+1}\)。
然後開始推 \(\operatorname{h}(n)\)：
const int N = 1e7+10;
int minprime[N];
int prime[N];
int withoutmin[N];
short pow_cnt[N];
short h[N];
void line(int max_size) {
const int half_size = max_size>>1;
register int i, upd, m_h;
register int *pm = prime;
for(i = 2;i <= half_size;++i) {
if(!minprime[i]) {
*++pm = minprime[i] = i;
pow_cnt[i] = withoutmin[i] = h[i] = 1;
}
if((upd = max_size/i) >= minprime[i]) {
register int *j;
for(j = prime+1;*j != minprime[i];++j) {
minprime[m_h = i*(*j)] = (*j);
pow_cnt[m_h] = 1;
withoutmin[m_h] = i;
if(pow_cnt[i] == 1)
h[m_h] = -h[i];
else
h[m_h] = 0;
}
minprime[m_h = i*(*j)] = (*j);
pow_cnt[m_h] = pow_cnt[i]+1;
withoutmin[m_h] = withoutmin[i];
if(withoutmin[m_h] == 1)
h[m_h] = h[i];
else if(pow_cnt[withoutmin[m_h]] == pow_cnt[m_h])
h[m_h] = -h[withoutmin[m_h]];
} else {
for(register int *j = prime+1;*j <= upd;++j) {
minprime[m_h = i*(*j)] = (*j);
pow_cnt[m_h] = 1;
withoutmin[m_h] = i;
if(pow_cnt[i] == 1)
h[m_h] = -h[i];
else
h[m_h] = 0;
}
}
}
for(;i <= max_size;++i)
if(!minprime[i])
h[i] = 1;
for(register int i = 2;i <= max_size;++i)
h[i] += h[i-1];
}
於是易得了：
#include <iostream>
const int N = 1e7+10;
int minprime[N];
int prime[N];
int withoutmin[N];
short pow_cnt[N];
short h[N];
void line(int max_size) {
const int half_size = max_size>>1;
register int i, upd, m_h;
register int *pm = prime;
for(i = 2;i <= half_size;++i) {
if(!minprime[i]) {
*++pm = minprime[i] = i;
pow_cnt[i] = withoutmin[i] = h[i] = 1;
}
if((upd = max_size/i) >= minprime[i]) {
register int *j;
for(j = prime+1;*j != minprime[i];++j) {
minprime[m_h = i*(*j)] = (*j);
pow_cnt[m_h] = 1;
withoutmin[m_h] = i;
if(pow_cnt[i] == 1)
h[m_h] = -h[i];
else
h[m_h] = 0;
}
minprime[m_h = i*(*j)] = (*j);
pow_cnt[m_h] = pow_cnt[i]+1;
withoutmin[m_h] = withoutmin[i];
if(withoutmin[m_h] == 1)
h[m_h] = h[i];
else if(pow_cnt[withoutmin[m_h]] == pow_cnt[m_h])
h[m_h] = -h[withoutmin[m_h]];
} else {
for(register int *j = prime+1;*j <= upd;++j) {
minprime[m_h = i*(*j)] = (*j);
pow_cnt[m_h] = 1;
withoutmin[m_h] = i;
if(pow_cnt[i] == 1)
h[m_h] = -h[i];
else
h[m_h] = 0;
}
}
}
for(;i <= max_size;++i)
if(!minprime[i])
h[i] = 1;
for(register int i = 2;i <= max_size;++i)
h[i] += h[i-1];
}
int n, m, T;
long long ans;
int main() {
line(10000000);
scanf("%d",&T);
while(T--) {
scanf("%d %d",&n,&m);
ans = 0;
if(n > m)
std :: swap(n,m);
for(register int l = 1, r;l <= n;l = r+1) {
r = std :: min(n/(n/l),m/(m/l));
ans += 1LL*(n/l)*(m/l)*1LL*(h[r]-h[l-1]);
}
printf("%lld\n",ans);
}
return 0;
}
简单的数学题
题解
先化式子：
考虑用杜教筛筛出 \(\operatorname{f}(n) = n^2 \times \sum\limits_{d \mid n} d \cdot \mu\left(\frac{n}{d}\right)\)。
考虑到 \(\sum\limits_{d \mid n} d \cdot \mu\left(\frac{n}{d}\right)\) 的意义是 \(\operatorname{id} * \mu\)，结果为 \(\varphi\)。所以 \(\operatorname{f}(n) = n^2 \cdot \varphi(n)\)
设 \(\operatorname{h} = \operatorname{g} * \operatorname{f}\)，则有 \(\operatorname{h}(n) = \sum\limits_{d \mid n} \operatorname{g}(d) \cdot \left(\frac{n}{d}\right)^2 \cdot \varphi(n)\)。
发现当 \(\operatorname{g}(n) = n^2\) 时，\(\operatorname{h}(n) = \sum\limits_{d \mid n} n^2 \cdot \varphi(d) = n^2 \times \sum\limits_{d \mid n}\varphi(n)\)。
想起 \(\sum\limits_{d \mid n} \varphi(n)\) 的意义是 \(\operatorname{1} * \varphi\)，结果为 \(\operatorname{id}\)。
所以 \(\operatorname{h}(n) = n^3\)。
代码
#include <iostream>
#include <map>
typedef long long llt;
int N = 8000000;
int H = 4000000;
int prime[8000005], minprime[8000005], *pm = prime;
int phi[8000005];
llt p, n;
int quick_pow(int _a,int _n) {
int _res = 1;
while(_n) {
if(_n&1)
_res = 1ll*_res*_a%p;
_a = 1ll*_a*_a%p;
_n >>= 1;
}
return _res;
}
void sieve() {
int i, *j, upd;
phi[1] = 1;
for(i = 2;i <= H;++i) {
if(!minprime[i]) {
minprime[i] = *++pm = i;
phi[i] = i-1;
}
if((upd = N/i) >= minprime[i]) {
for(j = prime+1;*j != minprime[i];++j) {
minprime[i*(*j)] = *j;
phi[i*(*j)] = 1ll*phi[i]*phi[*j]%p;
}
minprime[i*(*j)] = *j;
phi[i*(*j)] = 1ll*phi[i]*(*j)%p;
} else {
for(j = prime+1;*j <= upd;++j) {
minprime[i*(*j)] = *j;
phi[i*(*j)] = 1ll*phi[i]*phi[*j]%p;
if(j == pm)
break;
}
}
}
for(i = 2;i <= H;++i)
phi[i] = (phi[i-1]+1ll*i*i%p*phi[i])%p;
for(;i <= N;++i) {
if(!minprime[i])
phi[i] = i-1;
phi[i] = (phi[i-1]+1ll*i*i%p*phi[i])%p;
}
}
int inv_6, inv_2;
llt pre_2(llt x) {
x %= p;
return x*(x+1)%p*(2*x+1)%p*inv_6%p;
}
llt pre_3(llt x) {
x %= p;
llt res = 1ll*x*(x+1)%p*inv_2%p;
return 1ll*res*res%p;
}
std :: map<llt,llt> his;
llt dsv(llt x) {
if(x <= N)
return phi[x];
if(his.find(x) != his.end())
return his[x];
llt res = pre_3(x);
for(llt l = 2, r;l <= x;l = r+1) {
r = x/(x/l);
res = (res-(pre_2(r)-pre_2(l-1))*dsv(x/l))%p;
}
if(res < 0)
res += p;
return his[x] = res;
}
llt ans;
int main() {
scanf("%lld%lld",&p,&n);
if(n < N) {
N = n;
H = N>>1;
}
sieve();
inv_6 = quick_pow(6,p-2);
inv_2 = quick_pow(2,p-2);
for(llt l = 1, r;l <= n;l = r+1) {
r = n/(n/l);
ans = (ans+(dsv(r)-dsv(l-1))%p*pre_3(n/l))%p;
}
if(ans < 0)
ans += p;
printf("%lld\n",ans);
return 0;
}
欧拉反演
实际上莫比乌斯反演应用时最重要的就是 \(\sum\limits_{d \mid n}\mu(d) = \left[n = 1\right]\)，而欧拉反演最重要的是 \(\sum\limits_{d \mid n} \varphi(d) = n\)。
不过显然欧拉反演的应用范围相对要小，毕竟莫比乌斯反演有着 \(\mu * 1 = \epsilon\) 这样优美的性质。（尤其注意到 \(\epsilon\) 是狄利克雷卷积的单位元）
