题解：P4213 【模板】杜教筛

cogimyun · 2025-10-26 19:20:39 · 题解

前置知识

积性函数

顾名思义，积性函数是一类满足 f(ab)=f(a)\times f(b) 的函数，当然 f(ab)=f(a)\times f(b) 是有成立条件的，它的成立条件是 \gcd(a,b)=1。

线性筛

可以用 O(n) 的时间复杂度筛出积性函数 f(x),x\in[1,n] 的值。~~既然都来学习杜教筛了，一定会线性筛了吧~~。代码如下：

const int N=1e7;
bool pri[N];
int prime[N],mul[N],phi[N],cnt[N],idx;
inline void sieve(int x){
    mul[1]=1;
    phi[1]=1;
    memset(pri,true,sizeof pri);
    for(int i=2;i<=x;i++){
        if(pri[i]) prime[++idx]=i,mul[i]=-1,phi[i]=i-1,cnt[i]=1;
        for(int j=1;j<=idx&&prime[j]*i<=x;j++){
            pri[i*prime[j]]=false;
            if(i%prime[j]) mul[i*prime[j]]=-mul[i],phi[i*prime[j]]=phi[i]*(prime[j]-1),cnt[i*prime[j]]=cnt[i]+1;
            else{mul[i*prime[j]]=0,phi[i*prime[j]]=phi[i]*prime[j],cnt[i*prime[j]]=cnt[i];break;}
        }
    }
}

狄利克雷卷积

存在算术函数 f 和 g，其狄利克雷卷积记为 f*g，满足：

(f*g)(n)=\sum_{d|n}f(d)\times g(\frac{n}{d})

或

(f*g)(n)=\sum_{ab=n}f(a)\times g(b)

两种情况本质上是一致的，只是两种不同的写法。

这里给出几对常见的 f,g：

\mu*I=\epsilon$ 源于 $\sum\limits_{d|n} \mu(d)=[n=1]
\varphi*I=id$ 源于 $\sum\limits_{d|n}\varphi(d)=n
\mu*id=\varphi$ 源于 $\varphi(n)=\sum\limits_{i=1}^n [gcd(i,n)=1]=\sum\limits_{d|n}\frac{n}{d}\mu(d)
I*id=\sigma

其中 I(x)=1,\epsilon(x)=[x=1],id(x)=x。

整除分块

对于 \lfloor \frac{n}{l}\rfloor 明显存在一个区间 [l,r] 使得 \forall i\in[l,r],\lfloor \frac{n}{i}\rfloor=\lfloor \frac{n}{l}\rfloor，若给定 l，r 的最大值是 \lfloor\frac{n}{\lfloor\frac{n}{l}\rfloor}\rfloor，如果将 \lfloor \frac{n}{i}\rfloor 值相同的 i 分为一个块，块的数量小于 2\sqrt{n}，可以大大优化形似 \sum\limits_{i=1}^n f(i)\times g(\lfloor\frac{n}{i}\rfloor) 求和公式计算的时间复杂度，详见例题。

杜教筛

杜教筛是一种用于求数论函数前缀和的非线性算法，可以在 O(n^{\frac{2}{3}}) 的时间复杂度内求出 S(n)=\sum\limits_{i=1}^{n} f(i) 的值。

为了快速求 f(x) 的前缀和，我们可以为它找到一个函数 g(x) 使得 h(x)=(f*g)(x) 与 g(x) 的前缀和是好求的，这时就可以用狄利克雷卷积对式子进行一些变换了：

\begin{aligned}\sum\limits_{i=1}^{n}h(i)&=\sum\limits_{i=1}^{n}\sum_{d|i}g(d)f(\frac{n}{d})\\&=\sum\limits_{d=1}^{n}g(d)\sum\limits_{d|i}^nf(\frac{i}{d})\\&=\sum\limits_{d=1}^{n}g(d)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)\\&=\sum\limits_{d=1}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)\end{aligned}

其中记 S(x)=\sum\limits_{i=1}^{x}f(i)

于是我们可以得到：

\sum\limits_{i=1}^{n}h(i)=\sum\limits_{d=1}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)

\sum\limits_{i=1}^{n}h(i)=g(1)S(n)+\sum\limits_{d=2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)

g(1)S(n)=\sum\limits_{i=1}^{n}h(i)-\sum\limits_{d=2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)

根据这个式子我们可以写出如下伪代码：

long long get_fsum(int x){
    long long res=sumh(x);//h的前缀和
    for(int i=2;i<=x;i=r+1){
        int r=x/(x/i);//整除分块
        res-=(sumg(r)-sumg(i-1))*get_fsum(x/i);
    }
    return res;
}

我们假设 g,h 求前缀和的时间复杂度是 O(1) 的，那么这样的杜教筛时间复杂度是 O(n^{\frac{3}{4}}) 的，证明如下：

