杜教筛学习笔记

· · 个人记录

本来想放在 数论初步 里的，但是那篇文章实在太长了，编辑的时候标签页崩溃了。。。

上面的文章为前置知识。

参考资料：

OI Wiki

模板题题解

理论基础

杜教筛是一种可以在 \mathcal O(n^{\frac 23}) 的时间内算出一个积性函数的前缀和的算法。

数论基础

设有一个积性函数 \mathbf f，它的前缀和为 \mathbf S(n)=\sum_{i=1}^n\mathbf f(i)。再找一个积性函数 \mathbf g，考虑它们狄利克雷卷积的前缀和，可得：

\begin{aligned} \sum_{i=1}^n\mathbf{(f*g)}(i)&=\sum_{i=1}^n\sum_{d\mid i}\mathbf f(\frac id)\mathbf g(d)\\ &=\sum_{d=1}^n\mathbf g(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}\mathbf f(i)\\ &=\sum_{d=1}^n\mathbf g(d)\mathbf S(\lfloor\frac nd\rfloor) \end{aligned}

则有

\begin{aligned} \mathbf g(1)\mathbf S(n)&=\sum_{d=1}^n\mathbf g(d)\mathbf S(\lfloor\frac nd\rfloor)-\sum_{d=2}^n\mathbf g(d)\mathbf S(\lfloor\frac nd\rfloor)\\ &=\sum_{i=1}^n\mathbf{(f*g)}(i)-\sum_{d=2}^n\mathbf g(d)\mathbf S(\lfloor\frac nd\rfloor) \end{aligned}

由于 \mathbf g 是积性函数，所以 \mathbf g(1)=1。因此

\mathbf S(n)=\sum_{i=1}^n\mathbf{(f*g)}(i)-\sum_{d=2}^n\mathbf g(d)\mathbf S(\lfloor\frac nd\rfloor)

如果我们可以快速计算 \mathbf {f*g} 和 \mathbf g 的前缀和，就可以利用数论分块递归计算 \mathbf S。另外，为了优化时间复杂度，可以 \mathcal O(n^{\frac 23}) 预处理出前 n^{\frac 23} 个数的答案，之后再用杜教筛。然而我不会证明时间复杂度（不会微积分太不方便了）。 现在会了。见下文。

对于记录答案，可以用 map，但更好的做法是观察性质。我们知道

\left\lfloor\frac{\lfloor\frac xy\rfloor}z\right\rfloor=\left\lfloor\frac x{yz}\right\rfloor

证明：设 x=ay+b(b<y)，a=cz+d(d<z)，则左式显然为 c，而 x=(cz+d)y+b=cyz+dy+b，又因为 dy+b\le (z-1)y+y-1<yz，所以右式也等于 c，得证。

所以每次递归计算的一定是某个 \mathbf S(\lfloor\dfrac nx\rfloor)。由于当 x>n^{\frac 13} 时会直接获取预处理出的答案，所以需要记忆化的数一定少于 n^{\frac 13}，因此可以开一个哈希表记录答案。

所以我们要设计一个哈希函数 H，保证对于任意的 \lfloor\dfrac nx\rfloor\ne \lfloor\dfrac ny\rfloor\land x,y\le n^{\frac 13}，都有 H(\lfloor\dfrac nx\rfloor)\ne H(\lfloor\dfrac ny\rfloor)。可以取 H(i)=\lfloor\dfrac ni\rfloor。因为

k=H(\lfloor\dfrac nx\rfloor)=\left\lfloor\frac n{\lfloor\frac nx\rfloor}\right\rfloor

是使得取值为 \lfloor\dfrac nk\rfloor=\lfloor\dfrac nx\rfloor 的最大的 k（证明见前置知识中的数论分块部分），所以 \lfloor\dfrac nx\rfloor 和 H(\lfloor\dfrac nx\rfloor) 构成了一个双射。

时间复杂度证明

