浅(kou)谈(hu)杜教筛
引入
让我们从一道例题开始
这道题目需要我们求出\(\varphi(n)\)与\(\mu(n)\)这两个函数的前缀和
显然，这可以通过线性筛\(\mathcal O(n)\)完成，但这道题目的\(n\)的范围是\(\le 2^{31}\)，线性筛无法完成，我们需要一种更高效的算法————杜教筛
test
前置知识
1.线性筛
2.积性函数：如果函数\(f(n)\)满足当\(gcd(n,m)=1\)时\(f(nm)=f(n)f(m)\)，那么我们称\(f(n)\)为积性函数。
3.完全积性函数：就是积性函数去掉那个互质的条件。
4.艾弗森约定：对于布尔型变量\(x\)：\([x]\)表示\(x\)为真时取\(1\)，为假时取\(0\)。
5.一些常见积性函数
- 元函数\(\epsilon(n)=[n=1]\)，是完全积性的。
- 幂函数\(id_k(n)=n^k\)，是完全积性的。
- 单位函数\(I(n)=id_0(n)\)。
- 约数幂函数\(\sigma_k(n)=\sum_{d|n}d^k\)，是积性的。
- 欧拉函数\(\varphi(n)=\sum_{d≤n}[gcd(d,n)=1]\)，是积性的.
- 莫比乌斯函数\(μ(n)\)满足\(\sum_{d|n}μ(d)=[n=1]=\epsilon(n)\)是积性的
6.狄利克雷卷积
设\(f,g\)是两个数论函数，它们的狄利克雷卷积是：\((f \otimes g)(n)=\sum_{d|n} f(d) g(\frac nd)\)
性质：满足交换律，结合律
概述
杜教筛，因由圈内著名的杜老师（杜瑜皓）引进而得名，可以以\(\mathcal O(n^{\frac 23})\)的时间复杂度求解一类积性函数的前缀和
假设\(f(n)\)是我们需要求出的积性函数，令
我们再寻找一个积性函数\(g(n)\)，令\(h(n)=(f \otimes g)(n)\),那么
移项后得到
那么，如果我们能够较快地求出\(g(n)\)与\(h(n)\)的前缀和，那么减号右边的柿子可以整除分块后递归处理，我们就能够在低于线性的复杂度下完成对\(S(n)\)的求解
可以证明直接递归并使用\(unordered\_map\)或哈希表记忆化后该算法的复杂度是\(\mathcal O(n^{\frac 34})\)的
如果在此基础上，先线性筛求出前\(n^{\frac 23}\)项\(f\)的前缀和，那么我们就可以进一步优化到\(\mathcal O(n^{\frac 23})\)
这就是杜教筛的核心思想
例题
洛谷P4213【模板】杜教筛（Sum）
就是引入中的那道题。。。
sol.
根据经典卷积柿子\(\mu \otimes I=\epsilon,\varphi \otimes I=id_1\)
而\(\epsilon,I,id_1\)的前缀和都可以\(\mathcal O(1)\)求出
于是直接上模板就行了
#include<bits/stdc++.h>
using namespace std;
const int N=5000000;
typedef long long ll;
typedef unsigned int uint;
int T;
uint n;
ll phi[N+10],mu[N+10];
unordered_map<uint,ll> ans1,ans2;
inline ll sumphi(uint m){
if(m<=N) return phi[m];
if(ans1.count(m)) return ans1[m];
ll ret=1ll*m*(m+1)/2;uint r;
for(uint l=2;l<=m;l=r+1){
r=m/(m/l);
ret-=1ll*(r-l+1)*sumphi(m/l);
}
return ans1[m]=ret;
}
inline ll summu(uint m){
if(m<=N) return mu[m];
if(ans2.count(m)) return ans2[m];
ll ret=1;uint r;
for(uint l=2;l<=m;l=r+1){
r=m/(m/l);
ret-=1ll*(r-l+1)*summu(m/l);
}
return ans2[m]=ret;
}
int p[N+10],f[N+10],pcnt;
inline void init(){
mu[1]=phi[1]=1;
for(int i=2;i<=N;++i){
if(!f[i]) p[++pcnt]=i,phi[i]=i-1,mu[i]=-1;
for(int j=1;j<=pcnt&&i*p[j]<=N;++j){
int d=i*p[j];f[d]=1;
if(i%p[j]) mu[d]=-mu[i],phi[d]=phi[i]*phi[p[j]];
else{
mu[d]=0;phi[d]=phi[i]*p[j];
break;
}
}
}
for(int i=2;i<=N;++i) phi[i]+=phi[i-1],mu[i]+=mu[i-1];
}
int main(){
scanf("%d",&T);
init();
while(T--){
scanf("%u",&n);
printf("%lld %lld\n",sumphi(n),summu(n));
}
return 0;
}
洛谷P3768 简单的数学题
题目要求
sol.
因为
其中\((\sum_{i=1}^{\lceil \frac nd \rceil} i)^2\)可以数论分块后\(\mathcal O(\sqrt{n})\)求出
于是我们现在就只需要求\(f(x)=\varphi(x)x^2\)的前缀和了
因为
又因为\(d*\frac nd=n\)，于是
于是令\(g=n^2,h=n^3\)即可
显然\(g\)的前缀和\(=\frac {n(n+1)(2*n+1)}{6}\),\(h\)的前缀和\(=(\frac{n(n+1)}{2})^2\)，可以\(\mathcal O(1)\)求出
于是上杜教筛即可
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e6;
int phi[N+10],f[N+10],p[N+10],pcnt,mod;
ll n;
inline int add(int x,int y){return (x+y>=mod)?x+y-mod:x+y;}
inline int dec(int x,int y){return (x-y<0)?x-y+mod:x-y;}
inline void init(){
phi[1]=1;
for(int i=2;i<=N;++i){
if(!f[i]) p[++pcnt]=i,phi[i]=i-1;
for(int j=1;j<=pcnt&&i*p[j]<=N;++j){
f[i*p[j]]=1;
if(i%p[j]) phi[i*p[j]]=1ll*phi[i]*phi[p[j]]%mod;
else{
phi[i*p[j]]=1ll*phi[i]*p[j]%mod;
break;
}
}
}
for(int i=1;i<=N;++i) f[i]=add(1ll*phi[i]*i%mod*i%mod,f[i-1]);
}
unordered_map<ll,int> ans;
inline int ksm(int x,int y){
int ret=1;
for(;y;x=1ll*x*x%mod,y>>=1) if(y&1) ret=1ll*ret*x%mod;
return ret;
}
int iv2,iv6;
inline int sum(ll x){x%=mod;return 1ll*x*(x+1)%mod*iv2%mod;}
inline int sumpow(ll x){x%=mod;return 1ll*x*(x+1)%mod*(2*x+1)%mod*iv6%mod;}
inline int sumf(ll x){
if(x<=N) return f[x];
if(ans.count(x)) return ans[x];
ll ret=1ll*sum(x)*sum(x)%mod,r;
for(ll l=2;l<=x;l=r+1){
r=x/(x/l);
ret=dec(ret,1ll*dec(sumpow(r),sumpow(l-1))*sumf(x/l)%mod);
}
ans[x]=ret;
return ret;
}
int main(){
// freopen("lx.out","w",stdout);
scanf("%d%lld",&mod,&n);
init();
int ret=0;iv2=ksm(2,mod-2);iv6=ksm(6,mod-2);ll r;
for(ll l=1;l<=n;l=r+1){
r=n/(n/l);
ret=add(ret,1ll*dec(sumf(r),sumf(l-1))*sum(n/l)%mod*sum(n/l)%mod);
}
printf("%d\n",ret);
return 0;
}
想必你也注意到了，我们的杜教筛往往筛的是一些与\(\varphi,\mu\)等积性函数相关的函数，于是我们这里给出一些常见的卷积，其他的柿子大部分都能由这些导出
