杜教筛及相关知识小结

· · 算法·理论

制作不易，可以点赞收藏，以便复习。

upd on 2026/4/6： 杜教筛时间复杂度证明有误，已修改。

求通过。

切入正题

一切都来源于它。

杜教筛被用于处理一类数论函数的前缀和问题。对于数论函数 f，杜教筛可以在低于线性时间的复杂度内计算

——以上内容来自 [oi_wiki](https://oi-wiki.org/math/number-theory/du/)。 # 1. 前置知识 ## 1.1 积性函数 对于一个定义域为所有正整数的函数（数论函数），若满足 $\forall p,q\in N^*$ 且 $\gcd(p,q)=1$ 时，有 $f(pq)=f(p)f(q)$，则称此数论函数为积性函数。 特别的，若 $\forall p,q\in N^*$，都有 $f(pq)=f(p)f(q)$，则称此数论函数为**完全**积性函数。 这里举一些栗子，

比如 \varphi，\sigma，以及接下来1.3会讲的 \mu 都是积性函数。

单位函数 id(n)=n，元函数 \epsilon(n)=[n=1]，幂函数 I_k(n)=n^k（特殊的，I(n)=1）

其中 \left [ \right ] 是艾弗森括号，若 \left [ P \right ] 中 P 为真，则值为 1，否则为 0。（类似于布尔型变量）

1.2 数论分块（整除分块）

这是为了解决一类求和问题：\sum_{i=1}^{n} \lfloor \frac{n}{i} \rfloor。

这里如果我们打出表格（以 8 为例）

i 1 2 3 4 5 6 7 8
8/i 8 4 2 2 1 1 1 1

会发现有一些 \lfloor \frac{n}{i} \rfloor 会重复连续出现，具体表现为，若 l 为该块的左端点，则右端点为 \lfloor \frac{n}{ \lfloor \frac{n}{i} \rfloor } \rfloor。

证明如下

代码如下

int ans=0;
for(int l=1,r;l<=n;l=r+1){
    r=n/(n/l);
    ans+=(r-l+1)*(n/l);
}

这样我们分块的数量小于 2\sqrt{n} 种。时间复杂度为 O(\sqrt{n})

也简单证明一下

1.3 狄利克雷卷积

设 f，g 为两个数论函数，则定义他们的狄利克雷卷积为 (f*g)(n)=\sum_{d \mid n} f(d)g(\frac{n}{d})。

注意：

1.狄利克雷卷积满足交换律和结合律。

2.I*id=\sigma，\varphi*I=id。

3.特别的，对于任意数论函数 f，都有 f*\epsilon=\epsilon*f=f。

4.两个积性函数的狄利克雷卷积仍是积性函数。

1.4 莫比乌斯函数&莫比乌斯反演

先给出莫比乌斯函数的定义。

\mu(i)=\begin{cases} 1 & x=1 \\ (-1)^{k} & x= {\textstyle \prod_{i=1}^{k} p_i} \text{ 且 } p_i为互不相同的素数\\ 0 & \text{其他情况} \end{cases}

莫比乌斯函数也是一个积性函数，特别的

\sum_{d\mid n} \mu(d)=\epsilon(n)=\left [ n=1 \right ]

接下来就可以给出莫比乌斯反演公式了。

若 f 为一个算数函数（数论函数），F 为 f 的和函数，即 F(n)=\sum_{d\mid n} f(d)，则有

f(n)=\sum_{d\mid n} \mu(d)F(\frac{n}{d})

所以设 f 是算术函数，它的和函数 F(n)=\sum_{d\mid n} f(d)。如果 f 是积性函数，则 F 也是积性函数。如果 F 是极性函数，则 f 也是也是积性函数。

注意莫比乌斯反演不需要 f 是积性函数。

2.例题

2.1 P2261 余数求和

题目传送门。

题目大意

求 G(n,k)=\sum_{i=1}^{n} k\mod{i}。

思路

先给出一个有用的公式（用 a\%b 表示 a 除以 p 的余数）：

\begin{aligned} a \% p=a-p* \lfloor \frac{a}{p} \rfloor \end{aligned}

这个还是显然的。

则原式可化为：

\begin{aligned} G(n,k) &=\sum_{i=1}^{n} k\mod{i} \\ &=\sum_{i=1}^{n}k-i*\lfloor \frac{k}{i} \rfloor \\ &=\sum_{i=1}^{n} k-\sum_{i=1}^{n} i*\lfloor \frac{k}{i} \rfloor \\ &=nk-\sum_{i=1}^{n} i*\lfloor \frac{k}{i} \rfloor \end{aligned}

前一部分直接算，后一部分用整除分块。

注意 \lfloor \frac{k}{i} \rfloor 在 i 较大时会等于 0，细节看代码。

code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int main(){
    ll n,k,ans=0;
    cin>>n>>k;
    ans=n*k;
    for(ll l=1,r;l<=n;l=r+1){
        if(k/l) r=min(k/(k/l),n);    //不为0时
        else r=n;    //等于0时
        ans-=(k/l)*(l+r)*(r-l+1)/2;
    }
    cout<<ans;
    return 0;
}

2.2 P2424 约数和

题目传送门。

题目描述

定义 \sigma(n)=\sum_{d \mid n} d。

求 \sum_{i=x}^{y} \sigma(i)。

思路

还是推式子。

\begin{aligned} \sum_{i=x}^{y} \sigma(i) &= \sum_{i=1}^{y} \sigma(i)-\sum_{i=1}^{x-1} \sigma(i) \end{aligned}

记 g(n)=\sum_{i=1}^{n} \sigma(i)，则有 \sum_{i=x}^{y} \sigma(i)=g(y)-g(x-1)。

考虑 g 则有：

\begin{aligned} g(n) &=\sum_{i=1}^{n} \sigma(i) \\ &=\sum_{i=1}^{n} \sum_{d \mid i} d \\ &=\sum_{d=1}^{n} d\sum_{d\mid i}^{n} 1\\ &=\sum_{d=1}^{n} d\lfloor \frac{n}{d} \rfloor \end{aligned}

这里运用到一个很重要的思想，就是枚举一堆数的因数，和枚举因数，在这堆数中枚举倍数的效果是一样的。

code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll sum(ll n){
    ll ans=0;
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ans+=(l+r)*(r-l+1)/2*(n/l);
    }
    return ans;
}
int main(){
    ll x,y;
    cin>>x>>y;
    cout<<sum(y)-sum(x-1);
    return 0;
}

