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

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

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

Dirichlet 双曲线法

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

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

$$H(n)=\sum _{k=1}^{n}h(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(\lfloor \frac{n}{x}\rfloor \right)+\sum _{y=1}^{\lfloor y_0\rfloor }F\left(\lfloor \frac{n}{y}\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_0{y}_0=n$，所以由均值不等式可知，当 $x_0=y_0=\sqrt{n}$ 时，就得到最低时间复杂度 $O(\sqrt{n})$．

这并非新的结果．在 Dirichlet 双曲线法的表达式中，令 ，就得到

$$H(n)=\sum _{x=1}^{n}f(x)G\left(\lfloor \frac{n}{x}\rfloor \right).$$

利用 数论分块 的技巧，当 $F,G$ 的点值已知时，该式可以在 $O(\sqrt{n})$ 时间内计算．这实际上和本节得到的算法是几乎等价的：这两种方法需要的已知信息基本类似，计算的表达式也大致相同．

等价性的说明

细究数论分块的计算过程可以发现，实际计算的表达式为

$$H(n)=\sum _{y\in D(n)}(F\left(\lfloor \frac{n}{y}\rfloor \right)-F\left(\lfloor \frac{n}{y+1}\rfloor \right))G(y).$$

其中，$D(n)=\{\lfloor {\displaystyle \frac{n}{x}}\rfloor :1\le x\le n,x\in {\mathbf{N}}_{+}\}$ 是数论分块中的关键点集合．

根据数论分块的性质可知，对于 $x\le \sqrt{n}$，对应分块的高度 $y=\lfloor n/x\rfloor$ 各不相同，于是这些分块的长度都是 $1$，亦即

$$\lfloor \frac{n}{y+1}\rfloor +1=x=\lfloor \frac{n}{y}\rfloor .$$

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

$$I_1=\sum _{x=1}^{\lfloor \sqrt{n}\rfloor }(F(x)-F(x-1))G\left(\lfloor \frac{n}{x}\rfloor \right)=\sum _{x=1}^{\lfloor \sqrt{n}\rfloor }f(x)G\left(\lfloor \frac{n}{x}\rfloor \right).$$

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

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

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

  

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

  $$I_2=\sum _{y=1}^{\lfloor \sqrt{n}\rfloor }F\left(\lfloor \frac{n}{y}\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(\lfloor \frac{n}{x}\rfloor \right)+\sum _{y=1}^{\lfloor \sqrt{n}\rfloor }F\left(\lfloor \frac{n}{y}\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)=\{\lfloor \frac{n}{x}\rfloor :1\le x\le n,x\in {\mathbf{N}}_{+}\}$$

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

块筛及其卷积

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

$${\mathcal{S}}_h(n)=\{H(x):x\in D(n)\}.$$

本节讨论 块筛卷积 问题的计算方法：已知 $f,g$ 的块筛时，求它们的 Dirichlet 卷积 $h=f\ast g$ 的块筛．

朴素算法

朴素算法就是将块筛的计算看作是 $|D(n)|$ 次前缀和点值的计算．这样做的总时间复杂度为

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

但是，这一算法显然过于暴力．集合 $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)$ 且 ，通过 Dirichlet 双曲线法计算 $H$ 在 $x$ 处的点值．

对于一般情形，这样做的时间复杂度为

当 $z={\left({\displaystyle \frac{n}{\log n}}\right)}^{2/3}$ 时，总时间复杂度最小为 $O(n^{2/3}(\log n{)}^{1/3})$．

当然，Dirichlet 卷积 点值计算 的复杂度与 $f,g,h$ 的性质有关．对于 $f,g,h$ 有特殊性质的情形，最优分点和复杂度均略有不同：

• 如果 $f$ 或 $g$ 是积性的，那么当 $z={\left({\displaystyle \frac{n}{\log \log n}}\right)}^{2/3}$ 时，总时间复杂度最小为 $O(n^{2/3}(\log \log n{)}^{1/3})$．
• 如果 $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+\epsilon })$ 时间内计算卷积前缀和的全部必要信息．

快速块筛卷积

前置知识：快速傅里叶变换

注意

初学者可以跳过本节内容．

本节讨论周康阳在 2024 年集训队论文提出的快速块筛卷积算法．它可以在 $O(\sqrt{n}{\log }^{2}n)$ 时间内根据块筛 ${\mathcal{S}}_f$ 和 ${\mathcal{S}}_g$ 计算出它们卷积 ${\mathcal{S}}_h$ 的取值．这一算法不依赖于额外的点值信息和数论函数的积性，但实现较为复杂．

块筛卷积问题希望计算

$$h(z)=\sum _{xy=z}f(x)g(y)$$

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

首先，考虑 的所有点对于这一系列前缀和的贡献，那些 的贡献类似．因为所有贡献必须满足 $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)$ 时间内计算出这部分贡献．

然后，考虑 $\lfloor n/t\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}$ 且 的贡献．所有贡献必须满足 $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)$，可能会遗漏部分贡献．这只会发生在 时，所以有

这相当于

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

