OI-wiki 数论专题Dirichlet 双曲线法与杜教筛的亚线性前缀和算法全解析【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki本文基于 OI-wiki 仓库 hyperbola.md 编写完整继承并深化了该文档的全部数学推导、算法分析与例题解答并结合作者仓库中 hyperbola.cpp、fast_bs_conv.cpp、du_1.cpp、du_2.cpp 等参考实现从几何直观 → 算法推导 → 源码级实现 → 竞赛例题四个层次展开。导读在 OI / ICPC 竞赛中形如 $S(n)\sum_{i1}^{n}\mu(i)$、$\sum_{i1}^{n}\varphi(i)$ 的数论函数前缀和问题当 $n$ 高达 $10^{10}\sim 2^{31}$ 时无法用线性筛在可行时间内求解。Dirichlet 双曲线法Dirichlet hyperbola method是这一切亚线性算法的几何基石它把卷积前缀和看成双曲线 $xyn$ 下方整点的加权求和通过容斥把 $O(n)$ 的暴力求和压到 $O(\sqrt n)$。由它出发可以自然引出块筛block sieve及其快速卷积算法并最终推导出竞赛中最常用的杜教筛Du Jiao sieve。读完本文你将掌握双曲线法公式的几何含义与推导、$O(n^{3/4})$ 的朴素块筛卷积、基于点值信息的 $O(n^{2/3}(\log n)^{1/3})$ 优化、周康阳 2024 集训队论文中的 $O(\sqrt n\log^2 n)$ 快速块筛卷积算法以及杜教筛的推导、记忆化实现与三道经典例题的完整解法。前置知识Dirichlet 卷积含 Dirichlet 生成函数、数论分块含关键点集合 $D(n)$ 的性质。Dirichlet 双曲线法设 $f,g,h$ 是数论函数且 $h f\ast g$Dirichlet 卷积。那么利用卷积的定义$h$ 的前缀和$$ H(n) \sum_{k1}^nh(k) \sum_{k1}^n\sum_{xyk}f(x)g(y). $$几何解释求和式遍历的点集恰为第一象限不含坐标轴中双曲线 $xyn$ 下方的整点集合。设整点 $(x,y)$ 的权值为 $f(x)g(y)$那么 $H(n)$ 就是这一权值的和。从几何上看双曲线下方的整点区域是一个曲边三角形。直接枚举全部 $n^2$ 个候选点对显然不划算但若在双曲线上取一个分点 $(x_0,y_0)$就可以用容斥原理把该区域切成两块矩形可计算的部分$$ H(n) \sum_{x1}^{\lfloor x_0\rfloor}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right) \sum_{y1}^{\lfloor y_0\rfloor}F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)g(y) - F(\lfloor x_0\rfloor)G(\lfloor y_0\rfloor). $$其中$F,G$ 分别是 $f,g$ 的前缀和函数$(x_0,y_0)$ 是双曲线 $xyn$ 上任意一个点。表达式中第一项表示图中绿色区域$x\le x_0$ 且 $xy\le n$的权值和第二项表示图中橙色区域$y\le y_0$ 且 $xy\le n$的权值和第三项则是两个区域重叠部分$x\le x_0,y\le y_0$的权值和——这正是容斥原理中多加的要减掉。这个表达式仅含有 $\lfloor x_0\rfloor \lfloor y_0\rfloor 1$ 项注意第三项是标量修正不再需要求和。对于合理选择的 $(x_0,y_0)$它的计算复杂度显著优于暴力计算 $h(n)$ 的前缀和。这就是Dirichlet 双曲线法。卷积前缀和点值的计算Dirichlet 双曲线法最基本的应用就是计算前缀和函数的点值 $H(n)$。如果 $F,G$ 的点值已知或可以在 $O(1)$ 时间内计算进而 $f,g$ 的点值也已知那么双曲线法表达式中的每一项都可以在 $O(1)$ 时间内计算总时间复杂度就等于 $O(x_0y_0)$。因为 $x_0y_0n$所以由均值不等式可知当 $x_0y_0\sqrt{n}$ 时达到最低时间复杂度$$ O(\sqrt{n}). $$与数论分块的等价性这并非新的结果。在双曲线法的表达式中令 $x_0 n$就得到$$ H(n) \sum_{x1}^nf(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right). $$利用数论分块的技巧当 $F,G$ 的点值已知时该式同样可以在 $O(\sqrt{n})$ 时间内计算。细究数论分块的计算过程可以发现它实际计算的表达式为$$ H(n) \sum_{y\in D(n)}\left(F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)-F\left(\left\lfloor\dfrac{n}{y1}\right\rfloor\right)\right)G(y), $$其中 $D(n) \left{\left\lfloor\dfrac{n}{x}\right\rfloor : 1 \le x \le n,~x\in\mathbf N_\right}$ 是数论分块中的关键点集合同时是全体块高与全体块右端点的集合。数论分块的性质表明对 $x\le\sqrt{n}$对应分块高度 $y\lfloor n/x\rfloor$ 互不相同即这些分块长度均为 $1$。把剩余分块的和式做一次Abel 变换分部积分法的求和形式并结合 $\lfloor n/\lfloor\sqrt{n}\rfloor\rfloor \ge \lfloor\sqrt{n}\rfloor$ 分两种情形整理最终可得数论分块的计算过程除了一个 Abel 变换本质上就是在计算$$ H(n) \sum_{x1}^{\lfloor\sqrt{n}\rfloor}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right) \sum_{y1}^{\lfloor\sqrt{n}\rfloor}F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)g(y) - F(\lfloor\sqrt{n}\rfloor)G(\lfloor\sqrt{n}\rfloor). $$这正是 $(x_0,y_0)(\sqrt{n},\sqrt{n})$ 时双曲线法的表达式。因此可以说两种算法的计算过程几乎等价而且由于双曲线法利用了更多数论分块的性质避免了朴素数论分块中的不必要计算常数更小一些。关键观察——只需要稀疏点值在处理实际问题时已知 $F,G$ 的全部点值这一条件可能过强。但观察求和表达式可知其实只需要 $F$ 和 $G$ 在数论分块关键点集合 $D(n)$ 处的取值。根据数论分块的性质$D(n)$ 包含所有 $1\le x\le\sqrt{n}$ 的整数值因此已知 $F,G$ 在 $D(n)$ 处的取值就相当于已知 $f,g$ 在所有 $1\le x\le\sqrt{n}$ 处的取值。而 $|D(n)|\Theta(\sqrt{n})$所以计算 $H(n)$ 时只需要 $F,G$ 在一个稀疏集合处的点值信息。这个观察是优化数论函数前缀和计算的关键。块筛及其卷积有些时候$hf\ast g$ 并非最终需要计算前缀和的函数而只是中间步骤之一。根据前文分析为了后续计算需要求出前缀和函数 $H$ 在集合 $D(n)$ 处的取值。这就称为数论函数 $h$ 的块筛block sieve$$ \mathcal S_h(n) \left{H(x) : x \in D(n)\right}. $$本节讨论块筛卷积问题的计算方法已知 $f,g$ 的块筛时求它们的 Dirichlet 卷积 $h f\ast g$ 的块筛。朴素算法朴素算法就是把块筛的计算看作 $|D(n)|$ 次前缀和点值的计算。总时间复杂度为$$ \begin{aligned} O\left(\sum_{d\in D(n)}\sqrt{d}\right) O\left(\sum_{x1}^{\lfloor\sqrt{n}\rfloor}\sqrt{x} \sum_{x1}^{\lfloor\sqrt{n}\rfloor}\sqrt{\dfrac{n}{x}}\right) \ O\left(\int_1^{\sqrt{n}}\sqrt{x}\mathrm{d}x \int_1^{\sqrt{n}}\sqrt{\dfrac{n}{x}}\mathrm{d}x\right)\ O(n^{3/4}). \end{aligned} $$正是因为块筛是稀疏的只有 $\Theta(\sqrt n)$ 个点整个块筛可以在亚线性时间内求出。但这一算法显然过于暴力集合 $D(n)$ 中较小的那些元素取值相对稠密块筛中相邻两个前缀和相差并不大完全可以直接计算卷积 $h$ 的点值再求它的前缀和——这比对每个点都单独求一遍前缀和点值更快。例如对 $x 1,2,\cdots,\lfloor\sqrt{n}\rfloor$ 分别计算前缀和点值需要 $O\left(\sum_{x1}^{\lfloor\sqrt{n}\rfloor}\sqrt{x}\right) O(n^{3/4})$ 的时间而直接计算 $h$ 在这些点处的点值再累和只需要 $O(n^{1/2}\log n)$ 的时间。不过如果只知道 $f,g$ 的块筛而不知道更多信息这一思路无法继续优化复杂度块筛中只包含 $x\le\sqrt n$ 处的点值信息至多只能计算 $h$ 在 $1\le x\le\sqrt n$ 处的点值剩余前缀和点值的计算仍需要 $O(n^{3/4})$ 的时间。利用点值信息优化如果已知信息除 $f,g$ 的块筛外还包含它们的更多点值那么确实可以改进复杂度。实践中这一算法通常应用于 $f,g$ 的点值可以快速预处理的情形。选择 $z \ge \sqrt{n}$把卷积 $h$ 的块筛分为两部分小值部分计算 $f\ast g$ 的卷积 $h$ 在 $1\le x \le z$ 处的点值再直接累加求和得到 $H$ 在 $1\le x\le z$ 处的点值大值部分对 $x\in D(n)$ 且 $x z$ 的点通过 Dirichlet 双曲线法计算 $H$ 在 $x$ 处的点值。一般情形下时间复杂度为$$ \begin{aligned} O\left(z\log z \sum_{d\in D(n),~d\ge z}\sqrt{d}\right) O\left(z\log z \sum_{x1}^{n/z}\sqrt{\dfrac{n}{x}}\right)\ O\left(z\log z \int_1^{n/z}\sqrt{\dfrac{n}{x}}\mathrm{d}x\right)\ O\left(z\log z \dfrac{n}{\sqrt{z}}\right). \end{aligned} $$当 $z\left(\dfrac{n}{\log n}\right)^{2/3}$ 时总时间复杂度最小为$$ O\left(n^{2/3}(\log n)^{1/3}\right). $$当然Dirichlet 卷积点值计算的复杂度与 $f,g,h$ 的性质有关。对于 $f,g,h$ 有特殊性质的情形最优分点和复杂度均略有不同如果 $f$ 或 $g$ 是积性的那么当 $z\left(\dfrac{n}{\log\log n}\right)^{2/3}$ 时总时间复杂度最小为 $O\left(n^{2/3}(\log\log n)^{1/3}\right)$如果 $h$ 是积性的那么当 $zn^{2/3}$ 时总时间复杂度最小为 $O(n^{2/3})$。应用这一优化并不需要 $f,g$ 的全部点值而只需要它们在 $1\le x\le z$ 处的点值。因为算法同时得到了 $h$ 在 $1\le x\le z$ 处的点值所以当 $h$ 作为中间变量时同样可以利用 $h$ 的点值优化后续计算过程。数论函数的块筛再加上这些点值就构成一个增强版的块筛它们是在 $O(n^{2/3\varepsilon})$ 时间内计算卷积前缀和的全部必要信息。快速块筛卷积前置知识快速傅里叶变换。初学者可以跳过本节。本节讨论周康阳在 2024 年集训队论文《关于积性函数求和问题的一些进展》中提出的快速块筛卷积算法。它可以在 $O(\sqrt{n}\log^2n)$ 时间内根据块筛 $\mathcal S_f$ 和 $\mathcal S_g$ 计算出它们卷积 $\mathcal S_h$ 的取值。这一算法不依赖于额外的点值信息和数论函数的积性但实现较为复杂。块筛卷积问题希望计算 $h(z) \sum_{xyz} f(x)g(y)$ 在块筛 $D(n) {\left\lfloor n / t\right\rfloor : 1\le t \le n}$ 处的前缀和。对于这一问题单一贡献可以由 $(x,y,t)$ 标记——即将项 $f(x)g(y)$ 累加到 $\lfloor n / t\rfloor$ 处前缀和的过程。算法将这些贡献分成若干组处理第一步$x$ 或 $y$ 大于 $\sqrt n$ 的部分考虑 $x \sqrt{n}$ 的所有点对这一系列前缀和的贡献$y \sqrt{n}$ 的贡献类似。因为所有贡献必须满足 $xy \le \lfloor n / t\rfloor$即 $xyt \le n$所以只需要枚举所有可能的 $t,y$利用前缀和技巧以及块筛 $\mathcal S_f$ 中的信息就可以在$$ O\left(\sum_{t,y: ty\le\sqrt{n}} 1\right) O(\sqrt{n}\log n) $$时间内计算出这部分贡献。第二步$\lfloor n/t\rfloor \le \sqrt n$ 的部分考虑 $\lfloor n / t\rfloor\le\sqrt{n}$ 的这部分贡献。这一部分同样可以暴力枚举所有可能的 $x,y$ 完成时间复杂度仍然是$$ O\left(\sum_{x,y:xy\le\sqrt{n}}1\right) O(\sqrt{n}\log n). $$这一部分实际上得到了函数 $h$ 在 $D(n)$ 的前 $\lfloor\sqrt{n}\rfloor$ 个点值。第三步对数分桶近似考虑剩下的贡献即满足 $x,y\le\sqrt{n}$ 且 $\lfloor n / t\rfloor \sqrt{n}$ 的贡献。所有贡献必须满足 $xyt\le n$亦即 $\ln x \ln y \le \ln(n/t)$。取正数 $S$可以利用$$ \lceil S\ln x\rceil \lceil S\ln y\rceil \le S\ln(n/t) $$近似估计这一条件。定义多项式 $\sigma_f(u)$ 和 $\sigma_g(u)$使得系数 $[u^k]\sigma_f$ 与 $[u^k]\sigma_g$ 分别等于满足 $\lceil S\ln x\rceil k$ 时 $f(x)$ 的和、满足 $\lceil S\ln y\rceil k$ 时 $g(y)$ 的和只考虑 $x,y\le\sqrt{n}$ 的部分。利用快速傅里叶变换FFT得到乘积 $\sigma_f\sigma_g$其系数 $u^k$ 表示 $\lceil S\ln x\rceil \lceil S\ln y\rceil k$ 时 $f(x)g(y)$ 的和。由此只需对 $D(n)$ 中剩下每个 $\lfloor n/t\rfloor$ 找到满足 $k \le S\ln(n/t)$ 的最大 $k$ 值即可得到这一部分贡献的估计值。第四步误差修正近似条件可能遗漏部分贡献这只会发生在$$ S\ln x 1 S\ln y 1 \ge \lceil S\ln x\rceil \lceil S\ln y\rceil S\ln(n/t) \ge S\ln x S\ln y $$时等价于 $xyt \in (n\mathrm{e}^{-2/S},n]$这是一个长度为 $O(n/S)$ 的区间。枚举区间内所有可能的贡献 $(x,y,t)$逐个检验是否遗漏即可完成误差修正。为了快速枚举区间内所有贡献可以先筛出不超过 $\sqrt{n}$ 的全部素数用这些素数去除区间中的整数剩下的因子必然是大于 $\sqrt{n}$ 的素数由此得到区间内所有整数的素因数分解进而快速枚举所有可能的 $(x,y,t)$。复杂度分析估计贡献时需要对长度为 $S\log n$ 的多项式做乘法时间复杂度 $O(S\log n\log(S\log n))$误差修正时预处理素因数分解的时间复杂度为 $O(\sqrt{n}(n/S)\log\log n)$枚举区间内所有贡献的时间复杂度为 $O(\sum_{k\in(n\mathrm{e}^{-2/S},n]}d_3(k))$此处 $d_3(n)$ 表示将 $n$ 分解成三个有序整数乘积的方法数。解析数论的Piltz 除数问题结果指出$$ \sum_{k\le n}d_3(k) nP(\log n) O(n^{43/96\varepsilon}), $$其中 $P(\cdot)$ 是二次多项式。前两步时间复杂度已是 $O(\sqrt{n}\log n)$忽略所有 $o(\sqrt{n}\log n)$ 的项最后一部分贡献计算的时间复杂度为$$ O\left(S\log n\log(S\log n) \dfrac{n}{S} \log^2n\right). $$取 $S \sqrt{n}$就得到总时间复杂度$$ O(\sqrt{n}\log^2n). $$参考实现仓库 fast_bs_conv.cpp 提供了该算法的完整 C 实现。其核心结构BlockSieve用两个数组维护块筛struct BlockSieve { long long n, b; std::vectorint s1, s2; BlockSieve(long long _n) : n(_n), b(std::sqrt(_n 0.25l)), s1(b 1), s2(b 1) {} int operator[](long long x) { return x b ? s1[x] : s2[n / x]; } };operator[]把 $x\le b\lfloor\sqrt n\rfloor$ 的下标映射到稠密数组s1其余下标按 $\lfloor n/x\rfloor$ 映射到稀疏数组s2这正是 $D(n)$ 的两段式结构。主函数block_sieve_convolute依次实现前缀和差分求点值df/dg、$x$ 或 $y\sqrt n$ 部分的前缀和技巧、$xy\le b$ 的暴力枚举、按 $\lceil b\ln x\rceil$ 分桶后用 NTTntt_mul模 $998244353$原根 $3$做多项式乘法、以及利用素数筛分解区间 $(n\mathrm e^{-2/b},n]$ 内整数并 DFS 枚举 $(x,y,t)$ 的误差修正。文件末尾的test函数用随机数据对照双曲线法逐一验证 $D(n)$ 中每个点的取值输出Correct/Incorrect可以作为算法的正确性测试模板。杜教筛前文讨论了如何计算数论函数 Dirichlet 卷积的前缀和。本节考虑它的逆过程设 $f \ast g h$且 $f,h$ 已知计算 $g$ 的前缀和$$ G(n) \sum_{x1}^ng(x). $$换句话说本节考虑两个数论函数在 Dirichlet 卷积意义下的商的前缀和计算。本节总是假设 $f(1)\neq 0$以保证 $f$ 可逆Dirichlet 卷积逆元存在的充要条件见 dirichlet.md。为此在 Dirichlet 双曲线法表达式中令 $x_0 n$就得到$$ H(n) \sum_{x1}^{n}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right). $$直接解出 $G(n)$就得到杜教筛的表达式$$ G(n) \dfrac{1}{f(1)}\left(H(n)-\sum_{x2}^nf(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right)\right). $$实际上对于 $x_0 \ge 1$总有更一般的带分点形式$$ G(n) \dfrac{1}{f(1)}\left(H(n)-\sum_{x2}^{\lfloor x_0\rfloor}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right) - \sum_{y1}^{\lfloor y_0\rfloor}F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)g(y) F(\lfloor x_0\rfloor)G(\lfloor y_0\rfloor)\right). $$无论是哪种形式它都是一个关于 $G(n)$ 的递推关系式。为计算 $G(n)$ 的取值需要计算 $G$ 在 $D(n)\setminus{n}$ 处的取值。因为 $D(n)$ 具有递归结构——即对于 $m\in D(n)$总是有 $D(m)\subseteq D(n)$——所以在整个递归计算过程中只需要计算 $G$ 在 $D(n)$ 中元素处的取值各一次。换句话说计算 $G(n)$ 时实际上得到了 $g$ 的块筛 $\mathcal S_g(n)$。实现方式具体实现时可以采用递归方法并配合记忆化避免重复计算也可以采用迭代方法从小到大依次计算 $D(n)$ 中每个点处 $G$ 的取值。此时表达式中的求和式既可以用数论分块计算也可以用 Dirichlet 双曲线法计算。这些实现的复杂度是相同的。复杂度由于杜教筛总是得到块筛所以杜教筛的复杂度其实相当于计算块筛的复杂度。如果已知信息只有 $F,H$ 的块筛那么杜教筛的复杂度就是 $O(n^{3/4})$如果对某个 $z\ge\sqrt{n}$可以在 $T_0(z)$ 时间内预处理出 $g$ 在 $1\le x \le z$ 处的点值那么杜教筛的复杂度就是$$ O\left(T_0(z) \dfrac{n}{\sqrt{z}}\right). $$当 $g$ 是积性函数时可以应用线性筛即 $T_0(z)\Theta(z)$所以最优需要预处理到 $z n^{2/3}$ 处总时间复杂度为 $O(n^{2/3})$对于更一般的情形总时间复杂度则为 $O(n^{2/3}(\log n)^{1/3})$。这些都和块筛部分的分析完全一致。警告递归实现时不使用记忆化将导致复杂度错误。杜教筛表达式中计算 $G(n)$ 依赖 $D(n)\setminus{n}$ 中 $G$ 的取值保证复杂度的关键在于 $D(n)$ 的递归结构当 $m\in D(n)$ 时$D(m)\subseteq D(n)$。利用记忆化后算法复杂度为 $O(n^{3/4})$而不使用记忆化时设计算 $G(n)$ 的复杂度为 $T(n)$有$$ \begin{aligned} T(n) \Theta(\sqrt{n}) \sum_{d\in D(n),~d\neq n}T(d)\ \Theta(\sqrt{n}) \sum_{x1}^{\lfloor n/\lfloor\sqrt{n}\rfloor\rfloor - 1} T(x) \sum_{x 2}^{\lfloor\sqrt{n}\rfloor}T\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right). \end{aligned} $$利用类似主定理的证明思路可以说明最后一项主导该式的增长且 $T(n)\in\Theta(n^\alpha)$其中 $\alpha\approx 1.73$ 是 $\zeta(\alpha)2$ 的根——远远劣于记忆化版本的 $O(n^{3/4})$。应用要点应用杜教筛计算数论函数 $g$ 前缀和时关键在于找到合适的 $f,h$ 使得 $hf\ast g$ 且 $f,h$ 的块筛都容易计算。有些时候这样的 $f,h$ 是显然的如 $\varepsilon\mu\ast 1$、$\mathrm{id}\varphi\ast 1$另一些时候需要利用 Dirichlet 卷积的性质或通过计算相应的 Dirichlet 生成函数来找到相应分解。后文的例题展示了这些情形。例题例 1ARC116 C - Multiple Sequences给定正整数 $N$ 和 $M$计算有多少长度为 $N$ 的序列 $A$ 满足 $1 \le A_i \le M$ 且 $A_i$ 整除 $A_{i1}$。答案对 $998244353$ 取模。数据范围$1 \le N, M \le 2\times 10^5$。解答设长度为 $n$ 且 $A_nm$ 的序列数目为 $f_n(m)$那么最终答案就是 $\sum_{m1}^M f_N(m)$。考虑动态规划转移方程为$$ f_n(m) \sum_{k\mid m}f_{n-1}(k). $$利用 Dirichlet 卷积记号可记作 $f_n f_{n-1}\ast 1$其中 $1$ 是常值数论函数。注意到 $f_1 1$归纳可知 $f_n 1^{\ast n}$即 $f_n$ 是 $n$ 个常值函数的卷积。最终答案就是 $f_N$ 的前缀和。由于过程中只涉及积性函数利用前文介绍的 Dirichlet 卷积前缀和计算方法单次卷积前缀和的计算只需 $O(M^{2/3})$ 时间再利用快速幂的思想只需计算 $O(\log N)$ 次卷积前缀和即可得到 $f_N$ 的卷积前缀和的值。整体时间复杂度为 $O(M^{2/3}\log N)$。参考实现仓库 hyperbola.cpp 是该题的完整解法其结构恰好体现了块筛 点值信息 快速幂的组合dirichlet_convolute在假设 $h$ 是积性函数的前提下用线性筛框架计算卷积 $hf\ast g$ 的前缀点值——素数次幂处直接按定义求和rem[x] 1分支合数处利用积性拆成互素因子相乘h[x] h[rem[x]] * h[x / rem[x]]BlockSieve结构体维护 $n,z$ 以及三组数据 $f,F,F_2$其中 $F$ 保存 $1\le x\le z$ 处前缀和$F_2$ 保存 $xz$ 处按 $\lfloor n/x\rfloor$ 索引的前缀和sum(x)在两者间切换operator*实现两个块筛的卷积——对 $1\le i\le z$ 直接卷积累加对 $iz$ 用双曲线法公式枚举 $x\le\sqrt{k}$ 的 $f(x)G(k/x)F(k/x)g(x)$ 再减去重叠修正 $F(\sqrt{k})G(\sqrt{k})$main取 $zn^{2/3}$初始化常值函数 $1$ 的块筛po对卷积做二进制快速幂res res * po; po po * po;最终输出res.sum(n)。仓库还提供了该题的测试数据hyperbola.in内容为200000 200000与对应的标准答案 hyperbola.ans835917264可用于验证实现。例 2P4213【模板】杜教筛Sum设 $\mu$ 和 $\varphi$ 分别是莫比乌斯函数和欧拉函数。求 $S_1(n) \sum_{i1}^{n} \mu(i)$ 和 $S_2(n) \sum_{i1}^{n} \varphi(i)$ 的值。数据范围$1\leq n2^{31}$。解答注意到两个 Dirichlet 卷积关系式$$ \varepsilon \mu \ast 1,~ \operatorname{id} \varphi \ast 1. $$其中 $\varepsilon(n) [n1]$ 是 Dirichlet 卷积的单位元函数$\operatorname{id}(n) n$ 是恒等函数$1(n)1$ 是常值函数。这三个函数的前缀和都可以在 $O(1)$ 时间内计算且都是积性函数所以利用前文介绍的杜教筛方法可以在 $O(n^{2/3})$ 时间内计算。对于欧拉函数前缀和另一种方法是利用莫比乌斯反演$$ \begin{aligned} S_2(n) \sum_{i1}^n\varphi(i) \sum_{i1}^n\sum_{j1}^i[i\perp j] \ \sum_{i1}^n\sum_{j1}^i\sum_d\mu(d)[d\mid i][d\mid j] \ \sum_d\mu(d)\dfrac{1}{2}\left\lfloor\dfrac{n}{d}\right\rfloor\left(\left\lfloor\dfrac{n}{d}\right\rfloor1\right). \end{aligned} $$在数论分块过程中需要 $\mu(d)$ 的前缀和而这可以通过杜教筛预处理出来时间复杂度仍然是 $O(n^{2/3})$。参考实现仓库 du_1.cpp 给出模板实现其中杜教筛的记忆化递归写得非常清晰long long S_mu(long long x) { // 求mu的前缀和 if (x MAXN) return sum_mu[x]; if (mp_mu[x]) return mp_mu[x]; // 如果map中已有该大小的mu值则可直接返回 long long ret 1; for (long long i 2, j; i x; i j 1) { j x / (x / i); ret - S_mu(x / i) * (j - i 1); } return mp_mu[x] ret; // 路径压缩方便下次计算 }实现要点对 $x MAXN$ 直接查线性筛预处理的前缀和sum_mu线性筛部分同时求出了 $\mu$ 数组对更大的 $x$利用 $\varepsilon\mu\ast 1$ 即 $S_\mu(n)1-\sum_{i2}^{n}S_\mu(\lfloor n/i\rfloor)$配合数论分块j x / (x / i)递归求解并用map记忆化S_phi则基于 $\mathrm{id}\varphi\ast 1$ 或上述莫比乌斯反演式在数论分块中复用S_mu。main里 $T$ 组询问每组输出S_phi(n) S_mu(n)。例 3P3768 简单的数学题给定 $p,n$计算 $\sum_{i1}^n\sum_{j1}^nij\cdot\gcd(i,j)\pmod p$。数据范围$n\leq 10^{10}$$5\times 10^8\leq p\leq 1.1\times 10^9$ 且 $p$ 是质数。解答利用欧拉函数的性质做反演$$ \begin{aligned} T(n) \sum_{i1}^n\sum_{j1}^nij\cdot\gcd(i,j)\ \sum_{i1}^n\sum_{j1}^nij\sum_d\varphi(d)[d\mid i][d\mid j]\ \sum_d\varphi(d)\left(\sum_{i1}^{\lfloor n/d\rfloor}id\right)\left(\sum_{j1}^{\lfloor n/d\rfloor}jd\right)\ \sum_d d^2\varphi(d) F\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)^2. \end{aligned} $$其中 $F(n) \dfrac{1}{2}n(n1)$。该式可以通过数论分块计算但需要预处理 $d^2\varphi(d)$ 的前缀和。为此利用杜教筛。记 $f(n)(\operatorname{id}^2\varphi)(n)$ 和 $S(n)\sum_{i1}^n f(i)$。应用杜教筛的关键是构造函数 $g$使得 $f\ast g$ 与 $g$ 都可以快速求和。前文已讨论用杜教筛预处理 $\varphi$ 前缀和的方法它利用关系 $\mathrm{id} \varphi\ast 1$。相较于 $\varphi$这里的 $f$ 多了一个 $\operatorname{id}^2$由于 $\operatorname{id}$ 是完全积性函数利用 Dirichlet 卷积的性质只需将每一项都乘以 $\operatorname{id}^2$就得到$$ \operatorname{id}^3 f \ast \operatorname{id}^2. $$因为 $\operatorname{id}^2(n)n^2$ 与 $\operatorname{id}^3(n)n^3$ 的前缀和都可以在 $O(1)$ 时间内计算$f$ 的前缀和就可以在 $O(n^{2/3})$ 时间内预处理得到再加上数论分块整体时间复杂度仍是 $O(n^{2/3})$。另一种将类似积性函数表示为两函数之在 Dirichlet 卷积意义下的商的方法是利用 Dirichlet 生成函数。参考实现仓库 du_2.cpp 给出模板实现。它的线性筛到 $n^{2/3}$pn (long long)pow(n, 0.666667); prime_work(pn);正是 $zn^{2/3}$ 最优分点的代码体现prime_work用线性筛求出 $\varphi$ 并预计算 $s[i]i^2\varphi(i)$ 的前缀和模 $P$s2/s3给出平方和与立方和用费马小定理ksm求 $6$、$2$ 的逆元calc(k)利用 $\mathrm{id}^3f\ast\mathrm{id}^2$ 做记忆化递归s_map存储超过预处理上限的答案求和式用数论分块最终solve对外层和式 $\sum_d d^2\varphi(d)F(\lfloor n/d\rfloor)^2$ 做数论分块并调用calc。文件首行注释不要为了省什么内存把数组开小,会卡80也提醒了 $N5e6$ 这类数组规模的选型理由。习题AtCoder Xmas Contest 2019 D - Sum of (-1)^f(n)需要综合利用本文学到的块筛与卷积前缀和技术可以结合 dirichlet.md 与 mobius.md 阅读。参考资料与注释任之洲. 2016. 《积性函数求和的几种方法》. 2016 年信息学奥林匹克中国国家队候选队员论文.周康阳. 2024. 《关于积性函数求和问题的一些进展》. 2024 年信息学奥林匹克中国国家队候选队员论文. 快速块筛卷积算法出处数论分块的完整性质与 $D(n)$ 的递归结构见 sqrt-decomposition.mdDirichlet 卷积的代数性质、逆元存在条件与 Dirichlet 生成函数见 dirichlet.md快速块筛卷积的 NTT 部分可参考仓库 fast_bs_conv.cpp 中ntt_mul的 DIF/DIT 实现关于 Piltz 除数问题$d_3(n)$ 的均值公式可查阅除数求和函数相关的解析数论文献。【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考 
长河君 · 2026/9/13 19:49:52 · 阅读 1.2k OI-wiki 数论专题：Dirichlet 双曲线法与杜教筛的亚线性前缀和算法全解析
OI-wiki 数论专题Dirichlet 双曲线法与杜教筛的亚线性前缀和算法全解析【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki本文基于 OI-wiki 仓库 hyperbola.md 编写完整继承并深化了该文档的全部数学推导、算法分析与例题解答并结合作者仓库中 hyperbola.cpp、fast_bs_conv.cpp、du_1.cpp、du_2.cpp 等参考实现从几何直观 → 算法推导 → 源码级实现 → 竞赛例题四个层次展开。导读在 OI / ICPC 竞赛中形如 $S(n)\sum_{i1}^{n}\mu(i)$、$\sum_{i1}^{n}\varphi(i)$ 的数论函数前缀和问题当 $n$ 高达 $10^{10}\sim 2^{31}$ 时无法用线性筛在可行时间内求解。Dirichlet 双曲线法Dirichlet hyperbola method是这一切亚线性算法的几何基石它把卷积前缀和看成双曲线 $xyn$ 下方整点的加权求和通过容斥把 $O(n)$ 的暴力求和压到 $O(\sqrt n)$。由它出发可以自然引出块筛block sieve及其快速卷积算法并最终推导出竞赛中最常用的杜教筛Du Jiao sieve。读完本文你将掌握双曲线法公式的几何含义与推导、$O(n^{3/4})$ 的朴素块筛卷积、基于点值信息的 $O(n^{2/3}(\log n)^{1/3})$ 优化、周康阳 2024 集训队论文中的 $O(\sqrt n\log^2 n)$ 快速块筛卷积算法以及杜教筛的推导、记忆化实现与三道经典例题的完整解法。前置知识Dirichlet 卷积含 Dirichlet 生成函数、数论分块含关键点集合 $D(n)$ 的性质。Dirichlet 双曲线法设 $f,g,h$ 是数论函数且 $h f\ast g$Dirichlet 卷积。那么利用卷积的定义$h$ 的前缀和$$ H(n) \sum_{k1}^nh(k) \sum_{k1}^n\sum_{xyk}f(x)g(y). $$几何解释求和式遍历的点集恰为第一象限不含坐标轴中双曲线 $xyn$ 下方的整点集合。设整点 $(x,y)$ 的权值为 $f(x)g(y)$那么 $H(n)$ 就是这一权值的和。从几何上看双曲线下方的整点区域是一个曲边三角形。直接枚举全部 $n^2$ 个候选点对显然不划算但若在双曲线上取一个分点 $(x_0,y_0)$就可以用容斥原理把该区域切成两块矩形可计算的部分$$ H(n) \sum_{x1}^{\lfloor x_0\rfloor}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right) \sum_{y1}^{\lfloor y_0\rfloor}F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)g(y) - F(\lfloor x_0\rfloor)G(\lfloor y_0\rfloor). $$其中$F,G$ 分别是 $f,g$ 的前缀和函数$(x_0,y_0)$ 是双曲线 $xyn$ 上任意一个点。表达式中第一项表示图中绿色区域$x\le x_0$ 且 $xy\le n$的权值和第二项表示图中橙色区域$y\le y_0$ 且 $xy\le n$的权值和第三项则是两个区域重叠部分$x\le x_0,y\le y_0$的权值和——这正是容斥原理中多加的要减掉。这个表达式仅含有 $\lfloor x_0\rfloor \lfloor y_0\rfloor 1$ 项注意第三项是标量修正不再需要求和。对于合理选择的 $(x_0,y_0)$它的计算复杂度显著优于暴力计算 $h(n)$ 的前缀和。这就是Dirichlet 双曲线法。卷积前缀和点值的计算Dirichlet 双曲线法最基本的应用就是计算前缀和函数的点值 $H(n)$。如果 $F,G$ 的点值已知或可以在 $O(1)$ 时间内计算进而 $f,g$ 的点值也已知那么双曲线法表达式中的每一项都可以在 $O(1)$ 时间内计算总时间复杂度就等于 $O(x_0y_0)$。因为 $x_0y_0n$所以由均值不等式可知当 $x_0y_0\sqrt{n}$ 时达到最低时间复杂度$$ O(\sqrt{n}). $$与数论分块的等价性这并非新的结果。在双曲线法的表达式中令 $x_0 n$就得到$$ H(n) \sum_{x1}^nf(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right). $$利用数论分块的技巧当 $F,G$ 的点值已知时该式同样可以在 $O(\sqrt{n})$ 时间内计算。细究数论分块的计算过程可以发现它实际计算的表达式为$$ H(n) \sum_{y\in D(n)}\left(F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)-F\left(\left\lfloor\dfrac{n}{y1}\right\rfloor\right)\right)G(y), $$其中 $D(n) \left{\left\lfloor\dfrac{n}{x}\right\rfloor : 1 \le x \le n,~x\in\mathbf N_\right}$ 是数论分块中的关键点集合同时是全体块高与全体块右端点的集合。数论分块的性质表明对 $x\le\sqrt{n}$对应分块高度 $y\lfloor n/x\rfloor$ 互不相同即这些分块长度均为 $1$。把剩余分块的和式做一次Abel 变换分部积分法的求和形式并结合 $\lfloor n/\lfloor\sqrt{n}\rfloor\rfloor \ge \lfloor\sqrt{n}\rfloor$ 分两种情形整理最终可得数论分块的计算过程除了一个 Abel 变换本质上就是在计算$$ H(n) \sum_{x1}^{\lfloor\sqrt{n}\rfloor}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right) \sum_{y1}^{\lfloor\sqrt{n}\rfloor}F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)g(y) - F(\lfloor\sqrt{n}\rfloor)G(\lfloor\sqrt{n}\rfloor). $$这正是 $(x_0,y_0)(\sqrt{n},\sqrt{n})$ 时双曲线法的表达式。因此可以说两种算法的计算过程几乎等价而且由于双曲线法利用了更多数论分块的性质避免了朴素数论分块中的不必要计算常数更小一些。关键观察——只需要稀疏点值在处理实际问题时已知 $F,G$ 的全部点值这一条件可能过强。但观察求和表达式可知其实只需要 $F$ 和 $G$ 在数论分块关键点集合 $D(n)$ 处的取值。根据数论分块的性质$D(n)$ 包含所有 $1\le x\le\sqrt{n}$ 的整数值因此已知 $F,G$ 在 $D(n)$ 处的取值就相当于已知 $f,g$ 在所有 $1\le x\le\sqrt{n}$ 处的取值。而 $|D(n)|\Theta(\sqrt{n})$所以计算 $H(n)$ 时只需要 $F,G$ 在一个稀疏集合处的点值信息。这个观察是优化数论函数前缀和计算的关键。块筛及其卷积有些时候$hf\ast g$ 并非最终需要计算前缀和的函数而只是中间步骤之一。根据前文分析为了后续计算需要求出前缀和函数 $H$ 在集合 $D(n)$ 处的取值。这就称为数论函数 $h$ 的块筛block sieve$$ \mathcal S_h(n) \left{H(x) : x \in D(n)\right}. $$本节讨论块筛卷积问题的计算方法已知 $f,g$ 的块筛时求它们的 Dirichlet 卷积 $h f\ast g$ 的块筛。朴素算法朴素算法就是把块筛的计算看作 $|D(n)|$ 次前缀和点值的计算。总时间复杂度为$$ \begin{aligned} O\left(\sum_{d\in D(n)}\sqrt{d}\right) O\left(\sum_{x1}^{\lfloor\sqrt{n}\rfloor}\sqrt{x} \sum_{x1}^{\lfloor\sqrt{n}\rfloor}\sqrt{\dfrac{n}{x}}\right) \ O\left(\int_1^{\sqrt{n}}\sqrt{x}\mathrm{d}x \int_1^{\sqrt{n}}\sqrt{\dfrac{n}{x}}\mathrm{d}x\right)\ O(n^{3/4}). \end{aligned} $$正是因为块筛是稀疏的只有 $\Theta(\sqrt n)$ 个点整个块筛可以在亚线性时间内求出。但这一算法显然过于暴力集合 $D(n)$ 中较小的那些元素取值相对稠密块筛中相邻两个前缀和相差并不大完全可以直接计算卷积 $h$ 的点值再求它的前缀和——这比对每个点都单独求一遍前缀和点值更快。例如对 $x 1,2,\cdots,\lfloor\sqrt{n}\rfloor$ 分别计算前缀和点值需要 $O\left(\sum_{x1}^{\lfloor\sqrt{n}\rfloor}\sqrt{x}\right) O(n^{3/4})$ 的时间而直接计算 $h$ 在这些点处的点值再累和只需要 $O(n^{1/2}\log n)$ 的时间。不过如果只知道 $f,g$ 的块筛而不知道更多信息这一思路无法继续优化复杂度块筛中只包含 $x\le\sqrt n$ 处的点值信息至多只能计算 $h$ 在 $1\le x\le\sqrt n$ 处的点值剩余前缀和点值的计算仍需要 $O(n^{3/4})$ 的时间。利用点值信息优化如果已知信息除 $f,g$ 的块筛外还包含它们的更多点值那么确实可以改进复杂度。实践中这一算法通常应用于 $f,g$ 的点值可以快速预处理的情形。选择 $z \ge \sqrt{n}$把卷积 $h$ 的块筛分为两部分小值部分计算 $f\ast g$ 的卷积 $h$ 在 $1\le x \le z$ 处的点值再直接累加求和得到 $H$ 在 $1\le x\le z$ 处的点值大值部分对 $x\in D(n)$ 且 $x z$ 的点通过 Dirichlet 双曲线法计算 $H$ 在 $x$ 处的点值。一般情形下时间复杂度为$$ \begin{aligned} O\left(z\log z \sum_{d\in D(n),~d\ge z}\sqrt{d}\right) O\left(z\log z \sum_{x1}^{n/z}\sqrt{\dfrac{n}{x}}\right)\ O\left(z\log z \int_1^{n/z}\sqrt{\dfrac{n}{x}}\mathrm{d}x\right)\ O\left(z\log z \dfrac{n}{\sqrt{z}}\right). \end{aligned} $$当 $z\left(\dfrac{n}{\log n}\right)^{2/3}$ 时总时间复杂度最小为$$ O\left(n^{2/3}(\log n)^{1/3}\right). $$当然Dirichlet 卷积点值计算的复杂度与 $f,g,h$ 的性质有关。对于 $f,g,h$ 有特殊性质的情形最优分点和复杂度均略有不同如果 $f$ 或 $g$ 是积性的那么当 $z\left(\dfrac{n}{\log\log n}\right)^{2/3}$ 时总时间复杂度最小为 $O\left(n^{2/3}(\log\log n)^{1/3}\right)$如果 $h$ 是积性的那么当 $zn^{2/3}$ 时总时间复杂度最小为 $O(n^{2/3})$。应用这一优化并不需要 $f,g$ 的全部点值而只需要它们在 $1\le x\le z$ 处的点值。因为算法同时得到了 $h$ 在 $1\le x\le z$ 处的点值所以当 $h$ 作为中间变量时同样可以利用 $h$ 的点值优化后续计算过程。数论函数的块筛再加上这些点值就构成一个增强版的块筛它们是在 $O(n^{2/3\varepsilon})$ 时间内计算卷积前缀和的全部必要信息。快速块筛卷积前置知识快速傅里叶变换。初学者可以跳过本节。本节讨论周康阳在 2024 年集训队论文《关于积性函数求和问题的一些进展》中提出的快速块筛卷积算法。它可以在 $O(\sqrt{n}\log^2n)$ 时间内根据块筛 $\mathcal S_f$ 和 $\mathcal S_g$ 计算出它们卷积 $\mathcal S_h$ 的取值。这一算法不依赖于额外的点值信息和数论函数的积性但实现较为复杂。块筛卷积问题希望计算 $h(z) \sum_{xyz} f(x)g(y)$ 在块筛 $D(n) {\left\lfloor n / t\right\rfloor : 1\le t \le n}$ 处的前缀和。对于这一问题单一贡献可以由 $(x,y,t)$ 标记——即将项 $f(x)g(y)$ 累加到 $\lfloor n / t\rfloor$ 处前缀和的过程。算法将这些贡献分成若干组处理第一步$x$ 或 $y$ 大于 $\sqrt n$ 的部分考虑 $x \sqrt{n}$ 的所有点对这一系列前缀和的贡献$y \sqrt{n}$ 的贡献类似。因为所有贡献必须满足 $xy \le \lfloor n / t\rfloor$即 $xyt \le n$所以只需要枚举所有可能的 $t,y$利用前缀和技巧以及块筛 $\mathcal S_f$ 中的信息就可以在$$ O\left(\sum_{t,y: ty\le\sqrt{n}} 1\right) O(\sqrt{n}\log n) $$时间内计算出这部分贡献。第二步$\lfloor n/t\rfloor \le \sqrt n$ 的部分考虑 $\lfloor n / t\rfloor\le\sqrt{n}$ 的这部分贡献。这一部分同样可以暴力枚举所有可能的 $x,y$ 完成时间复杂度仍然是$$ O\left(\sum_{x,y:xy\le\sqrt{n}}1\right) O(\sqrt{n}\log n). $$这一部分实际上得到了函数 $h$ 在 $D(n)$ 的前 $\lfloor\sqrt{n}\rfloor$ 个点值。第三步对数分桶近似考虑剩下的贡献即满足 $x,y\le\sqrt{n}$ 且 $\lfloor n / t\rfloor \sqrt{n}$ 的贡献。所有贡献必须满足 $xyt\le n$亦即 $\ln x \ln y \le \ln(n/t)$。取正数 $S$可以利用$$ \lceil S\ln x\rceil \lceil S\ln y\rceil \le S\ln(n/t) $$近似估计这一条件。定义多项式 $\sigma_f(u)$ 和 $\sigma_g(u)$使得系数 $[u^k]\sigma_f$ 与 $[u^k]\sigma_g$ 分别等于满足 $\lceil S\ln x\rceil k$ 时 $f(x)$ 的和、满足 $\lceil S\ln y\rceil k$ 时 $g(y)$ 的和只考虑 $x,y\le\sqrt{n}$ 的部分。利用快速傅里叶变换FFT得到乘积 $\sigma_f\sigma_g$其系数 $u^k$ 表示 $\lceil S\ln x\rceil \lceil S\ln y\rceil k$ 时 $f(x)g(y)$ 的和。由此只需对 $D(n)$ 中剩下每个 $\lfloor n/t\rfloor$ 找到满足 $k \le S\ln(n/t)$ 的最大 $k$ 值即可得到这一部分贡献的估计值。第四步误差修正近似条件可能遗漏部分贡献这只会发生在$$ S\ln x 1 S\ln y 1 \ge \lceil S\ln x\rceil \lceil S\ln y\rceil S\ln(n/t) \ge S\ln x S\ln y $$时等价于 $xyt \in (n\mathrm{e}^{-2/S},n]$这是一个长度为 $O(n/S)$ 的区间。枚举区间内所有可能的贡献 $(x,y,t)$逐个检验是否遗漏即可完成误差修正。为了快速枚举区间内所有贡献可以先筛出不超过 $\sqrt{n}$ 的全部素数用这些素数去除区间中的整数剩下的因子必然是大于 $\sqrt{n}$ 的素数由此得到区间内所有整数的素因数分解进而快速枚举所有可能的 $(x,y,t)$。复杂度分析估计贡献时需要对长度为 $S\log n$ 的多项式做乘法时间复杂度 $O(S\log n\log(S\log n))$误差修正时预处理素因数分解的时间复杂度为 $O(\sqrt{n}(n/S)\log\log n)$枚举区间内所有贡献的时间复杂度为 $O(\sum_{k\in(n\mathrm{e}^{-2/S},n]}d_3(k))$此处 $d_3(n)$ 表示将 $n$ 分解成三个有序整数乘积的方法数。解析数论的Piltz 除数问题结果指出$$ \sum_{k\le n}d_3(k) nP(\log n) O(n^{43/96\varepsilon}), $$其中 $P(\cdot)$ 是二次多项式。前两步时间复杂度已是 $O(\sqrt{n}\log n)$忽略所有 $o(\sqrt{n}\log n)$ 的项最后一部分贡献计算的时间复杂度为$$ O\left(S\log n\log(S\log n) \dfrac{n}{S} \log^2n\right). $$取 $S \sqrt{n}$就得到总时间复杂度$$ O(\sqrt{n}\log^2n). $$参考实现仓库 fast_bs_conv.cpp 提供了该算法的完整 C 实现。其核心结构BlockSieve用两个数组维护块筛struct BlockSieve { long long n, b; std::vectorint s1, s2; BlockSieve(long long _n) : n(_n), b(std::sqrt(_n 0.25l)), s1(b 1), s2(b 1) {} int operator[](long long x) { return x b ? s1[x] : s2[n / x]; } };operator[]把 $x\le b\lfloor\sqrt n\rfloor$ 的下标映射到稠密数组s1其余下标按 $\lfloor n/x\rfloor$ 映射到稀疏数组s2这正是 $D(n)$ 的两段式结构。主函数block_sieve_convolute依次实现前缀和差分求点值df/dg、$x$ 或 $y\sqrt n$ 部分的前缀和技巧、$xy\le b$ 的暴力枚举、按 $\lceil b\ln x\rceil$ 分桶后用 NTTntt_mul模 $998244353$原根 $3$做多项式乘法、以及利用素数筛分解区间 $(n\mathrm e^{-2/b},n]$ 内整数并 DFS 枚举 $(x,y,t)$ 的误差修正。文件末尾的test函数用随机数据对照双曲线法逐一验证 $D(n)$ 中每个点的取值输出Correct/Incorrect可以作为算法的正确性测试模板。杜教筛前文讨论了如何计算数论函数 Dirichlet 卷积的前缀和。本节考虑它的逆过程设 $f \ast g h$且 $f,h$ 已知计算 $g$ 的前缀和$$ G(n) \sum_{x1}^ng(x). $$换句话说本节考虑两个数论函数在 Dirichlet 卷积意义下的商的前缀和计算。本节总是假设 $f(1)\neq 0$以保证 $f$ 可逆Dirichlet 卷积逆元存在的充要条件见 dirichlet.md。为此在 Dirichlet 双曲线法表达式中令 $x_0 n$就得到$$ H(n) \sum_{x1}^{n}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right). $$直接解出 $G(n)$就得到杜教筛的表达式$$ G(n) \dfrac{1}{f(1)}\left(H(n)-\sum_{x2}^nf(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right)\right). $$实际上对于 $x_0 \ge 1$总有更一般的带分点形式$$ G(n) \dfrac{1}{f(1)}\left(H(n)-\sum_{x2}^{\lfloor x_0\rfloor}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right) - \sum_{y1}^{\lfloor y_0\rfloor}F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)g(y) F(\lfloor x_0\rfloor)G(\lfloor y_0\rfloor)\right). $$无论是哪种形式它都是一个关于 $G(n)$ 的递推关系式。为计算 $G(n)$ 的取值需要计算 $G$ 在 $D(n)\setminus{n}$ 处的取值。因为 $D(n)$ 具有递归结构——即对于 $m\in D(n)$总是有 $D(m)\subseteq D(n)$——所以在整个递归计算过程中只需要计算 $G$ 在 $D(n)$ 中元素处的取值各一次。换句话说计算 $G(n)$ 时实际上得到了 $g$ 的块筛 $\mathcal S_g(n)$。实现方式具体实现时可以采用递归方法并配合记忆化避免重复计算也可以采用迭代方法从小到大依次计算 $D(n)$ 中每个点处 $G$ 的取值。此时表达式中的求和式既可以用数论分块计算也可以用 Dirichlet 双曲线法计算。这些实现的复杂度是相同的。复杂度由于杜教筛总是得到块筛所以杜教筛的复杂度其实相当于计算块筛的复杂度。如果已知信息只有 $F,H$ 的块筛那么杜教筛的复杂度就是 $O(n^{3/4})$如果对某个 $z\ge\sqrt{n}$可以在 $T_0(z)$ 时间内预处理出 $g$ 在 $1\le x \le z$ 处的点值那么杜教筛的复杂度就是$$ O\left(T_0(z) \dfrac{n}{\sqrt{z}}\right). $$当 $g$ 是积性函数时可以应用线性筛即 $T_0(z)\Theta(z)$所以最优需要预处理到 $z n^{2/3}$ 处总时间复杂度为 $O(n^{2/3})$对于更一般的情形总时间复杂度则为 $O(n^{2/3}(\log n)^{1/3})$。这些都和块筛部分的分析完全一致。警告递归实现时不使用记忆化将导致复杂度错误。杜教筛表达式中计算 $G(n)$ 依赖 $D(n)\setminus{n}$ 中 $G$ 的取值保证复杂度的关键在于 $D(n)$ 的递归结构当 $m\in D(n)$ 时$D(m)\subseteq D(n)$。利用记忆化后算法复杂度为 $O(n^{3/4})$而不使用记忆化时设计算 $G(n)$ 的复杂度为 $T(n)$有$$ \begin{aligned} T(n) \Theta(\sqrt{n}) \sum_{d\in D(n),~d\neq n}T(d)\ \Theta(\sqrt{n}) \sum_{x1}^{\lfloor n/\lfloor\sqrt{n}\rfloor\rfloor - 1} T(x) \sum_{x 2}^{\lfloor\sqrt{n}\rfloor}T\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right). \end{aligned} $$利用类似主定理的证明思路可以说明最后一项主导该式的增长且 $T(n)\in\Theta(n^\alpha)$其中 $\alpha\approx 1.73$ 是 $\zeta(\alpha)2$ 的根——远远劣于记忆化版本的 $O(n^{3/4})$。应用要点应用杜教筛计算数论函数 $g$ 前缀和时关键在于找到合适的 $f,h$ 使得 $hf\ast g$ 且 $f,h$ 的块筛都容易计算。有些时候这样的 $f,h$ 是显然的如 $\varepsilon\mu\ast 1$、$\mathrm{id}\varphi\ast 1$另一些时候需要利用 Dirichlet 卷积的性质或通过计算相应的 Dirichlet 生成函数来找到相应分解。后文的例题展示了这些情形。例题例 1ARC116 C - Multiple Sequences给定正整数 $N$ 和 $M$计算有多少长度为 $N$ 的序列 $A$ 满足 $1 \le A_i \le M$ 且 $A_i$ 整除 $A_{i1}$。答案对 $998244353$ 取模。数据范围$1 \le N, M \le 2\times 10^5$。解答设长度为 $n$ 且 $A_nm$ 的序列数目为 $f_n(m)$那么最终答案就是 $\sum_{m1}^M f_N(m)$。考虑动态规划转移方程为$$ f_n(m) \sum_{k\mid m}f_{n-1}(k). $$利用 Dirichlet 卷积记号可记作 $f_n f_{n-1}\ast 1$其中 $1$ 是常值数论函数。注意到 $f_1 1$归纳可知 $f_n 1^{\ast n}$即 $f_n$ 是 $n$ 个常值函数的卷积。最终答案就是 $f_N$ 的前缀和。由于过程中只涉及积性函数利用前文介绍的 Dirichlet 卷积前缀和计算方法单次卷积前缀和的计算只需 $O(M^{2/3})$ 时间再利用快速幂的思想只需计算 $O(\log N)$ 次卷积前缀和即可得到 $f_N$ 的卷积前缀和的值。整体时间复杂度为 $O(M^{2/3}\log N)$。参考实现仓库 hyperbola.cpp 是该题的完整解法其结构恰好体现了块筛 点值信息 快速幂的组合dirichlet_convolute在假设 $h$ 是积性函数的前提下用线性筛框架计算卷积 $hf\ast g$ 的前缀点值——素数次幂处直接按定义求和rem[x] 1分支合数处利用积性拆成互素因子相乘h[x] h[rem[x]] * h[x / rem[x]]BlockSieve结构体维护 $n,z$ 以及三组数据 $f,F,F_2$其中 $F$ 保存 $1\le x\le z$ 处前缀和$F_2$ 保存 $xz$ 处按 $\lfloor n/x\rfloor$ 索引的前缀和sum(x)在两者间切换operator*实现两个块筛的卷积——对 $1\le i\le z$ 直接卷积累加对 $iz$ 用双曲线法公式枚举 $x\le\sqrt{k}$ 的 $f(x)G(k/x)F(k/x)g(x)$ 再减去重叠修正 $F(\sqrt{k})G(\sqrt{k})$main取 $zn^{2/3}$初始化常值函数 $1$ 的块筛po对卷积做二进制快速幂res res * po; po po * po;最终输出res.sum(n)。仓库还提供了该题的测试数据hyperbola.in内容为200000 200000与对应的标准答案 hyperbola.ans835917264可用于验证实现。例 2P4213【模板】杜教筛Sum设 $\mu$ 和 $\varphi$ 分别是莫比乌斯函数和欧拉函数。求 $S_1(n) \sum_{i1}^{n} \mu(i)$ 和 $S_2(n) \sum_{i1}^{n} \varphi(i)$ 的值。数据范围$1\leq n2^{31}$。解答注意到两个 Dirichlet 卷积关系式$$ \varepsilon \mu \ast 1,~ \operatorname{id} \varphi \ast 1. $$其中 $\varepsilon(n) [n1]$ 是 Dirichlet 卷积的单位元函数$\operatorname{id}(n) n$ 是恒等函数$1(n)1$ 是常值函数。这三个函数的前缀和都可以在 $O(1)$ 时间内计算且都是积性函数所以利用前文介绍的杜教筛方法可以在 $O(n^{2/3})$ 时间内计算。对于欧拉函数前缀和另一种方法是利用莫比乌斯反演$$ \begin{aligned} S_2(n) \sum_{i1}^n\varphi(i) \sum_{i1}^n\sum_{j1}^i[i\perp j] \ \sum_{i1}^n\sum_{j1}^i\sum_d\mu(d)[d\mid i][d\mid j] \ \sum_d\mu(d)\dfrac{1}{2}\left\lfloor\dfrac{n}{d}\right\rfloor\left(\left\lfloor\dfrac{n}{d}\right\rfloor1\right). \end{aligned} $$在数论分块过程中需要 $\mu(d)$ 的前缀和而这可以通过杜教筛预处理出来时间复杂度仍然是 $O(n^{2/3})$。参考实现仓库 du_1.cpp 给出模板实现其中杜教筛的记忆化递归写得非常清晰long long S_mu(long long x) { // 求mu的前缀和 if (x MAXN) return sum_mu[x]; if (mp_mu[x]) return mp_mu[x]; // 如果map中已有该大小的mu值则可直接返回 long long ret 1; for (long long i 2, j; i x; i j 1) { j x / (x / i); ret - S_mu(x / i) * (j - i 1); } return mp_mu[x] ret; // 路径压缩方便下次计算 }实现要点对 $x MAXN$ 直接查线性筛预处理的前缀和sum_mu线性筛部分同时求出了 $\mu$ 数组对更大的 $x$利用 $\varepsilon\mu\ast 1$ 即 $S_\mu(n)1-\sum_{i2}^{n}S_\mu(\lfloor n/i\rfloor)$配合数论分块j x / (x / i)递归求解并用map记忆化S_phi则基于 $\mathrm{id}\varphi\ast 1$ 或上述莫比乌斯反演式在数论分块中复用S_mu。main里 $T$ 组询问每组输出S_phi(n) S_mu(n)。例 3P3768 简单的数学题给定 $p,n$计算 $\sum_{i1}^n\sum_{j1}^nij\cdot\gcd(i,j)\pmod p$。数据范围$n\leq 10^{10}$$5\times 10^8\leq p\leq 1.1\times 10^9$ 且 $p$ 是质数。解答利用欧拉函数的性质做反演$$ \begin{aligned} T(n) \sum_{i1}^n\sum_{j1}^nij\cdot\gcd(i,j)\ \sum_{i1}^n\sum_{j1}^nij\sum_d\varphi(d)[d\mid i][d\mid j]\ \sum_d\varphi(d)\left(\sum_{i1}^{\lfloor n/d\rfloor}id\right)\left(\sum_{j1}^{\lfloor n/d\rfloor}jd\right)\ \sum_d d^2\varphi(d) F\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)^2. \end{aligned} $$其中 $F(n) \dfrac{1}{2}n(n1)$。该式可以通过数论分块计算但需要预处理 $d^2\varphi(d)$ 的前缀和。为此利用杜教筛。记 $f(n)(\operatorname{id}^2\varphi)(n)$ 和 $S(n)\sum_{i1}^n f(i)$。应用杜教筛的关键是构造函数 $g$使得 $f\ast g$ 与 $g$ 都可以快速求和。前文已讨论用杜教筛预处理 $\varphi$ 前缀和的方法它利用关系 $\mathrm{id} \varphi\ast 1$。相较于 $\varphi$这里的 $f$ 多了一个 $\operatorname{id}^2$由于 $\operatorname{id}$ 是完全积性函数利用 Dirichlet 卷积的性质只需将每一项都乘以 $\operatorname{id}^2$就得到$$ \operatorname{id}^3 f \ast \operatorname{id}^2. $$因为 $\operatorname{id}^2(n)n^2$ 与 $\operatorname{id}^3(n)n^3$ 的前缀和都可以在 $O(1)$ 时间内计算$f$ 的前缀和就可以在 $O(n^{2/3})$ 时间内预处理得到再加上数论分块整体时间复杂度仍是 $O(n^{2/3})$。另一种将类似积性函数表示为两函数之在 Dirichlet 卷积意义下的商的方法是利用 Dirichlet 生成函数。参考实现仓库 du_2.cpp 给出模板实现。它的线性筛到 $n^{2/3}$pn (long long)pow(n, 0.666667); prime_work(pn);正是 $zn^{2/3}$ 最优分点的代码体现prime_work用线性筛求出 $\varphi$ 并预计算 $s[i]i^2\varphi(i)$ 的前缀和模 $P$s2/s3给出平方和与立方和用费马小定理ksm求 $6$、$2$ 的逆元calc(k)利用 $\mathrm{id}^3f\ast\mathrm{id}^2$ 做记忆化递归s_map存储超过预处理上限的答案求和式用数论分块最终solve对外层和式 $\sum_d d^2\varphi(d)F(\lfloor n/d\rfloor)^2$ 做数论分块并调用calc。文件首行注释不要为了省什么内存把数组开小,会卡80也提醒了 $N5e6$ 这类数组规模的选型理由。习题AtCoder Xmas Contest 2019 D - Sum of (-1)^f(n)需要综合利用本文学到的块筛与卷积前缀和技术可以结合 dirichlet.md 与 mobius.md 阅读。参考资料与注释任之洲. 2016. 《积性函数求和的几种方法》. 2016 年信息学奥林匹克中国国家队候选队员论文.周康阳. 2024. 《关于积性函数求和问题的一些进展》. 2024 年信息学奥林匹克中国国家队候选队员论文. 快速块筛卷积算法出处数论分块的完整性质与 $D(n)$ 的递归结构见 sqrt-decomposition.mdDirichlet 卷积的代数性质、逆元存在条件与 Dirichlet 生成函数见 dirichlet.md快速块筛卷积的 NTT 部分可参考仓库 fast_bs_conv.cpp 中ntt_mul的 DIF/DIT 实现关于 Piltz 除数问题$d_3(n)$ 的均值公式可查阅除数求和函数相关的解析数论文献。【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考