铃悬的数学小讲堂——杜教筛

铃悬的数学小讲堂——杜教筛

本篇算是【狄利克雷卷积及莫比乌斯反演】篇的后续 qwq。

本篇在读者掌握莫比乌斯反演与狄利克雷卷积的前提下，讲述了杜教筛算法，并伴有若干例题；另外还解释了如何利用贝尔级数更方便的使用杜教筛（要求具有一定生成函数知识）。

如果掌握莫比乌斯反演，也建议看一遍上篇博客；因为可能某些符号不太一样。

-2.符号及规定

同上篇。为方便此处再次写出：

3. 数论函数（见下）用小写粗体字母或普通希腊文字母（如 {\bf f,g,h},\mu,\epsilon ）表示。// QAQ希腊文字母不能粗体

-1.前备知识

至少要理解莫比乌斯反演篇中的知识（即使不知道证明）qwq。

0.数论分块

对于一个像 F(n)=\sum_{i=1}^n f(i)\lfloor\frac ni\rfloor 这样的式子，如果我们可以 O(1) 求出 f(i) 的区间和，那么则可以在 O(\sqrt n) 的复杂度内计算 F(n) 。这依赖于以下引理：

引理 0.1

不同的 \lfloor\frac ni\rfloor 只有 O(\sqrt n) 种。

证明实际上很简单：如果 i\geq\sqrt n，那么 1\leq\lfloor\frac ni\rfloor\leq\sqrt n ，又只取整数值，自然只有 O(\sqrt n) 种取值。否则 i<\sqrt n，但这样的 i 只有 \sqrt n 种，自然也只会有 O(\sqrt n) 种取值。

至于如何找出每一种取值，以及确定每个值对应的 i 的区间，可以参考以下代码：

int ans = 0;
for (int i = 1, j; i <= n; i = j + 1) {
  j = n / (n / i);
  // j 为最大的满足 floor(n/t)=floor(n/i) 的 t
  ans += (n / i) * S(i, j);
}

代码中 S(a,b) 指 \sum_{i=a}^n f(i)。这段代码即可求出 \sum_{i=1}^n f(i)\lfloor\frac ni\rfloor。

至于代码中的注释，则是因为

\left\lfloor\frac nt\right\rfloor\geq\left\lfloor\frac ni\right\rfloor\iff \frac nt\geq\left\lfloor\frac ni\right\rfloor\iff t\leq\frac n{\left\lfloor\frac ni\right\rfloor}\iff t\leq\left\lfloor\frac n{\left\lfloor\frac ni\right\rfloor}\right\rfloor

所以 t 的最大值就是 n/(n/i) 。

如果求和上界是 m 而不是 n （m\leq n），那就吧循环条件写成 i\leq m 并令 j=min(m,n/(n/i)) 即可。\sum_{i=1}^n f(i)g(\lfloor\frac ni\rfloor) 也可以用相同求法求。

同样的，如果求 \sum_{i=1}^{\min(n,m)}f(i)\lfloor\frac ni\rfloor\lfloor\frac mi\rfloor，也只需要令 j=min(n/(n/i),m/(m/i))，复杂度仍旧是 O(\sqrt{\max(n,m)})。

1.杜教筛

如果给定函数 \bf f,g，令 {\bf S}(n)=\sum_{i=1}^n{\bf f}(i)，则有

\sum_{i=1}^n({\bf f*g})(i)=\sum_{i=1}^n\sum_{xy=i}{\bf f}(x){\bf g}(y)=\sum_{y=1}^n{\bf g}(y)\sum_{xy\leq n}{\bf f}(x)=\sum_{y=1}^n{\bf g}(y){\bf S}\left(\left\lfloor\frac ny\right\rfloor\right)

换句话说，

{\bf g}(1){\bf S}(n)=\sum_{i=1}^n({\bf f*g})(i)-\sum_{y=2}^n{\bf g}(y){\bf S}\left(\left\lfloor\frac ny\right\rfloor\right)

这是由上面的等式右侧除了 y=1 一项之外都移项到左边得到的。

