把博客园图标替换成自己的图标
把博客园图标替换成自己的图标end

杜教筛入门

前言

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

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

一波分析

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

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

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

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

转化一下就是这个样子:

\[\sum_{i=1}^nh(i)=\sum_{d=1}^ng(d)·\sum_{i=1}^{\lfloor\frac nd\rfloor}f(i) \]

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

\[\sum_{i=1}^nh(i)=\sum_{d=1}^ng(d)·sum(\lfloor\frac nd\rfloor) \]

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

\[\sum_{i=1}^nh(i)=\sum_{d=2}^ng(d)·sum(\lfloor\frac nd\rfloor)+g(1)·sum(n) \]

移项可得:

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

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

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

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

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

实现方式

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

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

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

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

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

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

一道模板题

这里有一道模板题【洛谷4213】【模板】杜教筛(Sum)

大致题意就是让我们求\(\sum_{i=1}^n\phi(i)\)\(\sum_{i=1}^n\mu(i)\),而且\(N\le2^{31}-1\),不能用线性筛。

所以我们要用杜教筛。

首先考虑如何筛\(\sum_{i=1}^n\phi(i)\)

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

\(id\)函数是很好求的(\(id(n)=n\))。

于是可以得到:

\[\sum_{i=1}^n\phi(i)=\frac{\sum_{i=1}^nid(i)-\sum_{i=2}^nI(d)·(\sum_{i-1}^{\lfloor\frac nd\rfloor}\phi(i))}{I(1)} \]

\(I\)\(id\)两个函数的值代入,就可以得到一个无比简单的式子:

\[\sum_{i=1}^n\phi(i)=\frac{n(n+1)}2-\sum_{d=2}^n\sum_{i=1}^{\lfloor\frac nd\rfloor}\phi(i) \]

\(\sum_{i=1}^n\mu(i)\)也是同理的。

我们可以使用式子\(\mu*I=e\),化简得到:

\[\sum_{i=1}^n\mu(i)=1-\sum_{d=2}^n\sum_{i=1}^{\lfloor\frac nd\rfloor}\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('\n');
    return F.clear(),0;
}
posted @ 2018-10-30 15:05  TheLostWeak  阅读(588)  评论(0编辑  收藏  举报