题解 P4450 【双亲数】
fa_555
2019-07-26 09:33:10
**upd：256 天后的更新，修复了错误和公式格式，充实了内容**
**如果 $\LaTeX$ 渲染挂了请移步[洛谷博客](https://www.luogu.com.cn/blog/fa555/solution-p4450)或[个人博客](https://fa555.github.io/2019/%E9%A2%98%E8%A7%A3-P4450-%E3%80%90%E5%8F%8C%E4%BA%B2%E6%95%B0%E3%80%91/)**
under 题解 [P4450](https://www.luogu.org/problem/P4450)
*蒟蒻第一道莫反题*
*本题解适合新手学习基础，dalao 请跳过，不要嘲讽 qwq*
---
### 前置知识：
#### 单位函数（$\epsilon$）：
$$
\epsilon(i) = \begin{cases}
1 &, i = 1 \\
0 &, i \ne 1
\end{cases}
$$
即当且仅当 `i = 1` 时取值取值为 `1`，否则为 `0`。很简单对不对？
#### Dirichlet 卷积($\ast$)：
Dirichlet 卷积是数论函数之间的运算。即
$$
h = f \ast g
$$
其中 $f, g, h$ 均为数论函数。则
$$
h(n) = \sum_{d | n}{f(d) g \left(\frac nd \right)}
$$
也可以写作
$$
h(n) = \sum_{i · j = n}{f(i) g(j)}
$$
我们发现：$\epsilon$ 函数是 Dirichlet 卷积的单位元，即 $f = f \ast \epsilon$。
#### Mobius 函数($\mu$)：
$$
\mu(n) = \begin{cases}
1 &, n = 1 \\
(-1)^m &, n = p_1p_2 \cdots p_m \\
0 &, \mathrm{otherwise}
\end{cases}
$$
其中 $p_1,p_2, \cdots, p_m$ 是不同的质数。$\mu(n)$ 恰在 $n$ 无平方因子时取值非零。显然 $\mu$ 是积性函数。
$$
\begin{aligned}
\sum_{d | n}{\mu(d)} &= \epsilon(n) \\
\epsilon &= 1 \ast \mu
\end{aligned}
$$
#### Mobius 变换
设 $f$ 是数论函数，若数论函数 $g$ 满足
$$
\begin{aligned}
g &= f * \mathbf1 \\
\textrm{That is} \quad g(n) &= \sum_{d | n}{f(d)}
\end{aligned}
$$
则称 $g$ 是 $f$ 的 Mobius 变换。
推一下式子。
$$
g = f \ast 1 \Rightarrow f = f \ast \epsilon = f \ast 1 \ast \mu = g \ast \mu
$$
即
$$
g = f \ast 1 \Leftrightarrow f = g \ast \mu
$$
至此就可以来做这道题了。
---
回到了原题。题目要求的是
$$
\sum_{i = 1}^A{\sum_{j = 1}^B{[\gcd(i, j) = d]}}
$$
对式子进行套路性的变形
$$
\sum_{i = 1}^{\lfloor A / d \rfloor}{\sum_{j = 1}^{\lfloor B / d \rfloor}{[\gcd(i, j) = 1]}}
$$
出现了 $[ \gcd(i, j) = 1]$，也即 $\epsilon( \gcd(i, j))$。把 $\epsilon = 1 \ast \mu$ 代入，得到了
$$
\sum_{i = 1}^{\lfloor A / d \rfloor}{\sum_{j = 1}^{\lfloor B / d \rfloor}{\sum_{g | \gcd(i, j)}{\mu(g)}}}
$$
~~（似乎越来越丑了）~~此时可以把 $d | \gcd(i, j)$ 拆开：
$$
\begin{aligned}
& \sum_{i = 1}^{\lfloor A / d \rfloor}{\sum_{j = 1}^{\lfloor B / d \rfloor}{\sum_{g | i, g | j}{\mu(g)}}} \\
=& \sum_{g = 1}^{\min(A / d, B / d)}{\mu(g) \left\lfloor \frac{A}{dg} \right\rfloor \left\lfloor \frac{B}{dg} \right\rfloor}
\end{aligned}
$$
至此，计算答案可以 $O(n)$ 解决。
算法还可以继续优化。发现式子里还有两个形如 $\lfloor n / d \rfloor$ 的因子，因此可以整除分块优化成 $O(\sqrt{n})$ 的。配合 $O(n)$ 的线性筛 $\mu$ 函数，总复杂度 $O(n)$，可以轻松通过本题。
---
### code 1
``` cpp
std::bitset<1000003> v;
int p[78501], mu[1000003];
void sieve(int N) {
int m = 0;
v.set(1), mu[1] = 1;
for (int i = 2; i <= N; ++i) {
if (!v[i]) p[++m] = i, mu[i] = -1;
for (int j = 1; j <= m && i * p[j] <= N; ++j) {
v.set(i * p[j]);
if (i % p[j] == 0) break;
mu[i * p[j]] = -mu[i];
}
mu[i] += mu[i - 1];
}
}
int main() {
int A, B, d;
long long ans = 0;
std::cin >> A >> B >> d;
A /= d, B /= d;
if (A > B) std::swap(A, B);
sieve(A);
for (int l = 1, r, a, b; l <= A; l = r + 1) {
a = A / l, b = B / l;
r = std::min(A / a, B / b);
ans += 1ll * a * b * (mu[r] - mu[l - 1]);
}
std::cout << ans << '\n';
return 0;
}
```
### code 2
~~如果你像我一样闲得难受，可以学来杜教筛写写这题，复杂度好像还更高了。~~
不会杜教筛但想学的朋友可以来看我的[这篇题解](https://fa555.github.io/2020/杜教筛学习笔记/)。
这个多次询问的题目让我对杜教筛及空间优化有了更深入的理解（
``` cpp
std::bitset<10003> v, vm;
int MX, base, p[1231], mu[10003], mu_[10003];
inline int to(int x) {
return MX / x;
}
void sieve(int N); // 同 code 1
int getMu(int N) {
if (N <= base) return mu[N];
if (vm[to(N)]) return mu_[to(N)];
int ans = 0;
for (int l = 2, r, L; l <= N; l = r + 1) {
L = N / l, r = N / L;
ans += (r - l + 1) * getMu(L);
}
vm.set(to(N));
return mu_[to(N)] = 1 - ans;
}
int main() {
int A, B, d;
long long ans = 0;
std::cin >> A >> B >> d;
A /= d, B /= d;
if (A > B) std::swap(A, B);
MX = A, base = pow(A, 2. / 3);
sieve(base);
for (int l = 1, r, a, b, gR, gL; l <= A; l = r + 1) {
a = A / l, b = B / l;
r = std::min(A / a, B / b);
vm.reset(), gR = getMu(r);
vm.reset(), gL = getMu(l - 1);
ans += 1ll * a * b * (gR - gL);
}
std::cout << ans << '\n';
return 0;
}
```