这里还有一个类似的题目，只不过求的是因数个数函数。

题目传送门 P9611。

2.3 P3455 [POI 2007] ZAP-Queries

题目描述

求 \sum_{i=1}^{a} \sum_{j=1}^{b} \left[ \gcd(i,j)=d \right]。

思路

又又又推式子。

\begin{aligned} \sum_{i=1}^{a} \sum_{j=1}^{b} \left[ \gcd(i,j)=d \right]&=\sum_{i=1}^{\lfloor \frac{a}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{b}{d} \rfloor} \left[ \gcd(i,j)=1 \right]\\ \end{aligned}

记 n=\lfloor \frac{a}{d} \rfloor，m=\lfloor \frac{b}{d} \rfloor。

运用 \epsilon(n)=\sum_{d\mid n} \mu(d)，有

\begin{aligned} \sum_{i=1}^{n} \sum_{j=1}^{m} \left[ \gcd(i,j)=1 \right] &=\sum_{i=1}^{n}\sum_{j=1}^{m} \sum_{d\mid i,d\mid j} \mu(d) \\ &=\sum_{d=1}^{\operatorname{min(n,m)}} \mu(d) \lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor \end{aligned}

只要求出 \sum_{i=1}^{n} \mu(i)，就可以做整除分块了。

code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e4+10;
ll cnt,mu[N],vis[N],pri[N],sum[N];
void init()
{
    mu[1]=1;
    for(int i=2;i<N;i++){
        if(vis[i]==0){
            mu[i]=-1;
            vis[i]=1;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt;j++){
            if(i*pri[j]>=N){
                break;
            }
            vis[i*pri[j]]=1;
            if(i%pri[j]==0){
                mu[i*pri[j]]=0;
                break;
            }
            mu[i*pri[j]]=mu[i]*mu[pri[j]];
        }
    }
    sum[0]=0;
    sum[1]=mu[1];
    for(int i=2;i<N;i++){
        sum[i]=sum[i-1]+mu[i];
    }
}  //线性筛求莫比乌斯函数
ll kkk(ll a,ll b,ll d){
    a/=d;
    b/=d;
    ll ans=0,r=0;
    if(a<b) swap(a,b);
    for(ll l=1;l<=b;l=r+1){
        r=min(a/(a/l),b/(b/l));
        ans+=(sum[r]-sum[l-1])*(long long)(a/l)*(b/l);
    }
    return ans;
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL);
    cout.tie(NULL);
    init();
    int t;
    cin>>t;
    while(t--){
        ll a,b,d;
        cin>>a>>b>>d;
        cout<<kkk(a,b,d)<<endl;
    }
    return 0;
}

