「算法笔记」杜教筛

upd on 2022.1.13：小修一波，并加了复杂度证明。

一、前置概念

具体在 「算法笔记」莫比乌斯反演 写过，所以「前置概念」就简单写写。积性函数和完全积性函数就不写了。

狄利克雷卷积：对于两个数论函数 \(f,g\)，定义它们的狄利克雷卷积 \(h=f*g\) 为

\[\displaystyle h(x)=\sum_{d\mid x}f(d)g\left(\frac{x}{d}\right) \]

狄利克雷卷积满足交换律、结合律、对加法的分配律，有单位元 \(\varepsilon\)（其中 \(\varepsilon\) 为单位函数 \(\varepsilon(x)=[x=1]\)）。若 \(f,g\) 是积性函数，则 \(f*g\) 也是积性函数。

在狄利克雷卷积的意义下，\(\mu*1=\varepsilon\)，即 \(\mu\) 和 \(1\) 互为逆元。

常用卷积：

• \(\varepsilon=\mu*1\)，即 \(\varepsilon(n)=\sum_{d\mid n}\mu(d)\)。
• \(\varphi*1=\text{Id}\)，即 \(\sum_{d\mid n}\varphi(d)=n\)。两边同时卷一个 \(\mu\) 可得 \(\varphi=\text{Id}*\mu\)，即 \(\varphi(n)=\sum_{d\mid n}d\cdot \mu(\frac n d)\)。
• \(d=1*1 \Leftrightarrow d(n)=\sum_{d\mid n}1\)，\(\sigma_k=\text{Id}_k*1 \Leftrightarrow \sigma_k(n)=\sum_{d|n} d^k\)

提共因式：对于任意数论函数 \(f,g\) 和 完全积性函数 \(h\) 均有 \((f*g)\cdot h=(f\cdot h)*(g\cdot h)\)（注意中间是点乘，\((f\cdot g)(i)=f(i)g(i)\)），展开即可证。比如 \((\varphi\cdot \text{Id}_k)* (1\cdot \text{Id}_k)=(\varphi*1)\cdot \text{Id}_k=\text{Id}\cdot \text{Id}_k=\text{Id}_{k+1}\)。

二、杜教筛

1. 基本思想

对于数论函数 \(f\)，求 \(S(n)=\sum_{i=1}^nf(i)\)。

假如我们能快速计算 \(g\) 和 \(f*g\) 的前缀和，考虑怎么求 \(S\)。

\[\sum_{i=1}^n (f*g)(i)=\sum_{i=1}^n\sum_{d\mid i}f(\frac i d)g(d)=\sum_{d=1}^n g(d)\sum_{i=1}^{\lfloor\frac n i\rfloor}f(\frac{id}{d})=\sum_{d=1}^n g(d)\sum_{i=1}^{\lfloor\frac n i\rfloor}f(i)=\sum_{i=1}^n g(i)S(\lfloor\frac n i\rfloor) \]

则可以得到 \(S(n)\) 关于 \(S(\lfloor \frac{n}{i}\rfloor)\) 的递推式：

\[g(1)S(n)=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \]

如果我们可以快速求出 \(\sum_{i=1}^n(f*g)(i)\)，再用整除分块来求 \(\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)，就能求得 \(g(1)S(n)\)。也就是说，我们要找到一个合适的 \(g\)，使得 \(f*g\) 和 \(g\) 都能快速求前缀和。

\(S(\lfloor\frac{n}{i}\rfloor)\) 是一个子问题，对 \(S(n)\) 递归求解并记忆化即可。

线性筛预处理前 \(n^{\frac{2}{3}}\) 个 \(f(x)\) 的值，则时间复杂度为 \(\mathcal O(n^{\frac{2}{3}})\)。

（当然如果我们能快速求出 \(f,g\) 的前缀和，也可以求 \(f*g\) 的前缀和，根据最初那个式子整除分块右边推左边就好了）

2. 复杂度证明

一种不用积分的证法（不保证我没写错 qwq）。

先不考虑预处理。首先我们每次记忆化搜索用到的 \(S\) 可以写成 \(S(\lfloor\large\frac n k\normalsize\rfloor)\) 的形式，再加上每次整除分块用的 \(\sqrt x\)，总复杂度就是 \(\sum_{\large x,\exists k, \lfloor\frac n k\rfloor=x}\sqrt x\)。

