杜教筛：快速计算积性函数前缀和的利器

1. 这是用来解决什么问题的？

想象你是一个数学侦探，想统计从1到N（比如N=10^12，超级大）中，有多少个数的莫比乌斯函数值为1、-1或0，或者计算所有数的欧拉函数值加起来是多少。这种“累加”操作叫前缀和，记作：

S_f(n) = f(1)+f(2)+\dots+f(n)

如果N很小（比如1000），我们可以一个一个算。但当N像天文数字（比如10^12）时，连最快的线性筛法都要花很长时间（因为线性筛要枚举每个数）。杜教筛是一种“聪明”的数学+编程技巧，它利用数论卷积和整除分块，把计算时间降到大约 $O(N^{2/3})$ ，也就是比直接算快得多。

生活类比：你想统计全校所有班级的人数总和。每个班级的人数你可能要一个一个数（线性）。但如果你知道每个年级的总人数（前缀和），然后把不同年级的人数加起来，就可能快很多。杜教筛就像用已知的“年级总人数”去快速推算“全校总人数”，而不需要把每个班都数完。

2. 杜教筛的核心思想——用卷积“拆解”问题

2.1 什么是狄利克雷卷积？

两个函数 $f$ 和 $g$ 的狄利克雷卷积定义为：

(f * g)(n) = \sum_{d|n} f(d) \cdot g\left(\frac{n}{d}\right)

简单说，就是把 $n$ 拆成两个因子相乘，然后分别看两个函数在这些因子上的值，再加起来。例如：

如果 $f=1$ （常函数）， $g=1$ ，那么 $(1*1)(n)$ 就是 $n$ 的因数个数。
如果 $f=1$ ， $g=\mathrm{id}$ （恒等函数， $\mathrm{id}(n)=n$ ），那么 $(1*\mathrm{id})(n)$ 就是 $n$ 的所有因数之和。

2.2 核心公式

杜教筛的关键是：我们想求 $S_f(N)$ ，但直接求很难。于是我们构造一个简单易求的函数 $g$ ，让它们的卷积 $h = f * g$ 的前缀和很容易算出来。然后利用下面这个关系：

\sum_{i=1}^N h(i) = \sum_{i=1}^N g(i) \cdot S_f\left(\left\lfloor\frac{N}{i}\right\rfloor\right)

这个公式怎么来的？从卷积的定义交换求和顺序：

\begin{aligned} \sum_{i=1}^N h(i) &= \sum_{i=1}^N \sum_{d|i} f(d) g(i/d) \\ &= \sum_{d=1}^N g(d) \sum_{i=1}^{\lfloor N/d \rfloor} f(i) \quad (\text{先枚举 } d = i/d \text{ 的角色}) \\ &= \sum_{d=1}^N g(d) S_f\left(\left\lfloor\frac{N}{d}\right\rfloor\right) \end{aligned}

把 $d=1$ 的项单独拿出来，并假设 $g(1)=1$ （大多数积性函数都满足），得到：

\sum_{i=1}^N h(i) = S_f(N) + \sum_{d=2}^N g(d) S_f\left(\left\lfloor\frac{N}{d}\right\rfloor\right)

移项就得到杜教筛的核心公式：

\boxed{S_f(N) = S_h(N) - \sum_{d=2}^N g(d) \cdot S_f\left(\left\lfloor\frac{N}{d}\right\rfloor\right)}

这个公式右边递归地出现了 $S_f$ 本身，而且 $\lfloor N/d \rfloor$ 的值只有大约 $2\sqrt{N}$ 种（整除分块）。我们只需要递归计算这些值，并配合记忆化（缓存结果），就能大大降低时间复杂度。

3. 如何选择函数 $g$ ？

选 $g$ 就像“打配合”：要让 $h = f * g$ 的前缀和很简单，同时 $g$ 本身的前缀和也要容易算（因为公式里会用到 $g(d)$ 的累加，我们在整除分块时需要快速求一段 $g$ 的和）。

常用配对：

目标 $f$	选 $g$	卷积 $h$	$S_h(N)$	说明
$\mu$ （莫比乌斯函数）	1	$\mu * 1 = \varepsilon$ （单位函数，只有 $\varepsilon(1)=1$ ）	1	因为 $\sum_{i=1}^N \varepsilon(i) = 1$
$\varphi$ （欧拉函数）	1	$\varphi * 1 = \mathrm{id}$ （恒等函数）	$N(N+1)/2$	因为 $\sum_{i=1}^N i = \frac{N(N+1)}{2}$
$\mathrm{id}$	1	$\mathrm{id} * 1 = \sigma$ （除数之和）	需要另外求？	通常不用杜教筛，直接用公式

例1：求Möbius函数前缀和
取 $f=\mu, g=1$ ，则 $h=\mu*1=\varepsilon$ ，所以 $S_h(N)=1$ 。公式变为：

S_\mu(N) = 1 - \sum_{d=2}^N 1 \cdot S_\mu\left(\left\lfloor\frac{N}{d}\right\rfloor\right) = 1 - \sum_{d=2}^N S_\mu\left(\left\lfloor\frac{N}{d}\right\rfloor\right)

