杜教筛
如果您对下文中的某些步奏或名词尚未了解,可以参考我的前一篇文章 Here
杜教筛用于求积性函数前缀和
对于积性函数$ f$,令$ S(n)=\sum\limits_{i=1}^nf(i)$,求$ S(n)$ $( n<=10^9)$
我们找另一个积性函数$ g$, 对$ g$和$ f$做狄利克雷卷积,令结果为$ h$,
则有$ h(i)=\sum\limits_{d|i}g(d)f(\frac{d}{i})$
对卷积结果$ h$计算前缀和$ t$,得$ t(n)=\sum\limits_{i=1}^n \sum\limits_{d|i}g(d)f(\frac{i}{d})$
改为枚举$ d$,得$ t(n)=\sum\limits_{d=1}^ng(d) \sum\limits_{i=1}^\frac{n}{d}f(i)$
后半部分用$ S$代替得$ t(n)=\sum\limits_{d=1}^ng(d)s(\frac{n}{d})$
显然有$ g(1)S(n)=\sum\limits_{i=1}^ng(i)S(\frac{n}{i})-\sum\limits_{i=2}^ng(i)S(\left\lfloor\frac{n}{i}\right\rfloor)$
由于$ f$是积性函数,$ g(1)$恒等于1,化简可得
$ S(n)=t(n)-\sum\limits_{i=2}^ng(i)S(\left\lfloor\frac{n}{i}\right\rfloor)$
发现式子最右侧明显可以数论分块
那么如果能够找到一个合理的积性函数$ g$,使得$ g$的前缀和与狄利克雷卷积的前缀和$ t$都很好计算,那么就可以快速递归计算前缀和$ S$了
Samples
求 $ \sum\limits_{i=1}^n \mu(i)$ $ (i<=10^9)$
我们需要找到一个合适的积性函数$ g$使得可以快速计算出$ \sum\limits_{i=1}^n g(i)$ 以及$ \sum\limits_{i=1}^n (g*\mu)(i)$
由于$ \mu*1=e$,我们可以令$ g=1$
那么化简模型式可得$ S(n)=1-\sum\limits_{i=2}^nS(\frac{n}{i})$
我们可以线筛n<=500万的$ S(n)$,对于n>500万的进行记忆化搜索
总复杂度O(能过) $ O(n^\frac{2}{3})$
tips:记忆化的时候不需要使用map或者hash,因为搜索到的元素一定为$ n$的因数,所以对于$ 5000000<x<=10^9$的x我们只需要把答案存储在$ \frac{n}{x}$的位置上即可
求 $ \sum\limits_{i=1}^n \phi(i)$ $ (i<=10^9)$
同理,由于$ \phi*1=id$ 我们继续令$ g=1$
化简模型式可得$ S(n)=\frac{n(n+1)}{2}-\sum\limits_{i=2}^nS(\frac{n}{i})$
同理线筛+记忆化计算即可
模板题:BZOJ3944 Here
代码:
#include<cmath> #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #define rt register int #define l putchar('\n') #define ll long long #define r read() #define N 3000010 using namespace std; inline ll read() { register ll x = 0; char zf = 1; char ch; while (ch != '-' && !isdigit(ch)) ch = getchar(); if (ch == '-') zf = -1, ch = getchar(); while (isdigit(ch)) x = x * 10 + ch - '0', ch = getchar(); return x * zf; } void write(ll y) { if (y < 0) putchar('-'), y = -y; if (y > 9) write(y / 10); putchar(y % 10 + '0'); } inline void writeln(const ll y){write(y);putchar('\n');}int i,j,k,m,n,x,y,z,cnt; pair<ll,int>ret[10010];//记忆化 ll phi[N+10];int u[N+10]; pair<ll,int> anss(const int all,const unsigned x)//all用于记忆化,x表示求S(x),同时返回phi值和mu值 { if(x<N)return {phi[x],u[x]}; if(ret[all/x].first)return ret[all/x]; pair<ll,int> ans;ans={(ll)x*(x+1)>>1,1}; for(register int i=2;i<=x;)//数论分块递归计算 { const unsigned u=x/(x/i); const pair<ll,int>s=anss(all,x/i); ans.first-=(u-i+1)*s.first; ans.second-=(u-i+1)*s.second; i=u+1; } ret[all/x]=ans; return ans; } ll ss[800010];int top;bool pri[N+10]; int main() { phi[1]=u[1]=1; for(rt i=2;i<=N;i++)//线筛初始化 { if(!pri[i])ss[++top]=i,phi[i]=i-1,u[i]=-1; for(rt j=1;j<=top&&ss[j]*i<=N;j++) { pri[i*ss[j]]=1; phi[i*ss[j]]=phi[i]*(ss[j]-1); u[i*ss[j]]-=u[i]; if(i%ss[j]==0) { phi[i*ss[j]]=phi[i]*ss[j]; u[i*ss[j]]=0; break; } } } for(rt i=1;i<=N;i++)phi[i]+=phi[i-1];//计算前缀和 for(rt i=1;i<=N;i++)u[i]+=u[i-1]; for(rt t=r;t;t--) { x=r;memset(ret,0,sizeof(ret)); pair<ll,int>ans=anss(x,x); printf("%lld %d\n",ans.first,ans.second); } return 0; }