积性函数前缀和的处理——Min25筛

踩爆洲阁筛的算法终于出现了,向Min25大佬致敬


简介

为了解决一般积性函数前缀和问题,16年集训队论文中提出了洲阁筛

但是洲阁筛思维难度过高,难以理解

后来$Min25$在博客中提出了一种基于埃氏筛的扩展算法,大大降低了思维难度

此法对积性函数的唯一要求是$f(p)$是一个低次多项式,且$f(p^c)$可以快速计算

它的时间复杂度为:$O(\frac{n^{\frac{3}{4}}}{\log{n}})$,空间复杂度:$O(\sqrt{n})$


前置技能:埃氏筛

埃氏筛法每一轮找出第$i$个质数,并把这个质数的倍数筛掉

对所有$\sqrt{n}$以内的质数执行完之后就可以结束

第一次被$P_i$筛掉的$P_i$的倍数$x$必然满足$x\geq P_{i}^2$


第一部分:处理质数点的前缀和

对于每个$\left \lfloor \frac{n}{x} \right \rfloor$,我们先来计算下面的式子:

$$\sum_{i=2}^{\left \lfloor \frac{n}{x} \right \rfloor}[i是质数]f(i)$$

对于质数$p$,$f(p)$可以被分解成若干个$p^k$之和的形式,接下来只需考虑$f(p)=p^k$的情况

我们之所以这样拆分,是为了保证接下来的推理中$f(i)$是完全积性函数

值得一提的是,对于原本的$f(i)$,仅在质数点满足$f(p)=p^k$

但在第一部分的推导中,我们假设所有数的$f(i)$值都满足这个式子

这与埃氏筛的思想是一致的,即筛法开始之前,默认所有数都是质数

$$g(n,j)=\sum_{i=2}^{n}[i是质数或Min_p(i)>P_j]f(i)$$

其中$Min_p(i)$表示$i$的最小质因子

说人话,$g(n,j)$就是执行完第$j$轮埃氏筛后剩下的数的函数值之和

其中包括没被筛掉的合数与处理过的质数

那么我们要求的东西就是

$$\sum_{i=2}^{\left \lfloor \frac{n}{x} \right \rfloor}[i是质数]f(i)=g(n,+\infty)$$

我们尝试从$g(n,j-1)$推出$g(n,j)$的答案

1. 若$P_{j}^2>n$,那么第$j$轮埃氏筛没有筛掉任何数

2. 若$P_{j}^2\leq n$,那么第$j$轮埃氏筛会筛掉

$${x\leq \left\lfloor\frac{n}{P_j}\right\rfloor}且x不含小于P_j的质因子$$

这些$x$的$P_j$倍,而$g(\left\lfloor\frac{n}{P_j}\right\rfloor,j-1)$减去前$j-1$个质数的函数值就是这些$x$的函数值之和,同时我们利用完全积性函数的性质得到$f(xP_j)=f(x)f(P_j)$,提出$f(P_j)$就可以直接计算了

综合一下,我们得到了$g(n,j)$的递推式