这题有个类似的 P4450，还有一个扩展 P2257

简单讲一下扩展的思路：

令小于 \operatorname{min}(n,m) 的质数集为 P，且不妨令 n\le m，则原式为

\begin{aligned} \sum_{p\in P} \sum_{i=1}^{n}\sum_{j=1}^{m} \left [ \gcd(i,j)=p \right] &= \sum_{p\in P}\sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{p} \rfloor} \left[ \gcd(i,j)=1 \right] \\ &=\sum_{p\in P}\sum_{i=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{i=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{d\mid i,d\mid j} \mu(d)\\ &=\sum_{p\in P} \sum_{d=1}^{\lfloor \frac{n}{p} \rfloor} \mu(d) \lfloor \frac{n}{pd} \rfloor\lfloor \frac{m}{pd} \rfloor \end{aligned}

令 T=pd，则有

\begin{aligned} RHS &=\sum_{T=1}^{n} \lfloor \frac{n}{T} \rfloor\lfloor \frac{m}{T} \rfloor \sum_{p\mid T,p\in P} \mu(\lfloor \frac{T}{p} \rfloor) \end{aligned}

后面这一部分可以再 O(nloglogn) 的时间复杂度内预处理出来（和埃氏筛一样的时间复杂度）。然后就又可以整除分块了。

代码如下：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e7+10;
ll cnt,mu[N],vis[N],pri[N],f[N],sum[N];
void init()
{
    mu[1]=1;
    for(int i=2;i<N;i++){
        if(vis[i]==0){
            mu[i]=-1;
            vis[i]=1;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt;j++){
            if(i*pri[j]>=N){
                break;
            }
            vis[i*pri[j]]=1;
            if(i%pri[j]==0){
                mu[i*pri[j]]=0;
                break;
            }
            mu[i*pri[j]]=mu[i]*mu[pri[j]];
        }
    }
    for(int i=1;i<=cnt;i++){
        for(int j=1;j*pri[i]<N;j++){
            f[pri[i]*j]+=mu[j];
        }
    }
    sum[0]=0;
    sum[1]=f[1];
    for(int i=2;i<N;i++){
        sum[i]=sum[i-1]+f[i];
    }
}
ll kkk(ll a,ll b){
    ll ans=0,r=0;
    if(a<b) swap(a,b);
    for(ll l=1;l<=b;l=r+1){
        r=min(a/(a/l),b/(b/l));
        ans+=(sum[r]-sum[l-1])*(long long)(a/l)*(b/l);
    }
    return ans;
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL);
    cout.tie(NULL);
    init();
    int t;
    cin>>t;
    while(t--){
        ll a,b;
        cin>>a>>b;
        cout<<kkk(a,b)<<'\n';
    }
    return 0;
}

2.4 P5221 Product

题目传送门。

题目描述

求 \prod_{i=1}^{n}\prod_{j=1}^{n} \frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}。

浅浅退亿个式子：

\begin{aligned} \prod_{i=1}^{n}\prod_{j=1}^{n} \frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}&=\prod_{i=1}^{n}\prod_{j=1}^{n} \frac{ij}{\gcd(i,j)^2}\\ &=\prod_{i=1}^{n}\prod_{j=1}^{n} ij*(\prod_{i=1}^{n}\prod_{j=1}^{n} \gcd(i,j))^{(-2)} \end{aligned}

对于前半部分，有：

\begin{aligned} \prod_{i=1}^{n}\prod_{j=1}^{n}ij=\prod_{i=1}^{n}i^n(n!) &=(n!)^{2n} \end{aligned}

对于后半部分，有

\begin{aligned} \prod_{i=1}^{n}\prod_{j=1}^{n} gcd(i,j)=\prod_{d=1}^{n}d^{\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}\left[\gcd(i,j)=1 \right]} \end{aligned}

对于上面的式子就是 P2158 的故事

\sum_{i=1}^{x}\sum_{j=1}^{x} \left[gcd(i,j)=1 \right]=2\sum_{i=1}^{n}\varphi(i)-1

整合一下就 OK 了。

这里没有用莫比乌斯反演，只是帮大家感受一下连乘的推式子题。

注意这题的时空限制，所以对指数要对 mod-1 取模（欧拉定理）。

