<学习笔记> 杜教筛

合集 - 学习笔记(16)

1.<学习笔记> 组合数学2023-06-202.<学习笔记> 整除分块2023-08-133.<学习笔记> 线段树分治2023-09-204.<学习笔记> 拉格朗日插值2023-10-195.<学习笔记> 二分图2023-10-246.<学习笔记> 点分树2023-11-027.<学习笔记> 网络流2023-12-058.<学习笔记> 四边形不等式2023-12-209.<学习笔记> 后缀树(数)组2023-12-2610.<学习笔记> manacher 和 回文自动机2023-12-2611.<学习笔记> SAM2024-01-0112.<学习笔记> 斜率优化2024-01-0713.<学习笔记> 莫比乌斯反演2023-08-14

14.<学习笔记> 杜教筛2024-01-21

15.<学习笔记> 二项式反演2024-02-0216.<学习笔记> 笛卡尔树2024-02-06

收起

杜教筛

处理数论函数的前缀和问题，可以在低于线性的复杂度里求出 $S(n)=\sum _{i=1}^{n}f(i)$。

对于任意一个数论函数 $g$，必须满足 ：

$$\sum _{i=1}^n(f\ast g)(i)=\sum _{i=1}^n\sum _{d\mid i}g(d)\ast d(\frac{i}{d})$$

$$=\sum _{d=1}^{n}g(d)\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }f(j)=\sum _{d=1}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )$$

那么可以得到：

$$g(1)S(n)=\sum _{i=1}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor )-\sum _{i=2}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor )$$

$$=\sum _{i=1}^n(f\ast g)(i)-\sum _{i=2}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor )$$

时间复杂度最小为 $n^{\frac{2}{3}}$。

一般情况下，可以设出方便求前缀的 $g$ 和 $f\ast g$，然后进行递归求解。

例如根据 $I(n)=1$，$id(n)=n$，$ϵ(n)=[n=1]$

单位元 $ϵ$ 满足 $f\ast ϵ=f$

根据狄利克雷卷积得到几个性质：

1. $\mu \ast I=ϵ$

2. $\phi \ast I=id$

3. $\mu \ast id=\phi$

对于一，我们直接展开就是 $(\mu \ast I)(n)=\sum _{d\mid n}\mu (d)=ϵ$

对于二，还是展开 $(\phi \ast I)(n)=\sum _{d\mid n}\phi (d)$，然后可以证明 $\sum _{d\mid n}\phi (d)=n$。

证明

设 $f(x)$ 表示 $gcd(k,n)=x,(k<n)$ 的个数，那么有 $n=\sum _{i=1}^{n}f(i)$。

我们还知道 $f(x)=\phi (\frac{n}{x})$ 相当于就是 $gcd(k,n)=x$ 的个数相当于 $gcd(\frac{k}{x},\frac{n}{x})=1$ 的个数，那么就有 $n=\sum _{i\mid n}\phi (i)$。

对于三，我们可以另 $2$ 的右式卷上一个 $1$ 的左式，是不变的，消掉就可以证明。

一些题#

DZY Loves Math IV #

$$\mathrm{求}\sum _{i=1}^n\sum _{j=1}^m\phi (ij)$$

发现 $n$ 的范围很小，所以我们考虑枚举 $n$。

我们设 $n=\sum _{i=1}^{w(n)}{p}_i^{k_i}$， $p=\sum _{i=1}^{w(n)}{p}_i$， $q=\frac{n}{p}$，那么 $\phi (n)=q\phi (p)$，也就有 $\phi (ni)=q\phi (pi)$。以下有好几条性质是由于 $p$ 的质因子系数为一才会成立。

我们另

$$S(n,m)=q\sum _{i=1}^m\phi (p\ast j)$$

$$=q\sum _{i=1}^m\phi (\frac{p}{gcd(p,i)}gcd(p,i)\ast i)$$

