time: 2021-11-05 21:00:22
tags:
- 笔记 听课笔记 数论/积性函数 数论/莫比乌斯反演
数论函数
定义域为正整数，陪域为实数的函数。
积性函数
定义当 \((a,b)=1\) 时满足 \(f(ab)=f(a)f(b)\) 的函数为积性函数。而对于任意 \(a,b\)，\(f(ab)=f(a)f(b)\) 都成立的函数叫做完全积性函数。
常见的积性函数有
- 恒等函数 \(I(n)=1\)
- 幂函数 \(I_k(n)=n^k\)
- 单位函数 \(id(n)=n\)
- 元函数 \(\varepsilon (n)=\begin{cases}1,&n=1\\0,&n>1\end{cases}\)
- 因子和函数 \(\sigma(n)=\sum_{d\mid n}d\)
- 约数个数函数 \(d(n)=\sum_{d\mid n} 1\)
【例 关于积性】
如设 \(\sigma_k(n)=\sum_{d\mid n}d^k\)，则 \(\sigma_k(n)\) 也是积性函数。考虑设互质的两个数 \(a=\prod_{i=1}^x p_i^{q_i},b=\prod_{i=1}^yh_i^{s_i}\).
那么 \(\sigma_k(a)=\sum_{i=1}^x(1+p_i+\dots+p_i^{q_i})^k\)，\(\sigma_k(b)=\sum_{i=1}^y(1+h_i+\dots+h_i^{s_i})^k\)，所以显然有 \(\sigma_k(ab)=\sigma_k(a)\sigma_k(b)\).
线性筛积性函数
线性筛任何积性函数都要分两种情况考虑， 一是当前数 \(i\bmod p_j\ne0\)，二是 \(i\bmod p_j=0\). \(i\bmod p_j\ne 0\) 的情况意味着 \((i,p_j)=1\)，则 \(f(ip_j)=f(i)f(p_j)\) 即可。\(i\bmod p_j=0\) 的情况意味着 \(p_j\) 是 \(i\) 最小的质因数。设 \(i=p_j^q\times k\)，则我们可以记下 \(q\) 的值，如果能 \(O(1)\) 直接计算 \(f(p_j^{q+1})\)，或者 \(O(1)\) 从 \(f(p_j^q)\) 转移得到 \(f(p_j^{q+1})\) 的值，就可以算出 \(f(ip_j)=f(i)/f(p_j^q)\times f(p_j^{q+1})\) 或者 \(f(ip_j)=f(i/p_j^q)\times f(p_j^{q+1})\)，不管怎么算，都需要记下对于 \(i\) 来说 \(f(p_j^q)\) 的值。
不过对于一些特殊的积性函数，有特殊的计算技巧使得我们不用考虑那么多东西。
筛 \(\varphi\)
对于 \(p\mid i\) 的情况，有 \(\varphi(pi)=pi\prod_{j=1}^q\frac{s_i-1}{s_i}=p\times \varphi(i)\).
什么额外的东西都不用记。
筛 \(\mu\)
对于 \(p\mid i\) 的情况，发现 \(p\) 的幂次大于了 1，因此 \(\mu(pi)=0\).
什么额外的东西都不用记。
筛 \(d\)（约数个数）
对于 \(p\mid i\) 的情况，则 \(d(pi)=d(i)/d(p^q)\times (p^{q+1})=d(i)/(q+1)\times(q+2)\).
需要记下最小质因子的次数 \(q\)。
筛 \(\sigma\)（约数和）
对于 \(p\mid i\) 的情况，\(\sigma(pi)=\sigma(i)/\sigma(p^q)\times \sigma (p^{q+1})=\sigma(i)/\sigma(p^q)\times (\sigma(p^q)+p^{q+1})\).
需要记下 \(q\)，\(p^q\).
代码
什么？你说 \(\sigma\) 会爆 int？这与我有什么关系
#include <bits/stdc++.h>
using namespace std;
const int N = 1e7 + 5;
int flag[N], phi[N], mu[N], d[N], sig[N];
int q[N], pq[N], p[N], v, c;
void sieve(int n) {
flag[1] = true, phi[1] = 1, mu[1] = 1, d[1] = 1, sig[1] = 1;
for(int i = 2; i <= n; i++) {
if(!flag[i]) {
p[++c] = i, phi[i] = i - 1, mu[i] = -1, d[i] = 2, sig[i] = 1 + i;
q[i] = 1, pq[i] = i;
}
for(int j = 1; j <= c && (v = i * p[j]) <= n; j++) {
flag[v] = true;
if(i % p[j]) {
phi[v] = phi[p[j]] * phi[i];
mu[v] = -mu[i]; // mu[v] = mu[p[j]] * mu[i]
d[v] = d[p[j]] * d[i];
sig[v] = sig[p[j]] * sig[i];
q[v] = 1, pq[v] = p[j];
} else {
phi[v] = p[j] * phi[i];
mu[v] = 0;
d[v] = d[i] / (q[i] + 1) * (q[i] + 2);
sig[v] = sig[i] / sig[pq[i]] * (sig[pq[i]] + pq[i] * p[j]);
q[v] = q[i] + 1, pq[v] = pq[i] * p[j];
break;
}
}
}
}
int main() {
sieve(1e7);
for(int i = 1; i <= 20; i++) {
debug(i, phi[i], mu[i], d[i], sig[i], q[i], pq[i]);
}
return 0;
}
狄利克雷卷积
定义两个数论函数 \(f,g\) 的狄利克雷卷积 \(f\ast g=\sum_{d\mid n}f(d)g(\frac nd)\). 狄利克雷卷积具有交换律和结合律。两个积性函数在狄利克雷卷积后得到的函数仍然是积性函数。
\(\varphi\ast I=id\)
即 \(\sum_{d\mid n}\varphi(d)=n\). 考虑两个集合，\(S_1=\{(i,n)|1\le i\le n\}\)，\(S_2=\{(i,p)|1\le i\le p\le n,\gcd(i,p)=1,p\mid n\}\). 将 \(S_1\) 中的元素视作分数 \(\frac in\) 后构成分数集合 \(P_1\). 这些分数约分后一定在 \(S_2\) 构成的分数集合 \(P_2\) 中，所以 \(|S_1|\le|S_2|\).
而将 \(P_2\) 中的分数 \(\frac ip\) 化为 \(\frac {ig}n\)，其中 \(g=\frac np\) 后，由于 \(\frac ip\) 的值两两不同，所以 \(ig\) 也两两不同。而 \(1\le ig\le n\)，所以 \(P_2\) 中的元素也都在 \(P_1\) 中。
综上，\(|S_1|=|S_2|\)，即 \(n=\sum_{d\mid n}\varphi(d)\).
\(\mu\ast I=\varepsilon\)
当 \(n>1\) 时
当 \(n=1\) 时显然有 \(\mu\ast I=1\).
\(\mu\ast id=\varphi\)
第一种证法：
第二种证法：
考虑容斥计算 \(\varphi\)，比如 \(\varphi(60)=60-30-20-12+10+6+4-2=16\). 即有 0 个公共质因数的数的个数等于至少有 0 个的，减去至少有一个的，加上至少有两个的 ...
而右边的部分就是 \(\mu\ast id\).
整除分块（数论分块）
求 \(\sum_{i=1}^n\lfloor\frac ni\rfloor\)，如果直接枚举是 \(O(n)\) 的，使用整除分块的方法就是 \(O(\sqrt n)\).
考虑 \(\lfloor\frac ni\rfloor\) 在 \(i\le \sqrt n\) 的时候至多有 \(\sqrt n\) 种不同的取值，而在 \(i>\sqrt n\) 时，\(\lfloor\frac ni\rfloor< \sqrt n\)，至多有 \(\sqrt n\) 种不同的取值。因此 \(\lfloor\frac ni\rfloor\) 最多有 \(2\sqrt n\) 种不同的取值。
在算法实现时，先枚举 i=1 的情况，此时 \(\lfloor \frac ni\rfloor =n\)，我们试图 \(O(1)\) 算出 \(\lfloor \frac ni\rfloor\) 同样等于 \(n\) 的数的右端点。
\(\lfloor\frac ni\rfloor=k\) 实际上意味着 \(ki+p=n,0\le p<i\)，而如果 \(\lfloor\frac n{i+d}\rfloor=k\)，同样有 \(k(i+d)+p'=n,0\le p<i+d\). 于是有 \(p'=p-kd\)，所以 \(d\le\frac pk\)，\(d_{max}=\lfloor\frac pk\rfloor\).
所以
在算出这个右端点 \(i_{max}\) 后，下一个左端点就等于 \(i_{max}+1\).
莫比乌斯反演定理
【定理 莫比乌斯反演定理 形式 1】
若
则有
【证明】
【定理 莫比乌斯反演定理 形式 2】
这个其实在做题的时候更常用。
若
则
【证明】
![[files/Pasted image 20211110172030.png]]
莫比乌斯反演入门题
P2158 [SDOI2008] 仪仗队
作为体育委员，C 君负责这次运动会仪仗队的训练。仪仗队是由学生组成的 \(N \times N\) 的方阵，为了保证队伍在行进中整齐划一，C 君会跟在仪仗队的左后方，根据其视线所及的学生人数来判断队伍是否整齐（如下图）。
现在，C 君希望你告诉他队伍整齐时能看到的学生人数。
显然答案为 \(2\sum_{i=1}^{n-1}\varphi(i)+2\)，但是可以用一点高妙的方法做。
答案为
将 \([\gcd(i,j)=1]\) 转化为 \(\sum_{d\mid \gcd(i,j)}\mu(d)\) 是个常用技巧，后面也会用到。有些题解说这个就是莫比乌斯反演？？
注意 \(n=1\) 的时候答案为 0。。。
#include <bits/stdc++.h>
using namespace std;
const int N = 40005;
int n;
int flag[N], phi[N], mu[N], d[N], sig[N];
int q[N], pq[N], p[N], v, c;
void sieve(int n) {
flag[1] = true, phi[1] = 1, mu[1] = 1, d[1] = 1, sig[1] = 1;
for(int i = 2; i <= n; i++) {
if(!flag[i]) {
p[++c] = i, phi[i] = i - 1, mu[i] = -1, d[i] = 2, sig[i] = 1 + i;
q[i] = 1, pq[i] = i;
}
for(int j = 1; j <= c && (v = i * p[j]) <= n; j++) {
flag[v] = true;
if(i % p[j]) {
phi[v] = phi[p[j]] * phi[i];
mu[v] = -mu[i]; // mu[v] = mu[p[j]] * mu[i]
d[v] = d[p[j]] * d[i];
sig[v] = sig[p[j]] * sig[i];
q[v] = 1, pq[v] = p[j];
} else {
phi[v] = p[j] * phi[i];
mu[v] = 0;
d[v] = d[i] / (q[i] + 1) * (q[i] + 2);
sig[v] = sig[i] / sig[pq[i]] * (sig[pq[i]] + pq[i] * p[j]);
q[v] = q[i] + 1, pq[v] = pq[i] * p[j];
break;
}
}
}
}
int main() {
cin >> n;
sieve(n);
long long ans = 0;
for(int i = 1; i <= n - 1; i++)
ans += mu[i] * (long long)((n - 1) / i) * ((n - 1) / i);
if(n > 1) ans += 2;
cout << ans << '\n';
return 0;
}
P1829 [国家集训队]Crash的数字表格 / JZPTAB
今天的数学课上，Crash 小朋友学习了最小公倍数（Least Common Multiple）。对于两个正整数 \(a\) 和 \(b\)，\(\text{lcm}(a,b)\) 表示能同时整除 \(a\) 和 \(b\) 的最小正整数。例如，\(\text{lcm}(6, 8) = 24\)。
回到家后，Crash 还在想着课上学的东西，为了研究最小公倍数，他画了一张 \(n \times m\) 的表格。每个格子里写了一个数字，其中第 \(i\) 行第 \(j\) 列的那个格子里写着数为 \(\text{lcm}(i, j)\)。
看着这个表格，Crash 想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当 \(n\) 和 \(m\) 很大时，Crash 就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash 只想知道表格里所有数的和 \(\bmod 20101009\) 的值。
设 \(n\le m\)，答案为
现在将枚举 \(i,j\) 换为枚举 \(d\) 的倍数，得到
\(n/d\) 可以整除分块一下，乘上一个 \(\sqrt n\)，然后 \(n/dl\) 又可以整除分块一下，所以总复杂度是 \(O(n)\).
代码：
记得把减法的地方都加上 MD.
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e7 + 5, MD = 20101009;
int n, m;
int flag[N], mu[N], f[N], S[N], p[N], c;
void init(int n) {
flag[1] = true, mu[1] = 1, f[1] = 1, S[1] = 1;
for(int i = 2, v; i <= n; i++) {
if(!flag[i]) p[++c] = i, mu[i] = -1 + MD;
for(int j = 1; j <= c && (v = i * p[j]) <= n; j++) {
flag[v] = true;
if(i % p[j]) mu[v] = (-mu[i] + MD) % MD;
else {
mu[v] = 0;
break;
}
}
f[i] = (f[i - 1] + ((ll)i * i) % MD * mu[i] % MD) % MD;
S[i] = (S[i - 1] + i) % MD;
}
}
int sum(int x, int y) {
int res = 0;
for(int i = 1, j; i <= x; i = j + 1) {
j = min(x / (x / i), y / (y / i));
res += (ll)(f[j] - f[i - 1] + MD) * S[x / i] % MD * S[y / i] % MD;
res %= MD;
}
return res;
}
int calc(int n, int m) {
int res = 0;
for(int i = 1, j; i <= n; i = j + 1) {
j = min(n / (n / i), m / (m / i));
res += (ll)(S[j] - S[i - 1] + MD) * sum(n / i, m / i) % MD;
res %= MD;
}
return res;
}
int main() {
cin >> n >> m;
if(n > m) swap(n, m);
init(m);
cout << calc(n, m) << '\n';
return 0;
}
BZOJ2694 lcm
注意这一步相当于把 \(di,dj\) 换成了 \(i,j\)，相应的在后面就要乘上 \(d^2\). 然后用上面的技巧 \([\gcd(i,j)=1]=\sum_{l\mid\gcd(i,j)}\mu(l)\) 得到
令 \(T=dl,F(i)=\frac{x(x+1)}2\)，得到
如果此时先枚举 \(T\)，再枚举 \(d\)，那么 \(l\) 就可以直接确定，即 \(\frac Td\). 写出来就是
针对 \(F(A/T)F(B/T)\) 的不同取值，可以整除分块。后面的东西是个积性函数，但是线性筛不好筛，可以枚举 \(d\)，然后对 \(d\) 的倍数计算贡献。复杂度是 \(n(1+\frac 12+\frac 13+\dots +\frac 1n)=O(n\log n)\).
代码：
模数是 \(2^{30}\)，可以直接利用 unsigned int 类型的自动溢出。
#include <bits/stdc++.h>
#define int unsigned int
using namespace std;
const int N = 4e6 + 5;
int T, n, m;
int flag[N], f[N], mu[N], p[N], S[N], v, c;
void sieve1(int n) { // 线性筛 mu
mu[1] = 1;
for(int i = 2; i <= n; i++) {
if(!flag[i]) mu[i] = -1, p[++c] = i;
for(int j = 1; j <= n && (v = p[j] * i) <= n; j++) {
flag[v] = true;
if(i % p[j]) mu[v] = -mu[i];
else {
mu[v] = 0;
break;
}
}
}
}
void sieve2(int n) { // n log n 筛 f
S[1] = 1;
for(int i = 1; i <= n; i++) {
for(int j = i; j <= n; j += i) {
f[j] += i * mu[i] * mu[i] * mu[j / i] * (j / i) * (j / i);
}
S[i] = S[i - 1] + i;
}
for(int i = 1; i <= n; i++) f[i] += f[i - 1]; // 对 f 做前缀和，方便整除分块
}
int calc(int n, int m) {
int res = 0;
for(int i = 1, j; i <= n; i = j + 1) {
j = min(n / (n / i), m / (m / i));
res += (f[j] - f[i - 1]) * S[n / i] * S[m / i];
}
return res;
}
signed main() {
sieve1(4e6);
sieve2(4e6);
for(cin >> T; T; T--) {
cin >> n >> m;
if(n > m) swap(n, m);
cout << calc(n, m) % (1 << 30) << '\n';
}
return 0;
}
P3327 [SDOI2015]约数个数和
设 \(d(x)\) 为 \(x\) 的约数个数，给定 \(n,m\)，求
\[\sum_{i=1}^n\sum_{j=1}^md(ij) \]\(T,n,m\le 50000\)
【结论】
【证明】 设 \(ij\) 的一个质因子 \(p\) 的幂次为 \(c\)，那么 \(d(ij)\) 就是所有 \(c+1\) 相乘，意义是从 \(p^0,p^1,\dots,p^c\) 里面选出来一个。
那么设 \(i\) 的质因子 \(p\) 的幂次为 \(a\)，\(j\) 的质因子的幂次为 \(b\)。把「在 \(i\) 中选择 \(p^x\)」视作「在 \(ij\) 中选择 \(p^x\)」，把「在 \(j\) 中选择 \(p^y\)」视作「在 \(ij\) 中选择 \(p^{a+x}\)」. 这样就可以组合出 \(ij\) 所有的约数，而条件是要求 \(i\) 的约数和 \(j\) 的约数互质。
【结论】
【证明】
设 \(x=Ai+p,\ (0\le p<A)\)，则
由于 \(\frac pA<1\)，所以
原式得证。
现在开始做题了
然后这一步我直接用莫比乌斯函数 \(\sum_{d\mid \gcd(x,y)}\) 替换了 \([\gcd(x,y)=1]\)，弄出来
孩子懵逼了，虽然可以算，但是复杂度不对。
回到式 (1)，其实可以用一个技巧化简，把其中两个求和号搞掉。
又出现了式 (1) 中求和号底下有个整除条件的形式！于是可以确定枚举的东西的上下界，然后把条件丢到要求得东西里去。最后改变一下枚举的东西，就可以变得更简单。
这个技巧一直都在用，只是这里归纳一下。
把枚举 \(x,y\) 改成枚举 \(dx,dy\)，得到
令 \(S(x)=\sum_{i=1}^x\lfloor\frac xi\rfloor\)，考虑「枚举每个数的倍数，做出 1 的贡献」就等于「所有数的约数个数和」，有
所以 \(S(x)\) 可以 \(O(n)\) 线性筛。最后答案为
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 5e4 + 5;
int T, n, m;
int flag[N], p[N], mu[N], d[N], q[N], S[N], Smu[N], c;
void sieve(int n) {
flag[1] = true, mu[1] = 1, d[1] = 1, S[1] = 1, Smu[1] = 1;
for(int i = 2, v; i <= n; i++) {
if(!flag[i]) p[++c] = i, mu[i] = -1, d[i] = 2, q[i] = 1;
for(int j = 1; j <= c && (v = i * p[j]) <= n; j++) {
flag[v] = true;
if(i % p[j]) {
mu[v] = -mu[i];
d[v] = d[i] * d[p[j]];
q[v] = 1;
} else {
mu[v] = 0;
d[v] = d[i] / (q[i] + 1) * (q[i] + 2);
q[v] = q[i] + 1;
break;
}
}
S[i] = S[i - 1] + d[i];
Smu[i] = Smu[i - 1] + mu[i];
}
}
int calc(int n, int m) {
int res = 0;
for(int i = 1, j; i <= n; i = j + 1) {
j = min(n / (n / i), m / (m / i));
res += (Smu[j] - Smu[i - 1]) * S[n / i] * S[m / i];
}
return res;
}
signed main() {
ios::sync_with_stdio(false); cin.tie(nullptr);
sieve(5e4);
for(cin >> T; T; T--) {
cin >> n >> m;
if(n > m) swap(n, m);
cout << calc(n, m) << '\n';
}
return 0;
}