「算法笔记」杜教筛

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·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. 复杂度证明

提供一种不用积分的证法。

先不考虑预处理。首先我们每次记忆化搜索用到的 \(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:

裴老板指正:并不能按 \(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 \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;
} 

 

posted @ 2021-04-02 20:44  maoyiting  阅读(187)  评论(0编辑  收藏  举报