跳转至

狄利克雷双曲线法 & 杜教筛

前置知识：Dirichlet 卷积、数论分块

本文介绍 Dirichlet 双曲线法，它可用于在亚线性时间复杂度内计算两个数论函数的 Dirichlet 卷积的前缀和．通过对该算法的分析，可以引申出块筛的概念；据此，本文将进一步介绍块筛卷积的快速计算方法．最后，本文将介绍杜教筛，它可以用于计算两个数论函数在 Dirichlet 卷积意义下的商的前缀和．

Dirichlet 双曲线法

Dirichlet 双曲线法可以用于计算两个数论函数的 Dirichlet 卷积的前缀和．

设 \(f,g,h\) 是数论函数，且 \(h = f\ast g\)．那么，利用 Dirichlet 卷积的定义，\(h\) 的前缀和

\[ H(n) = \sum_{k=1}^nh(k) = \sum_{k=1}^n\sum_{xy=k}f(x)g(y). \]

求和式遍历的点集恰为第一象限（不含坐标轴）中双曲线 \(xy=n\) 下方的整点集合．设整点 \((x,y)\) 的权值为 \(f(x)g(y)\)，那么 \(H(n)\) 就是这一权值的和．

如图所示，这一权值和可以通过容斥原理计算：

\[ H(n) = \sum_{x=1}^{\lfloor x_0\rfloor}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right) + \sum_{y=1}^{\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)\) 是双曲线 \(xy=n\) 上任意一个点．表达式中，第一项表示图中绿色区域的权值和，第二项表示图中橙色区域的权值和，第三项则表示两个区域重叠部分的权值和．这个表达式仅含有 \(\lfloor x_0\rfloor + \lfloor y_0\rfloor + 1\) 项．对于合理选择的 \((x_0,y_0)\)，它的计算复杂度显著优于暴力计算 \(h(n)\) 的前缀和．这就是 Dirichlet 双曲线法（Dirichlet hyperbola method）．

卷积前缀和点值的计算

Dirichlet 双曲线法最基本的应用就是计算前缀和函数的点值 \(H(n)\)．

如果 \(F,G\) 的点值已知（或可以在 \(O(1)\) 时间内计算），进而 \(f,g\) 的点值也已知，那么 Dirichlet 双曲线法的表达式中每一项都可以在 \(O(1)\) 时间内计算，总时间复杂度就等于 \(O(x_0+y_0)\)．因为 \(x_0y_0=n\)，所以由均值不等式可知，当 \(x_0=y_0=\sqrt{n}\) 时，就得到最低时间复杂度 \(O(\sqrt{n})\)．

这并非新的结果．在 Dirichlet 双曲线法的表达式中，令 \(x_0 > n\)，就得到

\[ H(n) = \sum_{x=1}^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}{y+1}\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\)，亦即

\[ \left\lfloor\dfrac{n}{y+1}\right\rfloor + 1 = x = \left\lfloor\dfrac{n}{y}\right\rfloor. \]

因此，这些分块对应的求和为

\[ I_1 = \sum_{x=1}^{\lfloor\sqrt{n}\rfloor}\left(F(x)-F(x-1)\right)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right) = \sum_{x=1}^{\lfloor\sqrt{n}\rfloor}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right). \]

剩下的分块，高度 \(y\) 遍历 \(1,2,\cdots,y^*= \left\lfloor\dfrac{n}{\lfloor\sqrt{n}\rfloor + 1}\right\rfloor\)．对剩余的和式应用 Abel 变换（即分部积分法的求和形式），就得到

\[ \begin{aligned} I_2 &= \sum_{y=1}^{y^*}\left(F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)-F\left(\left\lfloor\dfrac{n}{y+1}\right\rfloor\right)\right)G(y) \\ &= \sum_{y=1}^{y^*}F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)g(y) - F\left(\left\lfloor\dfrac{n}{y^*+1}\right\rfloor\right)G(y^*). \end{aligned} \]