由此我们可以得到一个 {\bf S} 的计算方法（假设 ({\bf f*g})(i) 的前缀和与 {\bf g}(i) 的区间和都可以 O(1) 计算，以下计算复杂度时亦作此假设）：

int S(int n) {
  int ans = ???; // \sum_{i=1}^n(f*g)(i)
  for (int i = 2, j; i <= n; i = j + 1) {
    j = n / (n / i);
    ans -= Sg(i, j) * S(n / i);
  }
  ans /= g1;
  return ans;
}

其中 ans 初值为 \sum_{i=1}^n({\bf f*g})(i)，{\bf Sg}(i,j)=\sum_{k=i}^j{\bf g}(k)。

这段代码如何优化复杂度呢？首先估计计算 {\bf S}(N) 的过程中产生了多少递归计算。

以下，我们以 N 指最初给出的 n，以 n 指递归调用产生的 n。

引理1.1

对于任意正整数 x,y,z，有 \left\lfloor\frac{\left\lfloor\frac zx\right\rfloor}y\right\rfloor=\left\lfloor\frac z{xy}\right\rfloor

证明也很简单：假设 z=ax+b,a=cy+d，其中 0\leq b<x,0\leq d<y，那么显然 \left\lfloor\frac zx\right\rfloor=a,\left\lfloor\frac ay\right\rfloor=c，于是等号左边就等于 c。同样的，我们有 z=cxy+(dx+b)，其中 0\leq dx+b\leq(y-1)x+(x-1)=xy-1<xy，所以 \left\lfloor\frac z{xy}\right\rfloor=c。所以等号左右相等。

由此可以发现，上面的代码在计算 {\bf S}(N) 时，产生的（直接的或间接的）递归计算中的 n 都是某个 \left\lfloor\frac Nx\right\rfloor（因为如果这次的 n 是 \lfloor\frac nx\rfloor，下次一定是某个 \left\lfloor\frac{\lfloor\frac nx\rfloor}i\right\rfloor=\lfloor\frac n{xi}\rfloor)，从而根据引理 0.1 只有 O(\sqrt N) 次递归调用（如果 n 相等的只计算一次的话），记忆化即可。

如果略去递归调用的复杂度，显然计算一次 {\bf S}(n) 是 O(\sqrt n)。根据引理 0.1 的证明过程可知递归调用的 n 分别为 1,2,\dots,\sqrt N,\lfloor\frac N1\rfloor,\lfloor\frac N2\rfloor,\dots\lfloor\frac N{\sqrt N}\rfloor。于是复杂度为

O\left(\sum_{i=1}^{\sqrt N}\sqrt i+\sum_{i=1}^{\sqrt N}\sqrt{\left\lfloor\frac Ni\right\rfloor}\right)=O\left(\sum_{i=1}^{\sqrt N}\sqrt{\left\lfloor\frac Ni\right\rfloor}\right)=O\left(\int_1^{\sqrt N}\sqrt{\frac Nx}{\rm d}x\right)=O\left(N^{\frac34}\right)

第一步等号是因为后面的求和项显然比左边的大，所以在大 O 符号中左边的可以省略。第二步等号为积分近似。

但是我们不能满足于此。我们发现如果对于较小的 n 仍然 O(\sqrt n) 计算 {\bf S} 是不明智的，因为大部分情况下我们可以 O(n) 预处理出 {\bf S}(1),{\bf S}(2),\dots{\bf S}(n)（利用线性筛一类）。我们如果预处理出 1\dots N^c(c>\frac12) 的答案，那么需要递归计算的只剩下\lfloor\frac N1\rfloor,\lfloor\frac N2\rfloor,\dots\lfloor\frac N{N^{1-c}}\rfloor，于是复杂度变成

O\left(N^c+\sum_{i=1}^{N^{1-c}}\sqrt{\left\lfloor\frac Ni\right\rfloor}\right)=O\left(N^c+\int_1^{N^{1-c}}\sqrt{\frac Nx}{\rm d}x\right)=O\left(N^c+N^{1-\frac12 c}\right)

取 c=\frac23 时复杂度最优，为 O(N^{\frac23})。

