杜教筛学习笔记

数学 学习笔记

---

upd：内容比较老旧，但是不想更新了，到时候会写在另一个博客里。

主要是上一篇字数太多，太卡了，于是杜教筛就换了一篇写 多水了一篇博客

先讲一下前置知识：

前置知识

莫比乌斯函数

莫比乌斯函数，即$\mu(x)=\cases{1,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x=1\\(-1)^s,\ \ \ \ \ \ \ x=p_1\cdot p_2\cdots p_s(\forall_{1\leqslant i\leqslant x,1\leqslant j\leqslant x,i\ne j}p_i\ne p_j)\\0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ otherwise}$

狄利克雷卷积

狄利克雷卷积，可以表示为

$$设r=f\times g，则r(n)=\sum_{d\mid n}f(d)\times g(\frac{n}{d})$$

狄利克雷卷积满足交换律和结合律，且存在单位元$\epsilon(n)=[n==1]$

简证$\epsilon$为单位元：

$$有一数论函数f，则(f\times\epsilon)(n)=\sum_{d\mid n} f(d)\times\epsilon(\frac{n}{d})=f(n)\times\epsilon(1)+\sum_{d\mid n,d\ne 1}f(d)\times\epsilon(\frac{n}{d})\\=f(n)\times 1+\sum_{d\mid n,d\ne 1}f(d)\times0=f(n)+0=f(n)$$

例子：

$$\epsilon=\mu\times I，即\epsilon(n)=\sum_{d\mid n}\mu(d)\times I(\frac{n}{d})=\sum_{d\mid n}\mu(d)=\cases{1,\ &n=1\\0,\ &n>1}=[n==1]\\d=I\times I，即d(n)=\sum_{d\mid n}I(d)\times I(\frac{n}{d})=\sum_{d\mid n}1=n的因子个数$$

重要的数论函数中互相转化：

$$\mu\cdot I=\epsilon\\\varphi\cdot I=id\\id\cdot I=\sigma\\I\cdot I=d\\\mu\cdot id=\varphi$$
可以列一个表：

$\ $	$\epsilon$	$I$	$d$	$\mu$	$\varphi$	$\sigma$	$id$
$\epsilon$	$\epsilon$	$I$	$d$	$\mu$	$\varphi$	$\sigma$	$id$
$I$	$I$	$d$	$-$	$\epsilon$	$id$	$d\cdot id$	$\sigma$
$d$	$d$	$-$	$-$	$I$	$\sigma$	$-$	$I\cdot\sigma$
$\mu$	$\mu$	$\epsilon$	$I$	$-$	$-$	$id$	$\varphi$
$\varphi$	$\varphi$	$id$	$\sigma$	$-$	$-$	$id\cdot id$	$-$
$\sigma$	$\sigma$	$d\cdot id$	$-$	$id$	$id\cdot id$	$d\cdot id\cdot id$	$-$
$id$	$id$	$\sigma$	$I\cdot\sigma$	$\varphi$	$-$	$-$	$\varphi\cdot\sigma$

公式推导