例2：求欧拉函数前缀和
取 $f=\varphi, g=1$ ，则 $h=\varphi*1=\mathrm{id}$ ， $S_h(N)=\frac{N(N+1)}{2}$ 。公式变为：

S_\varphi(N) = \frac{N(N+1)}{2} - \sum_{d=2}^N S_\varphi\left(\left\lfloor\frac{N}{d}\right\rfloor\right)

4. 整除分块——把递归的层数压下来

公式里需要对 $d$ 从2到N求和。但直接循环 $N$ 次就太慢了。但是注意到： $\lfloor N/d \rfloor$ 随着 $d$ 变化，只有大约 $2\sqrt{N}$ 个不同的值。比如 $N=100$ ，当 $d=51$ 时 $\lfloor100/51\rfloor=1$ ，当 $d=100$ 时也是1，所以从51到100这一大段，商都是1。我们只需要把这一整段一起处理。

具体做法：对于当前 $d$ ，设 $q = \lfloor N/d \rfloor$ ，那么下一个不同商的起点是 $\lfloor N/q \rfloor + 1$ 。在这一段内， $g(d)$ 的和可以用前缀和快速求出（因为 $g$ 的前缀和也容易算）。

生活类比：假设你要统计班级里每个人的数学分数。直接问每个人（循环N次）太慢。但如果分数是重复的（比如很多人考了90分），你可以先按分数分组，统计每组的人数，然后用“人数×分数”快速算出总分。整除分块就是按“商”分组。

5. 编程实现——计算莫比乌斯函数前缀和（完整代码）

我们要实现 $S_\mu(N)$ 。按照以下步骤：

先用线性筛预处理前 $M$ 个数的 $\mu$ 和它的前缀和（ $M$ 通常取 $N^{2/3}$ 左右）。
对于大于 $M$ 的 $N$ ，使用杜教筛递归计算，并记忆化结果。
在递归中，利用整除分块快速求和。

5.1 C++ 完整代码（动态预处理上限）

#include <iostream>
#include <vector>
#include <unordered_map>
#include <cmath>
using namespace std;

// 全局变量，用于记忆化和前缀和
unordered_map<long long, long long> memo;  // 记忆化 S_mu 的结果
vector<long long> mu_sum;                  // 预处理的前缀和（下标从0开始，mu_sum[i] = S_mu(i)）
int limit;                                 // 预处理的上限

// 线性筛预处理 mu 和前缀和，直到 up_to
void precompute(int up_to) {
    vector<int> mu(up_to + 1, 0);
    vector<bool> is_composite(up_to + 1, false);
    vector<int> primes;
    mu[1] = 1;
    for (int i = 2; i <= up_to; i++) {
        if (!is_composite[i]) {
            primes.push_back(i);
            mu[i] = -1;
        }
        for (size_t j = 0; j < primes.size() && i * primes[j] <= up_to; j++) {
            int p = primes[j];
            int x = i * p;
            is_composite[x] = true;
            if (i % p == 0) {
                mu[x] = 0;      // 有平方因子
                break;
            } else {
                mu[x] = mu[i] * mu[p];
            }
        }
    }
    mu_sum.resize(up_to + 1);
    mu_sum[0] = 0;
    for (int i = 1; i <= up_to; i++) {
        mu_sum[i] = mu_sum[i-1] + mu[i];
    }
}

// 杜教筛求 S_mu(N)
long long S_mu(long long N) {
    if (N <= limit) return mu_sum[N];      // 已经预处理过了
    if (memo.count(N)) return memo[N];     // 记忆化
    long long res = 1;                     // S_h(N) = 1，因为 h = epsilon
    // 整除分块，d 从2开始
    long long d = 2;
    while (d <= N) {
        long long q = N / d;               // 商
        long long next_d = N / q + 1;      // 下一个不同的 d
        // 这一段的 d 范围是 [d, next_d-1]，它们的 floor(N/d) 都是 q
        // 需要求和 g(d) = 1，所以块长度 = next_d - d
        res -= (next_d - d) * S_mu(q);
        d = next_d;
    }
    memo[N] = res;
    return res;
}

int main() {
    long long N;
    cout << "请输入 N（例如 1000000 或 10000000000）: ";
    cin >> N;
    // 预处理上限取 N^(2/3)，但不要超过 1e7 以防内存爆炸
    limit = (int)min((long long)1e7, (long long)pow(N, 2.0/3.0));
    if (limit < 1) limit = 1;
    precompute(limit);
    cout << "S_mu(" << N << ") = " << S_mu(N) << endl;
    return 0;
}

代码解释：

mu_sum[i] 存储了从1到i的莫比乌斯函数值之和，方便直接查表。
memo 是一个哈希表，缓存大N的递归结果，避免重复计算。
在 S_mu 中，res 初始为1（对应 $S_h(N)=1$ ），然后减去所有 $d\ge 2$ 的贡献。
整除分块中，next_d - d 就是这一段的长度，乘以 S_mu(q)（因为 $g(d)=1$ ）。

5.2 Python 完整代码

import math
import sys
sys.setrecursionlimit(10**6)   # 防止递归深度太大