对于记忆化：可以使用 map，但会导致复杂度多一个 \log。

更好的做法是，由于递归调用的 n 总是某个 \lfloor\frac Nx\rfloor，并且 x>n^{\frac13} 的时候都会直接查预处理的表，所以可以直接记 Ans[x] 表示 \lfloor\frac Nx\rfloor 的答案。详细代码例题中会给出。

对于计算 \sum_{i=1}^n({\bf f*g})(i)，由于后面还有 O(\sqrt n) 的递归复杂度，所以这个东西只要在 O(\sqrt n) 时间算完就不影响杜教筛复杂度；或者因为 n 都是 \lfloor\frac Nx\rfloor，所以这个也可以套一层杜教筛。\sum_{i=1}^n{\bf g}(i) 也类似，因为要计算的 n 都是某个 \lfloor\frac Nx\rfloor，所以这个 \bf g 的前缀和也可以再套一层杜教筛qaq。

2.例题

1.P4213 【模板】杜教筛（Sum）

Description

给定 n\leq 2^{31}-1，求：{\bf S}_{\varphi}(n)=\sum_{i=1}^n\varphi(n)，以及 {\bf S}_{\mu}(n)=\sum_{i=1}^n\mu(n)。

Solution

首先两个答案一定在 long long 和 int 范围内，因为 0\leq\varphi(n)\leq n,-1\leq\mu(n)\leq1。

杜教筛的难点在于，给定一个 {\bf f}(n)，要找一个合适的函数 {\bf g}(n) 使 {\bf g},{\bf f*g} 的前缀和都可以快速计算。

对于 \varphi(n)，我们知道 \varphi=\mu*{\bf id},\varphi*{\bf 1}={\bf id}。所以选取 {\bf g}(n)=1 即可，此时 ({\bf f*g})(n)=n，很容易求和。

对于 \mu(n)，选取 {\bf g}(n)=1 也可，此时 ({\bf f*g})=\epsilon(n)=[n=1]，也很容易求和。

Code

下面给出求 {\bf S}_\varphi 的代码：

typedef long long LL;
const int maxn = 2147483647;
const int maxm = 2000000;

LL S1[maxm], S2[1300];
bool vis[1300];
int N;

LL S(int n) {
  if (n < maxm) return S1[n];
  int x = N / n; // 如果存在某个 x 使得 n = floor(N / x)，
                 // 选 x = floor(N / n) 一定可以。
  if (vis[x]) return S2[x];
  vis[x] = true;
  LL &ans = S2[x];
  ans = (LL)n * (n + 1) / 2; // 对 (f*g)(n) = n 求和
  for (int i = 2, j; i <= n; i = j + 1) {
    j = n / (n / i);
    ans -= (j - i + 1) * S(n / i);
  }
  return ans;
}

代码中 maxn, maxm, 1300 分别是 N,N^{\frac23},N^{\frac13} 的最大值，S1[i] 在执行之前应预处理，是当 n\leq N^{\frac23} 时的 {\bf S}_\varphi(n)。S2[x],vis[x] 用来记忆化 {\bf S}_\varphi(\lfloor\frac Nx\rfloor)。

至于代码中的第一句注释，则是因为 \lfloor\frac Nn\rfloor 是使得 xn\leq N\iff\lfloor\frac Nx\rfloor\geq n 的最大正整数 x 。如果这个 x 还不满足 \lfloor\frac Nx\rfloor=n 的话，就不可能存在其他 x 了。

2. P3768 简单的数学题

Description

给定 n，求

\sum_{i=1}^n\sum_{j=1}^n(ij\gcd(i,j))

对 p 取模。n,m\leq10^{10},5\times10^8\leq p\leq1.1\times10^9 且 p 是质数。时限 6s。

Solution

首先有