$$=q\sum _{i=1}^m\phi (\frac{p}{gcd(p,i)})\ast \phi (gcd(p,i)\ast i)$$

$$=q\sum _{i=1}^m\phi (\frac{p}{gcd(p,i)})\ast \phi (i)\ast gcd(p,i)$$

$$=q\sum _{i=1}^m\phi (\frac{p}{gcd(p,i)})\ast \phi (i)\ast \sum _{d\mid gcd(p,i)}\phi (d)$$

$$=q\sum _{i=1}^m\phi (i)\ast \sum _{d\mid gcd(p,i)}\phi (\frac{p\ast d}{gcd(p,i)})$$

$$=q\sum _{i=1}^m\phi (i)\ast \sum _{d\mid gcd(p,i)}\phi (\frac{p}{d})$$

$$=q\sum _{i=1}^m\phi (i)\ast \sum _{d\mid p,d\mid i}\phi (\frac{p}{d})$$

$$=q\sum _{i=1}^m\phi (i)\ast \sum _{d\mid p,d\mid i}\phi (\frac{p}{d})$$

尝试枚举 $d$

$$=q\sum _{d\mid p}\phi (\frac{p}{d})\ast \sum _{i=1}^{\lfloor \frac{m}{d}\rfloor }\phi (i\ast d)$$

$$=q\sum _{d\mid p}\phi (\frac{p}{d})\ast S(d,\lfloor \frac{m}{d}\rfloor )$$

当 $n$ 等于 $1$ 时我们直接暴力求 $m^{\frac{2}{3}}$，我们可以对 $S$ 进行记忆化。

code

#include<bits/stdc++.h>
const int mod=998244353;
#define int long long
using namespace std;
const int N=1e5;
int prime[N+5];
bool nprime[N+5];
int np[N+5],phi[N+5],sumphi[N+5];
void get_prime(){
    phi[1]=1;
    for(int i=2;i<=N;i++){
        if(!nprime[i]){
            prime[++prime[0]]=i;
            phi[i]=i-1;
            np[i]=i;
        }
        for(int j=1;j<=prime[0] && i*prime[j]<=N;j++){
            int k=i*prime[j];
            nprime[k]=1;
            if(i%prime[j]==0){
                np[k]=np[i];
                phi[k]=phi[i]*prime[j];
                break;
            }
            phi[k]=phi[i]*phi[prime[j]];
            np[k]=np[i]*np[prime[j]];
        }
    }
    for(int i=1;i<=N;i++) sumphi[i]=(sumphi[i-1]+phi[i])%mod;
}
unordered_map<int,int> Sphi,Ans[N+5];
int get_phi(int n){
    if(n<=N) return sumphi[n];
    if(Sphi[n]) return Sphi[n];
    int l=2,r;
    int ans=(n+1)*n/2;
    while(l<=n){
        int t=n/l;
        int r=n/t;
        ans=(ans-get_phi(n/l)*(r-l+1)%mod+mod)%mod;
        l=r+1;
    }
    return Sphi[n]=ans;
}
int S(int n,int m){
    if(m==0) return 0;
    if(Ans[n][m]) return Ans[n][m];
    if(n==1){
        Ans[n][m]=get_phi(m);
        return Ans[n][m];
    }
    int ans=n/np[n];
    int sum=0;
    int t=sqrt(np[n]);
    for(int d=1;d<=t;d++){
        if(np[n]%d==0){
            sum=(sum+phi[np[n]/d]*S(d,m/d)%mod)%mod;
            if(d!=np[n]/d){
                int tmp=np[n]/d;
                sum=(sum+phi[np[n]/tmp]*S(tmp,m/tmp)%mod)%mod;
            }
        }
    }
    return Ans[n][m]=ans*sum%mod;
}
signed main(){
    get_prime();
    int T;
    scanf("%lld",&T);
    while(T--){
        int n,m;
        scanf("%lld%lld",&n,&m);
        int ans=0;
        for(int i=1;i<=n;i++){
            ans=(ans+S(i,m))%mod;
        }
        printf("%lld\n",ans);
    }
}