$$g(n,j)=\left\{\begin{matrix}g(n,j-1),&P_{j}^2>n\\ g(n,j-1)-f(P_j)[g(\left \lfloor \frac{n}{P_j} \right \rfloor,j-1)-\sum_{i=1}^{j-1}f(P_j)],&P_{j}^2\leq n\end{matrix}\right.
$$

事实上,由于做完$j-1$轮埃氏筛后,区间$[1,P_{j-1}]$中就只剩下质数了,所以

$$\sum_{i=1}^{j-1}f(P_j)=g(P_{j-1},j-1)$$

这样可以简化运算,当然预处理也是可以的

初始值就是还没有执行埃氏筛的状态:

$$g(n,0)=\sum_{i=2}^{n}f(i)$$


第二部分:计算答案

现在我们对于每个$x=\left \lfloor \frac{n}{x} \right \rfloor$求出了$\sum_{i=2}^{\left \lfloor \frac{n}{x} \right \rfloor}[i是质数]f(i)$

$$S(n,j)=\sum_{i=2}^{n}[Min_p(i)\geq P_j]f(i)$$

质数点的答案为($g(n,+\infty)$减去小于$P_j$的质数的贡献)

$$g(n,+\infty)-\sum_{i=2}^{j-1}f(P_i)=g(n,+\infty)-g(P_{j-1},+\infty)$$

接下来考虑合数的贡献,枚举合数的最小质因子及其出现次数即可:

$$S(n,j)=g(n,+\infty)-g(P_{j-1},+\infty)+\sum_{k=j}^{P_{k}^2\leq n}\sum_{e=1}^{P_{k}^{e+1}\leq n}S(\left \lfloor \frac{n}{P_{k}^e} \right \rfloor,k+1)f(P_{k}^e)+f(P_{k}^{e+1})$$

$S(*,k+1)$保证了没有小于等于$P_k$的质因子,这样的话质数幂的情况需要单独考虑

最后的答案显然就是$S(n,1)+f(1)$


代码实现

我们预处理所有$x=\left \lfloor \frac{n}{i} \right \rfloor$的值,从小到大记为$Pos_1,Pos_2,\cdots,Pos_{m}$

1. 对于$1\leq x\leq\sqrt{n}$,显然$x=Pos_x$

2. 对于$x>\sqrt{n}$,观察下面的值:

$$\left\{\begin{matrix}Pos_m=n\\Pos_{m-1}=\left \lfloor \frac{n}{2} \right \rfloor\\ \cdots\\Pos_{\sqrt{n}}=\sqrt{n}\end{matrix}\right.$$

我们发现$Pos_x=\left \lfloor \frac{n}{m-x+1} \right \rfloor$,所以$x=m-\left \lfloor \frac{n}{Pos_x} \right \rfloor+1$

这样的话从值到编号就可以$O(1)$转化了

事实上,代码实现上还是有一些技巧的

第一部分预处理时,我们枚举质数$i$,然后分层从大到小递推

第二部分直接递归处理即可


【LOJ6235】区间素数个数

第一部分中,取$f(p)=1$即可,求出的$g(n,+\infty)$就是答案

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 316300
using namespace std;
typedef long long ll;
ll n,N,cnt,pos[MAXN<<1],g[MAXN<<1],prime[MAXN/10];
inline int Id(ll x){return x<=N?x:pos[0]-n/x+1;}
int main(){
    freopen(FILE".in","r",stdin);
    freopen(FILE".out","w",stdout);
    scanf("%lld",&n); N=sqrt(n);
    for(ll i=1;i<=n;i=pos[pos[0]]+1) pos[++pos[0]]=n/(n/i);
    for(int i=1;i<=pos[0];++i) g[i]=pos[i]-1;
    for(int i=2;i<=N;++i)if(g[i]!=g[i-1]){
        prime[++cnt]=i;
        for(int j=pos[0];1LL*i*i<=pos[j];--j) 
            g[j]-=g[Id(pos[j]/i)]-(cnt-1);
    }
    printf("%lld\n",g[pos[0]]);
    return 0;
}

 


【Luogu4213】杜教筛

我们用Min25筛去水杜教筛模板题,跑得挺快的

按上面的过程走就对了

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 46345
using namespace std;
typedef long long ll;
int n,N,cnt,pos[MAXN<<1];ll g[MAXN<<1],h[MAXN<<1],sum[MAXN/10],prime[MAXN/10];
inline int Id(int x){return x<=N?x:pos[0]-n/x+1;}
inline int read(){
    int x=0,f=1; char ch=getchar();
    while(ch>'9'||ch<'0') {if(ch=='-') f=-1; ch=getchar();}
    while(ch>='0'&&ch<='9') {x=x*10+ch-'0'; ch=getchar();}
    return x*f;
}
inline ll Sum_Phi(int n,int j){
    if(n<prime[j]) return 0;
    ll Ans=(h[Id(n)]-sum[j-1])-(g[Id(n)]-(j-1));
    for(int k=j;k<=cnt&&1LL*prime[k]*prime[k]<=n;++k){
        int fac=prime[k];
        for(int e=1;fac*prime[k]<=n;++e){
            Ans+=(Sum_Phi(n/fac,k+1)*(fac/prime[k])+fac)*(prime[k]-1);
            fac*=prime[k];
        }
    }return Ans;
}
inline ll Sum_Mu(int n,int j){
    if(n<prime[j]) return 0;
    ll Ans=(j-1)-g[Id(n)];
    for(int k=j;k<=cnt&&1LL*prime[k]*prime[k]<=n;++k)
        Ans-=Sum_Mu(n/prime[k],k+1);
    return Ans;
}
void solve(){
    n=read(); N=sqrt(n); pos[0]=cnt=0;
    for(int i=1;i<=n;i=pos[pos[0]]+1) pos[++pos[0]]=n/(n/i);
    for(int i=1;i<=pos[0];++i) g[i]=1LL*pos[i]-1;
    for(int i=1;i<=pos[0];++i) h[i]=1LL*pos[i]*(pos[i]+1)/2-1;
    for(int i=2;i<=N;++i)if(g[i]!=g[i-1]){
        prime[++cnt]=i; sum[cnt]=sum[cnt-1]+i;
        for(int j=pos[0];1LL*i*i<=pos[j];--j){
            g[j]=g[j]-(g[Id(pos[j]/i)]-(cnt-1));
            h[j]=h[j]-(h[Id(pos[j]/i)]-sum[cnt-1])*i;
        }
    }
    printf("%lld %lld\n",Sum_Phi(n,1)+1,Sum_Mu(n,1)+1);
}
int main(){
    freopen(FILE".in","r",stdin);
    freopen(FILE".out","w",stdout);
    int T=read(); while(T--) solve();
    return 0;
}

 


【LOJ6053】简单的函数

对于质数$P$的值

$$f(P)=\left\{\begin{matrix}P+1,P=2\\P-1,P\neq 2\end{matrix}\right.$$

2的情况单独处理,其余的就是套路了

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 100010
#define cadd(a,b) a=add(a,b)
#define csub(a,b) a=sub(a,b)
#define cmul(a,b) a=mul(a,b)
using namespace std;
typedef long long ll;
const int mod=1e9+7;
int N,cnt,g[MAXN<<1],h[MAXN<<1];ll n,pos[MAXN<<1],prime[MAXN/10];
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline int sub(int a,int b){return (a-=b)<0?a+mod:a;}
inline int mul(int a,int b){return 1LL*a*b%mod;}
inline int Pow(int a,int b){
    int ret=1;
    while(b){
        if(b&1) cmul(ret,a);
        cmul(a,a); b>>=1;
    }return ret;
}
inline int Inv(int x){return Pow(x,mod-2);}
inline int Id(ll x){return x<=N?x:pos[0]-n/x+1;}
inline int S(ll n,int j){
    if(n<prime[j]) return 0;
    int Ans=sub(h[Id(n)],h[prime[j-1]]);
    for(int k=j;k<=cnt&&prime[k]*prime[k]<=n;++k){
        ll fac=prime[k];
        for(int e=1;fac*prime[k]<=n;++e){
            cadd(Ans,add(mul(S(n/fac,k+1),prime[k]^e),prime[k]^(e+1)));
            fac*=prime[k];
        }
    }return Ans;
}
int main(){
    freopen(FILE".in","r",stdin);
    freopen(FILE".out","w",stdout);
    scanf("%lld",&n); N=sqrt(n);
    for(ll i=1;i<=n;i=pos[pos[0]]+1) pos[++pos[0]]=n/(n/i);
    for(int i=1;i<=pos[0];++i) g[i]=sub(pos[i]%mod,1);
    for(int i=1;i<=pos[0];++i) h[i]=mul(mul((pos[i]+2)%mod,(pos[i]-1)%mod),Inv(2));
    for(int i=2;i<=N;++i)if(g[i]!=g[i-1]){
        prime[++cnt]=i;
        for(int j=pos[0];1LL*i*i<=pos[j];--j){
            csub(g[j],sub(g[Id(pos[j]/i)],cnt-1));
            csub(h[j],mul(sub(h[Id(pos[j]/i)],h[prime[cnt-1]]),i));
        }
    }
    for(int i=1;i<=pos[0];++i) h[i]+=(i>1)*2-g[i];
    printf("%d\n",add(S(n,1),1));
    return 0;
}

 


【SPOJ DIVCNTK】Counting Divisors(general)

注意到:

$$f(P)=\sigma_0(P^k)=k+1$$

$$f(P^e)=\sigma_0(P^{ek})=ek+1$$

然后就无脑Min25筛了

至于对$2^{64}$取模就相当于开个unsigned long long让它自然溢出即可

顺便可以水掉下面的东西(四倍经验岂不美哉)

「SPOJ DIVCNT1」Counting Divisors

「SPOJ DIVCNT2」Counting Divisors (square)

「SPOJ DIVCNT3」Counting Divisors (cube)

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 100010
using namespace std;
typedef unsigned long long ll;
ll n,N,K,cnt,prime[MAXN/10],pos[MAXN<<1],g[MAXN<<1];
inline ll read(){
    ll x=0,f=1; char ch=getchar();
    while(ch>'9'||ch<'0') {if(ch=='-') f=-1; ch=getchar();}
    while(ch>='0'&&ch<='9') {x=x*10+ch-'0'; ch=getchar();}
    return x*f;
}
inline int Id(ll x){return x<=N?x:pos[0]-n/x+1;}
inline ll S(ll n,ll j){
    if(n<prime[j]) return 0;
    ll Ans=g[Id(n)]-g[prime[j-1]];
    for(int k=j;k<=cnt&&1LL*prime[k]*prime[k]<=n;++k){
        ll fac=prime[k];
        for(int e=1;1LL*fac*prime[k]<=n;++e){
            Ans+=S(n/fac,k+1)*(e*K+1)+(e+1)*K+1;
            fac=fac*prime[k];
        }
    }
    return Ans;
}
void solve(){
    n=read(); K=read(); N=sqrt(n); cnt=pos[0]=0;
    for(ll i=1;i<=n;i=pos[pos[0]]+1) pos[++pos[0]]=n/(n/i);
    for(int i=1;i<=pos[0];++i) g[i]=pos[i]-1;
    for(int i=2;i<=N;++i)if(g[i]!=g[i-1]){
        prime[++cnt]=i;
        for(int j=pos[0];1LL*i*i<=pos[j];--j)
            g[j]=g[j]-(g[Id(pos[j]/i)]-(cnt-1));
    }
    for(int i=1;i<=pos[0];++i) g[i]=g[i]*(K+1);
    printf("%llu\n",S(n,1)+1);
}
int main(){
    freopen(FILE".in","r",stdin);
    freopen(FILE".out","w",stdout);
    int T=read(); while(T--) solve();
    return 0;
}

 


【SPOJ APS2】Amazing Prime Sequence(hard)

我们知道对于质数的情况$f(p)=p$,接下来考虑合数的答案

事实上,埃氏筛每一轮中筛掉的数都是被最小质因子筛掉的,对答案的贡献是当前质数

我们只需要统计每一轮被筛掉的数个数,乘上当前质数就是贡献

这个东西在Min25筛第一部分可以完成

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 1200000
using namespace std;
typedef unsigned long long ll;
ll n,N,Ans,cnt,prime[MAXN/10],g[MAXN<<1],h[MAXN<<1],pos[MAXN<<1];
inline int Id(ll x){return x<=N?x:pos[0]-n/x+1;}
inline ll read(){
    ll x=0,f=1; char ch=getchar();
    while(ch>'9'||ch<'0') {if(ch=='-') f=-1; ch=getchar();}
    while(ch>='0'&&ch<='9') {x=x*10+ch-'0'; ch=getchar();}
    return x*f;
}
void solve(){
    n=read(); N=sqrt(n); Ans=cnt=pos[0]=0;
    for(ll i=1;i<=n;i=pos[pos[0]]+1) pos[++pos[0]]=n/(n/i);
    for(int i=1;i<=pos[0];++i) g[i]=pos[i]-1;
    for(int i=1;i<=pos[0];++i) h[i]=(pos[i]&1)?pos[i]*((pos[i]+1)/2)-1:pos[i]/2*(pos[i]+1)-1;
    for(int i=2;i<=N;++i)if(g[i]!=g[i-1]){ 
        prime[++cnt]=i;
        Ans+=(g[Id(n/i)]-(cnt-1))*i;//统计合数贡献
        for(int j=pos[0];1LL*i*i<=pos[j];--j){
            g[j]=g[j]-(g[Id(pos[j]/i)]-g[prime[cnt-1]]);
            h[j]=h[j]-(h[Id(pos[j]/i)]-h[prime[cnt-1]])*i;
        }
    }
    printf("%llu\n",Ans+h[pos[0]]);
}
int main(){
    freopen(FILE".in","r",stdin);
    freopen(FILE".out","w",stdout);
    int T=read(); while(T--) solve();
    return 0;
}

 


参考文献

《Min_25筛学习笔记》——CMXRYNP

《Min25筛》——租酥雨

《一些特殊的数论函数求和问题》——朱震霆

此条目发表在数论, 算法分类目录。将固定链接加入收藏夹。

发表评论

电子邮件地址不会被公开。 必填项已用*标注