杜教筛

摘要

复习杜教筛~

前言

莫比乌斯反演推完式子后,一般考虑线性筛和数论分块

当要求在低于线性时间的求解积性函数前缀和的问题的时候就会用到杜教筛

杜教筛

杜教筛被用来处理数论函数的前缀和问题。对于求解一个前缀和,杜教筛可以在低于线性时间的复杂度内求解

对于数论函数 ,要求我们计算 .

我们想办法构造一个 关于 的递推式

对于任意一个数论函数 ,必满足

略证:

就是对所有 的做贡献,因此变换枚举顺序,枚举 (分别对应新的

那么可以得到递推式

那么假如我们可以快速对 求和,并用数论分块求解 就可以在较短时间内求得 .

【例 1】

.

我们选取 函数作为 ,那么

于是我们就可以分块递归求解 ,复杂度 .

我们还可以预处理前 项,使递归部分复杂度降至 ,总复杂度 .

【例1例2】模板

LuoguP4213 【模板】杜教筛(Sum)

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e6+6,lim=5e6;
int T,n;

int pn[N],cnt;
bool vis[N];
ll phi[N];
int mu[N];
void sieve(int n){
vis[0]=vis[1]=1,mu[1]=1,phi[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i])pn[++cnt]=i,phi[i]=i-1,mu[i]=-1;
for(int j=1;j<=cnt;j++){const int x=i*pn[j];
if(x>n)break;
vis[x]=1;
if(i%pn[j])phi[x]=phi[i]*phi[pn[j]], mu[x]=mu[i]*mu[pn[j]];
else {phi[x]=phi[i]*pn[j], mu[x]=0; break;}
}
}
for(int i=2;i<=n;i++)phi[i]+=phi[i-1],mu[i]+=mu[i-1];
}
unordered_map<int,ll> mphi,mmu;
ll calc_phi(int n){
if(n<=lim)return phi[n];
if(mphi.find(n)!=mphi.end())return mphi[n];
ll res=(ll)n*(n+1)/2;
int p=2,q;//一定要开int不然卡常
while(p<=n){
q=n/(n/p);
res-=calc_phi(n/p)*(q-p+1);
p=q+1;
}
return mphi[n]=res;
}
int calc_mu(int n){
if(n<=lim)return mu[n];
if(mmu.find(n)!=mmu.end())return mmu[n];
int res=1,p=2,q;//p,q 一定开int不然卡常
while(p<=n){
q=n/(n/p);
res-=calc_mu(n/p)*(q-p+1);
p=q+1;
}
return mmu[n]=res;
}
int main(){
sieve(lim);
scanf("%d",&T);
while(T--){
scanf("%d",&n);
printf("%lld %d\n",calc_phi(n),calc_mu(n));
}
return 0;
}

【例 2】

.

同样的,

【例 3】带函数系数的

[LuoguP3768] 简单的数学题,求

是质数

利用 反演化为(见莫比乌斯反演 - 例 5

做数论分块, 的前缀和用杜教筛处理:

需要构造积性函数 ,使得 能快速求和

单纯的 的前缀和可以用 的杜教筛处理,但是这里的 多了一个 ,那么我们就卷一个 上去,让它变成常数:

化一下卷积

再化一下

非常友好的式子啊,分块求解即可

#include<cstdio>
#include<cmath>
#include<map>
#define int long long
using namespace std;
const signed N=5e6,NP=5e6,SZ=N;
int n,P,inv2,inv6,s[N];
signed phi[N],p[NP],cnt,pn;
bool bp[N];
map<int,int> s_map;
int ksm(int a,int m){// 求逆元用
int res=1;
while(m){
if(m&1)res=res*a%P;
a=a*a%P,m>>=1;
}
return res;
}
void prime_work(signed k){// 线性筛 phi,s
bp[0]=bp[1]=1,phi[1]=1;
for(signed i=2;i<=k;i++){
if(!bp[i])p[++cnt]=i,phi[i]=i-1;
for(signed j=1;j<=cnt&&i*p[j]<=k;j++){
bp[i*p[j]]=1;
if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
else phi[i*p[j]]=phi[i]*phi[p[j]];
}
}
for(signed i=1;i<=k;i++)s[i]=(1ll*i*i%P*phi[i]%P+s[i-1])%P;
}
int s3(int k){return k%=P,(k*(k+1)/2)%P*((k*(k+1)/2)%P)%P;}// 立方和
int s2(int k){return k%=P,k*(k+1)%P*(k*2+1)%P*inv6%P;}// 平方和
int calc(int k){// 计算 S(k)
if(k<=pn)return s[k];
if(s_map[k])return s_map[k];// 对于超过 pn 的用 map 离散存储
int res=s3(k),pre=1,cur;
for(int i=2,j;i<=k;i=j+1)
j=k/(k/i),cur=s2(j),res=(res-calc(k/i)*(cur-pre)%P)%P,pre=cur;
return s_map[k]=(res+P)%P;
}
int solve(){
int res=0,pre=0,cur;
for(int i=1,j;i<=n;i=j+1)
j=n/(n/i),cur=calc(j),res=(res+(s3(n/i)*(cur-pre))%P)%P,pre=cur;
return (res+P)%P;
}
signed main(){
scanf("%lld%lld",&P,&n);
inv2=ksm(2,P-2),inv6=ksm(6,P-2);
pn=(int)pow(n,0.666667);//n^(2/3)
prime_work(pn);
printf("%lld",solve());
return 0;
}// 不要为了省什么内存把数组开小...... 卡了好几次 80

一般化小结

也就是说,对于数论函数 的前缀和,我们需要找到另外一个数论函数 使得 的前缀和可以快速计算,然后套一个数论分块求解 .

参考文献

http://blog.leanote.com/post/richard_love_oi/%E8%8E%AB%E6%AF%94%E4%B9%8C%E6%96%AF%E5%8F%8D%E6%BC%94 , Richard’s Blog

https://tqyaaaaang.com/sieve-of-dujiao/ , tqyaaaaang’s Blog

任之洲 , 2016 , 积性函数求和的几种方法 , 国家集训队 2016 论文集