1. $$对于n\in\mathbb{Z^+}，\sum_{d\mid n}\mu(d)=\cases{1,\ &n=1\\0,\ &n>1}\\证明如下：\\\sum_{d\mid n}\mu(d)=\sum_{d\mid n,p^2\mid d(p\in\mathbb{P})}\mu(d)+\sum_{d\mid n,p^2\not\mid d(p\in\mathbb{P})}\mu(d)=0+\sum_{d\mid n,p^2\not\mid d(p\in\mathbb{P})}\mu(d)\\设n=\prod_{i=1}^s p_i^{k_i}，m=\prod_{i=1}^s p_i，则\sum_{d\mid n}\mu(d)=\sum_{d\mid n,p^2\not\mid d(p\in\mathbb{P})}\mu(d)=\sum_{d\mid m}\mu(d)\\设w\mid m，则有：\\当w有0个质因子，即w=1时，有1=C^0_s种选择，贡献为\mu(w)=1=(-1)^0\\当w有1个质因子时，有C^1_s种选择，贡献为\mu(w)=(-1)^1=-1\\当w有2个质因子时，有C^2_s种选择，贡献为\mu(w)=(-1)^2=1\\\vdots\\当w有i个质因子时，有C^i_s种选择，贡献为\mu(w)=(-1)^i\\\vdots\\当w有s个质因子时，有C^s_s种选择，贡献为\mu(w)=(-1)^n\\故\sum_{d\mid n}\mu(d)=\sum_{d\mid m}\mu(d)=\sum_{i=0}^sC^i_s\cdot (-1)^i=\sum_{i=0}^s C^i_s\cdot 1^{s-i}\cdot (-1)^i\\由二项式定理得\sum_{d\mid n}\mu(d)=\sum_{i=0}^s C^i_s\cdot 1^{s-i}\cdot (-1)^i=(1+(-1))^s=0^s=0$$
2. $$设I(x)=1,\epsilon(x)=[n==1]（即狄利克雷卷积的单位元），则有\mu为I的逆元，即\mu\times I=\epsilon，证明如下：\\由狄利克雷卷积得\mu\times I=\sum_{d\mid n}\mu(d)\times I(\frac{n}{d})=\sum_{d\mid n}\mu(d)\times 1=\sum_{d\mid n}\mu(d)\\由性质1得\mu\times I=\sum_{d\mid n}\mu(d)=\cases{1,\ &n=1\\0,\ &n>1}=[n==1]=\epsilon\\即\mu\times I=\epsilon，\mu为I的逆元$$
3. 莫比乌斯反演公式 $$设I(x)=1,\epsilon(x)=[n==1]，则对于数论函数f和F=f\times I，f=F\times\mu，证明如下：\\\because F=f\times I\\\therefore F\times\epsilon=f\times I\\由莫比乌斯函数性质2得\mu\times I=\epsilon，即F\times\mu\times I=f\times I\\\Rightarrow f=F\times\mu\\可以写成：当F(x)=\sum_{d\mid n}f(d)，则f(n)=\sum_{d\mid n}\mu(d)\times F(\frac{n}{d})$$
4. 莫比乌斯函数的积性证明 $$当\gcd(n,m)=1时，\mu(n)\times\mu(m)=\mu(n\times m)，证明如下：\\当n=1时，\mu(n)=1，即n\times m=m，\mu(n)\times\mu(m)=\mu(m)\\即\mu(n)\times\mu(m)=\mu(m)=\mu(1\times m)=\mu(n\times m)（m同理）\\当存在质数p使得p^2\mid n时，有\mu(n)=0，即\mu(n)\times\mu(m)=0\\此时肯定有p^2\mid n\times m，即\mu(n\times m)=0\\整理得\mu(n)\times\mu(m)=0=\mu(n\times m)（m同理）\\否则肯定有n=\prod_{i=1}^s p_i，且\forall_{i\ne j}p_i\ne p_j，m=\prod_{i=1}^t q_i，且\forall_{i\ne j}q_i\ne q_j\\\because\gcd(n,m)=1\\\therefore \forall_{i,j} p_i\ne q_j\\\Rightarrow\mu(n\times m)=(-1)^{s+t}\\\because\mu(n)=(-1)^s，\mu(m)=(-1)^t\\\therefore\mu(n)\times \mu(m)=(-1)^{s+t}=\mu(n\times m)$$
5. $$\sum_{d\mid\gcd(i,j)}\mu(d)=[\gcd(i,j)==1]\\由莫比乌斯函数性质1得\sum_{d\mid n}\mu(d)=[n==1]，将\gcd(i,j)带入就可以得到上式$$

整除分块

在莫比乌斯反演的应用题中，经常会有$\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor$的形式，但是$O(n)$的复杂度有时候仍然满足不了需求

为了解决这一问题，我们就要用到整除分块

很容易发现很多$\lfloor\frac{n}{i}\rfloor$的值是一样的，且所有的$\lfloor\frac{n}{i}\rfloor$为一个不下降子序列，呈块状分布

通过简单的计算可以得到，对于每个起点为$i$的块，它的值为$\lfloor\frac{n}{i}\rfloor$，终点为$\lfloor\frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor$，然后我们就可以用$O(\sqrt{n})$的算法计算$\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor$的值了

代码：

int l=1,r,ans=0;
if(n>m)
    swp(n,m);
while(l<=n){
    r=min(n/(n/l),mm/(m/l));
    ans+=(r-l+1)*(n/l);
    l=r+1;
}

杜教筛

我们要筛一个积性函数的前缀和，杜教筛可以帮我们做到$O(n^{\frac{3}{4}})$（当然可以记忆化常数优化）

我们设$S(n)=\sum_{i=1}^n f(i)$，其中$f(i)$为积性函数

我们找一个可以与$f$卷积的积性函数$g$，推一下它们的狄利克雷卷积前缀和：