• Lemma：\(\lfloor\sqrt x\rfloor=p\)，那么 \(x\in[p^2,(p+1)^2)\)，即 \(x\) 有 \((p+1)^2-p^2=\mathcal O(p)\) 个不同的取值（下面的根号为向下取整）。

• 若 \(x=\lfloor\large\frac n k\normalsize\rfloor\leq \sqrt n\)，对和式的贡献为 \(\sum\limits_{x=1}^{\sqrt x}\sqrt x\)，枚举 \(i=\sqrt x\)，根据 Lemma 有 \(\sum\limits_{i=1}^{n^{\frac 1 4}}i^2\)（有 \(\mathcal O(i)\) 个 \(\sqrt x=i\) 的 \(x\)，每个贡献 \(i\)），平方求和式个三次方的东西所以有 \(\mathcal O((n^\frac 1 4)^3)=\mathcal O(n^{\frac 3 4})\)。

• 若 \(x=\lfloor\large\frac n k\normalsize\rfloor>\sqrt n\)，那么对应的 \(k\leq \sqrt n\)，对和式的贡献为 \(\sum\limits_{k=1}^{\sqrt n}\sqrt{\large\frac n k\normalsize }=\sqrt n\sum\limits_{k=1}^{\sqrt n}\frac{1}{\sqrt k}\)，枚举 \(i=\sqrt k\)，根据 Lemma 有 \(\sqrt n\sum\limits_{i=1}^{n^{\frac 1 4}}\frac{i}{i}=n^{\frac 3 4}\)。

因此总复杂度也是 \(\mathcal O(n^{\frac 3 4})\) 级别的。还可以进一步优化杜教筛，即先线性筛出前 \(m\) 个答案，之后再用杜教筛。假设 \(m>\sqrt n\)，优化后复杂度变为 \(\mathcal O(m+\sum\limits_{k=1}^{\large\lfloor\frac n m\rfloor}\sqrt{\large\frac n k})=\mathcal O(m+\sqrt n\cdot \sqrt{\large\frac n m})\)，根据均值不等式当 \(m=\sqrt n\cdot \sqrt{\large\frac n m}\) 时复杂度最优，此时 \(m=n^{\frac 2 3}\)，时间复杂度达到最小值 \(\mathcal O(n^{\frac 2 3})\)。

实测 1s 杜教筛能跑 \(n=10^{10}\)。

upd on 2022.1.14：

pwj 指正：并不能按 \(m=\sqrt n\cdot \sqrt{\large\frac n m}=\large\frac{n}{\sqrt m}\) 算，因为 \(m+\large\frac{n}{\sqrt m}\normalsize\geq 2\sqrt{m\cdot \large\frac{n}{\sqrt m}}\) 会发现右边的 \(m\) 消不掉，而 \(m\) 是我们自己定的。但是数量级没错。

这里有个拆项凑系数用来证不等式的 trick：

\[m+\frac{n}{\sqrt m}=m+\frac{\frac n 2}{\sqrt m}+\frac{\frac n 2}{\sqrt m}\geq 3\sqrt[3]{m\cdot \frac{\frac n 2}{\sqrt m}\cdot \frac{\frac n 2}{\sqrt m}}=3\sqrt[3]{\frac{n^2}{4}}=\mathcal O(n^{\frac 2 3}) \]

当 \(m=\large\frac{\frac n 2}{\sqrt m}\) 时复杂度最优。

3. 例子

求 \(S(n)=\sum_{i=1}^n\varphi(i)\)。

\(\varphi*1=\text{Id}\)。显然 \(1\)、\(\text{Id}(x)=x\) 都可以快速求前缀和。取 \(g(x)=1\) 即可。

\(S(n)=\sum_{i=1}^n i-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)=\frac{n(n+1)}{2}-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)。

求 \(S(n)=\sum_{i=1}^n\mu(i)\)。

\(\mu*1=\varepsilon\)。显然 \(1\) 和 \(\varepsilon(x)=[x=1]\) 可以快速求前缀和。取 \(g(x)=1\) 即可。