数论分块的性质说明，\(\lfloor n/\lfloor \sqrt{n}\rfloor\rfloor \ge \lfloor \sqrt{n}\rfloor\)．由此，分两种情形：

  • 当 \(\lfloor n/\lfloor \sqrt{n}\rfloor\rfloor = \lfloor \sqrt{n}\rfloor\) 时，必然有 \(y^* < \lfloor \sqrt{n}\rfloor\)．由于它与 \(x = \lfloor \sqrt{n}\rfloor\) 所在块相邻，块高只能恰好相差一，即 \(y^* = \lfloor \sqrt{n}\rfloor - 1\)．此时，有

    \[ \begin{aligned} I_2 &= \sum_{y=1}^{\lfloor\sqrt{n}\rfloor-1}F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)g(y) - F(\lfloor\sqrt{n}\rfloor)G(\lfloor\sqrt{n}\rfloor-1)\\ &= \sum_{y=1}^{\lfloor\sqrt{n}\rfloor-1}F\left(\left\lfloor\dfrac{n}{y}\right\rfloor\right)g(y) - F(\lfloor\sqrt{n}\rfloor)\left(G(\lfloor\sqrt{n}\rfloor) - g(\lfloor\sqrt{n}\rfloor)\right) \\ &= \sum_{y=1}^{\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). \end{aligned} \]

  • 当 \(\lfloor n/\lfloor \sqrt{n}\rfloor\rfloor > \lfloor \sqrt{n}\rfloor\) 时，必然有 \(y^* = \lfloor \sqrt{n}\rfloor\)．直接代入求和式，就得到

    \[ I_2 = \sum_{y=1}^{\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). \]

综上，除了一个 Abel 变换，数论分块的计算过程实际上就是在计算如下算式：

\[ H(n) = \sum_{x=1}^{\lfloor\sqrt{n}\rfloor}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right) + \sum_{y=1}^{\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})\) 时 Dirichlet 双曲线法的表达式．因此可以说，两种算法的计算过程几乎等价．而且，由于 Dirichlet 双曲线法利用了更多数论分块的性质，避免了朴素数论分块中的不必要计算，常数更小一些．

在处理实际问题时，已知 \(F,G\) 点值这一条件可能过强．但是，Dirichlet 双曲线法（或对应的数论分块）其实并不需要 \(F,G\) 的全部点值信息．观察求和表达式可知，其实只需要 \(F\) 和 \(G\) 在数论分块关键点集合

\[ D(n) = \left\{\left\lfloor\dfrac{n}{x}\right\rfloor : 1 \le x \le n,~x\in\mathbf N_+\right\} \]

处的取值．它既是全体块高的集合，也是全体块的右端点集合．注意，根据数论分块的 性质 可知，这一集合包含所有 \(1\le x\le\sqrt{n}\) 的整数值，因此，已知 \(F,G\) 在 \(D(n)\) 处的取值就相当于已知 \(f,g\) 在所有 \(1\le x\le\sqrt{n}\) 处的取值．这一集合的大小是 \(\Theta(\sqrt{n})\) 的，所以，计算 \(H(n)\) 时，只需要知道 \(F,G\) 在一个稀疏集合处的点值信息．这个观察是优化数论函数前缀和计算的关键．

块筛及其卷积

有些时候，\(h=f\ast g\) 并非最终需要计算前缀和的函数，而只是中间步骤之一．根据前文分析，为了后续计算，需要求出前缀和函数 \(H\) 在集合 \(D(n)\) 处的取值．这就称为数论函数 \(h\) 的 块筛：

\[ \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_{x=1}^{\lfloor\sqrt{n}\rfloor}\sqrt{x} + \sum_{x=1}^{\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} \]

正是因为块筛是稀疏的，所以整个块筛可以在亚线性时间内求出．

