杜教筛——亚线性处理数论函数求和

问题引入

给定一个数论函数 \(f(x)\)，求 \(\sum\limits_{i=1}^nf(i)\)。

对 \(n \le 10^7\) 甚至 \(n \le 10^8\) 都是好做的，\(\mathcal O(n)\) 解决即可，但如果 \(n < 2^{31}\) 呢？

这就需要亚线性时间复杂度的算法，杜教筛 就是其一。

杜教筛

杜教筛是一种能在幂时间 \(\mathcal O(n^\frac23)\) 下对数论函数进行求和的算法。

记 \(S(n) = \sum\limits_{i=1}^nf(i)\)，所求即 \(S(n)\)。

考虑构造另一个函数 \(g(x)\)，然后做狄利克雷卷积：

\[(f * g)(n) = \sum_{d|n}f(d)g\left(\dfrac nd\right) \]

那么有：

\[\sum_{i=1}^n(f * g)(i) = \sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}f(i) = \sum_{d=1}^ng(d)S\left(\lfloor\dfrac nd\rfloor\right) \]

把 \(g(1)S(n)\) 拿出来，就有：

\[g(1)S(n) = \sum_{i=1}^n(f*g)(i) - \sum_{i=2}^ng(d)S\left(\lfloor\dfrac nd\rfloor\right) \]

然后关于 \(\lfloor\dfrac nd\rfloor\) 数论分块就行。

所以 \(g(x)\) 需要满足以下性质：

• \(\sum\limits_{i=1}^n(f*g)(i)\) 可以在亚线性时间复杂度内求得。
• \(g(x)\) 的前缀和可以在亚线性时间复杂度内求得。

关于时间复杂度，将求 \(\sum\limits_{i=1}^n(f*g)(i)\) 和求 \(g(x)\) 的前缀和的时间复杂度都视作 \(\mathcal O(1)\) 的话，可以有：

\[T(n) = \mathcal O(\sqrt n) + \mathcal O\left(\sum_{i=2}^{\sqrt n} T(i) + T\left(\lfloor\dfrac ni\rfloor\right)\right) \]

展开，有：

\[T(n) = \mathcal O(\sqrt n) + \mathcal O\left(\sum_{i=2}^\sqrt n\mathcal O(\sqrt i) + \sum_{j=2}^\sqrt i T(j) + \mathcal O\left(\sqrt{\lfloor\dfrac ni\rfloor}\right) + \sum_{j=2}^\sqrt{\lfloor\frac ni\rfloor} T(j)\right) \]

因为 \(\lim\limits_{x \to \infin} \dfrac{n^\frac12}{n^\frac14} = \infin\)，所以 \(\sum\limits_{j=1}^\sqrt i T(j)\) 和 \(\sum\limits_{j=1}^\sqrt{\lfloor\frac ni\rfloor} T(j)\) 可视作高阶无穷小，忽略不计；\(\sqrt i < \sqrt{\dfrac ni}\)，忽略不计：

\[T(n) = \mathcal O(\sqrt n) + \mathcal O\left(\sum_{i=2}^\sqrt n \sqrt{\lfloor\frac ni\rfloor} \right) < \mathcal O\left(\int_0^\sqrt n\sqrt{\frac nx}\,\mathrm dx \right) = \mathcal O(n^\frac34) \]

积分步骤戳这里 O.o

考虑优化：前 \(B > n^{0.5}\) 个前缀和暴力求。时间复杂度为：

\[\mathcal O(B) + \mathcal O\left(\int_0^\frac nb \sqrt\frac nx \,\mathrm dx \right) = \mathcal O(B) + \mathcal O\left(\dfrac n{\sqrt B}\right) \]

积分步骤戳这里 o.O

根据均值不等式，当 \(B\) 取 \(n^\frac23\) 时，时间复杂度最小，为 \(\mathcal O(n^\frac23)\)。

实际实现中，因为递归部分常数较大，\(B\) 可以取大一些。例如 \(n < 2^{31}\) 时，可以取 \(B = 10^7\)。

模板

\(\mathcal{PORTAL}\)

分别求：

\[\sum_{i=1}^n \varphi(i) \\ \sum_{i=1}^n \mu(i) \]

\(1 \le n \le 2^{31}\)。

显然杜教筛，分别考虑 \(g\) 就好。

\(f = \varphi\)

注意到 欧拉函数的性质 中有一条是 \(\sum\limits_{d | n} \varphi(d) = n\)，所以我们取 \(g = I\)（\(\forall x, I(x) = 1\)），就有 \((f * g)(n) = n\)，\(g\) 的前缀和也显然。

\(f = \mu\)

注意到 莫比乌斯函数的性质 中有一条是 \(\sum\limits_{d|n}\mu(d) = \begin{cases}1 & n = 1 \\ 0 & \text{otherwise} \end{cases}\)，所以我们同样取 \(g = I\)，就有 \((f * g)(n) = \begin{cases}1 & n = 1 \\ 0 & \text{otherwise} \end{cases}\)。

代码：

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

constexpr int N = 1e7 + 10;

int n, mu[N]; ll phi[N];
int prn, prime[N];
bool vis[N];

map<int, ll> S;

void init(int n) {
    phi[1] = mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) {
            phi[i] = i - 1, mu[i] = -1;
            prime[++prn] = i;
        }
        for (int j = 1; j <= prn && i * prime[j] <= n; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j]) phi[i * prime[j]] = phi[i] * phi[prime[j]], mu[i * prime[j]] = -mu[i];
            else {phi[i * prime[j]] = phi[i] * prime[j]; break;}
        }
    }
    for (int i = 2; i <= n; i++) phi[i] += phi[i - 1], mu[i] += mu[i - 1];
}

ll Sphi(int n) {
    if (n <= 1e7) return phi[n];
    if (S.find(n) != S.end()) return S[n];
    ll res = (n + 1ll) * n / 2;
    for (uint l = 2, r; l <= n; l = r + 1) r = n / (n / l), res -= (r - l + 1) * Sphi(n / l);
    return S[n] = res;
}

int Smu(int n) {
    if (n <= 1e7) return mu[n];
    if (S.find(n) != S.end()) return S[n];
    int res = 1;
    for (uint l = 2, r; l <= n; l = r + 1) r = n / (n / l), res -= (r - l + 1) * Smu(n / l);
    return S[n] = res;
}

void solve() {
    cin >> n; 
    S.clear();
    cout << Sphi(n) << ' ';
    S.clear();
    cout << Smu(n) << '\n';
}

int main() {
    ios::sync_with_stdio(0); cin.tie(nullptr), cout.tie(nullptr);
    init(1e7); int t; cin >> t;
    while (t--) solve();
    return 0;
}