Luogu P4213 【模板】杜教筛 题解

前置知识

积性函数

顾名思义，积性函数是一类满足 \(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\)：

1. \(\mu*I=\epsilon\) 源于 \(\sum\limits_{d|n} \mu(d)=[n=1]\)
2. \(\varphi*I=id\) 源于 \(\sum\limits_{d|n}\varphi(d)=n\)
3. \(\mu*id=\varphi\) 源于 \(\varphi(n)=\sum\limits_{i=1}^n [gcd(i,n)=1]=\sum\limits_{d|n}\frac{n}{d}\mu(d)\)
4. \(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;
}