杜教筛入门
前言
杜教筛真的是一个十分著名的筛法,它那玄学的\(O(n^{\frac23})\)时间复杂度真的是十分神奇。
它主要用途是求积性函数的前缀和(当然,根据差分思想,求一段区间内的值之和也是很简单的)。
一波分析
假设我们要求的是\(\sum_{i=1}^nf(i)\)(注意,\(f\)是一个积性函数)。
则第一步是寻找(或构造)两个函数\(h\)和\(g\)使得\(h=f*g\)(通常我都找不到两个这样的函数),其中\(*\)表示狄利克雷卷积。
然后我们可以进行一波分析:
转化一下就是这个样子:
如果用\(sum(n)\)来表示\(\sum_{i=1}^nf(i)\),那么就可以将式子转化成这个样子:
如果我们将右边式子中\(d=1\)的情况单独提出,就可以得到这样一个式子:
移项可得:
两边同时除以\(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\))。
于是可以得到:
将\(I\)和\(id\)两个函数的值代入,就可以得到一个无比简单的式子:
筛\(\sum_{i=1}^n\mu(i)\)也是同理的。
我们可以使用式子\(\mu*I=e\),化简得到:
这样转化之后就可以用杜教筛来筛了。
代码如下:
#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;
}