$$\sum_{i=1}^n(f\cdot g)(i)=\sum_{i=1}^n\sum_{j\mid i}f(j)\cdot g(\frac{i}{j})=\sum_{j=1}^n g(j)\sum_{i=1}^{\lfloor\frac{n}{j}\rfloor}f(i)=\sum_{j=1}^n g(j)\cdot S(\lfloor\frac{n}{j}\rfloor)=\sum_{i=1}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)$$

我们可以递归求解$S(\lfloor\frac{n}{j}\rfloor)$，而且可以设一个边界$maxn$，如果$n<maxn$则直接输出预处理的前缀和

然后我们考虑$g(1)\cdot S(n)$的值：

$$g(1)\cdot S(n)=\sum_{i=1}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)-\sum_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)=\sum_{i=1}^n(f\cdot g)(i)-\sum_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)$$

此时，如果找到一个合适的$g$使$g$与$f\cdot g$的前缀和可以很容易处理，那么就可以用数论分块递归求解

我们可以由这个表得出一些数论函数相应的函数：

$\ $	$\epsilon$	$I$	$d$	$\mu$	$\varphi$	$\sigma$	$id$
$\epsilon$	$\epsilon$	$I$	$d$	$\mu$	$\varphi$	$\sigma$	$id$
$I$	$I$	$d$	$-$	$\epsilon$	$id$	$d\cdot id$	$\sigma$
$d$	$d$	$-$	$-$	$I$	$\sigma$	$-$	$I\cdot\sigma$
$\mu$	$\mu$	$\epsilon$	$I$	$-$	$-$	$id$	$\varphi$
$\varphi$	$\varphi$	$id$	$\sigma$	$-$	$-$	$id\cdot id$	$-$
$\sigma$	$\sigma$	$d\cdot id$	$-$	$id$	$id\cdot id$	$d\cdot id\cdot id$	$-$
$id$	$id$	$\sigma$	$I\cdot\sigma$	$\varphi$	$-$	$-$	$\varphi\cdot\sigma$

即：

$f$	$g$	$f\cdot g$
$d$	$\varphi$	$\sigma$
$\mu$	$I$	$\epsilon$
$\varphi$	$I$	$id$
$\sigma$	$\mu$	$id$

杜教筛应用

1.例题1：P4213 【模板】杜教筛（Sum）
题意：共$T$组数据，每组数据给定$n$，求$\sum_{i=1}^n\varphi(i)$和$\sum_{i=1}^n\mu(i)$

我们根据上面的表很容易写出来，因此没有什么推导过程

代码：

#include<stdio.h>
#include<map>
using namespace std;
const long long maxn=3000005;
long long i,j,k,m,n,T,cnt;
long long a[maxn],p[maxn],varphi[maxn],mu[maxn],sumvarphi[maxn],summu[maxn];
map<long long,long long>ansvarphi,ansmu;
long long getvarphi(long long n){
    if(n<maxn)
        return sumvarphi[n];
    if(ansvarphi.count(n))
        return ansvarphi[n];
    long long l=2,r,res=n*(n+1)/2;
    while(l<=n){
        r=n/(n/l);
        res-=(r-l+1)*getvarphi(n/l);
        l=r+1;
    }
    return ansvarphi[n]=res;
}
long long getmu(long long n){
    if(n<maxn)
        return summu[n];
    if(ansmu.count(n))
        return ansmu[n];
    long long l=2,r,res=1;
    while(l<=n){
        r=n/(n/l);
        res-=(r-l+1)*getmu(n/l);
        l=r+1;
    }
    return ansmu[n]=res;
}
int main(){
    p[1]=varphi[1]=mu[1]=1;
    for(i=2;i<maxn;i++){
        if(p[i]==0)
            a[++cnt]=i,varphi[i]=i-1,mu[i]=-1;
        for(j=1;j<=cnt;j++){
            if(i*a[j]>=maxn)
                break;
            p[i*a[j]]=1;
            if(i%a[j]==0){
                varphi[i*a[j]]=varphi[i]*a[j],mu[i*a[j]]=0;
                break;
            }
            varphi[i*a[j]]=varphi[i]*(a[j]-1),mu[i*a[j]]=-mu[i];
        }
    }
    for(i=1;i<maxn;i++)
        sumvarphi[i]=sumvarphi[i-1]+varphi[i],summu[i]=summu[i-1]+mu[i];
    scanf("%lld",&T);
    while(T--){
        scanf("%lld",&n);
        printf("%lld %lld\n",getvarphi(n),getmu(n));
    }
    return 0;
}