Pyh 的求/毒瘤之神的考验 #

比上面那道题多一个多测，并且 $n,m$ 都为 $1e5$ 级别。

发现上面做法行不通，答案是没啥关系，唯一的关系是我都不会。

有一个结论：

$$\phi (i\ast j)=\frac{\phi (i)\ast \phi (j)\ast gcd(i,j)}{\phi (gcd(i,j))}$$

证明

$$\phi (i)\phi (j)=i\prod _{p\mid i}\frac{p-1}{p}j\prod _{p\mid j}\frac{p-1}{p}$$

$$=ij\prod _{p\mid i}\frac{p-1}{p}\prod _{p\mid j}\frac{p-1}{p}$$

$$=ij\prod _{p\mid ij}\frac{p-1}{p}\prod _{p\mid gcd(i,j)}\frac{p-1}{p}$$

如果在乘上一个 $gcd(i,j)$ 那么就可以变成 $\phi (ij)\ast \phi (gcd(i,j))$，得证。

那么我们那它推式子

$$\sum _{i=1}^n\sum _{i=1}^m\frac{\phi (i)\ast \phi (j)\ast gcd(i,j)}{\phi (gcd(i,j))}$$

$$\sum _{d=1}^n\sum _{i=1}^n\sum _{i=1}^m\frac{\phi (i)\ast \phi (j)\ast d}{\phi (d)}[gcd(i,j)=d]$$

$$\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{i=1}^{\lfloor \frac{m}{d}\rfloor }\frac{\phi (i)\ast \phi (j)\ast d}{\phi (d)}[gcd(i,j)=1]$$

$$\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{i=1}^{\lfloor \frac{m}{d}\rfloor }\frac{\phi (d\ast i)\ast \phi (d\ast j)\ast d}{\phi (d)}\sum _{k\mid i,k\mid j}\mu (k)$$

变为枚举 $k$

$$\sum _{d=1}^n\frac{d}{\phi (d)}\sum _{k=1}^{\lfloor \frac{n}{d}\rfloor }\mu (k)\sum _{i=1}^{\lfloor \frac{n}{kd}\rfloor }\phi (d\ast k\ast i)\sum _{i=1}^{\lfloor \frac{m}{kd}\rfloor }\phi (d\ast k\ast j)$$

另 $T=dk$

$$\sum _{T=1}^n\sum _{d\mid T}\mu (\frac{T}{d})\frac{d}{\phi (d)}\sum _{i=1}^{\lfloor \frac{n}{T}\rfloor }\phi (T\ast i)\sum _{i=1}^{\lfloor \frac{m}{T}\rfloor }\phi (T\ast j)$$

我们可以套路的另 $F(T)=\sum _{d\mid T}\mu (\frac{T}{d})\frac{d}{\phi (d)}$，这个也是很好预处理的。

另 $G(x,y)=\sum _{i=1}^x\phi (i\ast y)$，并且存在递推 $G(x,y)=G(x-1,y)+\phi (xy)$。

那么答案就变成了 $Ans=\sum _{T=1}^{n}F(T)G(\lfloor \frac{n}{T}\rfloor ,T)G(\lfloor \frac{m}{T}\rfloor ,T)$。

但是后面的那个东西好像无法直接整除分块，所以我们考虑根号分治，直接预处理一部分答案，我们另 $Ans(x,y,z)=\sum _{T=1}^{x}F(T)G(y,T)G(z,T)$，有转移 $Ans(x,y,z)=Ans(x-1,y,z)+F(x)\ast G(y,x)\ast G(z,x)$

