<学习笔记> 杜教筛

杜教筛

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

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

\[\sum_{i=1}^{n}(f*g)(i)=\sum_{i=1}^{n} \sum_{d \mid i} g(d)*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*g)(i)-\sum_{i=2}^{n} g(i)S(\lfloor \frac{n}{i} \rfloor) \]

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

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

例如根据 \(I(n)=1\)，\(id(n)=n\)，\(\epsilon(n)=[n=1]\)

单位元 \(\epsilon\) 满足 \(f * \epsilon= f\)

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

\(\mu * I =\epsilon\)
\(\varphi * I = id\)
\(\mu * id= \varphi\)

对于一，我们直接展开就是 \((\mu*I)(n)=\sum_{d \mid n} \mu(d)=\epsilon\)

对于二，还是展开 \((\varphi * I)(n)=\sum_{d \mid n} \varphi(d)\)，然后可以证明 \(\sum_{d \mid n} \varphi(d)=n\)。

证明

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

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

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

一些题

DZY Loves Math IV

\[求\sum_{i=1}^{n} \sum_{j=1}^{m} \varphi (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}\)，那么 \(\varphi(n)=q \varphi(p)\)，也就有 \(\varphi(ni)=q \varphi(pi)\)。以下有好几条性质是由于 \(p\) 的质因子系数为一才会成立。

我们另

\[S(n,m)= q \sum_{i=1}^{m} \varphi(p * j) \]

\[= q \sum_{i=1}^{m} \varphi(\frac{p}{\gcd(p,i)}\gcd(p,i)*i) \]

\[= q \sum_{i=1}^{m} \varphi(\frac{p}{\gcd(p,i)}) * \varphi(\gcd(p,i)*i) \]

\[= q \sum_{i=1}^{m} \varphi(\frac{p}{\gcd(p,i)}) * \varphi(i)* \gcd(p,i) \]

\[= q \sum_{i=1}^{m} \varphi(\frac{p}{\gcd(p,i)}) * \varphi(i)* \sum_{d \mid \gcd(p,i)} \varphi(d) \]

\[= q \sum_{i=1}^{m} \varphi(i)* \sum_{d \mid \gcd(p,i)} \varphi(\frac{p*d}{\gcd(p,i)}) \]

\[= q \sum_{i=1}^{m} \varphi(i)* \sum_{d \mid \gcd(p,i)} \varphi(\frac{p}{d}) \]

\[= q \sum_{i=1}^{m} \varphi(i)* \sum_{d \mid p,d \mid i} \varphi(\frac{p}{d}) \]

尝试枚举 \(d\)

\[= q \sum_{d \mid p} \varphi(\frac{p}{d})* \sum_{i=1}^{\lfloor \frac{m}{d} \rfloor} \varphi(i*d) \]

\[= q \sum_{d \mid p} \varphi(\frac{p}{d})* 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\) 级别。

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

有一个结论：

\[\varphi(i*j)=\frac{\varphi(i)*\varphi(j)*gcd(i,j)}{\varphi(gcd(i,j))} \]

证明

\[\varphi(i)\varphi(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)\) 那么就可以变成 \(\varphi(ij) *\varphi(\gcd(i,j))\)，得证。

那么我们那它推式子

\[\sum_{i=1}^{n} \sum_{i=1}^{m} \frac{\varphi(i)*\varphi(j)*gcd(i,j)}{\varphi(gcd(i,j))} \]

\[\sum_{d=1}^{n} \sum_{i=1}^{n} \sum_{i=1}^{m} \frac{\varphi(i)*\varphi(j)*d}{\varphi(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{\varphi(i)*\varphi(j)*d}{\varphi(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{\varphi(d*i)*\varphi(d*j)*d}{\varphi(d)} \sum_{k \mid i,k \mid j} \mu(k) \]

变为枚举 \(k\)

\[\sum_{d=1}^{n} \frac{d}{\varphi(d)} \sum_{k=1}^{\lfloor \frac{n}{d} \rfloor } \mu(k) \sum_{i=1}^{\lfloor \frac{n}{kd} \rfloor } \varphi(d*k*i) \sum_{i=1}^{\lfloor \frac{m}{kd} \rfloor } \varphi(d*k*j) \]

另 \(T=dk\)

\[\sum_{T=1}^{n} \sum_{d \mid T}\mu(\frac{T}{d}) \frac{d}{\varphi(d)}\sum_{i=1}^{\lfloor \frac{n}{T} \rfloor } \varphi(T*i) \sum_{i=1}^{\lfloor \frac{m}{T} \rfloor } \varphi(T*j) \]

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

另 \(G(x,y)=\sum_{i=1}^{x} \varphi(i*y)\)，并且存在递推 \(G(x,y)=G(x-1,y)+\varphi(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)*G(y,x)*G(z,x)\)

我们定义阈值为 \(B\)，首先我们预处理出来 \(G(x,y) (x \leq N, y \leq B)\)， \(Ans(x,y,z) (x \leq N, y \leq B,z \leq B)\)，对于 \(T \leq 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*i)\mu(i)= \lfloor \frac{m}{x}\rfloor ^ n\)

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

\[\sum_{T=1}^{m} g(T) \sum_{d \mid T} \mu(d) * \frac{T}{d} \]

后面的式子相当于 \(\mu * id\) ，也就等于 \(\varphi\) ，就变成

\[\sum_{T=1}^{m} g(T) \varphi(T) \]

直接用杜教筛就可以。

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

那么就有

\[\begin{aligned} &\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)} \varphi(d)\\ =&\sum_{d=1}^m \varphi(d) \left\lfloor\dfrac{m}{d}\right\rfloor^n \end{aligned}\]

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

posted @ 2024-01-21 19:45 _bloss 阅读(37) 评论(0) 收藏举报

刷新页面返回顶部

bloss

<学习笔记> 杜教筛

杜教筛

一些题

DZY Loves Math IV

Pyh 的求/毒瘤之神的考验

简单的最大公约数

公告