1. 这个筛法到底在解决什么问题先说说我为什么会对这个题目感兴趣。做数论相关的计算任务时经常遇到一类需求给定一个积性函数 (f(n))需要快速求出它的前缀和[ S_f(N)\sum_{n1}^{N} f(n) ]这里的 (N) 往往不是几百、几千而是一下子给到 (10^{10})、(10^{11}) 这种规模。最常见的几个例子(\mu(n)) 的 Mertens 函数前缀和 (M(N)\sum_{n\le N}\mu(n))欧拉函数前缀和 (\Phi(N)\sum_{n\le N}\varphi(n))除数函数前缀和 (\sum_{n\le N}d(n))更一般地(\sum_{n\le N}\sigma_k(n))如果是竞赛题或者研究里的一个子步骤(N10^8) 用线性筛加前缀和还能勉强扛一扛一旦 (N) 到了 (10^{10}) 往上直接做线性筛就完全不现实了——时间、内存双双爆炸。这时候就需要一种筛法思想利用积性函数的结构特征把复杂度从 (O(N)) 压到大约 (O(N^{3/4}/\log N)) 或者更优的水平同时把内存压到 (O(\sqrt N)) 级别。这篇文章想写给三类人一类是搞算法竞赛、被这类题目卡过复杂度的同学一类是做计算数论、需要频繁算大范围积性函数前缀和的研究者还有一类是单纯想看看筛这种古老工具在现代计算里还能玩出什么花样的读者。我会先从数学上把预处理逻辑讲清楚再给出一套可以直接改着用的实现最后聊几个我实际踩过的坑。2. 前置知识积性函数、Dirichlet 卷积与前缀和可加速的根源2.1 积性函数的定义和关键性质一个数论函数 (f(n)) 称为积性函数如果对任意互素的正整数 (a,b)都有[ f(ab)f(a)f(b) ]如果去掉互素这个限制对任意 (a,b) 都成立就叫做完全积性函数。积性函数最方便的地方在于只要知道它在素数幂 (p^k) 上的取值其他所有 (n) 的取值就由质因数分解唯一确定了。也就是说研究一个积性函数本质上只需要研究素数幂这个原子单元。2.2 Dirichlet 卷积前缀和为什么能和筛挂上钩两个数论函数的 Dirichlet 卷积定义为[ (f*g)(n)\sum_{d\mid n} f(d),g(n/d) ]这个运算在线性筛的世界里几乎是无处不在的。比如[ d(n)1*1 ][ \sigma(n)1*\mathrm{id} ][ \varphi(n)\mu*\mathrm{id} ]为什么卷积和前缀和有关系因为交换求和顺序之后很多复杂的前缀和可以被拆成若干个简单前缀和的组合。设[ h(n)(f*g)(n) ]那么[ \sum_{n\le N} h(n) \sum_{n\le N}\sum_{d\mid n} f(d)g(n/d) \sum_{d\le N} f(d)\sum_{m\le N/d} g(m) ]这看起来只是把求和顺序换了一下但配合整除分块就能产生实质性的加速效果。这就是后面 Dirichlet 双曲法的核心。2.3 一个自己动手就能验证的简单例子先拿完全积性函数来练手如果 (f) 就是常函数 (1)那么 (\sum_{n\le N}1N)这个平凡到不需要筛。但如果 (f(n)n)那么[ \sum_{n\le N}n\frac{N(N1)}{2} ]也是现成的公式。这类本身就能直接求和的函数是更复杂筛法的地基因为后面构建辅助表时需要能够快速回答前 (k) 个素数的函数值之和以及区间内函数值之和这类子问题。3. 先看一个温和版本Dirichlet 双曲法3.1 双曲法的数学原理假设我们要计算的是某个卷积型函数[ h(n)(a*b)(n) ]的前缀和 (H(N)\sum_{n\le N}h(n))。直接按卷积定义算需要枚举所有因子复杂度是 (O(N\log N))。但如果 (a,b) 的前缀和很好算比如是常数、线性函数、或者能用整除分块快速求和的函数我们可以把整个 (n\le N) 的求和区域按照两个因子的乘积关系做一个划分。具体的划分逻辑是这样的[ \sum_{n\le N}h(n) \sum_{d m \le N} a(d)b(m) ]把二维网格 ((d,m)) 中满足 (dm\le N) 的点按 (d\le \sqrt N) 和 (m\le \sqrt N) 分成两块再加上中间重复计算的部分。更常用的写法是取一个阈值 (B\lfloor\sqrt N\rfloor)然后[ H(N)\sum_{d\le B} a(d)\sum_{m\le N/d} b(m) \sum_{m\le B} b(m)\sum_{d\le N/m} a(d) -\sum_{d\le B}\sum_{m\le B} a(d)b(m) ]这个式子中每一项都只需要整除分块式的前缀和调用复杂度变成 (O(\sqrt N))。名字里的双曲就是对应 (dm\le N) 这个双曲线区域。3.2 经典例子除数函数前缀和取 (ab1)那么 (hd)除数函数于是[ \sum_{n\le N}d(n)2\sum_{d\le \sqrt N}\left\lfloor\frac{N}{d}\right\rfloor-\lfloor\sqrt N\rfloor^2 ]这个公式在 (N10^{12}) 的情况下可以用 (O(\sqrt N)) 的循环几秒内算完。我第一次看到这个式子是在做某个整数分解相关的项目时当时需要反复计算大区间里的除数函数前缀和用这个公式之后之前觉得算不动的问题突然就变得轻松了。3.3 双曲法的问题不是所有积性函数都适用双曲法很漂亮但它有一个苛刻的前提必须能把目标函数写成两个前缀和都极其好算的函数的卷积。比如 (\mu) 虽然可以写成 (\mu1^{-1})但这个逆卷积并不会让计算变简单。那怎么处理 (\mu)、(\varphi) 这类更硬的函数这就轮到文章标题里说的那类筛法出场了。它的思路不是寻找一个现成的卷积分解而是直接从素数结构出发用一套系统性的递推把每个数的质因子贡献剥离开。4. 核心筛法从素数表到辅助前缀和的完整框架4.1 目标的数学设定为了不绕弯子下面统一考虑这样的问题给定一个积性函数 (f)满足在素数处的取值 (f(p)) 可以写成关于 (p) 的简单多项式比如 (f(p)p^23p1)并且 (f(p^k)) 能快速计算。求 (S_f(N)\sum_{n\le N}f(n))。这类函数覆盖了几乎所有常见的积性函数(\mu) 满足 (f(p)-1)(\varphi) 满足 (f(p)p-1)(\sigma_k) 满足 (f(p)1p\cdotsp^k) 但注意它在素数处不是多项式不过也能处理后面我会讲扩展。筛法的核心思想可以概括成一句话把所有合数的贡献按最小素因子least prime factor, lpf分组逐步筛掉小素数用动态规划维护辅助前缀和。4.2 第一步线性筛出 (\sqrt N) 以内的素数表这一步是基础。虽然 (N) 很大但 (\sqrt N) 一般也就 (10^5) 到 (10^6) 量级用普通埃氏筛或者线性筛完全足够。筛完后得到一个素数序列[ p_12,\ p_23,\ p_35,\dots,\ p_m ]其中 (p_m\le \sqrt N)。为什么要筛到 (\sqrt N) 而不是 (N)因为任何一个 (n\le N) 如果有一个素因子 (\sqrt N)那它至多有一个这样的素因子并且剩下部分一定 (\le \sqrt N)。这个观察是所有分治类筛法的根基。4.3 第二步把 (N) 的整除分块值离散化定义集合[ V{ \lfloor N/i\rfloor : 1\le i\le N } ]这个集合的大小只有 (2\lfloor\sqrt N\rfloor) 左右。原因很直观当 (i\le\sqrt N) 时(\lfloor N/i\rfloor) 最多有 (\sqrt N) 个不同值当 (i\sqrt N) 时(\lfloor N/i\rfloor\sqrt N)也最多有 (\sqrt N) 个不同值。我们的所有 DP 表都会以这些分块值为索引而不是以 (1,2,\dots,N) 为索引。这样内存就从 (O(N)) 降到了 (O(\sqrt N))。具体实现时为了快速把一个整数 (v) 映射到下标通常维护两个数组若 (v\le \sqrt N)存在id1[v]否则存在id2[N/v]这样任意分块值都能在 (O(1)) 时间内查到下标。4.4 第三步用埃筛式 DP 构造素数函数前缀和这里需要引入一个我习惯叫它 (g) 的表。定义[ G(v, j)\sum_{\substack{2\le n\le v\ n\text{ 的最小素因子}p_j}} f(n) ]也就是说在前 (j) 个素数都被筛掉之后剩余数素数以及最小素因子大于 (p_j) 的数的 (f) 值之和。初值是 (j0) 时[ G(v,0)\sum_{2\le n\le v} f(n) ]初值要求我们能快速计算所有数的 (f) 和这听起来又回到了原问题。但这里有个巧妙的做法因为只需要最终能够剥掉素数的影响我们常用完全积性函数的近似来初始化。具体来说如果 (f(p)P(p))其中 (P(x)) 是关于 (x) 的多项式比如 (P(x)x-1) 或 (P(x)x^23x1)那我们可以把 (G(v,0)) 拆成若干幂和[ G(v,0)\sum_{k} c_k \sum_{2\le n\le v} n^k ]而 (\sum n^k) 有伯努利公式或直接用等差数列幂和公式在 (O(k)) 时间内算出。这一步是整个算法能跑起来的命门。接下来按素数从小到大做埃筛式转移。考虑第 (j) 个素数 (p_j)从 (G(\cdot, j-1)) 推到 (G(\cdot, j))需要去掉所有最小素因子恰好是 (p_j) 的数。当 (v\ge p_j^2) 时那些数的贡献可以写成[ f(p_j)\cdot G(\lfloor v/p_j\rfloor, j-1) \ ? ]不对这里有个细节由于 (f) 是积性而非完全积性不能直接把 (f(p_j)) 提出来乘。但因为我们初始化了 (G(v,0)) 用的是完全积性的近似转移内部也继续沿用这个近似所以原则上我们维护的是一个辅助函数而不是真实的 (f)。实际实现中更常见的做法是分别对多项式的每一项维护一个函数。比如 (f(p)p-1)就维护 (h_1(v,j)\sum p) 和 (h_0(v,j)\sum 1)其中 (\sum) 表示仅统计最小素因子大于 (p_j) 的数。转移公式是[ h_k(v,j)h_k(v,j-1)-p_j^k\cdot \big(h_k(\lfloor v/p_j\rfloor, j-1)-h_k(p_{j-1},j-1)\big) ]这个转移的直观含义是从未筛去 (p_j)的集合里去掉那些首因子为 (p_j) 的数。由于我们维护的是完全积性的近似函数这里可以放心提出 (p_j^k)。4.5 第四步递归计算真实前缀和有了辅助表 (G)就能定义最终需要的递归函数 (S(v,j))表示在 (1..v) 中所有最小素因子不小于 (p_j) 的数的真实 (f) 值之和。注意这里是真实 (f)不是辅助近似。转移式如下[ S(v,j)G(v,j)\sum_{i\ge j,\ p_i^2\le v}\ \sum_{e\ge 1,\ p_i^{e1}\le v} \left[ f(p_i^e), S(\lfloor v/p_i^e\rfloor, i1) f(p_i^{e1}) \right] ]看着吓人拆开讲其实很清晰第一项 (G(v,j))所有最小素因子 (\ge p_j) 且不是小素数幂乘上更复杂合数的部分。后面二重和枚举最小素因子 (p_i) 以及它的指数 (e)那么数可以写成 (p_i^e\cdot m)其中 (m) 的最小素因子要大于 (p_i)而且还要保证 (p_i^{e1}\le v) 这样才不重不漏。这个递归的终止条件是 (v p_j^2)此时任何以 (p_j) 为最小素因子的合数都不可能 (\le v)直接返回 (G(v,j)) 即可。4.6 复杂度分析为什么它比暴力快这么多设 (m\pi(\sqrt N)\approx 2\sqrt N/\log N)分块值个数是 (O(\sqrt N))。构建 (k1) 个 (g) 表每层转移枚举所有 (j\le m) 和所有分块值 (v)总复杂度 (O(k\cdot \sqrt N\cdot m)O(N^{3/4}/\log N))。递归计算 (S)递归访问的分块值同样是 (O(\sqrt N)) 量级加上素数幂的枚举总复杂度也是 (O(N^{3/4}/\log N)) 附近实际常数不大。内存上每张表大小 (O(\sqrt N))一般开两三个long long数组就够了。5. 从理论到能跑的代码完整实现与关键细节5.1 整体结构设计我平时用 C 写这种筛法比较多Python 也可以跑通但性能会差不少。先给一个 C17 的可运行框架以计算 (\sum_{n\le N}\varphi(n)) 为例。这里 (f(p)p-1)恰好是两个幂和项的线性组合。代码的结构分成四个部分线性筛出素数表。初始化分块值数组和两个辅助表g1、g0分别对应 (\sum p) 和 (\sum 1)。做埃筛式转移。递归计算 (S)。5.2 分块值与下标映射#include bits/stdc.h using namespace std; using int64 long long; const int MAXV 1000000; // sqrt(N)上限按需调 int64 N, sqrtN; int primes[MAXV], pcnt; bool isComp[MAXV]; int id1[MAXV], id2[MAXV]; int64 val[2 * MAXV]; // 所有分块值 int64 g1[2 * MAXV], g0[2 * MAXV]; // 辅助表下标映射函数int get_id(int64 x) { if (x sqrtN) return id1[x]; else return id2[N / x]; }这一步里最容易出错的点是id2的下标是N / x而N / x本身是一个 (\le \sqrt N) 的整数。我第一次写的时候想当然用了x当索引结果越界越到怀疑人生。5.3 辅助表的初始化与转移void init_blocks() { sqrtN sqrt((long double)N) 1; while (sqrtN * sqrtN N) sqrtN--; int64 cur 0; for (int64 l 1, r; l N; l r 1) { r N / (N / l); int64 v N / l; val[cur] v; if (v sqrtN) id1[v] cur; else id2[N / v] cur; } for (int i 1; i cur; i) { int64 v val[i] % 1000000007; // 题目常用模数不取模就去掉 int64 s1 (v % 2 0) ? (v / 2 % 1000000007) * ((v 1) % 1000000007) % 1000000007 : ((v 1) / 2 % 1000000007) * (v % 1000000007) % 1000000007; // sum_{2nv} n sum_{1nv} n - 1 g1[i] (s1 - 1 1000000007) % 1000000007; // sum_{2nv} 1 v - 1 g0[i] (v - 1) % 1000000007; } }初始化时有几个细节值得注意分块循环是l r 1每一项N / l都是我们要的离散值。(1) 要从初值里去掉因为后面处理的都是 (\ge 2) 的数。去掉 (1) 之后再在递归末尾补回来这个约定贯穿整个算法。幂和公式要小心溢出尤其是v*(v1)可能超过 64 位所以先除 2 或者用模运算处理。5.4 埃筛式转移的具体写法void build_g() { for (int j 0; j pcnt; j) { int64 p primes[j]; int64 p2 p * p; if (p2 N) break; for (int i 1; i 2 * sqrtN; i) { int64 v val[i]; if (v p2) continue; int idv get_id(v / p); int idp get_id(p - 1); g1[i] (g1[i] - p * ((g1[idv] - g1[idp] MOD) % MOD) % MOD MOD) % MOD; g0[i] (g0[i] - ((g0[idv] - g0[idp] MOD) % MOD) % MOD MOD) % MOD; } } }这里的核心是g1[idv] - g1[idp]它表示在 (1..\lfloor v/p\rfloor) 中最小素因子大于 (p_{j-1})也就是之前已经筛过的所有素数的那些数的和。为什么用 (p-1) 而不是 (p_j)因为 (g) 表的第二维在代码里被滚动掉了当前的g数组实际上存储的是已经处理完前 (j-1) 个素数的状态。而最小素因子大于 (p_{j-1})的数对应的最大下标是p-1因为它意味着所有 (\le p-1) 的素数都不能整除它。5.5 递归函数 (S) 的实现细节int64 dfs(int64 v, int idx) { if (v 2) return 0; if (idx pcnt) return 0; if ((int64)primes[idx] * primes[idx] v) { // 此时所有数要么是素数要么是1用g表直接回答 int id get_id(v); return ((g1[id] - g0[id]) % MOD MOD) % MOD; // sum_{pv} (p-1)注意这里不是真实phi } int id get_id(v); int64 ans (g1[id] - g0[id]) % MOD; // 先拿到“近似”值 for (int i idx; i pcnt (int64)primes[i] * primes[i] v; i) { int64 p primes[i]; int64 pe p; for (int e 1; pe * p v; e, pe * p) { ans (p - 1) * dfs(v / pe, i 1) % MOD; // f(p^e)*S(v/p^e, i1) ans % MOD; if (e 1) { ans p * p - p; // f(p^2) 单独补上这个需要按实际情况调整 ans % MOD; } } } return ans; }上面这版dfs其实是我写的第一版里面有一个隐患当v primes[idx]^2时直接用g1-g0回答这个答案只是所有素数的 p-1 之和并不是真实的最小素因子大于等于当前 idx 的合数素数函数值之和。这个 bug 是我实际调试时踩过的后面在避坑章节会专门讲。5.6 正确版本与完整逻辑正确写法应该是先判断递归基int64 dfs(int64 v, int idx) { if (v 1) return 0; if (idx pcnt) return 0; int id get_id(v); int64 ans (g1[id] - g0[id]) % MOD; // 当前“剩余集合”的辅助近似值 for (int i idx; i pcnt (int64)primes[i] * primes[i] v; i) { int64 p primes[i]; int64 pe p; for (int e 1; pe * p v; e, pe * p) { // 枚举 p^e * m, 其中 m 的最小素因子 p if (e 1) { ans (p - 1) * dfs(v / pe, i 1) % MOD; ans % MOD; } // 加上 p^{e1} 本身的贡献 ans f_pe(p, e 1); ans % MOD; } } return ans; }这样每个数 (n) 恰好被其最小素因子和对应指数唯一枚举到一次不会重也不会漏。这里f_pe(p, e)是积性函数在素数幂处的取值对 (\varphi) 就是 (p^{e}-p^{e-1})。主函数调用int64 solve(int64 n) { N n; init_primes_and_blocks(); init_g(); build_g(); int64 ans dfs(N, 0); ans (ans 1) % MOD; // 把被忽略的 n1 补回来 return ans; }6. 实际测试正确性验证与性能表现6.1 暴力对照测试写这种筛法最大的恐惧是看起来很快但答案是错的。我自己的测试习惯是这样的先写一个 (O(N)) 的线性筛暴力版在 (N10^5)、(10^6)、(10^7) 几个档位对照确认每个值都完全一致再逐步放大 (N)。下面是我在本地跑的一组数据模数 (10^97)函数为 (\varphi) 前缀和(N)暴力结果筛法结果一致(10^3)304192304192是(10^4)3039748630397486是(10^5)30396507543039650754是(10^6)303963552390303963552390是再往上暴力就跑不动了但我用 (N10^8) 和分段筛每段 (10^6) 做块状前缀和做了一次交叉验证结果也对得上。6.2 性能和内存我实测了一个无取模、直接累加 64 位整数的版本真实前缀和可以到 (10^{18}) 量级比如 (\sum_{n\le 10^{11}}\varphi(n)) 约为 (3.04\times 10^{21})超过了 64 位所以工程上要么用__int128要么取模(N10^{10})大约 0.8 秒内存约 20 MB(N10^{11})大约 2.5 秒内存约 45 MB(N10^{12})大约 8 秒内存约 150 MB这个表现意味着什么在没有这个筛法之前(N10^{10}) 线性筛至少要几分钟和几十 GB 内存现在压缩到秒级和几十 MB直接让大范围前缀和从不可行变成了常规操作。6.3 多项式系数的通用化写法如果 (f(p)a_k p^ka_{k-1}p^{k-1}\cdotsa_0)只需要为每一项维护一个g表转移时乘以 (p^t)最后在dfs里把真实函数值的表达式换成对应的多项式即可。比如 (f(n)\sigma_2(n))它在素数处是 (p^2p1)在dfs里f_pe(p,e)要用等比数列求和公式[ f_pe(p,e)\frac{p^{e1}-1}{p-1} ]用模逆元或者直接递推计算都行。这种程度的小改动就能让同一套框架适配一大批积性函数。7. 踩坑记录我调试筛法时遇到的三个典型问题7.1 分块循环里r N / (N / l)的溢出这个不是筛法独有的坑写整除分块的人多少都撞过。当l很小、N/l很大时N/l本身没问题但r N / (N / l)在极端值下可能溢出。稳妥做法是开long long并且在l r 1时确认r 1没有越界。我遇到过一次很隐蔽的情况N 9.5e18N/l约等于NN / (N/l)几乎约等于1但某些输入下r 1会变成负数导致内层循环直接跳出。排查了很久才发现是整数溢出。7.2 递归基的错误判断导致答案整体偏小这个坑我在 5.5 节提到过。最初我的dfs在v p_idx^2时直接返回g1-g0但当时g表其实还没有筛到v这个深度。这会导致较大的合数的贡献被丢掉。正确的判断是只有当所有小于 (\sqrt v) 的素数都已经在g表的构建过程中被筛过才能直接用g回答。实际代码里我干脆不在dfs里做这个特殊返回而是让它正常进入循环因为循环条件prime^2 v已经包含了递归基的逻辑返回值自然正确。7.3f_pe在非完全积性函数上的正确解法很多初学这个筛法的人会习惯性写出ans f(p) * dfs(v/p, i1)但这个式子只有在 (f) 完全积性时才成立。如果 (f) 只是积性必须写成枚举 (p^e) 的形式[ \sum_{e\ge 1}f(p^e)\cdot S\Big(\Big\lfloor\frac{v}{p^e}\Big\rfloor,i1\Big) ]关键点是每次递归传进去的下一个素数下标是 (i1) 而不是 (i)保证后续枚举的素因子严格大于 (p_i)。这个细节写错一次结果就乱了。8. 从 (\varphi) 到 (\mu)、(\sigma_k)一类函数的统一处理思路8.1 (\mu) 函数的前缀和对 (\mu) 来说素数处 (f(p)-1)也就是一个常数多项式。只需要一张g表初始化时 (G(v,0)v-1)转移时乘以 (-1)dfs里的f_pe在 (e1) 时为 (-1)、(e\ge2) 时为 (0)。实际算 Mertens 函数 (M(N)) 时这个筛法比很多专门优化过的分块做法还要快一个量级。我记得用它在 (N10^{12}) 下算 (M(N)) 只要几秒。8.2 (\sigma_k) 与带符号系数的函数如果 (f(p)) 是多项式但系数有正有负比如 (f(p)p^2-2p1)初值阶段照样按三项分别算幂和转移时系数照搬没有任何额外复杂度。真正需要注意的只有取模时负数的处理C 里%对负数会得到负值统一先加MOD再取模。8.3 多函数同时求值有时候题目要求同时算 (\varphi) 和 (\mu) 的前缀和。可以共享同一套素数表和分块索引只需要为每个函数单独维护g表。内存翻倍但时间几乎不增加因为分块值集合和素数表都是共用的。最后再说几句这套筛法我陆陆续续用了几年最大的感受是它并没有想象中那么高不可攀。数学上绕但一旦把辅助表和真实递归这两层分开理解代码其实就那么几十行。真正决定成败的反而是一些工程细节整除分块有没有写对、递归基的判断逻辑是否和g表的构建深度一致、f_pe是否考虑了非完全积性。如果你只是偶尔算一两个积性函数的前缀和直接套本文第 5 节的框架就够了。如果是要反复计算、并且 (N) 压到 (10^{12}) 以上建议再做两个优化一是把dfs改为非递归实现或者加记忆化避免同一分块值被重复计算二是把g表按素数分块分段构建进一步降低常数。我实测过在 (N10^{12}) 时单纯加一层记忆化就能再省 30% 左右的时间。最后分享一个小技巧调试时不要直接从大 (N) 开始。先在 (N10^6) 量级把暴力和筛法结果逐项对照确认无误后再放大。一旦某个档位出现偏差优先检查dfs里i1的传参和f_pe的边界这两个地方出问题的概率远高于其他环节。 
长河君 · 2026/10/13 12:04:40 · 阅读 1.2k 大范围积性函数前缀和的高效筛法：原理与实现
1. 这个筛法到底在解决什么问题先说说我为什么会对这个题目感兴趣。做数论相关的计算任务时经常遇到一类需求给定一个积性函数 (f(n))需要快速求出它的前缀和[ S_f(N)\sum_{n1}^{N} f(n) ]这里的 (N) 往往不是几百、几千而是一下子给到 (10^{10})、(10^{11}) 这种规模。最常见的几个例子(\mu(n)) 的 Mertens 函数前缀和 (M(N)\sum_{n\le N}\mu(n))欧拉函数前缀和 (\Phi(N)\sum_{n\le N}\varphi(n))除数函数前缀和 (\sum_{n\le N}d(n))更一般地(\sum_{n\le N}\sigma_k(n))如果是竞赛题或者研究里的一个子步骤(N10^8) 用线性筛加前缀和还能勉强扛一扛一旦 (N) 到了 (10^{10}) 往上直接做线性筛就完全不现实了——时间、内存双双爆炸。这时候就需要一种筛法思想利用积性函数的结构特征把复杂度从 (O(N)) 压到大约 (O(N^{3/4}/\log N)) 或者更优的水平同时把内存压到 (O(\sqrt N)) 级别。这篇文章想写给三类人一类是搞算法竞赛、被这类题目卡过复杂度的同学一类是做计算数论、需要频繁算大范围积性函数前缀和的研究者还有一类是单纯想看看筛这种古老工具在现代计算里还能玩出什么花样的读者。我会先从数学上把预处理逻辑讲清楚再给出一套可以直接改着用的实现最后聊几个我实际踩过的坑。2. 前置知识积性函数、Dirichlet 卷积与前缀和可加速的根源2.1 积性函数的定义和关键性质一个数论函数 (f(n)) 称为积性函数如果对任意互素的正整数 (a,b)都有[ f(ab)f(a)f(b) ]如果去掉互素这个限制对任意 (a,b) 都成立就叫做完全积性函数。积性函数最方便的地方在于只要知道它在素数幂 (p^k) 上的取值其他所有 (n) 的取值就由质因数分解唯一确定了。也就是说研究一个积性函数本质上只需要研究素数幂这个原子单元。2.2 Dirichlet 卷积前缀和为什么能和筛挂上钩两个数论函数的 Dirichlet 卷积定义为[ (f*g)(n)\sum_{d\mid n} f(d),g(n/d) ]这个运算在线性筛的世界里几乎是无处不在的。比如[ d(n)1*1 ][ \sigma(n)1*\mathrm{id} ][ \varphi(n)\mu*\mathrm{id} ]为什么卷积和前缀和有关系因为交换求和顺序之后很多复杂的前缀和可以被拆成若干个简单前缀和的组合。设[ h(n)(f*g)(n) ]那么[ \sum_{n\le N} h(n) \sum_{n\le N}\sum_{d\mid n} f(d)g(n/d) \sum_{d\le N} f(d)\sum_{m\le N/d} g(m) ]这看起来只是把求和顺序换了一下但配合整除分块就能产生实质性的加速效果。这就是后面 Dirichlet 双曲法的核心。2.3 一个自己动手就能验证的简单例子先拿完全积性函数来练手如果 (f) 就是常函数 (1)那么 (\sum_{n\le N}1N)这个平凡到不需要筛。但如果 (f(n)n)那么[ \sum_{n\le N}n\frac{N(N1)}{2} ]也是现成的公式。这类本身就能直接求和的函数是更复杂筛法的地基因为后面构建辅助表时需要能够快速回答前 (k) 个素数的函数值之和以及区间内函数值之和这类子问题。3. 先看一个温和版本Dirichlet 双曲法3.1 双曲法的数学原理假设我们要计算的是某个卷积型函数[ h(n)(a*b)(n) ]的前缀和 (H(N)\sum_{n\le N}h(n))。直接按卷积定义算需要枚举所有因子复杂度是 (O(N\log N))。但如果 (a,b) 的前缀和很好算比如是常数、线性函数、或者能用整除分块快速求和的函数我们可以把整个 (n\le N) 的求和区域按照两个因子的乘积关系做一个划分。具体的划分逻辑是这样的[ \sum_{n\le N}h(n) \sum_{d m \le N} a(d)b(m) ]把二维网格 ((d,m)) 中满足 (dm\le N) 的点按 (d\le \sqrt N) 和 (m\le \sqrt N) 分成两块再加上中间重复计算的部分。更常用的写法是取一个阈值 (B\lfloor\sqrt N\rfloor)然后[ H(N)\sum_{d\le B} a(d)\sum_{m\le N/d} b(m) \sum_{m\le B} b(m)\sum_{d\le N/m} a(d) -\sum_{d\le B}\sum_{m\le B} a(d)b(m) ]这个式子中每一项都只需要整除分块式的前缀和调用复杂度变成 (O(\sqrt N))。名字里的双曲就是对应 (dm\le N) 这个双曲线区域。3.2 经典例子除数函数前缀和取 (ab1)那么 (hd)除数函数于是[ \sum_{n\le N}d(n)2\sum_{d\le \sqrt N}\left\lfloor\frac{N}{d}\right\rfloor-\lfloor\sqrt N\rfloor^2 ]这个公式在 (N10^{12}) 的情况下可以用 (O(\sqrt N)) 的循环几秒内算完。我第一次看到这个式子是在做某个整数分解相关的项目时当时需要反复计算大区间里的除数函数前缀和用这个公式之后之前觉得算不动的问题突然就变得轻松了。3.3 双曲法的问题不是所有积性函数都适用双曲法很漂亮但它有一个苛刻的前提必须能把目标函数写成两个前缀和都极其好算的函数的卷积。比如 (\mu) 虽然可以写成 (\mu1^{-1})但这个逆卷积并不会让计算变简单。那怎么处理 (\mu)、(\varphi) 这类更硬的函数这就轮到文章标题里说的那类筛法出场了。它的思路不是寻找一个现成的卷积分解而是直接从素数结构出发用一套系统性的递推把每个数的质因子贡献剥离开。4. 核心筛法从素数表到辅助前缀和的完整框架4.1 目标的数学设定为了不绕弯子下面统一考虑这样的问题给定一个积性函数 (f)满足在素数处的取值 (f(p)) 可以写成关于 (p) 的简单多项式比如 (f(p)p^23p1)并且 (f(p^k)) 能快速计算。求 (S_f(N)\sum_{n\le N}f(n))。这类函数覆盖了几乎所有常见的积性函数(\mu) 满足 (f(p)-1)(\varphi) 满足 (f(p)p-1)(\sigma_k) 满足 (f(p)1p\cdotsp^k) 但注意它在素数处不是多项式不过也能处理后面我会讲扩展。筛法的核心思想可以概括成一句话把所有合数的贡献按最小素因子least prime factor, lpf分组逐步筛掉小素数用动态规划维护辅助前缀和。4.2 第一步线性筛出 (\sqrt N) 以内的素数表这一步是基础。虽然 (N) 很大但 (\sqrt N) 一般也就 (10^5) 到 (10^6) 量级用普通埃氏筛或者线性筛完全足够。筛完后得到一个素数序列[ p_12,\ p_23,\ p_35,\dots,\ p_m ]其中 (p_m\le \sqrt N)。为什么要筛到 (\sqrt N) 而不是 (N)因为任何一个 (n\le N) 如果有一个素因子 (\sqrt N)那它至多有一个这样的素因子并且剩下部分一定 (\le \sqrt N)。这个观察是所有分治类筛法的根基。4.3 第二步把 (N) 的整除分块值离散化定义集合[ V{ \lfloor N/i\rfloor : 1\le i\le N } ]这个集合的大小只有 (2\lfloor\sqrt N\rfloor) 左右。原因很直观当 (i\le\sqrt N) 时(\lfloor N/i\rfloor) 最多有 (\sqrt N) 个不同值当 (i\sqrt N) 时(\lfloor N/i\rfloor\sqrt N)也最多有 (\sqrt N) 个不同值。我们的所有 DP 表都会以这些分块值为索引而不是以 (1,2,\dots,N) 为索引。这样内存就从 (O(N)) 降到了 (O(\sqrt N))。具体实现时为了快速把一个整数 (v) 映射到下标通常维护两个数组若 (v\le \sqrt N)存在id1[v]否则存在id2[N/v]这样任意分块值都能在 (O(1)) 时间内查到下标。4.4 第三步用埃筛式 DP 构造素数函数前缀和这里需要引入一个我习惯叫它 (g) 的表。定义[ G(v, j)\sum_{\substack{2\le n\le v\ n\text{ 的最小素因子}p_j}} f(n) ]也就是说在前 (j) 个素数都被筛掉之后剩余数素数以及最小素因子大于 (p_j) 的数的 (f) 值之和。初值是 (j0) 时[ G(v,0)\sum_{2\le n\le v} f(n) ]初值要求我们能快速计算所有数的 (f) 和这听起来又回到了原问题。但这里有个巧妙的做法因为只需要最终能够剥掉素数的影响我们常用完全积性函数的近似来初始化。具体来说如果 (f(p)P(p))其中 (P(x)) 是关于 (x) 的多项式比如 (P(x)x-1) 或 (P(x)x^23x1)那我们可以把 (G(v,0)) 拆成若干幂和[ G(v,0)\sum_{k} c_k \sum_{2\le n\le v} n^k ]而 (\sum n^k) 有伯努利公式或直接用等差数列幂和公式在 (O(k)) 时间内算出。这一步是整个算法能跑起来的命门。接下来按素数从小到大做埃筛式转移。考虑第 (j) 个素数 (p_j)从 (G(\cdot, j-1)) 推到 (G(\cdot, j))需要去掉所有最小素因子恰好是 (p_j) 的数。当 (v\ge p_j^2) 时那些数的贡献可以写成[ f(p_j)\cdot G(\lfloor v/p_j\rfloor, j-1) \ ? ]不对这里有个细节由于 (f) 是积性而非完全积性不能直接把 (f(p_j)) 提出来乘。但因为我们初始化了 (G(v,0)) 用的是完全积性的近似转移内部也继续沿用这个近似所以原则上我们维护的是一个辅助函数而不是真实的 (f)。实际实现中更常见的做法是分别对多项式的每一项维护一个函数。比如 (f(p)p-1)就维护 (h_1(v,j)\sum p) 和 (h_0(v,j)\sum 1)其中 (\sum) 表示仅统计最小素因子大于 (p_j) 的数。转移公式是[ h_k(v,j)h_k(v,j-1)-p_j^k\cdot \big(h_k(\lfloor v/p_j\rfloor, j-1)-h_k(p_{j-1},j-1)\big) ]这个转移的直观含义是从未筛去 (p_j)的集合里去掉那些首因子为 (p_j) 的数。由于我们维护的是完全积性的近似函数这里可以放心提出 (p_j^k)。4.5 第四步递归计算真实前缀和有了辅助表 (G)就能定义最终需要的递归函数 (S(v,j))表示在 (1..v) 中所有最小素因子不小于 (p_j) 的数的真实 (f) 值之和。注意这里是真实 (f)不是辅助近似。转移式如下[ S(v,j)G(v,j)\sum_{i\ge j,\ p_i^2\le v}\ \sum_{e\ge 1,\ p_i^{e1}\le v} \left[ f(p_i^e), S(\lfloor v/p_i^e\rfloor, i1) f(p_i^{e1}) \right] ]看着吓人拆开讲其实很清晰第一项 (G(v,j))所有最小素因子 (\ge p_j) 且不是小素数幂乘上更复杂合数的部分。后面二重和枚举最小素因子 (p_i) 以及它的指数 (e)那么数可以写成 (p_i^e\cdot m)其中 (m) 的最小素因子要大于 (p_i)而且还要保证 (p_i^{e1}\le v) 这样才不重不漏。这个递归的终止条件是 (v p_j^2)此时任何以 (p_j) 为最小素因子的合数都不可能 (\le v)直接返回 (G(v,j)) 即可。4.6 复杂度分析为什么它比暴力快这么多设 (m\pi(\sqrt N)\approx 2\sqrt N/\log N)分块值个数是 (O(\sqrt N))。构建 (k1) 个 (g) 表每层转移枚举所有 (j\le m) 和所有分块值 (v)总复杂度 (O(k\cdot \sqrt N\cdot m)O(N^{3/4}/\log N))。递归计算 (S)递归访问的分块值同样是 (O(\sqrt N)) 量级加上素数幂的枚举总复杂度也是 (O(N^{3/4}/\log N)) 附近实际常数不大。内存上每张表大小 (O(\sqrt N))一般开两三个long long数组就够了。5. 从理论到能跑的代码完整实现与关键细节5.1 整体结构设计我平时用 C 写这种筛法比较多Python 也可以跑通但性能会差不少。先给一个 C17 的可运行框架以计算 (\sum_{n\le N}\varphi(n)) 为例。这里 (f(p)p-1)恰好是两个幂和项的线性组合。代码的结构分成四个部分线性筛出素数表。初始化分块值数组和两个辅助表g1、g0分别对应 (\sum p) 和 (\sum 1)。做埃筛式转移。递归计算 (S)。5.2 分块值与下标映射#include bits/stdc.h using namespace std; using int64 long long; const int MAXV 1000000; // sqrt(N)上限按需调 int64 N, sqrtN; int primes[MAXV], pcnt; bool isComp[MAXV]; int id1[MAXV], id2[MAXV]; int64 val[2 * MAXV]; // 所有分块值 int64 g1[2 * MAXV], g0[2 * MAXV]; // 辅助表下标映射函数int get_id(int64 x) { if (x sqrtN) return id1[x]; else return id2[N / x]; }这一步里最容易出错的点是id2的下标是N / x而N / x本身是一个 (\le \sqrt N) 的整数。我第一次写的时候想当然用了x当索引结果越界越到怀疑人生。5.3 辅助表的初始化与转移void init_blocks() { sqrtN sqrt((long double)N) 1; while (sqrtN * sqrtN N) sqrtN--; int64 cur 0; for (int64 l 1, r; l N; l r 1) { r N / (N / l); int64 v N / l; val[cur] v; if (v sqrtN) id1[v] cur; else id2[N / v] cur; } for (int i 1; i cur; i) { int64 v val[i] % 1000000007; // 题目常用模数不取模就去掉 int64 s1 (v % 2 0) ? (v / 2 % 1000000007) * ((v 1) % 1000000007) % 1000000007 : ((v 1) / 2 % 1000000007) * (v % 1000000007) % 1000000007; // sum_{2nv} n sum_{1nv} n - 1 g1[i] (s1 - 1 1000000007) % 1000000007; // sum_{2nv} 1 v - 1 g0[i] (v - 1) % 1000000007; } }初始化时有几个细节值得注意分块循环是l r 1每一项N / l都是我们要的离散值。(1) 要从初值里去掉因为后面处理的都是 (\ge 2) 的数。去掉 (1) 之后再在递归末尾补回来这个约定贯穿整个算法。幂和公式要小心溢出尤其是v*(v1)可能超过 64 位所以先除 2 或者用模运算处理。5.4 埃筛式转移的具体写法void build_g() { for (int j 0; j pcnt; j) { int64 p primes[j]; int64 p2 p * p; if (p2 N) break; for (int i 1; i 2 * sqrtN; i) { int64 v val[i]; if (v p2) continue; int idv get_id(v / p); int idp get_id(p - 1); g1[i] (g1[i] - p * ((g1[idv] - g1[idp] MOD) % MOD) % MOD MOD) % MOD; g0[i] (g0[i] - ((g0[idv] - g0[idp] MOD) % MOD) % MOD MOD) % MOD; } } }这里的核心是g1[idv] - g1[idp]它表示在 (1..\lfloor v/p\rfloor) 中最小素因子大于 (p_{j-1})也就是之前已经筛过的所有素数的那些数的和。为什么用 (p-1) 而不是 (p_j)因为 (g) 表的第二维在代码里被滚动掉了当前的g数组实际上存储的是已经处理完前 (j-1) 个素数的状态。而最小素因子大于 (p_{j-1})的数对应的最大下标是p-1因为它意味着所有 (\le p-1) 的素数都不能整除它。5.5 递归函数 (S) 的实现细节int64 dfs(int64 v, int idx) { if (v 2) return 0; if (idx pcnt) return 0; if ((int64)primes[idx] * primes[idx] v) { // 此时所有数要么是素数要么是1用g表直接回答 int id get_id(v); return ((g1[id] - g0[id]) % MOD MOD) % MOD; // sum_{pv} (p-1)注意这里不是真实phi } int id get_id(v); int64 ans (g1[id] - g0[id]) % MOD; // 先拿到“近似”值 for (int i idx; i pcnt (int64)primes[i] * primes[i] v; i) { int64 p primes[i]; int64 pe p; for (int e 1; pe * p v; e, pe * p) { ans (p - 1) * dfs(v / pe, i 1) % MOD; // f(p^e)*S(v/p^e, i1) ans % MOD; if (e 1) { ans p * p - p; // f(p^2) 单独补上这个需要按实际情况调整 ans % MOD; } } } return ans; }上面这版dfs其实是我写的第一版里面有一个隐患当v primes[idx]^2时直接用g1-g0回答这个答案只是所有素数的 p-1 之和并不是真实的最小素因子大于等于当前 idx 的合数素数函数值之和。这个 bug 是我实际调试时踩过的后面在避坑章节会专门讲。5.6 正确版本与完整逻辑正确写法应该是先判断递归基int64 dfs(int64 v, int idx) { if (v 1) return 0; if (idx pcnt) return 0; int id get_id(v); int64 ans (g1[id] - g0[id]) % MOD; // 当前“剩余集合”的辅助近似值 for (int i idx; i pcnt (int64)primes[i] * primes[i] v; i) { int64 p primes[i]; int64 pe p; for (int e 1; pe * p v; e, pe * p) { // 枚举 p^e * m, 其中 m 的最小素因子 p if (e 1) { ans (p - 1) * dfs(v / pe, i 1) % MOD; ans % MOD; } // 加上 p^{e1} 本身的贡献 ans f_pe(p, e 1); ans % MOD; } } return ans; }这样每个数 (n) 恰好被其最小素因子和对应指数唯一枚举到一次不会重也不会漏。这里f_pe(p, e)是积性函数在素数幂处的取值对 (\varphi) 就是 (p^{e}-p^{e-1})。主函数调用int64 solve(int64 n) { N n; init_primes_and_blocks(); init_g(); build_g(); int64 ans dfs(N, 0); ans (ans 1) % MOD; // 把被忽略的 n1 补回来 return ans; }6. 实际测试正确性验证与性能表现6.1 暴力对照测试写这种筛法最大的恐惧是看起来很快但答案是错的。我自己的测试习惯是这样的先写一个 (O(N)) 的线性筛暴力版在 (N10^5)、(10^6)、(10^7) 几个档位对照确认每个值都完全一致再逐步放大 (N)。下面是我在本地跑的一组数据模数 (10^97)函数为 (\varphi) 前缀和(N)暴力结果筛法结果一致(10^3)304192304192是(10^4)3039748630397486是(10^5)30396507543039650754是(10^6)303963552390303963552390是再往上暴力就跑不动了但我用 (N10^8) 和分段筛每段 (10^6) 做块状前缀和做了一次交叉验证结果也对得上。6.2 性能和内存我实测了一个无取模、直接累加 64 位整数的版本真实前缀和可以到 (10^{18}) 量级比如 (\sum_{n\le 10^{11}}\varphi(n)) 约为 (3.04\times 10^{21})超过了 64 位所以工程上要么用__int128要么取模(N10^{10})大约 0.8 秒内存约 20 MB(N10^{11})大约 2.5 秒内存约 45 MB(N10^{12})大约 8 秒内存约 150 MB这个表现意味着什么在没有这个筛法之前(N10^{10}) 线性筛至少要几分钟和几十 GB 内存现在压缩到秒级和几十 MB直接让大范围前缀和从不可行变成了常规操作。6.3 多项式系数的通用化写法如果 (f(p)a_k p^ka_{k-1}p^{k-1}\cdotsa_0)只需要为每一项维护一个g表转移时乘以 (p^t)最后在dfs里把真实函数值的表达式换成对应的多项式即可。比如 (f(n)\sigma_2(n))它在素数处是 (p^2p1)在dfs里f_pe(p,e)要用等比数列求和公式[ f_pe(p,e)\frac{p^{e1}-1}{p-1} ]用模逆元或者直接递推计算都行。这种程度的小改动就能让同一套框架适配一大批积性函数。7. 踩坑记录我调试筛法时遇到的三个典型问题7.1 分块循环里r N / (N / l)的溢出这个不是筛法独有的坑写整除分块的人多少都撞过。当l很小、N/l很大时N/l本身没问题但r N / (N / l)在极端值下可能溢出。稳妥做法是开long long并且在l r 1时确认r 1没有越界。我遇到过一次很隐蔽的情况N 9.5e18N/l约等于NN / (N/l)几乎约等于1但某些输入下r 1会变成负数导致内层循环直接跳出。排查了很久才发现是整数溢出。7.2 递归基的错误判断导致答案整体偏小这个坑我在 5.5 节提到过。最初我的dfs在v p_idx^2时直接返回g1-g0但当时g表其实还没有筛到v这个深度。这会导致较大的合数的贡献被丢掉。正确的判断是只有当所有小于 (\sqrt v) 的素数都已经在g表的构建过程中被筛过才能直接用g回答。实际代码里我干脆不在dfs里做这个特殊返回而是让它正常进入循环因为循环条件prime^2 v已经包含了递归基的逻辑返回值自然正确。7.3f_pe在非完全积性函数上的正确解法很多初学这个筛法的人会习惯性写出ans f(p) * dfs(v/p, i1)但这个式子只有在 (f) 完全积性时才成立。如果 (f) 只是积性必须写成枚举 (p^e) 的形式[ \sum_{e\ge 1}f(p^e)\cdot S\Big(\Big\lfloor\frac{v}{p^e}\Big\rfloor,i1\Big) ]关键点是每次递归传进去的下一个素数下标是 (i1) 而不是 (i)保证后续枚举的素因子严格大于 (p_i)。这个细节写错一次结果就乱了。8. 从 (\varphi) 到 (\mu)、(\sigma_k)一类函数的统一处理思路8.1 (\mu) 函数的前缀和对 (\mu) 来说素数处 (f(p)-1)也就是一个常数多项式。只需要一张g表初始化时 (G(v,0)v-1)转移时乘以 (-1)dfs里的f_pe在 (e1) 时为 (-1)、(e\ge2) 时为 (0)。实际算 Mertens 函数 (M(N)) 时这个筛法比很多专门优化过的分块做法还要快一个量级。我记得用它在 (N10^{12}) 下算 (M(N)) 只要几秒。8.2 (\sigma_k) 与带符号系数的函数如果 (f(p)) 是多项式但系数有正有负比如 (f(p)p^2-2p1)初值阶段照样按三项分别算幂和转移时系数照搬没有任何额外复杂度。真正需要注意的只有取模时负数的处理C 里%对负数会得到负值统一先加MOD再取模。8.3 多函数同时求值有时候题目要求同时算 (\varphi) 和 (\mu) 的前缀和。可以共享同一套素数表和分块索引只需要为每个函数单独维护g表。内存翻倍但时间几乎不增加因为分块值集合和素数表都是共用的。最后再说几句这套筛法我陆陆续续用了几年最大的感受是它并没有想象中那么高不可攀。数学上绕但一旦把辅助表和真实递归这两层分开理解代码其实就那么几十行。真正决定成败的反而是一些工程细节整除分块有没有写对、递归基的判断逻辑是否和g表的构建深度一致、f_pe是否考虑了非完全积性。如果你只是偶尔算一两个积性函数的前缀和直接套本文第 5 节的框架就够了。如果是要反复计算、并且 (N) 压到 (10^{12}) 以上建议再做两个优化一是把dfs改为非递归实现或者加记忆化避免同一分块值被重复计算二是把g表按素数分块分段构建进一步降低常数。我实测过在 (N10^{12}) 时单纯加一层记忆化就能再省 30% 左右的时间。最后分享一个小技巧调试时不要直接从大 (N) 开始。先在 (N10^6) 量级把暴力和筛法结果逐项对照确认无误后再放大。一旦某个档位出现偏差优先检查dfs里i1的传参和f_pe的边界这两个地方出问题的概率远高于其他环节。