我们定义阈值为 $B$，首先我们预处理出来 $G(x,y)(x\le N,y\le B)$， $Ans(x,y,z)(x\le N,y\le B,z\le B)$，对于 $T\le B$ 的我们直接暴力扫一遍，对于 $T>B$ 的，我们进行整除分块，直接查询 $Ans(r,\lfloor \frac{n}{T}\rfloor ,\lfloor \frac{m}{T}\rfloor )-Ans(l-1,\lfloor \frac{n}{T}\rfloor ,\lfloor \frac{m}{T}\rfloor )$ 就可以。

复杂度：预处理的话是 $n\ln n+n{B}^2$，单次查询的话 $\frac{m}{B}+\sqrt{n}$

code

#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
const int N=1e5;
const int B=50;
int prime[N+2];
bool nprime[N+2];
int phi[N+2],mul[N+2],F[N+2],inv[N+2],sumF[N+2];
vector<int> G[N+2];
vector<int> Ans[B+2][B+2];
inline int qpow(int x,int p){
    int ans=1;
    while(p){
        if(p&1) ans=(1ll*ans*x)%mod;
        x=(1ll*x*x)%mod;
        p>>=1;
    }
    return ans;
}
inline void get_prime(){
    phi[1]=mul[1]=1;
    for(int i=2;i<=N;i++){
        if(!nprime[i]){
            phi[i]=i-1;
            mul[i]=mod-1;
            prime[++prime[0]]=i;
        }
        for(int j=1;j<=prime[0] && i*prime[j]<=N;j++){
            int k=i*prime[j];
            nprime[k]=1;
            if(i%prime[j]==0){
                mul[k]=0;
                phi[k]=1ll*phi[i]*prime[j]%mod;
                break;
            }
            phi[k]=1ll*phi[i]*phi[prime[j]]%mod;
            mul[k]=1ll*mul[i]*mul[prime[j]]%mod;
        }
    }
}
int ask_G(int x,int y){
    if(y<=B) return G[x][y];
    int ans=0;
    for(int i=1;i<=x;i++) ans=(1ll*ans+phi[i*y])%mod;
    return ans;
}
inline int read(){
	int x(0);bool f(0);char ch=getchar();
	for(;ch<'0'||ch>'9';ch=getchar()) f^=ch=='-';
	for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+(ch^48);
	return f?x=-x:x;
}
inline void write(int x){
	x<0?x=-x,putchar('-'):0;static short Sta[50],top(0);
	do{Sta[++top]=x%10;x/=10;}while(x);
	while(top) putchar(Sta[top--]|48);
	putchar('\n');
}
signed main(){
    get_prime();
    for(int i=1;i<=N;i++) inv[i]=qpow(i,mod-2);
    for(int d=1;d<=N;d++){
        for(int k=1;d*k<=N;k++){
            int T=d*k;
            F[T]=(1ll*F[T]+1ll*d*inv[phi[d]]%mod*mul[k]%mod)%mod;
        }
    }
    for(int i=1;i<=N;i++) sumF[i]=(1ll*sumF[i-1]+F[i])%mod;
    G[0].resize(N+5);
    for(int x=1;x<=N;x++){
        G[x].resize(N/x+5);
        for(int y=1;x*y<=N;y++){
            G[x][y]=(1ll*G[x-1][y]+phi[x*y])%mod;
        }
    }
    for(int y=1;y<=B;++y){
        for(int z=1;z<=B;++z){
            Ans[y][z].resize(min(N/y,N/z)+5);
            for(int x=1;x*y<=N && x*z<=N;++x){
                Ans[y][z][x]=(1ll*Ans[y][z][x-1]+1ll*F[x]*ask_G(y,x)%mod*ask_G(z,x)%mod)%mod;
            }
        }
    }
    int T;
    T=read();
    while(T--){
        int n,m;
        n=read(),m=read();
        if(n>m) swap(n,m);
        long long ans=0;
        for(int t=1;t<=m/B;t++){
            ans=(ans+1ll*F[t]*G[n/t][t]%mod*G[m/t][t]%mod)%mod;
        }
        int l=m/B+1,r=0;
        while(l<=n){
            int ta=n/l,tb=m/l;
            if(!ta || !tb) break;
            r=min(n/ta,m/tb);
            ans=(ans+1ll*(Ans[ta][tb][r]-Ans[ta][tb][l-1]+mod)%mod)%mod;
            l=r+1;
        }
        write(ans);
    }
}

