浅谈数论函数求和法：杜教筛

XeCtera · 2024-02-28 01:29:10 · 题解

0. 前言

杜教筛是一种常用的数论函数求和方法。

本文尝试介绍如下内容：

• 基于块筛结构，更清晰地呈现杜教筛的核心思想与算法流程；

• 给出常数优秀的、非递归式的杜教筛实现；

• 给出杜教筛复杂度的正确证明；

• 网上已经有一种错误证法广泛流传，将指出它为什么错。

1. Dirichlet 双曲线法

对于数论函数 f，记 \displaystyle{\rm S}f(n)=\sum_{i\le n}f(i)。

我们常要解决的一类问题是：对于两个数论函数 f,g，计算其 Dirichlet 卷积 f*g 的前缀和。即：

可引入 Dirichlet 双曲线法对其进行快速计算：

若 x,y>0 且 xy=n，则

\displaystyle\sum_{ij\le n}f(i)g(j)=\sum_{i\le x}f(i){\rm S}g(n/i)+\sum_{j\le y}g(j){\rm S}f(n/j)-{\rm S}f(x){\rm S}g(y)

（这里用 / 表示整除）

证明是简单的：为计算左式，考虑分别计算 i\le x 和 j\le y 的部分，再减去公共部分，即得右式。

通常取 x=y=\sqrt n，这时能以最优的 O(\sqrt n) 复杂度计算右式。

2. 块筛

对于数论函数 f，将所有 {\rm S}f(i),{\rm S}f(n/i)~(1\le i\le\sqrt n) 的值称为 f（在 n 处）的块筛。

块筛包含了所有 {\rm S}f(n/i)~(1\le i\le n) 的值。这是因为 i>\sqrt n 时 (n/i)<\sqrt n。

回顾前文 Dirichlet 双曲线法的式子，可以发现，计算其中右式时，我们需要用到的就是 f,g 的块筛。

换言之，Dirichlet 双曲线法的作用在于：

已知 f,g 块筛的前提下，能以 O(\sqrt n) 的复杂度得到 {\rm S}(f*g)(n) 的值。

3. 块筛卷积

考虑如下问题：

假设 f*g=h，且 f,g 的块筛已知，试求 h 的块筛。

使用 Dirichlet 双曲线法，直接对 h 的块筛的每个位置分别求解。

注意到经典性质：对正整数 n,i,j，有 (n/i)/j=n/(ij)。

这保证了 f 在 n/i 处的块筛被 f 在 n 处的块筛包含。

也即 f,g 的块筛提供了所有需要的信息，因此做法是可行的。

4. 杜教筛

再考虑如下问题：

假设 f*g=h，且 g,h 的块筛已知，试求 {\rm S}f(n)。

考察 Dirichlet 双曲线法给出的式子：

这里 j=1 时出现了 {\rm S}f(n) 项，通过移项将其解出：

容易发现求解 {\rm S}f(n) 将依赖于 {\rm S}f(n/j) 的值。

为此，考虑如下递推流程：

1. 枚举 k=1,2,\cdots,\lfloor\sqrt n\rfloor，依次求解 {\rm S}f(k)；

2. 枚举 k=\lfloor\sqrt n\rfloor,\cdots,2,1，依次求解 {\rm S}f(n/k)。

这样就求出了 f 的块筛，也就得到所需的 {\rm S}f(n)。

5. 复杂度分析

上述两个问题的算法时间复杂度均为

实际上，对于规模较小的部分往往有更好的方法计算。

设定阈值 B~(B\ge \sqrt n)，则 1\sim B 的部分可能存在如下优化（因具体情况而异）：

• 暴力计算 Dirichlet 卷积，复杂度为 O(B\log B)；

• 使用 Dirichlet 前缀和技术，复杂度为 O(B\log\log B)；

• 直接用线性筛预处理，复杂度为 O(B)。