memo = {}
limit = 0
mu_sum = []

def precompute(up_to):
    """线性筛预处理 mu 和前缀和到 up_to"""
    mu = [0] * (up_to + 1)
    mu[1] = 1
    primes = []
    is_composite = [False] * (up_to + 1)
    for i in range(2, up_to + 1):
        if not is_composite[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            if i * p > up_to:
                break
            is_composite[i * p] = True
            if i % p == 0:
                mu[i * p] = 0
                break
            else:
                mu[i * p] = mu[i] * mu[p]
    # 计算前缀和
    mu_sum.append(0)   # 下标0占位
    for i in range(1, up_to + 1):
        mu_sum.append(mu_sum[-1] + mu[i])
    return mu_sum

def S_mu(N):
    """杜教筛求莫比乌斯函数前缀和"""
    if N <= limit:
        return mu_sum[N]
    if N in memo:
        return memo[N]
    res = 1   # S_h(N) = 1
    d = 2
    while d <= N:
        q = N // d
        next_d = N // q + 1
        # 这一段的长度乘以 S_mu(q)
        res -= (next_d - d) * S_mu(q)
        d = next_d
    memo[N] = res
    return res

if __name__ == "__main__":
    N = int(input("请输入 N："))
    # 预处理上限取 N^(2/3)，但不超过 10^7
    limit = min(int(1e7), int(N ** (2/3)))
    if limit < 1:
        limit = 1
    mu_sum = precompute(limit)
    print(f"S_mu({N}) = {S_mu(N)}")

6. 常见错误与避坑指南

遗忘记忆化：递归时如果不把结果存起来，会陷入指数级重复计算。必须用哈希表或数组缓存。
整型溢出：C++ 中如果 N 很大（如 10^12），N * (N+1) 可能超过 64 位？实际上 10^12 * 10^12 会溢出，但在计算 $S_h(N)$ 时通常用 N * (N+1) / 2，中间结果要用 long long，且要小心除以2的顺序。更好的是用 N / 2 * (N+1) 或强制类型转换。
预处理上限选择不当：上限太小会导致递归太多层，效率变差；上限太大可能内存超限。经验值是 $N^{2/3}$ 左右，不超过 1e7 通常安全。
整除分块边界错误：next_d = N / q + 1 要确保处理完 d 到最后 N+1 时跳出循环。注意当 q = 0 时（d > N），循环应该结束。
g(1) 不等于1：公式中假设 g(1)=1，如果 g(1)≠1，公式要稍作修改，但大多数常用积性函数 g(1)=1。
递归深度过大：Python 需要设置 sys.setrecursionlimit，C++ 使用循环或迭代，但递归深度一般不会太深（因为整除分块会快速缩小参数）。

7. 用杜教筛求欧拉函数前缀和

只要把代码中的 $f$ 和 $g$ 换一下就 OK。对于欧拉函数：

S_\varphi(N) = \frac{N(N+1)}{2} - \sum_{d=2}^N S_\varphi\left(\left\lfloor\frac{N}{d}\right\rfloor\right)

因此 S_phi(N) 的递归函数如下（C++ 风格）：

long long S_phi(long long N) {
    if (N <= limit) return phi_sum[N];
    if (memo.count(N)) return memo[N];
    long long res = N * (N + 1) / 2;   // S_h(N) = N(N+1)/2
    long long d = 2;
    while (d <= N) {
        long long q = N / d;
        long long next_d = N / q + 1;
        res -= (next_d - d) * S_phi(q);
        d = next_d;
    }
    memo[N] = res;
    return res;
}

只需要把线性筛部分改成筛 $\varphi$ 即可。

8. 总结与相关指引

杜教筛是竞赛中处理大范围积性函数前缀和的“杀手锏”。它的核心三步：

选一个好搭档 $g$ ，使得 $h = f * g$ 和 $g$ 的前缀和都容易算。
用核心公式递归： $S_f(N) = S_h(N) - \sum_{d=2}^N g(d) S_f(\lfloor N/d \rfloor)$ 。
整除分块 + 记忆化 加速递归。

掌握了杜教筛，你就能轻松应付 $N$ 高达 $10^{12}$ 的前缀和问题。如果想进一步学习，可以看看以下方向：

莫比乌斯反演：杜教筛的基石之一，用于推导很多等价关系。
狄利克雷卷积：更多有趣的搭配（如 $f = \mathrm{id}^k$ 与 $g = \mu$ ）。
数论分块（整除分块）的更多应用，比如快速求 $\sum \lfloor N/i \rfloor$ 。
欧拉函数与莫比乌斯函数的性质：理解它们的积性有助于构造 $g$ 。

试一试：请写出杜教筛求 $S_{id}(N) = 1+2+...+N$ 的代码？其实直接用公式，但你能用杜教筛验证公式吗？（提示：取 $f = \mathrm{id}, g = \mu$ ，则 $h = \mathrm{id} * \mu = \varphi$ ，然后 $S_\varphi$ 用杜教筛求，最后解出 $S_{id}$ ——有点绕，但很有趣。）

希望这篇文章能帮你打开杜教筛的大门！

杜教筛简介