<学习笔记> 杜教筛

杜教筛

处理数论函数的前缀和问题,可以在低于线性的复杂度里求出 S(n)=i=1nf(i)

对于任意一个数论函数 g,必须满足 :

i=1n(fg)(i)=i=1ndig(d)d(id)

=d=1ng(d)j=1ndf(j)=d=1ng(d)S(nd)

那么可以得到:

g(1)S(n)=i=1ng(i)S(ni)i=2ng(i)S(ni)

=i=1n(fg)(i)i=2ng(i)S(ni)

时间复杂度最小为 n23

一般情况下,可以设出方便求前缀的 gfg,然后进行递归求解。

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

单位元 ϵ 满足 fϵ=f

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

  1. μI=ϵ

  2. φI=id

  3. μid=φ

对于一,我们直接展开就是 (μI)(n)=dnμ(d)=ϵ

对于二,还是展开 (φI)(n)=dnφ(d),然后可以证明 dnφ(d)=n

证明

f(x) 表示 gcd(k,n)=x,(k<n) 的个数,那么有 n=i=1nf(i)

我们还知道 f(x)=φ(nx) 相当于就是 gcd(k,n)=x 的个数相当于 gcd(kx,nx)=1 的个数,那么就有 n=inφ(i)

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

一些题#

DZY Loves Math IV #

i=1nj=1mφ(ij)

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

我们设 n=i=1w(n)pikip=i=1w(n)piq=np,那么 φ(n)=qφ(p),也就有 φ(ni)=qφ(pi)。以下有好几条性质是由于 p 的质因子系数为一才会成立。

我们另

S(n,m)=qi=1mφ(pj)

=qi=1mφ(pgcd(p,i)gcd(p,i)i)

=qi=1mφ(pgcd(p,i))φ(gcd(p,i)i)

=qi=1mφ(pgcd(p,i))φ(i)gcd(p,i)

=qi=1mφ(pgcd(p,i))φ(i)dgcd(p,i)φ(d)

=qi=1mφ(i)dgcd(p,i)φ(pdgcd(p,i))

=qi=1mφ(i)dgcd(p,i)φ(pd)

=qi=1mφ(i)dp,diφ(pd)

=qi=1mφ(i)dp,diφ(pd)

尝试枚举 d

=qdpφ(pd)i=1mdφ(id)

=qdpφ(pd)S(d,md)

n 等于 1 时我们直接暴力求 m23,我们可以对 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 级别。

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

有一个结论:

φ(ij)=φ(i)φ(j)gcd(i,j)φ(gcd(i,j))

证明

φ(i)φ(j)=ipip1pjpjp1p

=ijpip1ppjp1p

=ijpijp1ppgcd(i,j)p1p

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

那么我们那它推式子

i=1ni=1mφ(i)φ(j)gcd(i,j)φ(gcd(i,j))

d=1ni=1ni=1mφ(i)φ(j)dφ(d)[gcd(i,j)=d]

d=1ni=1ndi=1mdφ(i)φ(j)dφ(d)[gcd(i,j)=1]

d=1ni=1ndi=1mdφ(di)φ(dj)dφ(d)ki,kjμ(k)

变为枚举 k

d=1ndφ(d)k=1ndμ(k)i=1nkdφ(dki)i=1mkdφ(dkj)

T=dk

T=1ndTμ(Td)dφ(d)i=1nTφ(Ti)i=1mTφ(Tj)

我们可以套路的另 F(T)=dTμ(Td)dφ(d),这个也是很好预处理的。

G(x,y)=i=1xφ(iy),并且存在递推 G(x,y)=G(x1,y)+φ(xy)

那么答案就变成了 Ans=T=1nF(T)G(nT,T)G(mT,T)

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

我们定义阈值为 B,首先我们预处理出来 G(x,y)(xN,yB)Ans(x,y,z)(xN,yB,zB),对于 TB 的我们直接暴力扫一遍,对于 T>B 的,我们进行整除分块,直接查询 Ans(r,nT,mT)Ans(l1,nT,mT) 就可以。

复杂度:预处理的话是 nlnn+nB2,单次查询的话 mB+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)gcdx 的个数,设 g(x) 表示 gcdx 倍数的个数,那么就有 f(x)=i=1mxg(xi)μ(i)=mxn

那么答案就是 x=1mxf(x),展开后另 T=xd 就有

T=1mg(T)dTμ(d)Td

后面的式子相当于 μid ,也就等于 φ ,就变成

T=1mg(T)φ(T)

直接用杜教筛就可以。

然后假如用欧拉反演的话,让 gcd(x,y)=dgcd(x,y)φ(x)

那么就有

i1=1mi2=1min=1mgcd(i1,i2,,in)=i1=1mi2=1min=1mdgcd(i1,i2,,in)φ(d)=d=1mφ(d)mdn

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;
}

作者:bloss

出处:https://www.cnblogs.com/jinjiaqioi/p/17978192

版权:本作品采用「署名-非商业性使用-相同方式共享 4.0 国际」许可协议进行许可。

posted @   _bloss  阅读(28)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 单线程的Redis速度为什么快?
· 展开说说关于C#中ORM框架的用法!
· Pantheons:用 TypeScript 打造主流大模型对话的一站式集成库
· SQL Server 2025 AI相关能力初探
· 为什么 退出登录 或 修改密码 无法使 token 失效
点击右上角即可分享
微信分享提示
more_horiz
keyboard_arrow_up dark_mode palette
选择主题
menu