一点数论函数前缀和喵
Pentiment
·
·
算法·理论
数论函数
数论函数就是定义域为 \mathbb Z_+ 的函数，通常值域为 \mathbb N。
常见的数论函数有：
若数论函数 f 满足，对于 x\perp y，有 f(xy)=f(x)f(y)，则称 f 是积性的。
一些常见的等式：
- $\sum\limits_{d|x}\varphi(x)=x$。对 $x$ 素因子分解，由 $\varphi$ 的积性即得。
- $\sum\limits_{d|x}\mu(x)=\varepsilon(x)$，注意到有重复素因子的数不会产生贡献，故转化成对集合的枚举。设 $k=\omega(x)$，答案为 $\sum\limits_{i=0}^{k}\dbinom ki(-1)^i=(1+(-1))^k=[k=0]$，则当且仅当 $n=1$ 时原式为 $1$，否则为 $0$。
### 整除分块
欲求 $\sum\limits_{i=1}^nf\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)$，需要 $\mathcal O(T(n)\sqrt n)$ 复杂度，其中 $T(n)$ 为计算 $f$ 的前缀和的复杂度。
$\left|\left\{\left.\left\lfloor\dfrac{n}{i}\right\rfloor\right|1\le i\le n\right\}\right|=\mathcal O(\sqrt n)$。分讨 $i\le\sqrt n$ 和 $i>\sqrt n$ 的情况立即得到。
设 $k=\left\lfloor\dfrac{n}{i}\right\rfloor$，则使得 $\left\lfloor\dfrac{n}{j}\right\rfloor=k$ 的最大 $j$ 为 $\left\lfloor\dfrac{n}{k}\right\rfloor$。证明的话，由 $\left\lfloor\dfrac{n}{k}\right\rfloor\ge\left\lfloor\dfrac{n}{\frac ni}\right\rfloor=\lfloor i\rfloor=i$ 即得。
故可以以 $\mathcal O(\sqrt n)$ 的复杂度枚举每一段，对每一段统计 $f$ 的前缀和即可。
#### [P2261 \[CQOI2007\] 余数求和](https://www.luogu.com.cn/problem/P2261)
题意：求：
$$
G(n,k)=\sum_{i=1}^nk\bmod i
$$
数据范围：$1\leq n,k\leq10^9$。
原式即为：
$$
\sum_{i=1}^nk\bmod i=\sum_{i=1}^n\left(k-i\left\lfloor\dfrac ki\right\rfloor\right)=nk-\sum_{i=1}^n i\left\lfloor\dfrac ki\right\rfloor
$$
整除分块即可。
```cpp
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int n, k;
int main() {
scanf("%d%d", &n, &k);
ll ans = (ll)n * k;
for (int l = 1, r; l <= n; l = r + 1) {
r = (k / l ? min(k / (k / l), n) : n);
ans -= (ll)(k / l) * (r - l + 1) * (l + r) / 2;
}
printf("%lld", ans);
}
```
### 欧拉反演
即上文中的 $\sum\limits_{d|x}\varphi(x)=x$。
例如要求 $\sum\limits_{i=1}^n\sum\limits_{j=1}^m\gcd(i,j)$，有：
$$
\begin{aligned}
\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)&=\sum_{i=1}^n\sum_{j=1}^m\sum_{d|\gcd(i,j)}\varphi(d)\\
&=\sum_{i=1}^n\sum_{j=1}^m\sum_{d|i,d|j}\varphi(d)\\
&=\sum_{d=1}^{\min(n,m)}\varphi(d)\left\lfloor\dfrac nd\right\rfloor\left\lfloor\dfrac md\right\rfloor
\end{aligned}
$$
直接整除分块即可。
### 莫比乌斯反演
若 $g(n)=\sum\limits_{d|n}f(d)$，则有 $f(n)=\sum\limits_{d|n}\mu(d)g\left(\dfrac nd\right)$。
证明：
$$
\sum_{d|n}\mu(d)g\left(\dfrac nd\right)=\sum_{d|n}\mu(d)\sum_{d'|\frac nd}f(d')=\sum_{d'|n}f(d')\sum_{d|\frac{n}{d'}}\mu(d)
$$
由于 $\sum\limits_{d|x}\mu(x)=\varepsilon(x)$，故原式当且仅当 $d'=n$ 时第二个和式不为 $0$，故结果为 $f(n)$。
以 $\gcd$ 的二重求和为例：
$$
\sum_{i=1}^n\sum_{j=1}^mf(\gcd(i,j))
$$
令 $f(n)=\sum\limits_{d|n}g(d)$，按照上面的套路，答案即为：
$$
\sum_{d=1}^{\min(n,m)}g(d)\left\lfloor\dfrac nd\right\rfloor\left\lfloor\dfrac md\right\rfloor
$$
例如欧拉反演里举的例子就是 $f=\text{id}$ 的情况，这时 $g=\varphi$。
#### [P4449 于神之怒加强版](https://www.luogu.com.cn/problem/P4449)
题意：求 $\sum\limits_{i=1}^n\sum\limits_{j=1}^m\gcd(i,j)^k$。
答案为：
$$
\sum_{i=1}^{\min(n,m)}\left\lfloor\dfrac ni\right\rfloor\left\lfloor\dfrac mi\right\rfloor\sum_{d|i}d^k\mu\left(\dfrac{i}{d}\right)
$$
考虑线性筛出后面那一坨东西。对于素数的幂，有 $f(p^x)=p^{kx}-p^{k(x-1)}=(p^k-1)p^{k(x-1)}$，那么当 $x>1$ 时有 $f(p^x)=p^xf(p^{x-1})$，当 $x=1$ 时有 $f(p)=p^k-1$。
故有：
$$
f(np)=\left\{
\begin{aligned}
&p^kf(n)&p\mid n\\
&(p^k-1)f(n)&p\nmid n
\end{aligned}
\right.
$$
线性筛后整除分块即可，时间复杂度为 $\mathcal O(n+T\sqrt n)$。
```cpp
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5e6 + 5, p = 1e9 + 7;
int T, n, m, k, tot, pri[N], f[N];
bool vis[N];
int qpow(int a, int b) {
int c = 1;
while (b) {
if (b & 1) c = (ll)c * a % p;
a = (ll)a * a % p, b >>= 1;
}
return c;
}
void sieve(int n) {
f[1] = 1;
for (int i = 2, x; i <= n; i++) {
if (!vis[i]) pri[++tot] = i, f[i] = qpow(i, k) - 1;
for (int j = 1; j <= tot && i * pri[j] <= n; j++) {
vis[x = i * pri[j]] = 1;
if (!(i % pri[j])) { f[x] = (ll)f[i] * (f[pri[j]] + 1) % p; break; }
f[x] = (ll)f[i] * f[pri[j]] % p;
}
}
for (int i = 1; i <= n; i++) f[i] = (f[i - 1] + f[i]) % p;
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> T >> k, sieve(N - 1);
while (T--) {
cin >> n >> m;
int u = min(n, m), sum = 0;
for (int l = 1, r; l <= u; l = r + 1) r = min(n / (n / l), m / (m / l)), sum = (sum + (ll)(n / l) * (m / l) % p * (f[r] - f[l - 1]) % p + p) % p;
cout << sum << '\n';
}
}
```
这里没有 Tricks，但是有一个警示后人：后面要用到前面的 $f$ 值时就不要边筛边做前缀和。
#### [P3312 [SDOI2014] 数表](https://www.luogu.com.cn/problem/P3312)
题意：定义 $f(i,j)=\sigma(\gcd(i,j))$，求：
$$
\sum_{i=1}^n\sum_{j=1}^m[f(i,j)\le a]f(i,j)
$$
原式即为：
$$
\sum_{i=1}^{\min(n,m)}\left\lfloor\dfrac ni\right\rfloor\left\lfloor\dfrac mi\right\rfloor\sum_{d|i}[\sigma(d)\le a]\sigma(d)\mu\left(\dfrac id\right)
$$
设后面那一坨是 $g(i)$，将询问按 $a$ 扫描线，每一次加入新的 $\sigma(d)$，会对其所有倍数处的 $g$ 产生贡献，查询前缀和，使用树状数组。
时间复杂度为 $\mathcal O(q\sqrt n\log n+n\log^2 n)$。
```cpp
#include <bits/stdc++.h>
using namespace std;
#define int unsigned
typedef long long ll;
const int N = 100005, M = 20005, qwq = 1e5;
int Q, tot, pri[N], mu[N], t[N], ans[N];
bool vis[N];
struct query { int id, n, m, a; } q[M];
bool cmp(query a, query b) { return a.a < b.a; }
pair<ll, int> p[N];
void sieve(int n) {
mu[1] = 1;
for (int i = 2, k; i <= n; i++) {
if (!vis[i]) pri[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot && i * pri[j] <= n; j++) {
vis[k = i * pri[j]] = 1;
if (!(i % pri[j])) { mu[k] = 0; break; }
mu[k] = -mu[i];
}
}
for (int i = 1; i <= n; i++) {
p[i].second = i;
for (int j = i; j <= n; j += i) p[j].first += i;
}
sort(p + 1, p + 1 + n);
}
void upd(int x, int d) { while (x <= qwq) t[x] += d, x += x & -x; }
int qry(int x) { int y = 0; while (x) y += t[x], x -= x & -x; return y; }
signed main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> Q;
for (int i = 1; i <= Q; i++) cin >> q[i].n >> q[i].m >> q[i].a, q[i].id = i;
sort(q + 1, q + 1 + Q, cmp);
sieve(qwq);
for (int i = 1, j = 1; i <= Q; i++) {
while (p[j].second && p[j].first <= q[i].a) {
int x = p[j].second;
for (int k = 1; x * k <= qwq; k++) upd(x * k, (int)p[j].first * mu[k]);
j++;
}
int n = q[i].n, m = q[i].m, u = min(n, m);
for (int l = 1, r; l <= u; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans[q[i].id] += (n / l) * (m / l) * (qry(r) - qry(l - 1));
}
}
for (int i = 1; i <= Q; i++) cout << (ans[i] & INT_MAX) << '\n';
}
```
#### [VLATTICE - Visible Lattice Points](https://www.luogu.com.cn/problem/SP7001)
题意：三维空间中有 $(n+1)^3$ 个点，坐标从 $(0,0,0)$ 到 $(n,n,n)$，问有多少个点可以从 $(0,0,0)$ 看到而不受遮挡。
分点的坐标讨论，答案为：
$$
3+3\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=1]+\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^n[\gcd(i,j,k)=1]
$$
考虑如何求解下式：
$$
\sum_{i_1=1}^n\sum_{i_2=1}^n\dots\sum_{i_k=1}^n[\gcd(i_1,i_2,\dots,i_k)=1]
$$
由于 $\sum\limits_{d|n}\mu(d)=[n=1]$，有：
$$
\sum_{i_1=1}^n\sum_{i_2=1}^n\dots\sum_{i_k=1}^n\sum_{d|i_1,\dots,d|i_k}\mu(d)
$$
交换求和符号：
$$
\sum_{d=1}^n\mu(d)\left\lfloor\dfrac{n}{d}\right\rfloor^k
$$
使用整除分块计算即可，时间复杂度为 $\mathcal O(T\sqrt n)$。
```cpp
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e6 + 5;
int T, n, tot, pri[N], mu[N];
bool vis[N];
void sieve(int n) {
mu[1] = 1;
for (int i = 2, k; i <= n; i++) {
if (!vis[i]) pri[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot && i * pri[j] <= n; j++) {
vis[k = i * pri[j]] = 1;
if (!(i % pri[j])) { mu[k] = 0; break; }
mu[k] = -mu[i];
}
}
for (int i = 1; i <= n; i++) mu[i] += mu[i - 1];
}
ll sm2(int n) {
ll sum = 0;
for (ll l = 1, r, x; l <= n; l = r + 1) {
r = n / (x = n / l);
sum += (mu[r] - mu[l - 1]) * x * x;
}
return sum;
}
ll sm3(int n) {
ll sum = 0;
for (ll l = 1, r, x; l <= n; l = r + 1) {
r = n / (x = n / l);
sum += (mu[r] - mu[l - 1]) * x * x * x;
}
return sum;
}
int main() {
sieve(1e6);
cin >> T;
while (T--) {
cin >> n;
cout << 3 + 3 * sm2(n) + sm3(n) << '\n';
}
}
```
### 狄利克雷卷积 & DGF
对于数论函数 $f,g$，定义狄利克雷卷积 $f*g$ 为：
$$
(f*g)(x)=\sum_{d\mid x}f(x)g\left(\dfrac nx\right)
$$
狄利克雷卷积显然满足交换律，结合律和对 $+$ 的分配律的证明也是容易的。
一些例子是 $1*1=\tau$，$1*\text{id}_k=\sigma_k$，$1*\varphi=\text{id}$，$1*\mu=\varepsilon$。
若 $f*g=h$，则记 $f=h/g$，下文中 $/$ 均指狄利克雷除法。
类比多项式乘法和 OGF，考虑对狄利克雷卷积构造一种生成函数。
由于这里是下标乘法，故生成函数的 $x$ 须在指数上，故有 DGF 的定义：
$$
\text{DGF}(f)(x)=\sum_{i=1}^{\infty}\dfrac{f_i}{i^x}
$$
容易验证两个数论函数做狄利克雷卷积等价于它们的 DGF 相乘。
例如 $\text{DGF}(1)=\zeta$，其中 $\zeta$ 是黎曼 zeta 函数，定义为：
$$
\zeta(x)=\sum_{i=1}^\infty\dfrac{1}{i^x}
$$
更多的例子：
- $\text{DGF}(1)(x)=\zeta(x)$；
- $\text{DGF}(\text{id}_k)(x)=\sum\limits_{i=1}^\infty\dfrac{x^k}{i^x}=\sum\limits_{i=1}^\infty\dfrac{1}{i^{x-k}}=\zeta(x-k)$；
- $\text{DGF}(\tau)(x)=\text{DGF}(1*1)(x)=\zeta^2(x)$；
- $\text{DGF}(\sigma_k)(x)=\text{DGF}(1*\text{id}_k)(x)=\zeta(x)\zeta(x-k)$；
- $\text{DGF}(\varepsilon)(x)=1$；
- $\text{DGF}(\varphi)(x)=\text{DGF}(\text{id}/1)(x)=\dfrac{\zeta(x-1)}{\zeta(x)}$；
- $\text{DGF}(\mu)(x)=\text{DGF}(\varepsilon/1)(x)=\dfrac{1}{\zeta(x)}$。
用 DGF 可以推出各种奇奇怪怪的式子，譬如莫比乌斯反演的本质就是把 $f*1=g$ 左边的 $1$ 除了过去变成 $f=g*\mu$。
### 杜教筛
现在欲快速求出 $\sum\limits_{i=1}^nf(i)$，一种方式是杜教筛，要求存在一数论函数 $g$，使得 $g$ 和 $f*g$ 的前缀和可以 $\mathcal O(1)$ 计算，如 $1*\varphi=\text{id}$。
推式子的过程比较 simple，设 $h=f*g$，先列出 $h$ 的前缀和的式子：
$$
\begin{aligned}
\sum_{i=1}^nh(i)&=\sum_{i=1}^n\sum_{d|i}g(d)f\left(\dfrac id\right)\\
&=\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}f(i)\\
&=g(1)\sum_{i=1}^nf(i)+\sum_{d=2}^ng(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}f(i)
\end{aligned}
$$
通常来说 $g(1)=1$，移项得：
$$
\sum_{i=1}^nf(i)=\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}f(i)
$$
最震撼人心的一步在于，只要我们线性筛出 $\mathcal O(n^{\frac23})$ 以内的 $f$ 的前缀和，并且加上记忆化，那么我们就可以直接递归求解 $\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}f(i)$，复杂度为 $\mathcal O(n^{\frac23})$。如何证明？
首先有引理：设 $S(n)=\left|\left\{\left.\left\lfloor\dfrac{n}{i}\right\rfloor\right|1\le i\le n\right\}\right|$，则 $\forall m\in S(n),S(m)\subseteq S(n)$。这一点等价于证 $\left\lfloor\dfrac{\lfloor\frac ab\rfloor}{c}\right\rfloor=\left\lfloor\dfrac{a}{bc}\right\rfloor$，设 $\dfrac ab-\left\lfloor\dfrac{a}{b}\right\rfloor=r$，由于 $\dfrac rc<1$ 即得。
故由于记忆化，每个 $\left\lfloor\dfrac{n}{i}\right\rfloor$ 只会被求一次。设预处理出 $m$ 以内的 $f$ 的前缀和，其中 $m>\sqrt n$，则有：
$$
\begin{aligned}
T(n)&=\mathcal O(m)+\sum_{i=1}^{\lfloor\frac nm\rfloor}\mathcal O\left(\sqrt{\dfrac ni}\right)\\
&=\mathcal O\left(m+\int_0^{\frac nm}\sqrt{\dfrac nx}\text{d}x\right)\\
&=\mathcal O\left(m+\dfrac{n}{\sqrt m}\right)
\end{aligned}
$$
由基本不等式，当 $m=\mathcal O(n^\frac{2}{3})$ 时有最优复杂度 $\mathcal O(n^\frac{2}{3})$。
例如要求 $\varphi$ 或 $\mu$ 的前缀和，就将 $1*\varphi=\text{id}$ 和 $1*\mu=\varepsilon$ 代进去就行了。
#### [P4318 完全平方数](https://www.luogu.com.cn/problem/P4318)
求第 $k$ 个不是完全平方数的倍数的数。
不是完全平方数等价于 $\mu(x)\not=0$，但这个不好求和，这里有一个 Trick，注意到 $\mu(x)=\{-1,0,1\}$，那么有 $\mu^2(x)=[\mu(x)\not=1]$。
故要求使得 $\sum\limits_{i=1}^n\mu^2(i)\ge k$ 的最小的 $n$，二分。
考虑找一个函数 $f$ 和 $\mu^2$ 贴贴，我们希望 $\mu^2*f=1$，比如平方判定函数 $f(x)=[x=m^2,m\in\mathbb Z_+]$，当且仅当 $d$ 取 $x$ 的最大平方因子时 $f(d)\mu^2\left(\dfrac nd\right)=1$，否则为 $0$，证明比较容易，分讨 $f,\mu^2$ 的取值情况即可。
根据杜教筛有：
$$
\sum_{i=1}^n\mu^2(i)=n-\sum_{d=2}^nf(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}\mu^2(i)
$$
**由于 $f$ 是判定函数，仅枚举平方数**：
$$
\sum_{i=1}^n\mu^2(i)=n-\sum_{d=2}^{\lfloor\sqrt n\rfloor}\sum_{i=1}^{\lfloor\frac{n}{d^2}\rfloor}\mu^2(i)
$$
直接做即可，时间复杂度为 $\mathcal O(Tn^{\frac23}\log n)$，可以过。
```cpp
#include <bits/stdc++.h>
#include <bits/extc++.h>
using namespace std;
using namespace __gnu_pbds;
typedef long long ll;
const int N = 1e6 + 5;
int T, k, tot, mu2[N], pri[N];
bool vis[N];
gp_hash_table<int, int> mp;
void sieve(int n) {
mu2[1] = 1;
for (int i = 2, k; i <= n; i++) {
if (!vis[i]) pri[++tot] = i, mu2[i] = 1;
for (int j = 1; j <= tot && i * pri[j] <= n; j++) {
vis[k = i * pri[j]] = 1;
if (!(i % pri[j])) { mu2[k] = 0; break; }
mu2[k] = mu2[i];
}
}
for (int i = 1; i <= n; i++) mu2[i] += mu2[i - 1];
}
int djs(int n) {
if (n < N) return mu2[n];
if (mp[n]) return mp[n];
int sum = n;
for (int d = 2; d * d <= n; d++) sum -= djs(n / (d * d));
return mp[n] = sum;
}
void solve() {
cin >> k;
ll L = 1, R = 1.7e9, mid, ans;
while (L <= R) {
mid = (L + R) >> 1;
if (djs(mid) >= k) ans = mid, R = mid - 1;
else L = mid + 1;
}
cout << ans << '\n';
}
int main() {
sieve(N - 1);
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> T;
while (T--) solve();
}
```
Trick：值域为 $\{-1,0,1\}$ 的函数 $f$，$f^2(x)=[f(x)\not=0]$；对于判定函数 $f$，可以仅枚举满足 $f(x)=1$ 的数以优化复杂度。
#### [AT_xmascon19_d Sum of (-1)^f(n)](https://www.luogu.com.cn/problem/AT_xmascon19_d)
题意：定义 $\lambda(n)=(-1)^{\Omega(n)}$，求 $\sum\limits_{i=1}^n\lambda(i)$。
刘维尔函数 $\lambda$ 是完全积性的，考虑构造 $\lambda*g=h$。尝试令 $g=1$，则：
$$
h(p^k)=\sum_{i=0}^k(-1)^i=[2|k]
$$
由于 $\lambda,1$ 都是积性的，则 $h$ 也是积性的，实际上有 $h(n)=[n=m^2,m\in\mathbb Z_+]$，即平方和判定函数。这个函数在 P4318 中也有出现。
得到结论 $\lambda*1=h$。事实上有 $h$ 的 DGF 为 $\zeta(2x)$，那么 $\lambda$ 的 DGF 就是 $\dfrac{\zeta(2x)}{\zeta(x)}$，结合 P4318 的结论可以得出 $\mu^2$ 的 DGF 是 $\dfrac{\zeta(x)}{\zeta(2x)}$，故有 $\lambda*\mu^2=\varepsilon$。
$1$ 和 $h$ 都方便求前缀和，直接使用杜教筛即可。
```cpp
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 3e7 + 5;
int tot, pri[N], lv[N];
bool vis[N];
ll n;
unordered_map<ll, ll> mp;
ll djs(ll n) {
if (n < N) return lv[n];
if (mp[n]) return mp[n];
ll sum = sqrtl(n);
for (ll l = 2, r, k; l <= n; l = r + 1) r = n / (k = n / l), sum -= djs(k) * (r - l + 1);
return mp[n] = sum;
}
void sieve(int n) {
lv[1] = 1;
for (int i = 2; i <= n; i++) {
if (!vis[i]) pri[++tot] = i, lv[i] = -1;
for (int j = 1; j <= tot && i * pri[j] <= n; j++) {
vis[i * pri[j]] = 1, lv[i * pri[j]] = -lv[i];
if (!(i % pri[j])) break;
}
lv[i] += lv[i - 1];
}
}
int main() {
sieve(N - 1);
cin >> n, cout << djs(n);
}
```
Trick：将刘维尔函数和平方和判定函数、$\mu^2$ 等一起考虑。
---
附：[这里](https://www.luogu.com.cn/discuss/1034959)提到的性质，杜教筛求的是 $f$ 的块筛（一次杜教筛筛出了所有 $x=\left\lfloor\dfrac{n}{i}\right\rfloor$ 处的点值），故若后面只用到了块筛值，则只有一开始预处理时有 $\mathcal O(n^\frac{2}{3})$ 的复杂度。
### 贝尔级数
对于**积性函数** $f$，我们一般只需要关心 $f$ 在 $p^k$ 处的取值，考虑构造关于 $p^k$ 的级数。
定义积性函数 $f$，其在模 $p$ 意义下的贝尔级数为：
$$
F_p(x)=\sum_{i=0}^\infty x^if(p^i)
$$
容易验证，在模 $p$ 意义下，两个数论函数做狄利克雷卷积等价于它们的贝尔级数相乘。
贝尔级数相对于 DGF 的优点在于，前者是多项式的形式，故可以用研究普通生成函数的方式研究贝尔级数。
一些常见数论函数的贝尔级数：
- $\varepsilon\Rightarrow \text{E}_p(x)=1$；
- $1\Rightarrow \text{I}_p(x)=\dfrac{1}{1-x}$；
- $\text{id}_k\Rightarrow (\text{Id}_k)_p(x)=\dfrac{1}{1-p^kx}$；
- $\tau\Rightarrow\text{T}_p(x)=\dfrac{1}{(1-x)^2}$；
- $\sigma_k\Rightarrow\Sigma_k(x)=\dfrac{1}{(1-x)(1-p^kx)}$；
- $\varphi\Rightarrow (\Phi_k)_p(x)=\dfrac{1-x}{1-px}$；
- $\mu\Rightarrow \text{M}_p(x)=1-x$；
- $\mu^2\Rightarrow(\text{M}^2)_p(x)=1+x$；
- $\lambda\Rightarrow\Lambda_p(x)=\dfrac{1}{1+x}$；
- $k$ 次方判定函数 $\text{chk}_k\Rightarrow(\text{Chk}_k)_p(x)=\dfrac{1}{1-x^k}$。
可以类比 DGF，使用贝尔级数研究数论函数之间做卷积的结果，或者当 $k$ 上界不大时可以直接使用多项式乘法快速得到两个数论函数卷积后在素数幂处的取值。
定义数论函数点积 $f\cdot g$ 为：
$$
(f\cdot g)(x)=f(x)g(x)
$$
对于数论函数 $f$，其贝尔级数为 $\text{F}_p(x)$，则 $f\cdot\varepsilon$ 的贝尔级数为 $\text{E}_p(x)$，$f\cdot\text{id}_k$ 的贝尔级数为 $\text{F}_p(p^kx)$。
### Powerful Number 筛
如果正整数 $n$ 的素因子分解式中每个素因子的次数都大于 $1$，则称 $n$ 是 Powerful Number。
每个 Powerful Number 都可以被写成 $a^2b^3$ 的形式。证明：对于 $p^k$，若 $2|k$ 则将 $a$ 乘上 $p^{\frac k2}$，否则将 $a$ 乘上 $p^{\frac{k-3}{2}}$，将 $b$ 乘上 $p$。
定理：$n$ 以内的 Powerful Number 有 $\mathcal O(\sqrt n)$ 个。
证明：枚举 $a$：
$$
\mathcal O\left(\sum_{i=1}^{\sqrt n}\sqrt[3]{\dfrac{n}{x^2}}\right)=\mathcal O\left(\int_1^{\sqrt n}\sqrt[3]{\dfrac{n}{x^2}}\text{d}x\right)=\mathcal O(\sqrt n)
$$
找 Powerful Number 的过程直接暴力搜索出各素数的指数即可，时间复杂度为 $\mathcal O(\sqrt n)$。
顺带一提，Powerful Number 的补集叫 Square Free Number，就是 P4318 中计数的那个东西。
Powerful Number 筛用于求出一些积性函数的前缀和，如果要求 $f$ 的前缀和，需要找到一个积性函数 $g$，使得其在素数处的取值和 $f$ 相等，并且能快速求出 $g$ 的块筛。
令 $h=f/g$，由于 $f,g$ 都是积性函数，则 $h$ 也是积性函数，则 $h(1)=1$。
对于素数 $p$，$f(p)=g(1)h(p)+g(p)h(1)=g(p)+h(p)$，由 $f(p)=g(p)$ 推出 $h(p)=0$。由 $h$ 是积性函数推出 $h(x)$ 仅当 $x$ 为 Powerful Number 时不为 $0$，这样的 $x$ 仅有 $\mathcal O(\sqrt n)$ 个。
由于 $\sum\limits_{i=1}^nf(i)=\sum\limits_{d=1}^nh(d)\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}g(i)$，有值的 $h(d)$ 仅有 $\mathcal O(\sqrt n)$ 个直接枚举，如果 $g$ 的前缀和可以直接 $\mathcal O(1)$ 求那复杂度就是 $\mathcal O(\sqrt n)$ 的，否则需要求出 $g$ 的块筛，比如用杜教筛就是 $\mathcal O(n^{\frac23})$ 的。
#### [P5325 【模板】Min_25 筛](https://www.luogu.com.cn/problem/P5325)
~~【模板】Powerful Number 筛~~
其实 PN 筛和 Min_25 筛对 $f$ 的要求差不多是一样的，所以很多时候都可以在 Min_25 的题里看到 PN 筛的身影。~~但是跑不过 Min25 啊 qwq~~
这个题基本上直接把 $g$ 给我们了，不难看出有 $g(n)=n\varphi(n)$。
为了求出 $g$ 的块筛，使用杜教筛，$g$ 的 DGF 为 $\dfrac{\zeta(x-1)}{\zeta(x-2)}$，故有 $g*\text{id}_2=\text{id}$，直接套板子即可。
现在要做的就是快速求出 $h$ 在 Powerful Number 处的点值。由于 $h$ 是积性函数，只需要考虑 $h(p^k)$ 的值即可。
由于：
$$
f(p^k)=\sum_{i=0}^kh(p^i)g(p^{k-i})
$$
移项：
$$
h(p^k)=f(p^k)-\sum_{i=0}^{k-1}h(p^i)g(p^{k-i})
$$
这部分直接计算复杂度大概是 $\mathcal O(\sqrt n)$ 级别的，不太会证。
求出 $h(p^k)$ 后 dfs 求 Powerful Number 的时候利用积性直接计算即可。
时间复杂度瓶颈在求出 $g$ 的块筛，复杂度为杜教筛的 $\mathcal O(n^{\frac23})$。
顺便学了一下杜教筛不用 `unordered_map` 的写法，注意到要用到它当且仅当当前的 $x>n^{\frac 23}$，这时候开另外一个数组，存到 $\left\lfloor\dfrac nx\right\rfloor$ 下标即可，避免了 `unordered_map` 的大常数。~~但还是跑不过 Min25 啊 qwq~~
小心爆 `long long`。
```cpp
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5e6 + 5, p = 1e9 + 7, inv2 = (p + 1) / 2, inv6 = (p + 1) / 6;
int ans, tot, pri[N], g[N], G[N];
bool vis[N];
ll n;
void sieve(int n) {
g[1] = 1;
for (int i = 2, k; i <= n; i++) {
if (!vis[i]) pri[++tot] = i, g[i] = i - 1;
for (int j = 1; j <= tot && i * pri[j] <= n; j++) {
vis[k = i * pri[j]] = 1;
if (!(i % pri[j])) { g[k] = g[i] * pri[j]; break; }
g[k] = g[i] * (pri[j] - 1);
}
g[i] = (g[i - 1] + (ll)i * g[i]) % p;
}
}
int djs(ll x) {
if (x < N) return g[x];
if (G[n / x]) return G[n / x];
int X = x % p, sum = (ll)X * (X + 1) % p * (X * 2 + 1) % p * inv6 % p;
for (ll l = 2, r; l <= x; l = r + 1) {
r = x / (x / l);
sum = (sum - (l + r) % p * (r - l + 1) % p * inv2 % p * djs(x / l) % p + p) % p;
}
return G[n / x] = sum;
}
void pns(int x, ll m, int h) {
if (x > tot || m > n / ((ll)pri[x] * pri[x])) {
ll k = n / m;
if (k < N) ans = (ans + (ll)h * g[k]) % p;
else ans = (ans + (ll)h * G[n / k]) % p;
return;
}
pns(x + 1, m, h);
ll k = (ll)pri[x] * pri[x];
for (int i = 2; m * k <= n; i++, k *= pri[x]) pns(x + 1, m * k, h * (k % p) % p * (pri[x] - 1) % p * (i - 1) % p);
}
int main() {
cin >> n;
sieve(N - 1), djs(n), pns(1, 1, 1);
cout << ans;
}
```
#### [DIVCNT3 - Counting Divisors (cube)](https://www.luogu.com.cn/problem/SP20174)
题意：求 $\sum\limits_{i=1}^n\tau(i^3)$。
设积性函数 $f(p^k)=\tau(p^{3k})=3k+1$，现在要求 $f$ 的前缀和。使用 Powerful Number 筛，$f$ 在素数处的取值为 $4$，构造 $f=g*h$。为了好筛，令 $g=\tau*\tau$，其 DGF 为 $\zeta^2(x)$，设法求出它的块筛：
$$
\sum_{i=1}^ng(i)=\sum_{d=1}^n\tau(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}\tau(i)
$$
先求出 $\tau$ 的块筛，然后整除分块，需要 $g$ 在 $n^{\frac23}$ 内的线性筛才能保证复杂度。
有 $g(p^k)=\sum\limits_{i=0}^k(i+1)(k-i+1)$，线筛的时候需要把每个 $i$ 的最小素因子除干净，复杂度仍然是线性的，证明参见 [here](https://www.luogu.com/article/ja3e33rb)。（sto syp11 orz）
这样可以在 $\mathcal O(n^{\frac{2}{3}})$ 的时间复杂度内求出 $g$ 的块筛。要求出 $h$ 在 Powerful Number 处的点值，直接使用 PN 筛模板里的做法即可。
由于有多组数据，需要调整阈值。当 $B=\mathcal O((\sum n)^{\frac 23})$ 时有最优复杂度 $O((\sum n)^{\frac 23})$。实际上 SPOJ 上的题面保证了 $\sum n\le2\times10^{11}$，可以通过。
代码咕咕咕。
### Meissel-Lehmer
Meissel-Lehmer 是一种用于快速计算 $\pi(n)$ 的算法，其中 $\pi(n)$ 表示 $n$ 以内的素数个数。
先线性筛出 $\sqrt n$ 内的素数。然后阈值分治，算出 $>n^{\frac13}$ 的数中不会被 $\le n^{\frac13}$ 的素数筛掉的数，设 $dp_{i,j}$ 为 $i$ 以内不会被 $j$ 以内的素数筛掉的数。没有被筛掉的数如果是合数，则一定是两个 $>n^{\frac13}$ 的素数之积，直接枚举小的素数去计算大素数的方案数即可。有：
$$
\pi(n)=\pi(n^{\frac13})+dp_{n,\pi(n^{\frac13})}-1-\sum_{p=n^{\frac13}+1}^{\sqrt n}\left(\pi\left(\dfrac np\right)-\pi(p)+1\right)
$$
对于 dp 数组的转移，第 $j$ 个素数会筛掉所有没有被前面的素数筛掉的 $p_j$ 的倍数，故 $dp_{i,j}=dp_{i,j-1}-dp_{\frac{i}{p_j},j-1}$。
dp 的边界条件为 $dp_{i,0}=i$ 以及 $dp_{i,j}=\max(\pi(i)-j+1,0)$（当 $p_j^2\ge i$ 时）。
递推出一部分 dp 的值，将沿途求出的 $\pi$ 值存进桶里，这样可以在 $\mathcal O(n^{\frac23})$ 的复杂度内求出 $\pi$ 的块筛。
#### [AT_xmascon19_e Sum of f(n)](https://www.luogu.com.cn/problem/AT_xmascon19_e)
题意：求 $\sum\limits_{i=1}^n\Omega(i)$。
对每个素数分别统计答案：
$$
\sum_{p\in\mathbb P}^n\sum_{p^k\le n}\left\lfloor\dfrac{n}{p^k}\right\rfloor
$$
容易发现对于 $>\sqrt n$ 的素数其只会产生一次贡献。则对于 $k>1$，直接枚举范围内的素数即可；对于 $k=1$，要求：
$$
\sum_{p\in\mathbb P}^n\left\lfloor\dfrac{n}{p}\right\rfloor
$$
整除分块，使用带记忆化的 Meissel-Lehmer 预处理出 $\pi$ 的块筛即可。
```cpp
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e6 + 5, I = 1.8e6, J = 60;
ll n, P[N], ans;
int f[I + 5][J + 5], tot, pri[80005], pi[N];
bool vis[N];
void init() {
for (int i = 2; i < N; i++) {
if (!vis[i]) pri[++tot] = i;
for (int j = 1; j <= tot && i * pri[j] < N; j++) {
vis[i * pri[j]] = 1;
if (!(i % pri[j])) break;
}
}
for (int i = 1; i < N; i++) pi[i] = pi[i - 1] + !vis[i];
for (int i = 1; i <= I; i++) f[i][0] = i;
for (int i = 1; i <= I; i++) {
for (int j = 1; j <= J; j++) f[i][j] = f[i][j - 1] - f[i / pri[j]][j - 1];
}
}
ll dp(ll i, ll j) {
if (i <= I && j <= J) return f[i][j];
if (!i) return 0; if (!j) return i;
if (i < N && (ll)pri[j] * pri[j] >= i) return max(0ll, pi[i] - j);
return dp(i, j - 1) - dp(i / pri[j], j - 1);
}
ll Pi(ll x) {
if (!x) return 0;
if (x < N) return pi[x] - 1;
if (P[n / x]) return P[n / x];
ll m = pi[(int)cbrt(x)] - 1, ans = dp(x, m) + m - 1; m++;
while ((ll)pri[m] * pri[m] <= x) ans -= Pi(x / pri[m]) - m + 1, m++;
return P[n / x] = ans;
}
signed main() {
cin >> n, init(), Pi(n);
for (int i = 1; (ll)pri[i] * pri[i] <= n; i++) {
ll m = (ll)pri[i] * pri[i];
while (m <= n) ans += n / m, m *= pri[i];
}
for (ll l = 1, r, k; l <= n; l = r + 1) {
r = n / (k = n / l);
ans += k * (Pi(r) - Pi(l - 1));
}
cout << ans;
}
```
#### [LOJ6181 某个套路求和题](https://loj.ac/p/6181)
题意：定义 $f(n)=\prod\limits_{d|n}\mu(d)$，求 $\sum\limits_{i=1}^nf(i)$。
若 $n$ 中有平方因子则 $f(n)=0$，否则 $f(n)$ 仅和 $\omega(n)$ 有关。设 $k=\omega(n)$，有：
$$
f(n)=\prod_{i=0}^k(-1)^i\binom ki=(-1)^{\sum\limits_{i=0}^ki\binom ki}=(-1)^{k2^{k-1}}
$$
最后一步使用的组合恒等式可以通过对 $(1+x)^n$ 求导得到。
当且仅当 $k=1$ 时 $k2^{k-1}$ 是奇数，此时 $f(n)=-1$；否则 $f(n)=1$。
故有：
$$
\sum_{i=1}^nf(i)=\sum_{i=1}^n\mu^2(i)-2\pi(n)
$$
求 $\mu^2$ 的前缀和的方法已经提过很多次了，这里不再赘述。$\pi(n)$ 直接使用 Meissel-Lehmer 即可。
有点卡空间。
```cpp
#include <bits/stdc++.h>
#include <bits/extc++.h>
using namespace std;
using namespace __gnu_pbds;
#define int long long
const int N = 8e6 + 5, I = 3e5, J = 35, mod = 998244353;
int n;
signed tot, mu2[N], pri[N], pi[N], f[I + 5][J + 5];
bool vis[N];
gp_hash_table<int, int> mp;
void sieve() {
mu2[1] = 1;
for (int i = 2, k; i < N; i++) {
if (!vis[i]) pri[++tot] = i, mu2[i] = 1;
for (int j = 1; j <= tot && i * pri[j] < N; j++) {
vis[k = i * pri[j]] = 1;
if (!(i % pri[j])) { mu2[k] = 0; break; }
mu2[k] = mu2[i];
}
}
for (int i = 1; i < N; i++) mu2[i] += mu2[i - 1];
for (int i = 1; i < N; i++) pi[i] = pi[i - 1] + !vis[i];
for (int i = 1; i <= I; i++) f[i][0] = i;
for (int i = 1; i <= I; i++) {
for (int j = 1; j <= J; j++) f[i][j] = f[i][j - 1] - f[i / pri[j]][j - 1];
}
// cerr << "qwq\n";
}
int djs(int n) {
if (n < N) return mu2[n];
if (mp[n]) return mp[n];
int sum = n;
for (int d = 2; d * d <= n; d++) sum -= djs(n / (d * d));
return mp[n] = sum;
}
int dp(int i, int j) {
if (i <= I && j <= J) return f[i][j];
if (!i) return 0; if (!j) return i;
if (i < N && pri[j] * pri[j] >= i) return max(0ll, pi[i] - j);
return dp(i, j - 1) - dp(i / pri[j], j - 1);
}
int Pi(int n) {
// cerr << n << '\n';
if (n < N) return pi[n] - 1;
int m = pi[(int)cbrt(n)] - 1, ans = dp(n, m) + m - 1; m++;
// cerr << m << '\n';
while ((int)pri[m] * pri[m] <= n) /*cerr << n << ' ' << pri[m] << '\n', */ans -= Pi(n / pri[m]) - Pi(pri[m]) + 1, m++;
return ans;
}
signed main() {
sieve(), cin >> n;
// cerr << djs(n) << ' ' << "qwq\n";
// cerr << Pi(n) << ' ' << "qwq\n";
cout << ((djs(n) - 2 * Pi(n)) % mod + mod) % mod;
}
```
---
最后推首歌：[Doping Dance](https://www.bilibili.com/video/BV1BT411E7oP)，搭配 PV 食用更佳 ♪(^∇^*)