简单的最大公约数#

首先莫反做法，设 $f(x)$ 为 $gcd$ 为 $x$ 的个数，设 $g(x)$ 表示 $gcd$ 为 $x$ 倍数的个数，那么就有 $f(x)=\sum _{i=1}^{\frac{m}{x}}g(x\ast i)\mu (i)=\lfloor \frac{m}{x}{\rfloor }^n$

那么答案就是 $\sum _{x=1}^{m}x\ast f(x)$，展开后另 $T=x\ast d$ 就有

$$\sum _{T=1}^{m}g(T)\sum _{d\mid T}\mu (d)\ast \frac{T}{d}$$

后面的式子相当于 $\mu \ast id$ ，也就等于 $\phi$ ，就变成

$$\sum _{T=1}^{m}g(T)\phi (T)$$

直接用杜教筛就可以。

然后假如用欧拉反演的话，让 $gcd(x,y)=\sum _{d\mid gcd(x,y)}\phi (x)$

那么就有

$$\begin{array}{rl}& \sum _{i_1=1}^m\sum _{i_2=1}^m\cdots \sum _{i_n=1}^{m}gcd(i_1,i_2,\cdots ,i_n)\\ =& \sum _{i_1=1}^m\sum _{i_2=1}^m\cdots \sum _{i_n=1}^m\sum _{d\mid gcd(i_1,i_2,\cdots ,i_n)}\phi (d)\\ =& \sum _{d=1}^m\phi (d){\lfloor {\displaystyle \frac{m}{d}}\rfloor }^n\end{array}$$

code

//简单的最大公约数
#include<bits/stdc++.h>
#define int unsigned long long
using namespace std;
const int N=3*1e6;
int prime[N+5];
bool nprime[N+5];
int sumphi[N+5],phi[N+5];
unordered_map<int,int> Sumphi;
inline void get_prime(){
    phi[1]=1;
    for(register int i=2;i<=N;++i){
        if(!nprime[i]){
            prime[++prime[0]]=i;
            phi[i]=i-1;
        }
        for(register int j=1;j<=prime[0] && i*prime[j]<=N;j++){
            int k=prime[j]*i;
            nprime[k]=1;
            if(i%prime[j]==0){
                phi[k]=phi[i]*prime[j];
                break;
            }
            phi[k]=phi[i]*phi[prime[j]];
        }
    }
    for(register int i=1;i<=N;++i){
        sumphi[i]=sumphi[i-1]+phi[i];
    }
}
inline int get_phi(int n){
    if(n<=N) return sumphi[n];
    if(Sumphi[n]) return Sumphi[n];
    int ans=0;
    if(n&1)
        ans=(n+1)/2*n;
    else
        ans=n/2*(n+1);
    int l=2;
    while(l<=n){
        int t=n/l;
        int r=n/t;
        ans-=(r-l+1)*get_phi(n/l);
        l=r+1;
    }
    return Sumphi[n]=ans;
}
int qpow(int x,int p){
    int ans=1;
    while(p){
        if(p&1) ans=(ans*x);
        x=(x*x);
        p>>=1;
    }
    return ans;
}
signed main(){
    get_prime();
    int n,m;
    cin>>n>>m;
    int ans=0;
    int l=1;
    while(l<=m){
        int t=m/l;
        int r=min(m/t,m);
        ans=(ans+qpow(t,n)*(get_phi(r)-get_phi(l-1)));
        l=r+1;
    }
    cout<<ans<<endl;
}