杜教筛入门

Posted 2021-01-16 chenxiaoran666

前言

杜教筛真的是一个十分著名的筛法，它那玄学的(O(n^{frac23}))时间复杂度真的是十分神奇。

它主要用途是求积性函数的前缀和（当然，根据差分思想，求一段区间内的值之和也是很简单的）。

---

一波分析

假设我们要求的是(sum_{i=1}^nf(i))（注意，(f)是一个积性函数）。

则第一步是寻找（或构造）两个函数(h)和(g)使得(h=f*g)（通常我都找不到两个这样的函数），其中(*)表示狄利克雷卷积。

然后我们可以进行一波分析：

[sum_{i=1}^nh(i)=sum_{i=1}^nsum_{d|i}g(d)·f(frac id)]

转化一下就是这个样子：

[sum_{i=1}^nh(i)=sum_{d=1}^ng(d)·sum_{i=1}^{lfloorfrac nd floor}f(i)]

如果用(sum(n))来表示(sum_{i=1}^nf(i))，那么就可以将式子转化成这个样子：

[sum_{i=1}^nh(i)=sum_{d=1}^ng(d)·sum(lfloorfrac nd floor)]

如果我们将右边式子中(d=1)的情况单独提出，就可以得到这样一个式子：

[sum_{i=1}^nh(i)=sum_{d=2}^ng(d)·sum(lfloorfrac nd floor)+g(1)·sum(n)]

移项可得：

[g(1)sum(n)=sum_{i=1}^nh(i)-sum_{d=2}^ng(d)·sum(lfloorfrac nd floor)]

两边同时除以(g(1))可以得到：

[sum(n)=frac{sum_{i=1}^nh(i)-sum_{d=2}^ng(d)·sum(lfloorfrac nd floor)}{g(1)}]

而(sum(n))就是我们所要求的。

所以就不难发现，在(h)函数能够快速求出时（如(e,I,id)等函数就是非常好的选择），我们就能用杜教筛，在(O(n^{frac23}))的时间内求出(sum(n))。

---

实现方式

杜教筛采取了递归的实现方式，而且还使用记忆化加以优化。

根据上面那个式子，我们可以先预处理出(g)函数的前缀和，然后对(sum(lfloorfrac nd floor))进行除法分块。

而除法分块的过程中，我们需要再一次调用(sum)函数，这就是一个递归的过程。

在(n)比较小的时候，我们可以用线性筛预处理出答案。

于是，就可以通过这种方式来实现我们复杂度为(O(n^{frac23}))的杜教筛了。

注意，记忆化的过程可以使用(map)，当然也能使用传说中无比神奇(O(1))查询的(unordered\\_map)，但是据(hl666)大佬的亲身实践，貌似会跑得更慢... ...

---

一道模板题

这里有一道模板题：【洛谷4213】【模板】杜教筛（Sum）。

大致题意就是让我们求(sum_{i=1}^nphi(i))和(sum_{i=1}^nmu(i))，而且(Nle2^{31}-1)，不能用线性筛。

所以我们要用杜教筛。

首先考虑如何筛(sum_{i=1}^nphi(i))。

关于(phi)函数，有这样一个式子(phi*I=id)。

而(id)函数是很好求的（(id(n)=n)）。

于是可以得到：

[sum_{i=1}^nphi(i)=frac{sum_{i=1}^nid(i)-sum_{i=2}^nI(d)·(sum_{i-1}^{lfloorfrac nd floor}phi(i))}{I(1)}]

将(I)和(id)两个函数的值代入，就可以得到一个无比简单的式子：

[sum_{i=1}^nphi(i)=frac{n(n+1)}2-sum_{d=2}^nsum_{i=1}^{lfloorfrac nd floor}phi(i)]

筛(sum_{i=1}^nmu(i))也是同理的。

我们可以使用式子(mu*I=e)，化简得到：

