杜教筛

摘要

莫比乌斯反演推完式子后,一般考虑线性筛和数论分块。当要求在低于线性时间的求解积性函数前缀和的问题的时候就会用到杜教筛。

符号说明

简单解释一下本文可能用到的记号的含义。

  1. 狄利克雷卷积符号
  2. 表示数论函数 的前缀和。即

杜教筛

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

对于数论函数 ,计算

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

首先有一个小小的结论:对于任意一个数论函数 ,必满足

写成卷积的形式就是

略证:

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

通过上面的结论,我们可以得到一个递推式

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

【例 1】莫比乌斯函数前缀和

我们选取 函数作为 ,那么

于是我们就可以直接递归求解 ,复杂度 。我们还可以预处理前 项,使递归部分复杂度降至 ,总复杂度 .

【例 2】欧拉函数前缀和

同样的,利用 得到

【例 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;
}

【例 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 论文集


  转载请注明: Sshwy's Blog 杜教筛

 上一篇
[SDOI2009] 晨跑 [SDOI2009] 晨跑
题意在加权有向图中求一路径集合 满足 起点为 ,终点为 ,即 . 集合内的路径除了起点终
2019.01.12
下一篇 
莫比乌斯反演小结 莫比乌斯反演小结
摘要 复习莫比乌斯反演~ 前言OI 的数论涉及求和的部分,一般采用 暴力计算;但是当上界过大的时候就需要考虑数论求和法。 常用的技巧有前缀和、差分、组合计数、等差数列求和、矩阵快速幂等,这些技巧都建立在数学原理的基础上 如果
2019.01.10
  目录