但是，这一算法显然过于暴力．集合 \(D(n)\) 中较小的那些元素取值相对稠密，所以块筛中相邻两个前缀和相差并不大，完全可以直接计算卷积 \(h\) 的点值，再求它的前缀和．这比起对每个点都单独求一遍前缀和点值要更快．例如，对 \(x = 1,2,\cdots,\lfloor\sqrt{n}\rfloor\) 分别计算前缀和点值需要

\[ O\left(\sum_{x=1}^{\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_{x=1}^{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\) 是积性的，那么当 \(z=n^{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_{xy=z} f(x)g(y) \]

在块筛 \(D(n) = \{\left\lfloor n / t\right\rfloor : 1\le t \le n\}\) 处的前缀和．对于这一问题，单一贡献可以由 \((x,y,t)\) 标记，即将项 \(f(x)g(y)\) 累加到 \(\left\lfloor n / t\right\rfloor\) 处前缀和的过程．算法将这些贡献分成若干组处理．

首先，考虑 \(x > \sqrt{n}\) 的所有点对于这一系列前缀和的贡献，那些 \(y > \sqrt{n}\) 的贡献类似．因为所有贡献必须满足 \(xy \le \lfloor n / t\rfloor\)，即 \(xyt \le n\)，所以只需要枚举所有可能的 \(t,y\)，利用前缀和技巧以及块筛 \(\mathcal S_f\) 中的信息，就可以在 \(O(\sum_{t,y: ty\le\sqrt{n}} 1) = O(\sqrt{n}\log n)\) 时间内计算出这部分贡献．

然后，考虑 \(\left\lfloor n / t\right\rfloor\le\sqrt{n}\) 的这部分贡献．这一部分同样可以暴力枚举所有可能的 \(x,y\) 完成，时间复杂度仍然是 \(O(\sum_{x,y:xy\le\sqrt{n}}1) = O(\sqrt{n}\log n)\)．这一部分实际上得到了函数 \(h\) 在 \(D(n)\) 的前 \(\lfloor\sqrt{n}\rfloor\) 个点值．

考虑剩下的贡献，即满足 \(x,y\le\sqrt{n}\) 且 \(\left\lfloor n / t\right\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](\sigma_f\sigma_g)\) 就表示 \(\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\) 值，就可以得到这一部分贡献的估计值．

最后，对于这一部分贡献进行误差修正．前文试图利用更强的条件 \(\lceil S\ln x\rceil + \lceil S\ln y\rceil \le S\ln(n/t)\) 去近似 \(\ln x + \ln y \le \ln(n/t)\)，可能会遗漏部分贡献．这只会发生在 \(\lceil S\ln x\rceil + \lceil S\ln y\rceil > S\ln(n/t)\) 时，所以有

\[ 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\) 分解成三个有序整数乘积的方法数．解析数论的结果1指出：

\[ \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)\)．这也是该算法的总时间复杂度．

参考实现 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
struct BlockSieve {
  long long n, b;
  std::vector<int> 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]; }
};