$$\sum _{k\le n}{d}_3(k)=nP(\log n)+O(n^{43/96+\epsilon }),$$

其中，$P(\cdot )$ 是二次多项式．因为前两步时间复杂度已经是 $O(\sqrt{n}\log n)$，忽略所有 $o(\sqrt{n}\log n)$ 的项，最后这一部分贡献计算的时间复杂度为

$$O(S\log n\log (S\log n)+\frac{n}{S}{\log }^{2}n).$$

取 $S=\sqrt{n}$，就得到时间复杂度为 $O(\sqrt{n}{\log }^{2}n)$．这也是该算法的总时间复杂度．

参考实现

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}^{n}g(x).$$

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

为此，在 Dirichlet 双曲线法表达式中令 ，就得到

$$H(n)=\sum _{x=1}^{n}f(x)G\left(\lfloor \frac{n}{x}\rfloor \right).$$

直接解出 $G(n)$，就得到

$$G(n)=\frac{1}{f(1)}(H(n)-\sum _{x=2}^{n}f(x)G\left(\lfloor \frac{n}{x}\rfloor \right)).$$

这就是杜教筛的表达式．实际上对于 $x_0\ge 1$，总有

$$G(n)=\frac{1}{f(1)}(H(n)-\sum _{x=2}^{\lfloor x_0\rfloor }f(x)G\left(\lfloor \frac{n}{x}\rfloor \right)-\sum _{y=1}^{\lfloor y_0\rfloor }F\left(\lfloor \frac{n}{y}\rfloor \right)g(y)+F(\lfloor x_0\rfloor )G(\lfloor y_0\rfloor )).$$

无论是哪种形式，它都是一个关于 $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(T_0(z)+\frac{n}{\sqrt{z}}).$$

当 $g$ 是积性函数时，可以应用线性筛，即 $T_0(z)=\mathrm{\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)$．那么，有

利用类似 主定理 的证明思路，可以说明，最后一项主导了该式的增长，且 $T(n)\in \mathrm{\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)$．

参考实现

#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$ 和 $\phi$ 分别是莫比乌斯函数和欧拉函数．求 $S_1(n)=\sum _{i=1}^n\mu (i)$ 和 $S_2(n)=\sum _{i=1}^n\phi (i)$ 的值．数据范围：.

解答

注意到 Dirichlet 卷积关系式：

$$\epsilon =\mu \ast 1,\mathrm{id}=\phi \ast 1.$$

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

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

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

参考实现

#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}^{n}ij\cdot gcd(i,j)(\mathrm{mod}p).$$

数据范围：$n\le {10}^{10}$，$5\times {10}^8\le p\le 1.1\times {10}^9$ 且 $p$ 是质数．

解答

利用 欧拉函数 的性质，做如下反演：

其中，$F(n)={\displaystyle \frac{1}{2}}n(n+1)$．该式可以通过数论分块计算，但需要预处理出 $d^2\phi (d)$ 的前缀和．

为此，可以利用杜教筛．记 $f(n)=({\mathrm{id}}^2\phi )(n)$ 和 $S(n)=\sum _{i=1}^{n}f(i)$．应用杜教筛的关键是构造出函数 $g$，使得 $f\ast g$ 和 $g$ 都可以快速求和．前文已经讨论了用杜教筛预处理 $\phi$ 前缀和的方法，这只需要利用关系：$\mathrm{id}=\phi \ast 1$．相较于 $\phi$，这里的 $f$ 多了一个 ${\mathrm{id}}^2$．由于 $\mathrm{id}$ 是完全积性函数，利用 Dirichlet 卷积的 性质，只需要类似地将每一项都乘以 ${\mathrm{id}}^2$，就得到

$${\mathrm{id}}^3=f\ast {\mathrm{id}}^2.$$

因为 ${\mathrm{id}}^2(n)=n^2$ 和 ${\mathrm{id}}^3(n)=n^3$ 的前缀和都可以在 $O(1)$ 时间内计算，$f$ 的前缀和就可以在 $O(n^{2/3})$ 时间内预处理得到．再加上数论分块，整体时间复杂度仍然是 $O(n^{2/3})$．

另一种将类似积性函数表示为两函数之（在 Dirichlet 卷积意义下的）商的方法是利用 Dirichlet 生成函数．

参考实现

// 不要为了省什么内存把数组开小,会卡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;
}

习题

• AtCoder Xmas Contest 2019 D - Sum of (-1)^f(n)

参考资料与注释

• 任之洲．2016．《积性函数求和的几种方法》．2016 年信息学奥林匹克中国国家队候选队员论文．
• 周康阳．2024．《关于积性函数求和问题的一些进展》．2024 年信息学奥林匹克中国国家队候选队员论文．
• Dirichlet hyperbola method - Wikipedia
• 杜教筛的时空复杂度分析 - riteme.site
• OI 中常用数论函数求和法的简化陈述 by negiizhao - UOJ
• Dirichlet 積と、数論関数の累積和 by maspy
• Dirichlet convolution. Part 1: Fast prefix sum computations by adamant - Codeforces

---

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

---

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