问题背景
在数论问题中，经常需要计算积性函数的前缀和。例如，求前 $n$ 个正整数的欧拉函数之和 $\sum_{i=1}^{n} \varphi(i)$，或莫比乌斯函数之和 $\sum_{i=1}^{n} \mu(i)$。当 $n$ 达到 $10^{11}$ 级别时，线性筛法 $O(n)$ 的时间和空间复杂度都不可接受。
杜教筛（Dujiao Sieve）由杜瑜皓提出，能够在 $O(n^{2/3})$ 或 $O(n^{3/4})$ 的时间复杂度内计算积性函数前缀和，空间复杂度仅为 $O(\sqrt{n})$。其核心思想是利用狄利克雷卷积将难求的前缀和转化为易求的前缀和组合。
核心原理
狄利克雷卷积
两个数论函数 $f$ 和 $g$ 的狄利克雷卷积定义为：
$$(f * g)(n) = \sum_{d|n} f(d) \cdot g(n/d)$$
关键推导
设积性函数 $f$，我们希望求 $S(n) = \sum_{i=1}^{n} f(i)$。若存在数论函数 $g$ 使得 $f * g = h$，且 $h$ 和 $g$ 的前缀和都容易计算，则有：
$$\sum_{i=1}^{n} (f * g)(i) = \sum_{i=1}^{n} \sum_{d|i} f(d) \cdot g(i/d) = \sum_{d=1}^{n} g(d) \sum_{k=1}^{\lfloor n/d \rfloor} f(k) = \sum_{d=1}^{n} g(d) \cdot S(\lfloor n/d \rfloor)$$
将 $d=1$ 的项分离出来：
$$h(1) \cdot S(n) + \sum_{d=2}^{n} g(d) \cdot S(\lfloor n/d \rfloor) = \sum_{i=1}^{n} h(i)$$
由于 $\lfloor n/d \rfloor$ 只有 $O(\sqrt{n})$ 种不同的取值，我们可以递归+记忆化计算所有需要的 $S$ 值。
常用组合
| 目标函数 $f$ | 辅助函数 $g$ | 卷积结果 $h = f * g$ |
|---|---|---|
| $\varphi$（欧拉函数） | $I(n)=1$ | $id(n)=n$ |
| $\mu$（莫比乌斯函数） | $I(n)=1$ | $\varepsilon(n)=[n=1]$ |
对于欧拉函数：$\sum_{d|n} \varphi(d) = n$，因此 $\varphi * I = id$。
代入公式：
$$\sum_{i=1}^{n} i = \sum_{d=1}^{n} S_{\varphi}(\lfloor n/d \rfloor)$$
$$\frac{n(n+1)}{2} = S_{\varphi}(n) + \sum_{d=2}^{n} S_{\varphi}(\lfloor n/d \rfloor)$$
整理得：
$$S_{\varphi}(n) = \frac{n(n+1)}{2} – \sum_{d=2}^{n} S_{\varphi}(\lfloor n/d \rfloor)$$
Java实现
import java.util.*;
/**
* 杜教筛：计算积性函数前缀和
* 时间复杂度：O(n^(2/3)) 预处理 + O(n^(2/3)) 查询
* 空间复杂度：O(n^(2/3))
*/
public class DujiaoSieve {
// 预处理的上限，取 n^(2/3) 级别
private static final int PRE_LIMIT = 5000000;
// 线性筛相关数组
private int[] primes;
private boolean[] isComposite;
private int primeCount;
// 预处理的前缀和数组（用于小范围直接查表）
private long[] prefixPhi; // 欧拉函数前缀和
private long[] prefixMu; // 莫比乌斯函数前缀和
// 记忆化哈希表：存储大范围的查询结果
private HashMap<Long, Long> memoPhi;
private HashMap<Long, Long> memoMu;
public DujiaoSieve() {
initLinearSieve();
memoPhi = new HashMap<>();
memoMu = new HashMap<>();
}
/**
* 线性筛预处理欧拉函数和莫比乌斯函数
* 时间复杂度：O(PRE_LIMIT)
*/
private void initLinearSieve() {
int n = PRE_LIMIT;
isComposite = new boolean[n + 1];
primes = new int[n + 1];
primeCount = 0;
// phi[i] 表示 i 的欧拉函数值
int[] phi = new int[n + 1];
// mu[i] 表示 i 的莫比乌斯函数值
int[] mu = new int[n + 1];
phi[1] = 1;
mu[1] = 1;
for (int i = 2; i <= n; i++) {
if (!isComposite[i]) {
primes[primeCount++] = i;
phi[i] = i - 1; // 质数p的欧拉函数值为 p-1
mu[i] = -1; // 质数p的莫比乌斯函数值为 -1
}
for (int j = 0; j < primeCount && i * primes[j] <= n; j++) {
isComposite[i * primes[j]] = true;
if (i % primes[j] == 0) {
// i 是 primes[j] 的倍数
phi[i * primes[j]] = phi[i] * primes[j];
mu[i * primes[j]] = 0; // 含有平方因子
break;
} else {
phi[i * primes[j]] = phi[i] * (primes[j] - 1);
mu[i * primes[j]] = -mu[i];
}
}
}
// 计算前缀和
prefixPhi = new long[n + 1];
prefixMu = new long[n + 1];
for (int i = 1; i <= n; i++) {
prefixPhi[i] = prefixPhi[i - 1] + phi[i];
prefixMu[i] = prefixMu[i - 1] + mu[i];
}
}
/**
* 计算欧拉函数前缀和：sum_{i=1}^{n} phi(i)
* 利用公式：phi * I = id
* S_phi(n) = n*(n+1)/2 - sum_{d=2}^{n} S_phi(n/d)
*/
public long sumPhi(long n) {
if (n <= 0) return 0;
// 小范围直接查表
if (n <= PRE_LIMIT) {
return prefixPhi[(int) n];
}
// 大范围查记忆化表
if (memoPhi.containsKey(n)) {
return memoPhi.get(n);
}
// 核心公式：S_phi(n) = n*(n+1)/2 - sum_{d=2}^{n} S_phi(n/d)
// 利用整除分块优化：n/d 在一段区间内值相同
long result = n * (n + 1) / 2;
long l = 2;
while (l <= n) {
long q = n / l; // 当前块的值
long r = n / q; // 当前块的右边界
// sum_{d=l}^{r} S_phi(n/d) = (r - l + 1) * S_phi(q)
result -= (r - l + 1) * sumPhi(q);
l = r + 1;
}
memoPhi.put(n, result);
return result;
}
/**
* 计算莫比乌斯函数前缀和：sum_{i=1}^{n} mu(i)
* 利用公式：mu * I = epsilon
* S_mu(n) = 1 - sum_{d=2}^{n} S_mu(n/d)
*/
public long sumMu(long n) {
if (n <= 0) return 0;
if (n <= PRE_LIMIT) {
return prefixMu[(int) n];
}
if (memoMu.containsKey(n)) {
return memoMu.get(n);
}
// 核心公式：S_mu(n) = 1 - sum_{d=2}^{n} S_mu(n/d)
long result = 1;
long l = 2;
while (l <= n) {
long q = n / l;
long r = n / q;
result -= (r - l + 1) * sumMu(q);
l = r + 1;
}
memoMu.put(n, result);
return result;
}
/**
* 通用杜教筛框架
* 对于满足 f * g = h 的积性函数，计算 sum_{i=1}^{n} f(i)
* @param n 上限
* @param sumH 函数 h 的前缀和计算器 h_sum(n) = sum_{i=1}^{n} h(i)
* @param sumG 函数 g 的前缀和计算器 g_sum(n) = sum_{i=1}^{n} g(i)
* @param smallSum 小范围直接查表的前缀和数组
* @return sum_{i=1}^{n} f(i)
*/
public long dujiaoSieve(long n,
java.util.function.LongUnaryOperator sumH,
java.util.function.LongUnaryOperator sumG,
long[] smallSum) {
if (n <= 0) return 0;
if (n < smallSum.length) {
return smallSum[(int) n];
}
// 使用通用记忆化（简化版，实际使用时需要为不同函数单独存储）
long result = sumH.applyAsLong(n);
long l = 2;
while (l <= n) {
long q = n / l;
long r = n / q;
// 减去 g(d) * S_f(n/d) 的贡献
// 这里假设 g(1) = 1，否则需要调整
long gSum = sumG.applyAsLong(r) - sumG.applyAsLong(l - 1);
// 注意：通用版本需要递归调用，但为了避免复杂度爆炸
// 实际使用时建议针对具体函数特化
l = r + 1;
}
return result;
}
/**
* 主函数：测试杜教筛性能
*/
public static void main(String[] args) {
DujiaoSieve sieve = new DujiaoSieve();
// 测试小规模数据，与线性筛结果对比验证正确性
System.out.println("=== 正确性验证 ===");
long[] testValues = {1, 10, 100, 1000, 10000, 100000, 1000000};
for (long n : testValues) {
long phiSum = sieve.sumPhi(n);
long muSum = sieve.sumMu(n);
System.out.printf("n = %d: sum_phi = %d, sum_mu = %d%n", n, phiSum, muSum);
}
// 大规模测试，展示杜教筛的优势
System.out.println("\n=== 大规模测试 ===");
long[] largeTests = {
1_000_000_000L, // 10^9
10_000_000_000L, // 10^10
100_000_000_000L, // 10^11
1_000_000_000_000L // 10^12
};
for (long n : largeTests) {
long start = System.currentTimeMillis();
long phiSum = sieve.sumPhi(n);
long end = System.currentTimeMillis();
System.out.printf("n = %d: sum_phi = %d, time = %d ms%n",
n, phiSum, (end - start));
}
// 验证欧拉函数前缀和的一个数论性质：
// sum_{i=1}^{n} phi(i) 约等于 3*n^2 / PI^2
System.out.println("\n=== 渐近性质验证 ===");
long n = 1_000_000_000L;
long phiSum = sieve.sumPhi(n);
double approx = 3.0 * n * n / (Math.PI * Math.PI);
System.out.printf("n = %d%n", n);
System.out.printf("sum_phi(n) = %d%n", phiSum);
System.out.printf("3*n^2/PI^2 = %.2f%n", approx);
System.out.printf("相对误差 = %.4f%%%n",
Math.abs(phiSum - approx) / approx * 100);
}
}
整除分块技巧详解
杜教筛效率的关键在于整除分块。对于固定的 $n$，当 $d$ 在某个区间 $[l, r]$ 内变化时，$\lfloor n/d \rfloor$ 的值保持不变。
具体地，若当前 $l$ 已知，则 $q = \lfloor n/l \rfloor$，而最大的 $r$ 满足 $\lfloor n/r \rfloor = q$，即 $r = \lfloor n/q \rfloor$。
这意味着内层求和 $\sum_{d=l}^{r} S(\lfloor n/d \rfloor)$ 中，$\lfloor n/d \rfloor$ 都是同一个值 $q$，可以直接用 $(r – l + 1) \times S(q)$ 计算。
对于 $n$ 的每一种不同取值，$\lfloor n/d \rfloor$ 只有 $O(\sqrt{n})$ 种，因此单次查询的时间复杂度为 $O(\sqrt{n})$。但由于递归调用的存在，总复杂度需要更精细的分析。
复杂度分析
时间复杂度
设预处理上限为 $N = n^{2/3}$。
- 预处理阶段：线性筛 $O(N)$，计算前缀和 $O(N)$。
- 查询阶段：对于 $m > N$ 的查询，递归调用的参数都是 $\lfloor n/d \rfloor$ 的形式，这些值只有 $O(\sqrt{n})$ 种。但由于多层递归，实际访问的不同参数个数约为 $O(n / N) = O(n^{1/3})$。每个参数的计算涉及整除分块 $O(\sqrt{n})$，但利用记忆化后，总复杂度为 $O(n^{2/3})$。
综合最优取 $N = n^{2/3}$，总时间复杂度为 $O(n^{2/3})$。
空间复杂度
- 预处理数组：$O(N) = O(n^{2/3})$
- 记忆化哈希表：$O(n^{2/3})$
应用场景
- 数论函数求和：计算欧拉函数、莫比乌斯函数、约数函数等前缀和
- 莫比乌斯反演：$F(n) = \sum_{d|n} f(d) \Leftrightarrow f(n) = \sum_{d|n} \mu(d) F(n/d)$
- 组合计数：需要计算与互质相关的计数问题
- 信息竞赛：处理 $n \leq 10^{11}$ 的数论问题
扩展阅读
杜教筛可以与洲阁筛、Min_25筛结合使用，处理更一般的积性函数求和问题。其中 Min_25筛 在 $O(n^{3/4} / \log n)$ 或 $O(n^{2/3})$ 的复杂度内可以处理几乎所有积性函数前缀和，是信息学竞赛中的终极数论工具。