\begin{aligned}&\sum_{i=1}^n\sum_{j=1}^n(ij\gcd(i,j))\\=&\sum_{i=1}^n\sum_{j=1}^n\left(ij\sum_{d\mid i,d\mid j}\varphi(d)\right)\\=&\sum_{d=1}^{n}\varphi(d)\left(\sum_{1\leq i\leq n,d\mid i}i\right)\left(\sum_{1\leq j\leq n,d\mid j}j\right)\\=&\sum_{d=1}^{n}d^2\varphi(d)S\left(\left\lfloor\frac nd\right\rfloor\right)^2\end{aligned}

式中 S(k)=\sum_{i=1}^k i=\frac{k(k+1)}2。第一步等式是因为 \sum_{d\mid k}\varphi(d)=k。第二步等式是交换了求和顺序，把对 d 的求和放到最外层。第三步等式是把两个括号里的公因数 d 提出来。

根据老套路，我们数论分块，并对每一块的 d^2\varphi(d) 求区间和。

考虑使用杜教筛计算两个区间和相减。但是一个问题是时间复杂度。

假设我们已经有了函数 LL S(LL n) 用来求 d^2\varphi(d) 的前缀和，那么接下来应该这样：

for (int i = 1, j; i <= n && i <= m; i = j + 1) {
  j = n / (n / i);
  LL t = Sum(n / i);
  ans = (ans + (S(j) - S(i - 1)) * t % p * t % p) % p;
}

可以发现 j 始终是某个 \lfloor\frac nx\rfloor ，而 i-1 就是上一次的 j ，自然也是某个 \lfloor\frac nx\rfloor 。从而恰好用到了杜教筛所有的递归计算，仍然只有 O(n^{\frac23}) 的复杂度。

第二个问题是对于 {\bf f}(n)=n^2\varphi(n)，如何找到合适的 {\bf g}。

考虑在数论函数的狄利克雷卷积之外定义点乘 {\bf f\cdot g}，指逐项相乘，即 ({\bf f\cdot g})(n)={\bf f}(n){\bf g}(n)；定义函数 {\bf f}(n) 是完全积性函数，当且仅当对于所有的 n,m 都有 {\bf f}(nm)={\bf f}(n){\bf f}(m)。像 \epsilon(n)=[n=1],{\bf id}^k(n)=n^k 都在此列。

引理 2.1

若 {\bf f} 是完全积性函数，{\bf g, h} 是数论函数，则 ({\bf f\cdot g})*({\bf f\cdot h})={\bf f}\cdot({\bf g*h})。

证明是很显然的：

\begin{aligned}&\left(({\bf f\cdot g})*({\bf f\cdot h})\right)(n)\\=&\sum_{d\mid n}{\bf f}(d){\bf g}(d){\bf f}\left(\frac nd\right){\bf h}\left(\frac nd\right)\\=&{\bf f}(n)\sum_{d\mid n}{\bf g}(d){\bf h}\left(\frac nd\right)\\=&\left({\bf f}\cdot({\bf g*h})\right)(n)\end{aligned}

由此，对于 {\bf f}(n)=n^2\varphi(n)，即 {\bf f}={\bf id}^2\cdot\varphi，可取 {\bf g}={\bf id}^2\cdot{\bf 1}，则根据上述引理有 {\bf f*g}={\bf id}^2\cdot(\varphi*{\bf 1})，即 {\bf f*g}={\bf id}^3。

所以令 {\bf g}(n)=n^2\cdot1=n^2，则立得

{\bf S}(n)=\sum_{i=1}^ni^3-\sum_{i=2}^ni^2{\bf S}\left(\left\lfloor\frac ni\right\rfloor\right)

至于如何对 i^2,i^3 求和，我们有 \sum_{i=1}^ni^3=\frac{n^2(n+1)^2}4，\sum_{i=1}^ni^2=\frac{n(n+1)(2n+1)}6。

具体代码作为练习（雾），其实和上面的 \varphi(n) 的杜教筛是类似的。

3. BZOJ4176 Lucas的数论

这是权限题。

Description

给定 n，求 \sum_{i=1}^n\sum_{j=1}^n{\bf d}(ij)。{\bf d}(k) 是 k 的约数个数。