code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+10;
const ll mod=104857601;
int cnt,phi[N],pri[N];
bool vis[N];
void sieve(){
    phi[1]=1;vis[1]=1;
    for(int i=2;i<N;i++){
        if(!vis[i]){
            phi[i]=i-1;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt&&i*pri[j]<N;j++){
            vis[i*pri[j]]=1;
            if(i%pri[j]==0){
                phi[i*pri[j]]=phi[i]*pri[j];
                break; 
            }
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
    }
    for(int i=1;i<N;i++){
        phi[i]=(phi[i]+phi[i-1])%(mod-1);
    }
}
inline ll qp(ll a,ll b){
    ll sum=1;
    while(b){
        if(b&1) sum=sum*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return sum;
} 
int main(){
    //freopen(".in","r",stdin);
    //freopen(".out","w",stdout);
    ios::sync_with_stdio(false);
    cin.tie(NULL);
    cout.tie(NULL);
    sieve();
    ll n;
    cin>>n;
    ll ans1=1,ans2=1;
    for(int i=1;i<=n;i++){
        ans1=ans1*i%mod;
        ans2=ans2*qp(i,phi[n/i]*2-1)%mod;
    }
    cout<<(qp(ans1,2*n)*qp(ans2,(mod-2)*2)%mod);
    return 0;
}

P7486 「Stoi2031」彩虹

题目传送门。

非常毒瘤的一道题。

题目描述

求 \prod_{i=l}^{r}\prod_{j=l}^{r} \operatorname{lcm}(i,j)^{\operatorname{lcm}(i,j)}。

思路

准备好纸笔，一起让暴风雨来得更猛烈些吧。

令 S(n,m)=\prod_{i=1}^{n} \prod_{j=1}^{n} \operatorname{lcm}(i,j)^{\operatorname{lcm}(i,j)}，且我们不妨令 n\le m。

所以原式 LHS=\frac{S(r,r)S(l,l)}{S(l,r)^2}。

接下来就是比较毒瘤的过程，这有一篇视频很好。

还有一篇好不容易才找到的正确题解。

这里直接给出 code（注意本题卡常）。

code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+10;
const ll mod=32465177;
ll t,n;
template < typename T >
void read(T &x){
    x=0;int f=1;char ch=getchar();
    while(!isdigit(ch)){
        if(ch=='-') f*=-1;
        ch=getchar();
    }
    while(isdigit(ch)){
        x=(x<<1)+(x<<3)+(ch&15);
        ch=getchar(); 
    }
    x*=f;
}
void write(int x)
{
    if(x<0)
        putchar('-'),x=-x;
    if(x>9)
        write(x/10);
    putchar(x%10+'0');
    return;
}
inline ll qp(ll a,ll b){
    ll sum=1;
    while(b){
        if(b&1) sum=sum*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return sum;
}
ll cnt,mu[N],vis[N],pri[N],f[N],finv[N],fprod[N],sum[N],x1[N],x2[N],g[N],ginv[N],sumx1[N];
inline void pre(){
    mu[1]=1;
    vis[1]=1;
    for(int i=2;i<=n;i++){
        if(!vis[i]){
            mu[i]=-1;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=n;j++){
            vis[i*pri[j]]=1;
            if(i%pri[j]==0){
                mu[i*pri[j]]=0;
                break;
            }
            mu[i*pri[j]]=-mu[i];
        }
    }

    fprod[0]=1;
    for(int i=1;i<=n;i++){
        f[i]=qp(i,i);
        finv[i]=qp(f[i],mod-2);
        fprod[i]=fprod[i-1]*f[i]%mod;
        sum[i]=1ll*i*(i+1)/2;
        x2[i]=1;
    } 
    for(ll i=1;i<=n;i++){
        for(ll j=1;j*i<=n;j++){
            x1[i*j]=((x1[i*j]+i*mu[i])%(mod-1)+mod-1)%(mod-1);
            if(mu[i]>0){
                x2[i*j]=x2[i*j]*f[i]%mod;
            }
            else if(mu[i]<0){
                x2[i*j]=x2[i*j]*finv[i]%mod;
            }
        }
    }
    g[0]=ginv[0]=1;
    for(int i=1;i<=n;i++){
        sumx1[i]=(sumx1[i-1]+i*x1[i])%(mod-1);
        g[i]=g[i-1]*qp(qp(i,x1[i])*x2[i]%mod,i)%mod;
        ginv[i]=qp(g[i],mod-2);
    }
}
inline ll yy(ll n,ll m){
    return qp(fprod[n],sum[m]%(mod-1))*qp(fprod[m],sum[n]%(mod-1))%mod;
}
inline ll s(ll n,ll m){
    ll ans=1;
    if(n>=m) swap(n,m);
    for(ll l=1,r;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans=ans*(qp(yy(n/l,m/l),((sumx1[r]-sumx1[l-1])%(mod-1)+mod-1)%(mod-1))*qp(g[r]*ginv[l-1]%mod,(__int128)sum[n/l]*sum[m/l]%(mod-1))%mod)%mod;
    }
    return ans;
}
inline ll solve(ll l,ll r){
    ll ans1=s(r,r);
    ll ans2=s(l-1,l-1);
    ll ans3=s(l-1,r);
    return (ans1*ans2%mod)*qp(ans3,2*(mod-2))%mod;
}
int main(){
    //freopen(".in","r",stdin);
    //freopen(".out","w",stdout);
    ios::sync_with_stdio(false);
    cin.tie(NULL);
    cout.tie(NULL);
    read(t);read(n);
    pre(); 
    while(t--){
        int l,r;
        read(l);read(r);
        write(solve(l,r));
        putchar('\n');
    }
    return 0;
}

还有一道毒瘤 P5518，大家可以踊跃挑战。

3. 杜教筛介绍

3.1 推导

接下来就是 无聊 有趣的推式子。

问题的原型是这样的，设 f 为一个积性函数，在低于线性时间复杂度内计算 S(n)=\sum_{i=1}^{n} f(i)。

为了使其低于线性时间复杂度，尝试构造形如 \sum_{i=1}^{n} \lfloor \frac{n}{i} \rfloor 的整除分块。

构造两个积性函数 h 和 g，满足狄利克雷卷积（以下简称卷积） h=g*f。

由卷积的定义得 h(i)=\sum_{d\mid i} g(d)f(\frac{i}{d})，对 h(i) 求和，有

\sum_{i=1}^{n} h(i)&=\sum_{i=1}^{n} \sum_{d\mid i}g(d)f(\frac{i}{d})\\ &=\sum_{d=1}^{n}g(d)\sum_{d|i}^{n}f(\frac{i}{d})\text{(“杜教筛变换”)}\\ &=\sum_{d=1}^{n}g(d)\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} f(i)\\ &=\sum_{d=1}^{n}g(d)S(\lfloor \frac{n}{d} \rfloor) \end{aligned}

最后一步得到 \sum_{i=1}^{n} h(i)=\sum_{d=1}^{n}f(d)S(\lfloor \frac{n}{d} \rfloor)。则

\begin{aligned} \sum_{i=1}^{n} h(i)&=g(1)S(n)+\sum_{d=2}^{n}g(d)S(\lfloor \frac{n}{d} \rfloor)\\ g(1)S(n)&=\sum_{i=1}^{n}h(i)-\sum_{i=2}^{n}g(i)S(\lfloor \frac{n}{i} \rfloor) \end{aligned}

这就是杜教筛公式。

3.2 实现

我们可以用递归来实现一个朴素版本（伪代码）：

typedef long long ll;
ll getsum(ll n){
    ll ans=sumh(n);   //h 的前缀和
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        ans-=(sumg(r)-sumg(l-1))*(getsum(n/l));   //整除分块
    }
    return ans;
}

我们会发现这版代码时间复杂度为 O(n^{\frac{3}{4}})。

但我们发现其实还可以优化，比如用线性筛预处理一部分，还有在递归时加入记忆化。

先以此题为例，有特殊性质 n=\sum_{d|n} \varphi(d)，[n=1]=\sum_{d|n} \mu(d)。

对于 \varphi，有 id=\varphi*I，所以有

\begin{aligned} S(n)&=\sum_{i=1}^{n}i-\sum_{i=2}^{n}S(\lfloor \frac{n}{i} \rfloor)\\ &=\frac{n(n+1)}{2}-\sum_{i=2}^{n}S(\lfloor \frac{n}{i} \rfloor) \end{aligned}

对于 \mu，有 \epsilon=\mu*I，所以有

\begin{aligned} S(n)&=\sum_{i=1}^{n}\epsilon(i)-\sum_{i=2}^{n}S(\lfloor \frac{n}{i} \rfloor)\\ &=1-\sum_{i=2}^{n}S(\lfloor \frac{n}{i} \rfloor) \end{aligned}

代码用 unordered_map 实现如下

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e6+10;
ll mu[N],phi[N],vis[N],pri[N];
unordered_map<ll,ll> sumphi;
unordered_map<ll,ll> summu;
void init()
{
    vis[0]=vis[1]=1;
    mu[1]=phi[1]=1;
    int cnt=0;
    for(int i=2;i<N;i++){
        if(vis[i]==0){
            phi[i]=i-1;
            mu[i]=-1;
            vis[i]=1;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt;j++){
            if(i*pri[j]>=N){
                break;
            }
            vis[i*pri[j]]=1;
            if(i%pri[j]==0){
                phi[i*pri[j]]=phi[i]*pri[j];
                mu[i*pri[j]]=0;
                break;
            }
            phi[i*pri[j]]=phi[i]*phi[pri[j]];
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<N;i++){
        phi[i]+=phi[i-1];
        mu[i]+=mu[i-1];
    }
}
ll getsphi(ll x){
    if(x<N){
        return phi[x];
    }
    else if(sumphi.count(x)){
        return sumphi[x];
    }
    ll ans=x*(x+1)/2;
    for(ll l=2,r;l<=x;l=r+1){
        r=x/(x/l);
        ans-=(r-l+1)*getsphi(x/l);
    }
    sumphi[x]=ans;
    return ans;
}//杜教筛
ll getsmu(ll x){
    if(x<N){
        return mu[x];
    }
    else if(summu.count(x)){
        return summu[x];
    }
    ll ans=1;
    for(ll l=2,r;l<=x;l=r+1){
        r=x/(x/l);
        ans-=(r-l+1)*getsmu(x/l);
    }
    summu[x]=ans;
    return ans;
}//杜教筛
int main(){
    init();
    int t;
    cin>>t;
    while(t--){
        ll n;
        cin>>n;
        cout<<getsphi(n)<<" "<<getsmu(n)<<endl;
    }
    return 0;
}

不过这版代码有个小问题，就是 unordered_map 会有一点常数。所以我们可以直接建一个长度为 2\sqrt{n} 的数组。

改进如下：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e6+10,M=1<<17;
ll n,mu[N],phi[N],vis[N],pri[N],sumphi[M],summu[M];
ll ind(ll x){
    ll t=sqrt(n);
    if(x<=t) return x;
    if(x>t) return t+n/x;
}
void init()
{
    vis[0]=vis[1]=1;
    mu[1]=phi[1]=1;
    int cnt=0;
    for(int i=2;i<N;i++){
        if(vis[i]==0){
            phi[i]=i-1;
            mu[i]=-1;
            vis[i]=1;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt;j++){
            if(i*pri[j]>=N){
                break;
            }
            vis[i*pri[j]]=1;
            if(i%pri[j]==0){
                phi[i*pri[j]]=phi[i]*pri[j];
                mu[i*pri[j]]=0;
                break;
            }
            phi[i*pri[j]]=phi[i]*phi[pri[j]];
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<N;i++){
        phi[i]+=phi[i-1];
        mu[i]+=mu[i-1];
    }
}
ll getsphi(ll x){
    if(x<N){
        return phi[x];
    }
    else if(sumphi[ind(x)]!=-1){
        return sumphi[ind(x)];
    }
    \\杜教筛略去
}
ll getsmu(ll x){
    if(x<N){
        return mu[x];
    }
    else if(summu[ind(x)]!=-1){
        return summu[ind(x)];
    }
    \\杜教筛略去
}
int main(){
    init();
    int t;
    cin>>t;
    while(t--){
        memset(summu,-1,sizeof(summu));
        memset(sumphi,-1,sizeof(sumphi));    //记得清空
        cin>>n;
        cout<<getsphi(n)<<" "<<getsmu(n)<<endl;
    }
    return 0;
}

这样时间复杂度是 O(n^{\frac{2}{3}})。

3.3 时间复杂度证明（过于比较复杂，可略过）

回到公式 g(1)S(n)=\sum_{i=1}^{n}h(i)-\sum_{i=2}^{n}g(i)S(\lfloor \frac{n}{i} \rfloor)。

如果函数 g 和 h 选的好的话，可以在 O(1) 时间复杂度内求出前缀和。

设计算 \sum_{i=1}^{n} (f * g)(i) 和 \sum_{i=1}^{n} g(i) 的时间复杂度均为 O(1)。设计算 S(n) 的时间复杂度为 T(n)，则： ::::error[错误证法]

展开第一层，有

\begin{aligned} T(n)&=O(\sqrt{n})+O \left[ \sum_{i=2}^{\sqrt{n}} T(\frac{n}{i}) \right ] \end{aligned}

在展开一层，有

\begin{aligned} T(\frac{n}{i})&=O(\sqrt{\frac{n}{i}})+O\left[ \sum_{k=2}^{\sqrt{\frac{n}{i}}} T(\frac{\frac{n}{i}}{k}) \right]\\ &=O(\sqrt{\frac{n}{i}}) \end{aligned}

第二项是高阶小量，可以直接省略。

回带可以得到

\begin{aligned} T(n)&=O(\sqrt{n})+O \left[ \sum_{i=2}^{\sqrt{n}} T(\frac{n}{i}) \right ]\\ &=O(\sqrt{n})+O \left[ \sum_{i=2}^{\sqrt{n}} (\sqrt{\frac{n}{i}}) \right ] \end{aligned}

可以看出最后一项最大，只看它就行了。

于是就有

\begin{aligned} T(n)&=O \left[ \sum_{i=2}^{\sqrt{n}} (\sqrt{\frac{n}{i}}) \right ] \approx O(\int_{0}^{\sqrt n} \sqrt{\frac{n}{x}}\mathop{} \! \mathrm{d} x)\\&=O(n^{\frac{3}{4}}) \end{aligned}

:::warning[错误原因（引自 OI Wiki）] 在“视为高阶无穷小量”时，这里的论断过于突然，我们仔细分析一下，将 T\left(\left\lfloor \dfrac{n}{i} \right\rfloor\right) 代入 T(n) 的式子里，有：

\begin{aligned} T(n) & = \Theta\left(\sqrt{n}\right) + O\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} \sqrt{\frac{n}{i}}\right) + O\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} \sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right)\right) \\ & = O\left(\sqrt{n} + \int_{0}^{\sqrt{n}} \sqrt{\frac{n}{x}} \mathrm{d}x\right) + O\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} \sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right)\right) \\ & = O\left(n^{3/4}\right) + O\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} \sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right)\right) \end{aligned}