\(S(n)=1-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)。

//Luogu P4213
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e6+5;
int t,n=2e6,cnt,p[N],phi[N],mu[N];
bool vis[N];
map<int,int>Smu,Sphi;
int S_mu(int n){
    if(n<=2e6) return mu[n];
    if(Smu[n]) return Smu[n];
    int ans=0;
    for(int l=2,r=0;l<=n;l=r+1)    //注意从 2 开始
        r=n/(n/l),ans+=(r-l+1)*S_mu(n/l); 
    return Smu[n]=1-ans;
}
int S_phi(int n){
    if(n<=2e6) return phi[n];
    if(Sphi[n]) return Sphi[n];
    int ans=0;
    for(int l=2,r=0;l<=n;l=r+1)
        r=n/(n/l),ans+=(r-l+1)*S_phi(n/l);
    return Sphi[n]=n*(n+1)/2-ans; 
}
signed main(){
    scanf("%lld",&t);
    vis[0]=vis[1]=1,phi[1]=mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!vis[i]) p[++cnt]=i,phi[i]=i-1,mu[i]=-1; 
        for(int j=1;j<=cnt&&i*p[j]<=n;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j],mu[i*p[j]]=0;break;}
            phi[i*p[j]]=phi[i]*phi[p[j]],mu[i*p[j]]=-mu[i]; 
        }
    }
    for(int i=1;i<=n;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
    while(t--){
        scanf("%lld",&n);
        printf("%lld %lld\n",S_phi(n),S_mu(n));
    }
    return 0;
} 

再比如说求 \(S(n)=\sum_{i=1}^n i\cdot \varphi(i)\)。

很明显 \(f(x)=x\cdot\varphi(x)\) 是积性函数。我们尝试构造出合适的 \(g(x)\)。

\[(f*g)(x)=\sum_{d\mid x}d\cdot \varphi(d)\cdot g(\frac x d) \]

考虑将 \(d\) 消掉，这样可以利用 \(\sum_{d\mid x}\varphi(d)=x\) 推出一些东西。尝试 \(g=\text{Id}\)。

\[(f*g)(x)=\sum_{d\mid x}d\cdot \varphi(d)\cdot \frac{x}{d}=x\sum_{d\mid x}\varphi(d)=x^2 \]

取 \(g=\text{Id}\) 即可。

求 \(S(n)=\sum_{i=1}^n\varphi(i)\cdot i^k\)。

按照提公因式卷上一个 \(1\cdot \text{Id}_k\)，那么 \((\varphi\cdot \text{Id}_k)*(1\cdot\text{Id}_k)=(\varphi*1)\cdot \text{Id}_k=\text{Id}\cdot\text{Id}_k=\text{Id}_{k+1}\)。

取 \(g=\text{Id}_k\) 即可。

求 \(S(n)=\sum_{i=1}^n \mu(i)\cdot i^k\)。

按照提公因式卷上一个 \(1\cdot \text{Id}_k\)，那么 \((\mu\cdot \text{Id}_k)*(1\cdot \text{Id}_k)=(\mu*1)\cdot \text{Id}_k=\varepsilon\cdot \text{Id}_k=\varepsilon\)。取 \(g=\text{Id}_k\) 即可。

三、狄利克雷生成函数

有些时候不容易看出 \(g(x)\) 取什么，比如：求 \(S(n)=\sum_{i=1}^n i\cdot \varphi(i)\)。凑是一种方法；而使用 狄利克雷生成函数（DGF） 可以从另一个角度求出 \(g(x)\)。

对于无穷序列 \(f_1,f_2,\cdots\)，定义其狄利克雷生成函数为

\[\tilde F(x)=\sum_{i\geq 1}\frac{f_i}{i^x} \]

与普通生成函数以及指数生成函数的对比：

普通生成函数：\(F(x)=\sum_{i\geq 0}f_ix^i\)。

指数生成函数：\(\hat F(x)=\sum_{i\geq 0}f_i\frac{x^i}{i!}\)。

对于两个序列 \(f,g\)，其 DGF 的乘积对应 \(f,g\) 的狄利克雷卷积序列：