// Fast block sieve convolution by ZKY.
BlockSieve block_sieve_convolute(const BlockSieve& f, const BlockSieve& g) {
  // Special cases for n <= 3 are not implemented.
  // The upper bound is a loose bound for NTT to work for M = 998244353.
  assert(f.n == g.n && f.n > 3 && f.n < 1e11);
  long long n = f.n, b = f.b;
  BlockSieve h(n);
  std::vector<int> df(b + 1), dg(b + 1);
  for (int i = 1; i <= b; ++i) {
    df[i] = (f.s1[i] + M - f.s1[i - 1]) % M;
    dg[i] = (g.s1[i] + M - g.s1[i - 1]) % M;
  }
  // Prefix sum trick for x or y > sqrt(n).
  for (int t = 1; t <= n / (b + 1); ++t) {
    long long x;
    for (int y = 1; (x = n / (t * y)) > b; ++y) {
      (h.s2[t] += (long long)dg[y] * (f.s2[n / x] + M - f.s1[b]) % M) %= M;
      (h.s2[t] += (long long)df[y] * (g.s2[n / x] + M - g.s1[b]) % M) %= M;
    }
  }
  // Brutal force for small x*y.
  for (int x = 1; x <= b; ++x) {
    for (int y = 1; x * y <= b; ++y) {
      (h.s1[x * y] += (long long)df[x] * dg[y] % M) %= M;
    }
  }
  for (int i = 1; i <= b; ++i) {
    (h.s1[i] += h.s1[i - 1]) %= M;
  }
  // Approximation.
  std::vector<int> id1(b + 1);
  for (int i = 1; i <= b; ++i) {
    id1[i] = std::ceil(b * std::log(i));
  }
  std::vector<long double> id2(b + 1);
  for (int i = 1; i <= b; ++i) {
    id2[i] = b * std::log(1.0l * n / i);
  }
  std::vector<int> sf(id1[b] + 1), sg(id1[b] + 1);
  for (int i = 1; i <= b; ++i) {
    (sf[id1[i]] += df[i]) %= M;
    (sg[id1[i]] += dg[i]) %= M;
  }
  auto sh = ntt_mul(sf, sg);
  for (int i = 1; i < (int)sh.size(); ++i) {
    (sh[i] += sh[i - 1]) %= M;
  }
  for (int t = 1; t <= n / (b + 1); ++t) {
    int k = std::min((int)std::floor(id2[t]), (int)sh.size() - 1);
    (h.s2[t] += sh[k]) %= M;
  }
  // Correction.
  auto primes = get_primes(b);  // Obtain primes not exceeding b.
  long long l = n * std::exp(-2.0l / b);
  std::vector<long long> co(n - l + 1);
  std::iota(co.begin(), co.end(), l);
  std::vector<std::vector<std::pair<long long, int>>> fac(n - l + 1);
  for (int p : primes) {
    for (int i = (l + p - 1) / p * p - l; i <= n - l; i += p) {
      int e = 0;
      for (; co[i] % p == 0; co[i] /= p) ++e;
      fac[i].emplace_back(p, e);
    }
  }
  for (int i = 0; i <= n - l; ++i) {
    if (co[i] != 1) {
      fac[i].emplace_back(co[i], 1);
    }
  }
  long long num = 0;
  std::vector<std::pair<long long, int>> cur;
  auto dfs = [&](auto&& dfs, int k, int x, int y) {
    if (k == (int)cur.size()) {
      auto t = num / x / y;
      if (t && t <= n / (b + 1) && id1[x] + id1[y] > id2[t]) {
        (h.s2[t] += (long long)df[x] * dg[y] % M) %= M;
      }
      return;
    }
    auto p = cur[k].first;
    auto e = cur[k].second;
    long long nx = x;
    for (int e1 = 0; e1 <= e; ++e1) {
      if (e1) {
        nx *= p;
        if (nx > b) break;
      }
      long long ny = y;
      for (int e2 = 0; e1 + e2 <= e; ++e2) {
        if (e2) {
          ny *= p;
          if (ny > b) break;
        }
        dfs(dfs, k + 1, nx, ny);
      }
    }
  };
  for (int i = 0; i <= n - l; ++i) {
    num = l + i;
    cur = fac[i];
    dfs(dfs, 0, 1, 1);
  }
  return h;
}

杜教筛

前文讨论了如何计算数论函数 Dirichlet 卷积的前缀和．本节考虑它的逆过程：设 \(f \ast g = h\)，且 \(f,h\) 已知，计算 \(g\) 的前缀和

\[ G(n) = \sum_{x=1}^ng(x). \]

换句话说，本节考虑两个数论函数 Dirichlet 卷积意义的商的前缀和计算．本节总是假设 \(f(1)\neq 0\)，以保证 \(f\) 可逆．

为此，在 Dirichlet 双曲线法表达式中令 \(x_0 > n\)，就得到