剩余部分仍用杜教筛计算，其复杂度变为

假设 1\sim B 部分能用线性筛预处理，则总时间复杂度为

取 B=n^{\frac23} 得到最优时间复杂度 O(n^{\frac23})。

此时的空间复杂度也为 O(n^{\frac23})。

6. 实现

以洛谷P4213 【模板】杜教筛为例：要求 {\rm S}\varphi(n),{\rm S}\mu(n)。

构造卷积式 \mu*1=\varepsilon，则可用杜教筛求出 \mu 的块筛。

再注意到 \varphi=\mu*{\rm id}，用 Dirichlet 双曲线法即可求出 {\rm S}\varphi(n)。

具体实现可以参考下方代码。

#include<bits/stdc++.h>
#define rep(i,l,r) for (int i=l; i<=r; ++i)
using namespace std;
typedef long long i64;

const int N=2e6,M=5e4;
int T,mx,k,n[12],p[160000];
bitset<N>v;
int mu[N],smu[N],Smu[M];

void init(int n) { //线性筛预处理
    mu[1]=1;
    rep(i,2,n) {
        if (!v[i]) mu[i]=-1,p[++k]=i;
        for (int j=1,t; j<=k&&(t=i*p[j])<=n; ++j) {
            v[t]=1;
            if (i%p[j]==0) break;
            mu[t]=-mu[i];
        }
    }
    rep(i,1,n) smu[i]=smu[i-1]+mu[i];
}
inline i64 S(int x) { return x*(x+1ll)>>1; }

void solve(int n) {
    int sq=sqrt(n);
    for (int i=n/mx+1; i<=sq; ++i) Smu[i]=smu[n/i];
    for (int i=n/mx; i; --i) { //杜教筛
        int m=n/i,B=sqrt(m);
        int s=1+smu[B]*B-m,t,t2;
        rep(j,2,B)
            t=i*j,t2=m/j,
            s-=(t>sq?smu[t2]:Smu[t])+t2*mu[j];
        Smu[i]=s;
    }
    i64 Sphi=-S(sq)*smu[sq];
    rep(i,1,sq) Sphi+=1ll*i*Smu[i]+mu[i]*S(n/i);
    cout<<Sphi<<' '<<Smu[1]<<'\n';
}

int main() {
    ios::sync_with_stdio(0);
    cin.tie(0); cin>>T;
    rep(i,1,T) cin>>n[i];
    mx=pow(*max_element(n+1,n+T+1),2./3);
    init(mx);
    rep(i,1,T) solve(n[i]);
    return 0;
}

7. 另一种复杂度分析的勘误

此证法认为，杜教筛的复杂度由下式给出：

并由此得到：

这里有一步所谓“舍去高阶小量”的处理，实际上是错误的。

下面用简单的放缩来说明 T(n)=\Omega(n)，从而否定此证法。

当 n\le 100 时，有 T(n)\ge 1\ge \dfrac n{100} 。

假设 n\le k 时均有 100\cdot T(n)\ge n（其中 k\ge 100 ），则当 n=k+1 时，有

从而可归纳证得 100\cdot T(n)\ge n 对任意正整数 n 成立。

即得 T(n)=\Omega(n)，而非我们希望的 O(n^{\frac34})。

对于杜教筛的递归式实现的复杂度，正确的分析是：

• 应当注意到这一实现是带记忆化的。
  这保证了递归过程中每个 {\rm S}f(n/i) 只被计算一次，与递推式实现等价。
  因此两种实现具有相同的复杂度，即 O(n^{\frac34})。

• 若不带记忆化，我们已经证明这时的复杂度 T(n)=\Omega(n)，不符预期。

相比之下，递推式实现的复杂度更为直观，并且常数也较优秀，可以认为是更好的实现方法。

8. 后记

参考文献：negiizhao - OI中常用数论函数求和法的简化陈述

感谢 @渐变色 提供的帮助和指导。