设在 x 时我们计算 S(x) 附加计算的 S(y),y=\lfloor\frac{x}{i}\rfloor\ i=2,3,...,n 其中 y 构成的集合为 T(x)=\{y|y=\lfloor\frac{x}{i}\rfloor,i=2,3,...,x\}，根据整除分块的结论，我们可以得到 \forall z\in T(x) 必然有 T(z)\subseteq T(x)，我们可以记忆化记录 \forall z\in T(n) 的 S(z) 即可，而这样的 z 个数仅有 \sqrt{n} 个，计算 S(z) 的时间复杂度是基于整除分块的 O(\sqrt{z}）的，记 T'(n)=T(n)\cup\{n\}，时间复杂度即为：

O(\sum_{z\in T'(n)}\sqrt{z})

我们可以将 T'(n) 拆为 A=\{x|x\in T'(n),x>\sqrt{n}\} 与 B=\{x|x\in T'(n),x\le\sqrt{n}\}

对于 A 部分：

\sum_{z\in A}\sqrt{z}\approx\sum\limits_{i=1}^{\sqrt{n}}\sqrt{\frac{n}{i}}=\sqrt{n}\sum\limits_{i=1}^{\sqrt{n}}\sqrt{\frac{1}{i}}

已知 \sum\limits_{i=1}^{k}\sqrt{\frac{1}{i}}\sim2\sqrt{k}，所以：

\sum_{z\in A}\sqrt{z}\approx\sqrt{n}\cdot2\sqrt{\sqrt{n}}=2n^{\frac{3}{4}}

对于 B 部分

\sum_{z\in B}\sqrt{z}=\int_{1}^{\sqrt{n}}\sqrt{z}\ dz=\left.\frac{2}{3}z^{\frac{3}{2}}\right|_1^{\sqrt{n}}

分别带入上下界得：

\sum_{z\in B}\sqrt{z}=\left.\frac{2}{3}z^{\frac{3}{2}}\right|_1^{\sqrt{n}}=\frac{2}{3}(z^{\frac{3}{4}}-1)

于是总时间复杂度为：

O(\sum_{z\in T'(n)}\sqrt{z})=O(2n^{\frac{3}{4}}+\frac{2}{3}(z^{\frac{3}{4}}-1))=O(n^{\frac{3}{4}})

我们还可以继续优化杜教筛，我们先用线性筛筛出 i\le n^{\frac{2}{3}} 的 S(i)，此时线性筛与杜教筛的时间复杂度均为 O(n^{\frac{2}{3}})，达到最优，证明过程与上面类似。

实际应用

接下来我们回归题目，对于求 \mu 的前缀和就可以用前置知识中给出的 \mu*I=\epsilon，这非常容易实现。同理求 \varphi 的前缀和就可以用 \varphi*I=id，依旧容易实现，接下来给出代码。

CODE

#include<bits/stdc++.h>
using namespace std;
#define int long long 
int t,pri[1000005],prime[1000005],idx,mul[1000005],phi[1000005];
unordered_map<int,int>mull,phii,book_mul,book_phi;
const int len=1000000;
void init(){
    mul[1]=1;
    phi[1]=1;
    memset(pri,true,sizeof pri);
    for(int i=2;i<=len;i++){
        if(pri[i]) {prime[++idx]=i;mul[i]=-1;phi[i]=i-1;}
        for(int j=1;j<=idx&&i*prime[j]<=len;j++){
            pri[i*prime[j]]=false;
            if(i%prime[j]){mul[i*prime[j]]=-mul[i];phi[i*prime[j]]=phi[i]*(prime[j]-1);}
            else{mul[i*prime[j]]=0;phi[i*prime[j]]=phi[i]*prime[j];break;}
        }
    }
    for(int i=1;i<=len;i++)mul[i]+=mul[i-1],phi[i]+=phi[i-1];
}
int get_mul(int x){
    if(x<=len)return mul[x];
    else if(book_mul[x]) return mull[x];
    int r,res=1;
    for(int i=2;i<=x;i=r+1){
        r=x/(x/i);
        res-=(r-i+1)*get_mul(x/i);
    }
    book_mul[x]=1,mull[x]=res;
    return res;
}
int get_phi(int x){
    if(x<=len)return phi[x];
    else if(book_phi[x]) return phii[x];
    int r,res=(x+1)*x/2;
    for(int i=2;i<=x;i=r+1){
        r=x/(x/i);
        res-=(r-i+1)*get_phi(x/i);
    }
    book_phi[x]=1,phii[x]=res;
    return res;
}
signed main(){
    init();
    cin>>t;
    while(t--){
        int n;
        cin>>n;
        cout<<get_phi(n)<<" "<<get_mul(n)<<endl;
    }
    return 0;
}