2.例题2：简单的数学题
题意：求$(\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot\gcd(i,j))\bmod p$
数据范围：$n\leqslant 10^{10},p\leqslant 1.1\cdot 10^9$
相信大家都看了之前莫比乌斯反演的讲解，因此这里过程简单一点：

$$\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot\gcd(i,j)=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[\gcd(i,j)==1](i\cdot d)\cdot(j\cdot d)\cdot d\\=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}i\cdot j\cdot d^3\cdot(\sum_{k\mid\gcd(i,j)}\mu(k))=\sum_{d=1}^n d^3\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum_{i=1}^{\lfloor\frac{n}{d\cdot k}\rfloor}(i\cdot k)\sum_{j=1}^{\lfloor\frac{n}{d\cdot k}\rfloor}(j\cdot k)\\=\sum_{d=1}^n d^3\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\cdot k^2\sum_{i=1}^{\lfloor\frac{n}{d\cdot k}\rfloor}i\sum_{j=1}^{\lfloor\frac{n}{d\cdot k}\rfloor}j$$
设$sum(i)=\frac{i\cdot(i+1)}{2}$，即$1$加到$i$的和
因此： $$\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot\gcd(i,j)=\sum_{d=1}^n\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}d^3\cdot\mu(k)\cdot k^2\cdot sum(\lfloor\frac{n}{d\cdot k}\rfloor)^2\\=\sum_{t=1}^n t^2\cdot sum(\lfloor\frac{n}{t}\rfloor)^2\cdot\sum_{d\mid t} d\cdot\mu(\frac{t}{d})$$
用狄利克雷卷积推一下： $$\varphi\cdot I=id\Rightarrow\varphi\cdot I\cdot\mu=id\cdot\mu\Rightarrow id\cdot\mu=\varphi$$
因此继续化简：
$$\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot\gcd(i,j)=\sum_{t=1}^n t^2\cdot sum(\lfloor\frac{n}{t}\rfloor)\cdot\varphi(t)$$
设$f(n)=n^2\cdot\varphi(n)$
我们很容易看出$f$是积性函数，因为：
由于$\varphi$为积性函数，若$\gcd(n,m)=1$，则$f(n)$与$f(m)$相乘得：
$$f(n)\cdot f(m)=n^2\cdot\varphi(n)\cdot m^2\cdot\varphi(m)=(n\cdot m)^2\cdot\varphi(n\cdot m)=f(n\cdot m)$$
因此我们可以用杜教筛来处理$f$的前缀和，原式化为：$\sum_{i=1}^n sum(\lfloor\frac{n}{t}\rfloor)\cdot f(t)$，此时原式可以数论分块来做
考虑如何筛$f$，我们设另一个数论函数为$g$，则它们的卷积为：$(f\cdot g)(x)=\sum_{d\mid x}f(d)\cdot g(\frac{x}{d})=\sum_{d\mid x}d^2\cdot\varphi(d)\cdot g(\frac{x}{d})$
我们尝试用$g$消去$d^2$这个项，设$g(x)=x^2$，然后可以推出来： $$(f\cdot g)(x)=\sum_{d\mid x} d^2\cdot\varphi(d)\cdot(\frac{x}{d})^2=\sum_{d\mid x}x^2\cdot\varphi(d)=x^2\cdot\sum_{d\mid x}\varphi(d)$$
由欧拉反演$\sum_{d\mid n}\varphi(d)=n$得： $$(f\cdot g)(x)=x^3$$
因此可以用杜教筛的式子： $$g(1)\cdot S(n)=\sum_{i=1}^n(f\cdot g)(i)-\sum_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\\Rightarrow 1^2\cdot S(n)=\sum_{i=1}^n i^3-\sum_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\\Rightarrow S(n)=\sum_{i=1}^n i^3-\sum_{i=2}^n i^2\cdot S(\lfloor\frac{n}{i}\rfloor)$$
我们可以用小学奥数的思路推一下$\sum_{i=1}i^3$
首先$i^3=i^3-i+i=i(i^2-1)+i=(i-1)\cdot i\cdot(i+1)+i$，因此： $$\sum_{i=1}i^3=\sum_{i=1}((i-1)\cdot i\cdot(i+1)+i)=\sum_{i=1}^n (i-1)\cdot i\cdot(i+1)+\sum_{i=1}^n i\\=\frac{1}{4}\sum_{i=1}^n(i-1)\cdot i\cdot(i+1)\cdot 4+\sum_{i=1}^ni=\frac{1}{4}\sum_{i=1}^n(i-1)\cdot i\cdot(i+1)\cdot((i+2)-(i-2))+\sum_{i=1}^ni\\=\frac{1}{4}\sum_{i=1}^n(-(i-2)\cdot(i-1)\cdot i\cdot(i+1)+(i-1)\cdot i\cdot(i+1)\cdot(i+2))+\sum_{i=1}^n i\\=\frac{1}{4}(-(-1)\cdot 0\cdot 1\cdot 2+0\cdot1\cdot2\cdot3-0\cdot1\cdot2\cdot3+1\cdot2\cdot3\cdot4-1\cdot2\cdot3\cdot4+\cdots+(n-1)\cdot n\cdot(n+1)\cdot(n+2))+\sum_{i=1}^n i\\=\frac{1}{4}(n-1)\cdot n\cdot(n+1)\cdot(n+2)-\frac{n\cdot(n+1)}{2}\\=\frac{n\cdot(n+1)\cdot((n-1)\cdot(n+2))}{4}+\frac{n\cdot(n+1)\cdot 2}{4}\\=\frac{n\cdot(n+1)\cdot(n^2+n-2)+n\cdot(n+1)\cdot 2}{4}\\=\frac{n\cdot(n+1)(n^2+n-2+2)}{4}\\=\frac{n\cdot(n+1)\cdot(n\cdot(n+1))}{4}\\=\frac{n^2\cdot(n+1)^2}{4}$$
那么$\sum_{i=1}^ni^2$肯定也很显然了，就不推了，最终的式子为$\frac{n\cdot(n+1)\cdot(2n+1)}{6}$
代入杜教筛的式子得： $$S(n)=\frac{n^2\cdot(n+1)^2}{4}-\sum_{i=2}^n i^2\cdot S(\lfloor\frac{n}{i}\rfloor)=sum(n)^2-\sum_{i=2}^n i^2\cdot S(\lfloor\frac{n}{i}\rfloor)$$
用数论分块和递归可以在很快的时间内用$\sum_{i=1}^ni^2=\frac{n\cdot(n+1)\cdot(2n+1)}{6}$计算出$\sum_{i=2}^n i^2\cdot S(\lfloor\frac{n}{i}\rfloor)$
然后预处理一些$S(n)$就可以了
最后代入原式： $$\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot\gcd(i,j)=\sum_{i=1}^n sum(\lfloor\frac{n}{t}\rfloor)\cdot f(t)$$
然后数论分块就结束了