我们考虑 \displaystyle\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} \sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right) 这部分，不难发现：

\begin{aligned} \sum_{i=2}^{\lfloor\sqrt{n}\rfloor} \sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right) & = \Omega\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} T\left(\left\lfloor\frac{n}{i} \cdot \left\lfloor\sqrt{\frac{n}{i}}\right\rfloor^{-1}\right\rfloor\right)\right) \\ & = \Omega\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} T\left(\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor\right)\right) \end{aligned}

由于没有引入记忆化，因此上式中的 T\left(\left\lfloor\sqrt{\dfrac{n}{i}}\right\rfloor\right) 仍然是 \Omega\left(\left(\dfrac{n}{i}\right)^{1/4}\right) 的，进而所谓的「高阶无穷小」部分是不可以舍去的。

实际上杜教筛的亚线性时间复杂度是由记忆化保证的。只有使用了记忆化之后才能保证不会出现那个多重求和的项。

:::
::::success[正确证明（引用自 OI Wiki）]
令 R(n) = \left\{\left\lfloor \dfrac{n}{k} \right\rfloor : k = 2, 3, \ldots, n\right\}。利用数论分块的性质可知，对任意的 m \in R(n)，都有 R(m) \subseteq R(n)（证明见此处）。也就是说，使用记忆化之后，只需要对所有 k \in R(n) 计算一次 S(k) 就可以得到 R(n) 的值。而这些点的数目 $
\begin{aligned} T(n) & = \sum_{k \in R(n)} T(k) \\ & = \Theta(\sqrt{n}) + \sum_{k=1}^{\lfloor\sqrt{n}\rfloor} O(\sqrt{k}) + \sum_{k=2}^{\lfloor\sqrt{n}\rfloor} O\left(\sqrt{\frac{n}{k}}\right) \\ & = O\left(\int_{0}^{\sqrt{n}} \left(\sqrt{x} + \sqrt{\frac{n}{x}}\right) \mathrm{d}x\right) \\ & = O\left(n^{3/4}\right). \end{aligned}