[sum_{i=1}^nmu(i)=1-sum_{d=2}^nsum_{i=1}^{lfloorfrac nd floor}mu(i)]

这样转化之后就可以用杜教筛来筛了。

代码如下：

#include<bits/stdc++.h>
#define ten(x) (((x)<<3)+((x)<<1))
#define LL long long
using namespace std;
int n;
class FIO
{
    private:
        #define Fsize 100000
        #define tc() (A==B&&(B=(A=Fin)+fread(Fin,1,100000,stdin),A==B)?EOF:*A++)
        #define pc(ch) (FoutSize<Fsize?Fout[FoutSize++]=ch:(fwrite(Fout,1,Fsize,stdout),Fout[(FoutSize=0)++]=ch))
        int Top,FoutSize;char ch,*A,*B,Fin[Fsize],Fout[Fsize],Stack[Fsize];
    public:
        inline void read(int &x) {x=0;while(!isdigit(ch=tc()));while(x=ten(x)+(ch&15),isdigit(ch=tc()));}
        inline void write(LL x) {if(!x) return (void)(pc('0'));if(x<0) pc('-'),x=-x;while(x) Stack[++Top]=x%10+48,x/=10;while(Top) pc(Stack[Top--]);}
        inline void write_char(char x) {pc(x);}
        inline void clear() {fwrite(Fout,1,FoutSize,stdout);}
}F;
class Class_DuSieve//杜教筛
{
    private:
        #define Size 5000000
        int Prime_cnt,Prime[Size+5],phi[Size+5],mu[Size+5],SumMu[Size+5];LL SumPhi[Size+5];bool IsNotPrime[Size+5];
        map<LL,LL> MemmoryPhi;map<int,int> MemmoryMu;//用map存储答案
    public:
        Class_DuSieve()//预处理
        {
            register int i,j;
            for(phi[1]=mu[1]=1,i=2;i<=Size;++i)//先用一遍线性筛预处理
            {
                if(!IsNotPrime[i]) Prime[++Prime_cnt]=i,phi[i]=i-1,mu[i]=-1;
                for(j=1;j<=Prime_cnt&&1LL*i*Prime[j]<=Size;++j) 
                    if(IsNotPrime[i*Prime[j]]=1,i%Prime[j]) phi[i*Prime[j]]=phi[i]*(Prime[j]-1),mu[i*Prime[j]]=-mu[i];else {phi[i*Prime[j]]=phi[i]*Prime[j];break;}
            }
            for(i=1;i<=Size;++i) SumPhi[i]=SumPhi[i-1]+phi[i],SumMu[i]=SumMu[i-1]+mu[i];//求前缀和
        }
        inline LL GetPhi(int x)//求phi的前缀和
        {
            if(x<=Size) return SumPhi[x];//如果已经线性筛求出过答案，就返回已经求出的答案
            if(MemmoryPhi[x]) return MemmoryPhi[x];//如果记忆化过，返回答案
            register int l,r;register LL ans=1LL*x*(x+1)>>1;
            for(l=2;l<=x;l=r+1) r=x/(x/l),ans-=(r-l+1)*GetPhi(x/l);//除法分块，递归调用
            return MemmoryPhi[x]=ans;//记忆化
        }
        inline int GetMu(int x)//求mu的前缀和（和上面的函数基本上一样）
        {
            if(x<=Size) return SumMu[x];
            if(MemmoryMu[x]) return MemmoryMu[x];
            register int l,r,ans=1;
            for(l=2;l<=x;l=r+1) r=x/(x/l),ans-=(r-l+1)*GetMu(x/l);
            return MemmoryMu[x]=ans;
        }
}DuSieve;
int main()
{
    register int T;F.read(T);
    while(T--) F.read(n),F.write(DuSieve.GetPhi(n)),F.write_char(' '),F.write(DuSieve.GetMu(n)),F.write_char('
');
    return F.clear(),0;
}