莫比乌斯反演
数论分块(整除分块)
给定整数 \(n, k\),计算 \(ans = \sum\limits_{i = 1}^ n k \% i\)。
\(ans = kn - \sum\limits_{i = 1}^n k \lfloor \frac{k}{i} \rfloor\),看起来做不了。
实际上至多只有 \(O(\sqrt k)种\lfloor \frac{k}{i} \rfloor\),且每一种都连续,证明如下:
- 对于 \(i \le \sqrt k\),只有 \(\sqrt k\) 种。
- 对于 \(i > sqrt(k)\),\(\lfloor \frac{k}{i} \rfloor \le k\),只有 \(\sqrt k\) 种。
#include <bits/stdc++.h>
using namespace std;
int n, k;
long long ans;
int main() {
ios::sync_with_stdio(0), cin.tie(0);
cin >> n >> k, ans = n * k;
for (int l = 1, r; l <= n; l = r + 1) {
int x = k / l; r = min(n, (x ? k / x : n));
ans -= (l + r) * (r - l + 1ll) * x / 2;
// l ~ r 的 k / i = x
}
cout << ans;
return 0;
}
Mobius 函数
Mobius 函数 \(\mu(n) = \begin{cases} 1, n = 1 \\ (-1)^k, n = p_1p_2\dots p_k \\ 0, otherwise \end{cases}\)。\(p_1, p_2, \dots p_k\) 为质数。
不难发现 Mobius 函数为积性函数。\(\mu(n)\) 恰在 \(n\) 有平方因子时为 \(0\)。
性质
\(\sum\limits_{d | n} \mu(d) = \epsilon(n)\),记作 \(\mu * 1 = \epsilon\)。
证明
设 \(n\) 有 \(k\) 个质因子。\(\mu(d) \ne 0\) 当且仅当 \(d\) 无平方因子,故 \(d\) 的每个质因子指数只能为 \(0/1\)。
\(\sum\limits_{d | n} \mu(d) = \sum\limits_{i = 0}^k (-1)^k C_{k}^i = (1 - 1)^k\)。
当 \(k = 0\),即 \(n = 1\) 时 \(\sum\limits_{d / n} \mu(d) = 1\),否则为 \(0\)。
求法
因为 Mobius 函数为积性函数,所以可以 \(O(n)\) 用线性筛算。
void getmu() {
mu[1] = 1;
for (int i = 2; i < MAXN; i++) {
if (!vis[i]) p[++k] = i, mu[i] = -1;
for (int j = 1; j <= k && p[j] * i < MAXN; j++) {
vis[p[j] * i] = 1;
if (i % p[j] == 0) break;
mu[p[j] * i] = -mu[i];
}
}
}
Mobius 反演
设 \(f, g\) 为数论函数,\(g(n) = \sum\limits_{d \mid n} f(d) \ \ \ \ \ \ (1)\),即 \(g = f * 1\)。
\((1)\) 的充要条件为 \(f(n) = \sum\limits_{d \mid n} g(d)\mu(\frac{n}{d})\)。
证明:\(g = f * 1, f = f * \epsilon = f * (1 * \mu) = f * 1 * \mu = g * \mu\)。
另一种形式
\(F(d) = \sum\limits_{n | d} f(n) \iff f(n) = \sum\limits_{n | d} F(d)\mu(\frac{d}{n})\)。
使用容斥的方式更好理解,\(f(n)\) 为 \(x = n\) 的答案,\(F(n)\) 为 \(x\) 为 \(n\) 的倍数的答案。
以 \(f(1)\) 为例,\(f(1) = F(1)\) 减去 \(x = 2n, 3n, 5n, \dots\) 的答案,即 \(F(1) - \sum\limits_{p \in P} F(p)\)。但 \(F(6 = 2 \times 3), F(10 = 2 \times 5), \dots\) 被多减了,于是要加上 \(\sum\limits_{p1, p2 \in P} F(p1p2)\),但 \(F(2 \times 3 \times 5) \dots\) 又被多加了。由此可以发现规律 \(F(i)\) 的系数为 \(\mu(i)\),所以 \(f(1) = \sum\limits_{i = 1} \mu(i)F(i)\)。其他 \(f(x)\) 同理。
应用
利用 Dirichlet 卷积可以解决一系列求和问题。常见做法是使用 一个 Dirichlet 卷积替换求和式中的一部分,然后调换求和顺序, 最终降低时间复杂度。
经常利用的卷积有 \(\mu ∗ 1 = ϵ\) 和 \(Id = φ ∗ 1\)。有时会利用 \(\mu * Id = \varphi\)(证明:\(\mu * Id = \mu * (1 * \varphi) = \mu * 1 * \varphi = \epsilon * \varphi = \varphi\))。
常用的结论:\([gcd(i, j) == 1] = \sum\limits_{k \mid gcd(i, j)} \mu(k)\)(证明:\(\mu * 1 = \epsilon\))。
题目1:洛谷 P2257
多组数据,给定 \(n, m\),求 \(1 \le x ≤ n, 1 \le y \le m\) 且 \(gcd(x,y)\) 为质数的 \((x,y)\) 有多少对。
令 \(f(T) = \sum\limits_{g \in P, g \mid T} g\mu(\frac{T}{g})\),可 \(O(n)\) 预处理出 \(f(T)\),\(ans = \sum\limits_{T = 1} f(T)\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor\),求出 \(f\) 的前缀和然后整除分块即可。
时间复杂度:\(O(N + \sum \sqrt n)\)。
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int MAXN = 1e7 + 1;
bool vis[MAXN];
ll s[MAXN];
int T, n, m, k, mu[MAXN], p[MAXN];
void build() { // 预处理
mu[1] = 1;
for (int i = 2; i < MAXN; i++) {
if (!vis[i]) p[++k] = i, mu[i] = -1;
for (int j = 1; j <= k && p[j] * i < MAXN; j++) {
vis[p[j] * i] = 1;
if (i % p[j] == 0) break;
mu[p[j] * i] = -mu[i];
}
}
for (int i = 1; i <= k; i++) {
for (int j = 1; p[i] * j < MAXN; j++) {
s[p[i] * j] += mu[j];
}
}
for (int i = 1; i < MAXN; i++) {
s[i] += s[i - 1];
}
}
void Solve() {
cin >> n >> m;
if (n > m) swap(n, m);
ll ans = 0;
for (int l = 1, r; l <= n; l = r + 1) { // 整除分块
int d1 = n / l, r1 = n / d1;
int d2 = m / l, r2 = m / d2;
r = min(r1, r2);
ans += 1ll * d1 * d2 * (s[r] - s[l - 1]);
}
cout << ans << '\n';
}
int main() {
ios::sync_with_stdio(0), cin.tie(0);
for (cin >> T, build(); T--; ) {
Solve();
}
return 0;
}
题目2:洛谷 P1829
稍微加强:洛谷 P3911。
给定 \(n, m\),求 \(\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m \operatorname{lcm}(i, j)\)。
因为 \(\text{lcm(i, j)} = \frac{ij}{\gcd(i, j)}\),所以分开处理分子分母即可。
时间复杂度:\(O(n)\)。
题目3:洛谷 P3312
多组数据。给定 \(n, m, a\),求 \(\sum\limits_{i = 1}^n\sum\limits_{j = 1}^md(\gcd(i, j))[d(\gcd(i, j)) \le a]\),其中 \(d(i)\) 表示 \(d_i\) 的约数和。
在 前置知识 中我们提到 \(d(i)\) 为积性函数,也是可以线性筛求的。具体来说:
- 若 \(i\) 为质数,\(d(i) = i + 1\)。
- 否则设 \(i\) 的最小质因子为 \(p\),\(x = \frac{i}{p}\)。
- 若 \(p \mid x\),即 \(p, x\) 互质,\(d(i) = d(p)d(x)\)。
- 否则,不妨设 \(x = p^k r\),\(d(i) = d(p) + d(r) \cdot p^{k + 1}\)。
- 同时维护 \(p^k, d\) 即可。
// f[i]: i 的约数之和,g[i]: i 的最小质因子 p, 指数 k, g[i] = p^k;
void build() {
mu[1] = 1, f[1] = 1, a[1] = 1;
for (int i = 2; i < MAXN; i++) {
a[i] = i;
if (!vis[i]) {
p[++k] = i, mu[i] = -1;
g[i] = i, f[i] = i + 1;
}
//if (i <= 50) cout << i << ' ' << f[i] << '\n';
for (int j = 1; j <= k && p[j] * i < MAXN; j++) {
int u = p[j] * i; vis[u] = 1;
if (i % p[j] == 0) {
g[u] = g[i] * p[j];
f[u] = (f[i] + 1ll * f[i / g[i]] * g[u]) % Mod;
break;
}
f[u] = 1ll * f[i] * f[p[j]] % Mod;
g[u] = p[j], mu[u] = -mu[i];
}
}
预处理好 \(d(i)\) 后,开始推式子。
令 \(F(T) = \sum\limits_{g \mid T, d(g) \le a} \mu(\frac{T}{g})\),可以将 \(g\) 按 \(d(g)\) 排序,询问按 \(a\) 排序,每次加入 \(g\) 然后更新 \(F\),询问时整除分块即可。
注意,整除分块需要求 \(F\) 前缀和,更新需要树状数组。
时间复杂度:\(O(N \log^2 N + \sum\limits \sqrt n)\)。
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int MAXN = 1e5 + 3;
const ll Mod = (1ll << 31);
struct Query {
int n, m, a, id;
bool operator <(const Query &i) const {
return a < i.a;
}
} q[MAXN];
bool vis[MAXN];
ll tr[MAXN]; // F 的前缀和
int T, k, mu[MAXN], p[MAXN], f[MAXN], g[MAXN], a[MAXN], ans[MAXN];
// f[i]: i 的约数之和,g[i]: i 的最小质因子 p, 指数 k, g[i] = ik;
void build() {
mu[1] = 1, f[1] = 1, a[1] = 1;
for (int i = 2; i < MAXN; i++) {
a[i] = i;
if (!vis[i]) {
p[++k] = i, mu[i] = -1;
g[i] = i, f[i] = i + 1;
}
for (int j = 1; j <= k && p[j] * i < MAXN; j++) {
int u = p[j] * i; vis[u] = 1;
if (i % p[j] == 0) {
g[u] = g[i] * p[j];
f[u] = (f[i] + 1ll * f[i / g[i]] * g[u]) % Mod;
break;
}
f[u] = 1ll * f[i] * f[p[j]] % Mod;
g[u] = p[j], mu[u] = -mu[i];
}
}
sort(a + 1, a + MAXN, [](int x, int y) { // 按约数和排序
return f[x] < f[y];
});
}
inline int lowbit(int x) {
return x & (-x);
}
void add(int x, int v) {
for (; x < MAXN; x += lowbit(x)) {
tr[x] += v;
}
}
ll query(int x) {
ll ret = 0;
for (; x; x -= lowbit(x)) {
ret += tr[x];
}
return ret % Mod;
}
int Solve(int n, int m) {
ll ans = 0;
if (n > m) swap(n, m);
for (int l = 1, r; l <= n; l = r + 1) {
int d1 = n / l, d2 = m / l;
r = min(n / d1, m / d2);
ans += (query(r) - query(l - 1)) * d1 * d2;
}
return ans % Mod;
}
int main() {
ios::sync_with_stdio(0), cin.tie(0);
cin >> T, build();
for (int i = 1; i <= T; q[i].id = i, i++) {
cin >> q[i].n >> q[i].m >> q[i].a;
}
sort(q + 1, q + T + 1);
for (int i = 1, j = 1; i <= T; i++) {
for (; j < MAXN && f[a[j]] <= q[i].a; j++) {
int g = a[j];
for (int x = 1; x * g < MAXN; x++) {
add(x * g, f[g] * mu[x]);
}
}
ans[q[i].id] = Solve(q[i].n, q[i].m);
}
for (int i = 1; i <= T; i++) {
cout << (0ll + ans[i] + Mod) % Mod << '\n';
}
return 0;
}
题目4:洛谷 P3327
多组数据。给定 \(n, m\),求 \(\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m d(ij)\),\(d(i)\) 为 \(i\) 的约数个数。
看到 \(d(ij)\) 不知如何是好,这是一个小结论:\(d(ij) = \sum\limits_{x | i}\sum\limits_{y | j} [\gcd(x, y) == 1]\)。
证明:设 \(ij = p_1^{i_1 + j_1}p_2^{i_2 + j_2}\dots p_k^{i_k + j_k}\),约数个数为 \((i_1 + j_1 + 1)(i_2 + j_2 + 1) \cdots(i_k + j_k + 1)\)。对于每个质因子 \(p_u\) 都是独立算贡献的。
设 \(x\) 的唯一分解中 \(p_u\) 的指数为 \(a \in [0, i_u]\),\(y\) 的唯一分解中 \(p_u\) 的指数为 \(b \in [0, j_u]\),为了使 \(\gcd(x, y) = 1\),则 \(\min(a,b) = 0\),所以共有 \(i_k + j_k + 1\) 种选择方式,刚好对应约数个数中 \(p_u\) 的贡献。
设 \(f(i) = \sum\limits_{x = 1}^i \lfloor \frac{i}{x} \rfloor\),可有整出分块在 \(O(N \sqrt N)\) 的时间内预处理出来。
\(ans = \sum\limits_{k = 1} \mu(k) f(\lfloor \frac{n}{k} \rfloor)f(\lfloor \frac{m}{k} \rfloor)\),然后整除分块即可。
时间复杂度:\(O(\sum \sqrt n + N \sqrt N)\)。
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 5e4 + 1;
bool vis[MAXN];
int T, n, m, k, f[MAXN], p[MAXN], mu[MAXN], s[MAXN];
void build() {
s[1] = mu[1] = 1;
for (int i = 1; i < MAXN; i++) {
for (int l = 1, r; l <= i; l = r + 1) {
int d = i / l; r = i / d;
f[i] += (r - l + 1) * d;
}
}
for (int i = 2; i < MAXN; i++) {
if (!vis[i]) {
p[++k] = i, mu[i] = -1;
}
s[i] = s[i - 1] + mu[i];
for (int j = 1; j <= k && p[j] * i < MAXN; j++) {
vis[p[j] * i] = 1;
if (i % p[j] == 0) break;
mu[p[j] * i] = -mu[i];
}
}
}
void Solve() {
cin >> n >> m;
if (n > m) swap(n, m);
long long ans = 0;
for (int l = 1, r; l <= n; l = r + 1) {
int d1 = n / l, d2 = m / l;
r = min(n / d1, m / d2);
ans += 1ll * (s[r] - s[l - 1]) * f[d1] * f[d2];
}
cout << ans << '\n';
}
int main() {
ios::sync_with_stdio(0), cin.tie(0);
for (cin >> T, build(); T--; ) {
Solve();
}
return 0;
}
浙公网安备 33010602011771号