代码：

#include<stdio.h>
#include<map>
using namespace std;
const long long maxn=5000005;
long long i,j,k,m,n,mod,cnt,inv6;
long long a[maxn],p[maxn],varphi[maxn],S[maxn];
map<long long,long long>ans;
inline long long getsum(long long n){
    return n*(n+1)/2%mod;
}
inline long long calc(long long n){
    return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
long long ksm(long long a,long long b,long long mod){
    long long res=1;
    while(b){
        if(b&1)
            res=(res*a)%mod;
        a=(a*a)%mod,b>>=1;
    }
    return res;
}
long long getS(long long n){
    if(n<maxn)
        return S[n];
    if(ans.count(n))
        return ans[n];
    long long l=2,r,res=getsum(n)*getsum(n)%mod;
    while(l<=n){
        r=n/(n/l);
        res=((res-(calc(r)-calc(l-1)+mod)%mod*getS(n/l)%mod)%mod+mod)%mod;
        l=r+1;
    }
    res=(res%mod+mod)%mod;
    return ans[n]=res;
}
long long getans(long long n){
    long long l=1,r,res=0;
    while(l<=n){
        r=n/(n/l);
        res=(res+(getS(r)-getS(l-1)+mod)%mod*getsum(n/l)%mod*getsum(n/l)%mod)%mod;
        l=r+1;
    }
    res=(res%mod+mod)%mod;
    return res;
}
int main(){
    scanf("%lld%lld",&mod,&n);
    inv6=ksm(6,mod-2,mod);
    p[1]=varphi[1]=1;
    for(i=2;i<maxn;i++){
        if(p[i]==0)
            a[++cnt]=i,varphi[i]=i-1;
        for(j=1;j<=cnt;j++){
            if(i*a[j]>=maxn)
                break;
            p[i*a[j]]=1;
            if(i%a[j]==0){
                varphi[i*a[j]]=varphi[i]*a[j];
                break;
            }
            varphi[i*a[j]]=varphi[i]*(a[j]-1);
        }
    }
    for(i=1;i<maxn;i++)
        S[i]=(S[i-1]+varphi[i]%mod*i%mod*i%mod)%mod;
    printf("%lld\n",getans(n));
    return 0;
}