#### Solution 对于 ${\bf d}(ij)$，考虑某个 $d\mid ij$ 来说，令 $a=\frac i{\gcd(i,d)},b=\frac d{\gcd(i,d)}$，显然 $a\perp b,a\mid i,b\mid j$ （前两个显然，最后一个是因为$b\gcd(i,d)\mid ij\Rightarrow b\mid aj\Rightarrow b\mid j$ 。最后一步是因为 $a\perp b$ ）。另外，任意给定两个数 $a,b$ 使得 $a\perp b,a\mid i,b\mid j$ 都有 $d=\frac iab\mid ij$ ，且 $\frac i{\gcd(i,d)}=a,\frac d{\gcd(i,d)}=b$。 由此可知 $$\sum_{d\mid ij}1=\sum_{a\mid i,b\mid j,a\perp b}1=\sum_{a\mid i,b\mid j}[a\perp b]$$ 所以 $$\begin{aligned}ans=&\sum_{i=1}^n\sum_{j=1}^n{\bf d}(ij)\\=&\sum_{i=1}^n\sum_{j=1}^n\sum_{a\mid i}\sum_{b\mid j}[a\perp b]\\=&\sum_{i=1}^n\sum_{j=1}^n\sum_{a\mid i}\sum_{b\mid j}\sum_{d\mid a,d\mid b}\mu(d)\\=&\sum_{d=1}^n\mu(d)\left(\sum_{d\mid a}\sum_{a\mid i\leq n}1\right)\left(\sum_{d\mid b}\sum_{b\mid j\leq n}1\right)\\=&\sum_{d=1}^n\mu(d)\left(\sum_{a=1}^{\lfloor\frac nd\rfloor}\left\lfloor\frac {\lfloor\frac nd\rfloor}a\right\rfloor\right)^2\end{aligned}$$ 令 $F(n)=\sum_{i=1}^n\left\lfloor\frac ni\right\rfloor=\sum_{i=1}^n\sum_{ij\leq n}1=\sum_{t=1}^n{\bf d}(t)$， $$ans=\sum_{d=1}^n\mu(d)F\left(\left\lfloor\frac nd\right\rfloor\right)^2$$ 对 $d$ 数论分块，则需要快速求出 $\mu(d)$ 的前缀和与 $F\left(\left\lfloor\frac nd\right\rfloor\right)$。前者可以直接杜教筛；后者可以借鉴杜教筛的思想。 考虑到单次计算 $F(n)$ 复杂度为 $O(\sqrt n)$（数论分块），而预处理 $F(1)\dots F(n)$ 复杂度是 $O(n)$（因为可以线性筛出 ${\bf d}$，则可以预处理出 $F(1)\dots F(n^{\frac23})$，后面的 $F$ 依次 $O(\sqrt n)$ 算，复杂度仍然是 $O(n^{\frac23})$。 虽然实际上这道题不预处理的话 $O(n^{\frac34})$ 也能过，但是无所谓了qaq。 ### 4. 奇怪的题目 #### Description 给定 $n$ ，求 $\sum_{i=1}^n(\mu^2*({\bf id}\cdot\mu))(i)$。 其中 $n\leq10^9$，$\mu^2(n)=(\mu(n))^2$。 #### Solution 令 ${\bf f}=\mu^2*({\bf id}\cdot\mu)$，则 ${\bf id}*f=({\bf id}*({\bf id}\cdot\mu))*\mu^2=\mu^2$。所以取 ${\bf g}(n)=n$，因此 $${\bf S}(n)=\sum_{i=1}^n\mu^2(i)-\sum_{i=2}^ni{\bf S}\left(\left\lfloor\frac ni\right\rfloor\right)$$ 剩下的就是如何计算 $\sum_{i=1}^n\mu^2(i)$。 可以发现 $\mu^2(i)=[i\text{不含平方因子}]$。所以定义 ${\bf s}(i)=\max_{d^2\mid i}d$，则 $\mu^2(i)=[{\bf s}(i)=1]$。 容易发现 $d^2\mid i\iff d\mid{\bf s}(i)$， 所以 $$\sum_{i=1}^n\mu^2(i)=\sum_{i=1}^n[{\bf s}(i)=1]=\sum_{i=1}^n\sum_{d^2\mid i}\mu(d)=\sum_{d=1}^{\lfloor\sqrt n\rfloor}\mu(d)\left\lfloor\frac n{d^2}\right\rfloor$$ 可以暴力枚举 $d$，在 $O(\sqrt n)$ 时间内求出。由于杜教筛本来就有 $O(\sqrt n)$ 的计算，这并不会增加复杂度。 ## 3.贝尔级数 有的时候我们不太容易找到杜教筛的合适的 ${\bf g}$，这时候有一个很厉害的东西：贝尔级数。 （Warning: 以下要求大致了解生成函数相关） ### 3.1 定义 对于积性函数 ${\bf f}(n)$ ，定义 ${\bf f}$ 模 $p$ 的贝尔级数为 $${\bf f}_p(x)=\sum_{i=0}^{\infty}{\bf f}(p^i)x^i$$ 例如 $\epsilon_p(x)=1,{\bf 1}_p(x)=\frac1{1-x},{\bf id}^k_p(x)=\frac1{1-p^kx},\mu_p(x)=1-x,{\bf d}_p(x)=\frac1{(1-x)^2},\varphi_p(x)=\frac{1-x}{1-px}$。 一个很好的性质是（实际上就是把狄利克雷卷积变成了一般卷积） $$ ({\bf f*g})_p(x)={\bf f}_p(x){\bf g}_p(x) $$ 接下来我们就要生拼硬凑出一些奇怪的例题了： ### 3.2 例题 令积性函数 ${\bf f}$ 满足 ${\bf f}(p^c)=p^c+[c>0](-1)^c$，求 ${\bf f}$ 的前缀和。 可以发现 ${\bf f}_p(x)=-1+\sum_{i=0}^{\infty}\left(p^ix^i+(-1)^ix^i\right)=\frac1{1-px}+\frac1{1+x}-1$， 所以取 ${\bf g}_p(x)=(1-px)(1+x)$，即得 $$({\bf f*g})_p(x)=1+x+1-px-(1-px)(1+x)=1+px^2$$ 所以问题转化成了如何求 ${\bf g}$ 的前缀和，以及如何求 ${\bf h}_p(x)=1+px^2$ 的前缀和。 显然可知 $(1-px)$ 对应 ${\bf id}\cdot\mu$，而 $(1+x)$ 对应 $\mu^2$。（指 $\mu^2(n)=(\mu(n))^2=|\mu(n)|$）。 所以 ${\bf g}=({\bf id}\cdot\mu)*\mu^2$，其前缀和可以通过之前的那道“奇怪的题目”的做法求出。 而对于 ${\bf h}(n)$ 的前缀和，可以发现只有 $n$ 是完全平方数，且每个质因数的指数都是 2 （等价于 $\mu(\sqrt n)\neq0$）时 ${\bf h}(n)=\sqrt n$，其他时候 ${\bf h}(n)=0$。所以很显然的， $$\sum_{i=1}^n{\bf h}(i)=\sum_{i=1}^{\lfloor\sqrt n\rfloor}i\mu^2(i)$$ 枚举 $i$ 即可在 $O(\sqrt n)$ 时间内算出。所以本题解决了，复杂度是 $O(n^{\frac23})$。 ## 4.总结 ~~本篇文章，是一位在 OI 中数学方面小有名气的选手无私的馈赠。五道杜教筛例题，涵盖了杜教筛从入门到进阶的大部分套路。你可以通过这篇文章，对杜教筛进行较为深入的了解。详细的复杂度证明、精心挑选的例题和各种不同的套路与 trick，能让你对杜教筛有一个较为全面的掌握。作者相信，这篇漂亮的博客，可以给拼搏于 OI 的逐梦之路上的你，提供一个有力的援助。~~ 本文在读者已经掌握莫比乌斯反演的基础上讲述了杜教筛，并提出了若干例题，涵盖了部分杜教筛题目的套路；既可以作为杜教筛的入门文章，又可以让部分选手更进一步地了解杜教筛，更是一篇不错的杜教筛复习博客。