「算法笔记」杜教筛

一、前置概念

具体在 「算法笔记」莫比乌斯反演 写过，所以「前置概念」就简单写写。积性函数和完全积性函数就不写了。

狄利克雷卷积：对于两个数论函数 (f,g)，定义它们的狄利克雷卷积 (h=f*g) 为

[displaystyle h(x)=sum_{dmid x}f(d)gleft(frac{x}{d} ight) ]

狄利克雷卷积满足交换律、结合律、对加法的分配律，有单位元 (varepsilon)（其中 (varepsilon) 为单位函数 (varepsilon(x)=[x=1])）。若 (f,g) 是积性函数，则 (f*g) 也是积性函数。

在狄利克雷卷积的意义下，(mu*1=varepsilon)，即 (mu) 和 (1) 互为逆元。

常用卷积：(mu*1=varepsilon)，(varphi*1= ext{Id}Leftrightarrow ext{Id}*mu=varphi)（因为 (mu*1=varepsilon)，所以两边同时卷 (mu) 即可）。

二、杜教筛

1. 基本思想

对于数论函数 (f)，求 (S(n)=sum_{i=1}^nf(i))。

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

[sum_{i=1}^nsum_{dmid i}f(d)gleft(frac{i}{d} ight)=sum_{i=1}^ng(i)Sleft(leftlfloorfrac{n}{i} ight floor ight)Leftrightarrowsum_{i=1}^n(f*g)(i)=sum_{i=1}^ng(i)Sleft(leftlfloorfrac{n}{i} ight floor ight) ]

略证：

[sum_{i=1}^nsum_{dmid i}f(d)gleft(frac{i}{d} ight)=sum_{i=1}^nsum_{dmid i}fleft(frac{i}{d} ight)g(d)=sum_{d=1}^n g(d)sum_{i=1}^{lfloorfrac{n}{d} floor}f(frac{di}{d})=sum_{i=1}^ng(i)Sleft(leftlfloorfrac{n}{i} ight floor ight) ]

则可以得到 (S(n)) 关于 (S(lfloor frac{n}{i} floor)) 的递推式：

[g(1)S(n)=sum_{i=1}^n(f*g)(i)-sum_{i=2}^ng(i)Sleft(leftlfloorfrac{n}{i} ight floor ight) ]

如果我们可以快速求出 (sum_{i=1}^n(f*g)(i))，再用整除分块来求 (sum_{i=2}^ng(i)Sleft(leftlfloorfrac{n}{i} ight floor ight))，就能求得 (g(1)S(n))。也就是说，我们要找到一个合适的 (g)，使得 (f*g) 和 (g) 都能快速求前缀和。

(S(lfloorfrac{n}{i} floor)) 是一个子问题，对 (S(n)) 递归求解并记忆化即可。

线性筛预处理前 (n^{frac{2}{3}}) 个 (f(x)) 的值，则时间复杂度为 (mathcal O(n^{frac{2}{3}}))，证明略。

2. 例子

求 (S(n)=sum_{i=1}^nvarphi(i))。

(varphi*1= ext{Id})。显然 ( ext{Id}(x)=x)、(1) 可以快速求前缀和。取 (g(x)=1) 即可。

(S(n)=sum_{i=1}^n i-sum_{i=2}^nSleft(leftlfloorfrac{n}{i} ight floor ight)=frac{n(n+1)}{2}-sum_{i=2}^nSleft(leftlfloorfrac{n}{i} ight floor ight))。

求 (S(n)=sum_{i=1}^nmu(i))。

(mu*1=varepsilon)。显然 (varepsilon(x)=[x=1])、(1) 可以快速求前缀和。取 (g(x)=1) 即可。

(S(n)=1-sum_{i=2}^nSleft(leftlfloorfrac{n}{i} ight floor ight))。

//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
",S_phi(n),S_mu(n));
    }
    return 0;
} 

三、狄利克雷生成函数

p>有些时候不容易看出 (g(x)) 取什么，比如：求 (S(n)=sum_{i=1}^n icdot varphi(i))。凑是一种方法；而使用 狄利克雷生成函数（DGF） 可以从另一个角度求出 (g(x))。

对于无穷序列 (f_1,f_2,cdots)，定义其狄利克雷生成函数为

[ ilde F(x)=sum_{igeq 1}frac{f_i}{i^x} ]

与普通生成函数以及指数生成函数的对比：

普通生成函数：(F(x)=sum_{igeq 0}f_ix^i)。

指数生成函数：(hat F(x)=sum_{igeq 0}f_ifrac{x^i}{i!})。

对于两个序列 (f,g)，其 DGF 的乘积对应 (f,g) 的狄利克雷卷积序列：

[ ilde F(x) ilde G(x)=sum_{igeq 1}sum_{jgeq 1}frac{f_i}{i^x}frac{g_j}{j^x}=sum_{igeq 1}sum_{jgeq 1}frac{f_ig_j}{(ij)^x}=sum_dsum_{imid d}frac{f_ig_{frac{d}{i}}}{d^x} ]

如果序列 (f) 满足积性（积性函数），其 DGF 可以由质数幂处的取值表示：

[ ilde F=prod_{pin ext{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}^nsum_{j=1}^n icdot jcdot gcd(i,j) ight)mod p ]

(nleq 10^{10},5 imes 10^8leq pleq 1.1 imes 10^9)，(p) 是质数。

Solution：可以用 (varphi*1= ext{Id}) 推式子，即 ( ext{Id}(x)=sum_{dmid x}varphi(d))。

[sum_{i=1}^nsum_{j=1}^n icdot jcdot gcd(i,j)=sum_{i=1}^nsum_{j=1}^n icdot jsum_{dmid gcd(i,j)}varphi(d) =sum_{d}varphi(d)sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{n}{d} floor}icdot jcdot d^2\ =sum_d varphi(d)d^2sum_{i=1}^{lfloorfrac{n}{d} floor}icdot frac{(1+lfloorfrac{n}{d} floor) imes lfloorfrac{n}{d} floor}{2}=sum_d varphi(d)d^2left(frac{lfloorfrac{n}{d} floor}{2} ight)^2 ]

用杜教筛求 (S(n)=sum_{i=1}^nvarphi(i)i^2)，然后整除分块即可。

现在要用杜教筛来筛 (f(x)=varphi(x)x^2)，那么要找到一个合适的 (g)，使得 (f*g) 和 (g) 都可以快速求前缀和。可以用 狄利克雷生成函数（DGF） 推出。取 (g= ext{Id}_2) 即可，(f*g= ext{Id}_3)（幂函数 ( ext{Id}_k(n)=n^k)）。

#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
",(ans+mod)%mod);
    return 0;
}