若我们可以预处理出一部分 S(k)，其中 k = 1, 2, \ldots, m，m \geq \lfloor\sqrt{n}\rfloor。设预处理的时间复杂度为 T_0(m)，则此时的 T(n) 为：

\begin{aligned} T(n) & = T_0(m) + \sum_{k \in R(n); k > m} T(k) \\ & = T_0(m) + \sum_{k=1}^{\lfloor n/m \rfloor} O\left(\sqrt{\frac{n}{k}}\right) \\ & = O\left(T_0(m) + \int_{0}^{n/m} \sqrt{\frac{n}{x}} \mathrm{d}x\right) \\ & = O\left(T_0(m) + \frac{n}{\sqrt{m}}\right). \end{aligned}
若 T_0(m) = O(m)（如线性筛），由均值不等式可知：当 m = \Theta\left(n^{2/3}\right) 时，T(n) 取得最小值 O\left(n^{2/3}\right)。

所以该算法的关键就是要找到合适的 h 和 g，使得它们可以快速（大多是情况为 O(1)）求出前缀和。

比如求解的函数由形如 \sum_{d\mid n} f(d) 的性质，就容易找到 g 和 h。

3.4 应用

建议大家自己推一遍式子，感受其中的“杜教筛变换”。

然后大家可以完成 P3768。