\[\tilde F(x)\tilde G(x)=\sum_{i\geq 1}\sum_{j\geq 1}\frac{f_i}{i^x}\frac{g_j}{j^x}=\sum_{i\geq 1}\sum_{j\geq 1}\frac{f_ig_j}{(ij)^x}=\sum_d\sum_{i\mid d}\frac{f_ig_{\frac{d}{i}}}{d^x} \]

如果序列 \(f\) 满足积性（积性函数），其 DGF 可以由质数幂处的取值表示：

\[\tilde F=\prod_{p\in\text{prime}}(1+\frac{f(p)}{p^x}+\frac{f(p^2)}{p^{2x}}+\frac{f(p^3)}{p^{3x}}+\cdots) \]

可以考虑 \(i=\prod p_j^{e_j}\)，\(f(i)=\prod f(p_j^{e_j})\)。那么有：

\[\frac{f_i}{i^x}=\prod \frac{f(p_j^{e_j})}{(p_j^{e_j})^x} \]

咕咕咕……

四、例题

Luogu P3768 简单的数学题

Problem：给定 \(n,p\)，求：

\[\left(\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot \gcd(i,j)\right)\bmod p \]

\(n\leq 10^{10},5\times 10^8\leq p\leq 1.1\times 10^9\)，\(p\) 是质数。

Solution：

可以用 \(\varphi*1=\text{Id}\) 推式子，即 \(\sum_{d\mid x}\varphi(d)=x\)。

\[ans=\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot \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)^2 \]

设 \(s(n)=(\sum_{i=1}^n i)^2=(\frac{(1+n)n}{2})^2\)，则 \(ans=\sum_{d=1}^n\varphi(d)d^2s(\lfloor\frac{n}{d}\rfloor)\)，如果能快速求出 \(f(x)=\varphi(x)x^2\) 的前缀和，那么整除分块就做完了。

考虑用杜教筛来筛 \(f(x)=\varphi(x)x^2\)，要找到一个合适的 \(g\)，使得 \(f*g\) 和 \(g\) 都可以快速求前缀和。

\[(f*g)(x)=\sum_{d\mid x}\varphi(d)d^2g(\frac x d) \]

考虑把这里的 \(d^2\) 消掉，我们尝试 \(g=\text{Id}_2\)（幂函数 \(\text{Id}_k(n)=n^k\)），那么：

\[(f*g)(x)=\sum_{d\mid x}\varphi(d)d^2(\frac x d)^2=\left(\sum_{d\mid x}\varphi(d)\right)x^2=x\cdot x^2=x^3 \]

取 \(g=\text{Id}_2\) 即可，\(f*g=\text{Id}_3\)。当然也可以用 狄利克雷生成函数（DGF） 推出。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e6+5;
int n,m=2e6,mod,cnt,p[N],phi[N],inv2,inv6,ans;
bool vis[N];
map<int,int>s;
int mul(int x,int n,int mod){
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans;
}
int s2(int n){    //平方和 
    n%=mod;
    return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
int s3(int n){    //立方和 
    n%=mod; 
    int x=n%mod*(n+1)%mod*inv2%mod;
    return x*x%mod;
}
int S(int n){
    if(n<=m) return phi[n];
    if(s[n]) return s[n];
    int ans=0;
    for(int l=2,r=0;l<=n;l=r+1)
        r=n/(n/l),ans=(ans+(s2(r)-s2(l-1))%mod*S(n/l)%mod)%mod;
    return s[n]=((s3(n)-ans)%mod+mod)%mod;
}
signed main(){
    scanf("%lld%lld",&mod,&n);
    vis[0]=vis[1]=1,phi[1]=1;
    for(int i=2;i<=m;i++){
        if(!vis[i]) p[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*p[j]<=m;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
            phi[i*p[j]]=phi[i]*phi[p[j]]; 
        }
    }
    for(int i=1;i<=m;i++) phi[i]=(phi[i-1]+i*i%mod*phi[i]%mod)%mod;
    inv2=mul(2,mod-2,mod),inv6=mul(6,mod-2,mod);
    for(int l=1,r=0;l<=n;l=r+1)
        r=n/(n/l),ans=(ans+s3(n/l)*(S(r)-S(l-1))%mod)%mod;
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}