我们已经证明了 \mathbf S 的自变量取值一定是某个 \lfloor\dfrac nx\rfloor，所以只有 O(\sqrt n) 种状态。最多的时候状态集合为 \{x\mid x\in\mathbb N\cap[1,\sqrt n]\}\cap\{\lfloor\dfrac nx\rfloor\mid x\in\mathbb N\cap[1,\sqrt n]\}。而计算一个状态 x 的时间复杂度为 O(\sqrt x)。所以如果不预处理，总的时间复杂度为

\begin{aligned} O(\sum_{i=1}^{\sqrt n}\sqrt i+\sum_{i=1}^{\sqrt n}\sqrt{\lfloor\frac ni\rfloor})&=O(\sum_{i=1}^{\sqrt n}\sqrt{\lfloor\frac ni\rfloor})\\ &=O(\int_1^{\sqrt n}\sqrt{\frac nx}\mathrm dx)\\ &=O(\sqrt n\int_1^{\sqrt n}\frac 1{\sqrt x}\mathrm dx)\\ &=O(n^{\frac 12}\int_1^{\sqrt n}x^{-\frac 12}\mathrm dx)\\ &=O(n^{\frac 34}) \end{aligned}

如果线性预处理出 [1,n^a] 的答案且 a>\dfrac 12，则时间复杂度为

\begin{aligned} O(n^a+\sum_{i=1}^{n^{1-a}}\sqrt{\lfloor\frac ni\rfloor})&=O(n^a+\int_1^{n^{1-a}}\sqrt{\frac nx}\mathrm dx)\\ &=O(n^a+n^{\frac 12}\int_1^{n^{1-a}}x^{-\frac 12}\mathrm dx)\\ &=O(n^a+n^{1-\frac{a}{2}}) \end{aligned}

可以看出，a 取 \dfrac 23 时最优，时间复杂度为 O(n^{\frac 23})。

代码实现

模板题，提交记录。

知道杜教筛的理论后，只要知道 \varphi*\mathbf 1=\mathbf{id}，\mu*\mathbf 1=\varepsilon，以及数据的毒瘤程度，就可以写出代码了：

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
using namespace std;
typedef long long ll;
const int N=1.7e6,sN=1300;
int _n,v[N],phi[N],mu[N];
ll preP[N],preM[N],sumP[sN],sumM[sN];
vector<int>pr;
inline void init(){
    preP[1]=phi[1]=1;preM[1]=mu[1]=1;
    for(int i=2;i<N;i++){
        if(!v[i])pr.push_back(v[i]=i),mu[i]=-1,phi[i]=i-1;
        for(int p:pr){
            if(p>v[i]||p>(N-1)/i)break;
            v[i*p]=p;if(v[i]==p)phi[i*p]=phi[i]*p,mu[i*p]=0;
            else phi[i*p]=phi[i]*phi[p],mu[i*p]=mu[i]*mu[p];
        }
        preP[i]=preP[i-1]+phi[i];preM[i]=preM[i-1]+mu[i];
    }
}
ll getSumPhi(int n){
    if(n<N)return preP[n];
    if(sumP[_n/n])return sumP[_n/n];
    ll ret=((ll)n*(n+1))>>1;for(int l=2,r=1;l<=n;l=r+1){
        r=n/(n/l);
        ret-=getSumPhi(n/l)*(r-l+1);
    }
    return sumP[_n/n]=ret;
}
ll getSumMu(int n){
    if(n<N)return preM[n];
    if(sumM[_n/n])return sumM[_n/n];
    ll ret=1;for(int l=2,r=1;l<=n;l=r+1){
        r=n/(n/l);
        ret-=getSumMu(n/l)*(r-l+1);
    }
    return sumM[_n/n]=ret;
}
inline void ct(){
    scanf("%d",&_n);
    if(_n==2147483647){// 当 n=(1<<31)-1 时，int 会爆，所以要特判
        printf("%lld %d\n",1401784457568941916,9569);
        return;
    }
    memset(sumP,0,sizeof(sumP));
    memset(sumM,0,sizeof(sumM));
    printf("%lld %lld\n",getSumPhi(_n),getSumMu(_n));
}
int main(){
    init();int T;scanf("%d",&T);
    while(T--)ct();
    return 0;
}