有些莫比乌斯反演的推式子题可能会要求在低于线性时间复杂度内求出莫比乌斯函数前缀和，这时就可以用杜教筛。

3.4 例题

P3768：

\begin{aligned} \sum_{i=1}^{n}\sum_{j=1}^{n} ij\gcd(i,j) &=\sum_{i=1}^{n}\sum_{j=1}^{n} ij\sum_{d\mid \gcd(i,j)} \varphi(d)\\ &=\sum_{d=1}^{n}\varphi(d) d^2 \sum_{i=1}^{\lfloor \frac{n}{d}\rfloor} i\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor} j\\ \end{aligned}

这里令 sum(x)=\sum_{i=1}^{x}i，则有

\begin{aligned} RHS&=\sum_{d=1}^{n} \varphi(d)d^2 sum(\lfloor \frac{n}{d} \rfloor)\\ \text{令} f(d)&=\varphi(d) d^2\\ g(d)&=d^2\\ S(d)&=\sum_{i=1}^{d} f(i)\\ \text{则} S(n)&=\sum_{i=1}^{n} (f*g)(i)-\sum_{i=2}^{n} g(i)S(\lfloor \frac{n}{i} \rfloor)\\ &=\sum_{i=1}^{n} i^3-\sum_{i=2}^{n} i^2S(\lfloor \frac{n}{i} \rfloor)\\ &=sum(n)^2-\sum_{i=2}^{n} i^2S(\lfloor \frac{n}{i} \rfloor) \end{aligned} 代码如下： ```cpp line-numbers #include<bits/stdc++.h> using namespace std; typedef long long ll; ll p,n,inv6,inv2; ll fp(ll k,ll x) { ll sum=1; while(x) { if(x&1) { sum=sum*k%p; } k=k*k%p; x>>=1; } return sum; }//快速幂 const int N=5e6+10; ll phi[N+1],vis[N+1],pri[N+1],sum[N+1]; void init() { phi[1]=1; int cnt=0; for(ll i=2;i<N;i++){ if(vis[i]==0){ phi[i]=i-1; vis[i]=1; pri[++cnt]=i; } for(ll j=1;j<=cnt;j++){ if(i*pri[j]>=N){ break; } vis[i*pri[j]]=1; if(i%pri[j]==0){ phi[i*pri[j]]=(phi[i]*pri[j])%p; break; } phi[i*pri[j]]=(phi[i]*phi[pri[j]])%p; } } sum[1]=1*1*phi[1]; for(ll i=2;i<N;i++){ sum[i]=(sum[i-1]+((i*i)%p*phi[i])%p)%p; } }//预处理 ll sumy(ll x){x%=p;return x*(x+1)%p*inv2%p;} //一次方求和 ll sumq(ll x){ x%=p; return x*(x+1)%p*(2*x+1)%p*inv6%p; } //二次方求和 unordered_map<ll,ll> sumsq; ll ds(ll x){ if(x<N){ return sum[x]; } if(sumsq.count(x)){ return sumsq[x]; } else{ ll ans=0; ll tmp=sumy(x); ans=tmp*tmp%p; for(ll i=2,j;i<=x;i=j+1){ j=x/(x/i); ll stmp=(sumq(j)-sumq(i-1)+p)%p; ans=(ans-stmp*ds(x/i)%p+p)%p; } sumsq[x]=(ans+p)%p; return (ans+p)%p; } } //杜教筛 int main(){ cin>>p>>n; inv2=fp(2,p-2); inv6=fp(6,p-2); //求逆元 init(); ll ans=0; for(ll i=1,j;i<=n;i=j+1){ j=n/(n/i); ll tmp=sumy(n/i); tmp=tmp*tmp%p; ll cnt=(ds(j)-ds(i-1)+p)%p; ans+=tmp*cnt%p; ans%=p; } //数论分块 cout<<(ans+p)%p; return 0; } ``` # 4. 文献来源 1. 杜教筛的介绍：Oi-wiki 2. 前置知识、杜教筛的推导和时间复杂度证明：《算法竞赛》P454~P458，罗勇军等著