\[ H(n) = \sum_{x=1}^{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_{x=2}^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_{x=2}^{\lfloor x_0\rfloor}f(x)G\left(\left\lfloor\dfrac{n}{x}\right\rfloor\right) - \sum_{y=1}^{\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_{x=1}^{\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\) 的根．

应用杜教筛计算数论函数 \(g\) 前缀和时，关键在于找到合适的 \(f,h\) 使得 \(h=f\ast g\) 且 \(f,h\) 的块筛都容易计算．在有些时候，这样的 \(f,h\) 是显然的；在另一些时候，需要利用 Dirichlet 卷积的性质或通过计算相应的 Dirichlet 生成函数来找到相应的分解．后文的例题展示了这些情形．

例题

本节介绍几道计算数论函数前缀和的题目．

AtCoder Regular Contest 116 C - Multiple Sequences

给定正整数 \(N\) 和 \(M\)．需要计算有多少长度为 \(N\) 的序列 \(A\) 满足 \(1 \le A_i \le M\) 且 \(A_i\) 整除 \(A_{i+1}\)．答案对 \(998244353\) 取模．数据范围：\(1 \le N, M \le 2\times 10^5\)．

解答

设长度为 \(n\) 且 \(A_n=m\) 的序列数目为 \(f_n(m)\)．那么，最终的答案就是 \(\sum_{m=1}^M f_N(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)\)．

参考实现 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#include <cassert>
#include <cmath>
#include <iostream>
#include <vector>

constexpr int M = 998244353;

// Compute the Dirichlet convolution h = f * g.
// Assume that h is multiplicative.
auto dirichlet_convolute(const std::vector<int>& f, const std::vector<int>& g) {
  int n = f.size() - 1;
  std::vector<int> h(n + 1), primes, rem(n + 1), lpf(n + 1);
  std::vector<bool> vis(n + 1);
  h[1] = 1;
  for (int x = 2; x <= n; ++x) {
    if (!vis[x]) {
      primes.push_back(x);
      rem[x] = 1;
      lpf[x] = x;
    }
    for (int p : primes) {
      if (x * p > n) break;
      vis[x * p] = true;
      rem[x * p] = x % p ? x : rem[x];
      lpf[x * p] = p;
      if (x % p == 0) break;
    }
    if (rem[x] == 1) {  // prime powers.
      for (int k = x; k; k /= lpf[x]) {
        (h[x] += (long long)f[k] * g[x / k] % M) %= M;
      }
    } else {  // other cases.
      h[x] = (long long)h[rem[x]] * h[x / rem[x]] % M;
    }
  }
  return h;
}

// Block sieve convolution for multiplicative functions.
struct BlockSieve {
  int n, z;
  std::vector<int> f, F, F2;

  BlockSieve(int _n, int _z)
      : n(_n), z(_z), f(_z + 1), F(_z + 1), F2(_n / _z + 1) {}

  int sum(int x) const { return x <= z ? F[x] : F2[n / x]; }

  BlockSieve operator*(const BlockSieve& oth) const {
    assert(n == oth.n && z == oth.z);
    BlockSieve res(n, z);
    res.f = dirichlet_convolute(f, oth.f);
    for (int i = 1; i <= z; ++i) {
      res.F[i] = (res.F[i - 1] + res.f[i]) % M;
    }
    for (int i = 1; i <= n / z; ++i) {
      int k = n / i;
      int sqr = std::sqrt(k + 0.25l);
      for (int x = 1; x <= sqr; ++x) {
        (res.F2[i] += (long long)f[x] * oth.sum(k / x) % M) %= M;
        (res.F2[i] += (long long)sum(k / x) * oth.f[x] % M) %= M;
      }
      (res.F2[i] += (M - (long long)sum(sqr) * oth.sum(sqr) % M)) %= M;
    }
    return res;
  }
};

int main() {
  int n, k;
  std::cin >> k >> n;
  int z = std::pow(n, 2.0l / 3);
  // Initialization of inputs.
  BlockSieve po(n, z);
  for (int i = 1; i <= z; ++i) {
    po.f[i] = 1;
    po.F[i] = i;
  }
  for (int i = 1; i <= n / z; ++i) {
    po.F2[i] = n / i;
  }
  BlockSieve res(n, z);
  res.f[1] = 1;
  for (int i = 1; i <= z; ++i) {
    res.F[i] = 1;
  }
  for (int i = 1; i <= n / z; ++i) {
    res.F2[i] = 1;
  }
  // Binary exponentiation.
  for (; k; k >>= 1) {
    if (k & 1) res = res * po;
    po = po * po;
  }
  std::cout << res.sum(n) << std::endl;
  return 0;
}
P4213【模板】杜教筛（Sum）

设 \(\mu\) 和 \(\varphi\) 分别是莫比乌斯函数和欧拉函数．求 \(S_1(n)= \sum_{i=1}^{n} \mu(i)\) 和 \(S_2(n)= \sum_{i=1}^{n} \varphi(i)\) 的值．数据范围：\(1\leq n<2^{31}\).

解答

注意到 Dirichlet 卷积关系式：

\[ \varepsilon = \mu \ast 1,~ \operatorname{id} = \varphi \ast 1. \]

其中，\(\varepsilon(n) = [n=1]\) 是 Dirichlet 卷积的单位元函数，\(\operatorname{id}(n) = n\) 是恒等函数，\(1(n) = 1\) 是常值函数．因为这三个函数的前缀和都可以在 \(O(1)\) 时间内计算且都是积性函数，所以，利用前文介绍的杜教筛方法就可以在 \(O(n^{2/3})\) 时间内计算．

对于欧拉函数的前缀和的计算，另一种方法是利用 莫比乌斯反演．注意到

\[ \begin{aligned} S_2(n) &= \sum_{i=1}^n\varphi(i) = \sum_{i=1}^n\sum_{j=1}^i[i\perp j] \\ &= \sum_{i=1}^n\sum_{j=1}^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\rfloor+1\right). \end{aligned} \]

在数论分块的过程中，需要 \(\mu(d)\) 的前缀和，而这可以通过杜教筛预处理出来．时间复杂度仍然是 \(O(n^{2/3})\)．

参考实现 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include <cstring>
#include <iostream>
#include <map>
using namespace std;
constexpr int MAXN = 2000010;
long long T, n, pri[MAXN], cur, mu[MAXN], sum_mu[MAXN];
bool vis[MAXN];
map<long long, long long> mp_mu;

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 = (long long)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;  // 路径压缩，方便下次计算
}

long long S_phi(long long x) {  // 求phi的前缀和
  long long ret = (long long)0;
  long long j;
  for (long long i = 1; i <= x; i = j + 1) {
    j = x / (x / i);
    ret += (S_mu(j) - S_mu(i - 1)) * (x / i) * (x / i);
  }
  return (ret - 1) / 2 + 1;
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  cin >> T;
  mu[1] = 1;
  for (int i = 2; i < MAXN; i++) {  // 线性筛预处理mu数组
    if (!vis[i]) {
      pri[++cur] = i;
      mu[i] = -1;
    }
    for (int j = 1; j <= cur && i * pri[j] < MAXN; j++) {
      vis[i * pri[j]] = true;
      if (i % pri[j])
        mu[i * pri[j]] = -mu[i];
      else {
        mu[i * pri[j]] = 0;
        break;
      }
    }
  }
  for (int i = 1; i < MAXN; i++)
    sum_mu[i] = sum_mu[i - 1] + mu[i];  // 求mu数组前缀和
  while (T--) {
    cin >> n;
    cout << S_phi(n) << ' ' << S_mu(n) << '\n';
  }
  return 0;
}
「LuoguP3768」简单的数学题

给定 \(p,n\)，计算

\[ \sum_{i=1}^n\sum_{j=1}^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_{i=1}^n\sum_{j=1}^nij\cdot\gcd(i,j)\\ &= \sum_{i=1}^n\sum_{j=1}^nij\sum_d\varphi(d)[d\mid i][d\mid j]\\ &= \sum_d\varphi(d)\left(\sum_{i=1}^{\lfloor n/d\rfloor}id\right)\left(\sum_{j=1}^{\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(n+1)\)．该式可以通过数论分块计算，但需要预处理出 \(d^2\varphi(d)\) 的前缀和．

为此，可以利用杜教筛．记 \(f(n)=(\operatorname{id}^2\varphi)(n)\) 和 \(S(n)=\sum_{i=1}^n f(i)\)．应用杜教筛的关键是构造出函数 \(g\)，使得 \(f\ast g\) 和 \(g\) 都可以快速求和．前文已经讨论了用杜教筛预处理 \(\varphi\) 前缀和的方法，这只需要利用关系：\(\operatorname{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 生成函数．

参考实现 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
// 不要为了省什么内存把数组开小,会卡80
#include <cmath>
#include <iostream>
#include <map>
using namespace std;
constexpr int N = 5e6, NP = 5e6, SZ = N;
long long n, P, inv2, inv6, s[N];
int phi[N], p[NP], cnt, pn;
bool bp[N];
map<long long, long long> s_map;

long long ksm(long long a, long long m) {  // 求逆元用
  long long res = 1;
  while (m) {
    if (m & 1) res = res * a % P;
    a = a * a % P, m >>= 1;
  }
  return res;
}

void prime_work(int k) {  // 线性筛phi，s
  bp[0] = bp[1] = true, phi[1] = 1;
  for (int i = 2; i <= k; i++) {
    if (!bp[i]) p[++cnt] = i, phi[i] = i - 1;
    for (int j = 1; j <= cnt && i * p[j] <= k; j++) {
      bp[i * p[j]] = true;
      if (i % p[j] == 0) {
        phi[i * p[j]] = phi[i] * p[j];
        break;
      } else
        phi[i * p[j]] = phi[i] * phi[p[j]];
    }
  }
  for (int i = 1; i <= k; i++)
    s[i] = (1ll * i * i % P * phi[i] % P + s[i - 1]) % P;
}

long long s3(long long k) {  // 立方和
  return k %= P, (k * (k + 1) / 2) % P * ((k * (k + 1) / 2) % P) % P;
}

long long s2(long long k) {  // 平方和
  return k %= P, k * (k + 1) % P * (k * 2 + 1) % P * inv6 % P;
}

long long calc(long long k) {  // 计算S(k)
  if (k <= pn) return s[k];
  if (s_map[k]) return s_map[k];  // 对于超过pn的用map离散存储
  long long res = s3(k), pre = 1, cur;
  for (long long i = 2, j; i <= k; i = j + 1)
    j = k / (k / i), cur = s2(j),
    res = (res - calc(k / i) * (cur - pre) % P) % P, pre = cur;
  return s_map[k] = (res + P) % P;
}

long long solve() {
  long long res = 0, pre = 0, cur;
  for (long long i = 1, j; i <= n; i = j + 1) {
    j = n / (n / i);
    cur = calc(j);
    res = (res + (s3(n / i) * (cur - pre)) % P) % P;
    pre = cur;
  }
  return (res + P) % P;
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  cin >> P >> n;
  inv2 = ksm(2, P - 2), inv6 = ksm(6, P - 2);
  pn = (long long)pow(n, 0.666667);  // n^(2/3)
  prime_work(pn);
  cout << solve();
  return 0;
}

习题

参考资料与注释

  1. Piltz divisor problem - Divisor summatory function - Wikipedia ↩

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：Backl1ght, c-forrest, Enter-tainer, Great-designer, Henry-ZHR, huayucaiji, hsfzLZH1, Ir1d, kenlig, ksyx, Marcythm, MegaOwIer, Menci, Nanarikom, nanmenyangde, ouuan, purple-vine, shawlleyw, sshwy, StudyingFather, Tiphereth-A, Xeonacid
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用