min25筛学习笔记&模板详解

min25筛算法推导&模板详解

问题引入

引入函数${\sigma }_0(n)$为n的正因数的数量，求$S(n,k)={\sum }_{i=1}^n{\sigma }_0(i^k)$,$n,k\le 10^{10}$

简化分析与引入

对于${\sigma }_0(n)$函数，显然有：
$$当p为素数，n为p^k{\sigma }_0(n)=k+1$$
且显然${\sigma }_0$为积性函数

当$n\le 10^7$时可以通过

#include<bits/stdc++.h> 
const int size=1e7+5;
bool prime[size];
int p[size];
int sigma[size];
int tot=0;
int k;
void init()
{
	sigma[1]=1;
    for(int i=2;i<size;i++) prime[i]=true;
    for(int i=2;i<size;i++)
    {
        if(prime[i]) 
        {
            p[tot++]=i;
            sigma[i]=k+1;
        }
        for(int j=0;j<tot&&p[j]*i<size;j++)
        {
            prime[i*p[j]]=false;
            if(i%p[j]==0)
            {
                sigma[i*p[j]]=sigma[i]*2-sigma[i/p[j]];
                break;
            }
            else 
            {
                sigma[i*p[j]]=sigma[i]*(k+1);
            }
        }
    }
}
int main()
{
	int n;
	scanf("%d%d",&n,&k);
	init();
	int ans=0;
	for(int i=1;i<=n;i++)
	{
		ans+=sigma[i];
	}
	printf("%d\n",ans);
}

这样对于每个n，k我们都可以线性求解出答案。但是当n最大可以达到$10^{10}$时，这种做法显然是不行的为此就需要引入一种亚线性筛：“min25筛”

虽然min25筛乃至杜教筛名字上都带有一个筛字但是实际上并不能称作一种筛，准确地说，其应该是一种可以亚线性地复杂度下求出积性函数前缀和的算法。其中min25筛求解积性函数前缀和时需要满足条件

• $f(x)$当x为质数时有一个多项式表示
• $f(x^c)$在x为质数时可以快速计算

这两条条件都将在下面的推导种起到重要的作用。

前置知识

埃拉托斯特尼筛法

在线性求解一些数论问题时常常需要用到线性筛，但是有事也可以使用复杂度稍差一点，但是依然有效的埃拉托斯特尼筛法（下面简称埃氏筛法）。其主要的想法是每筛出来一个素数，就将其倍数也删去，这样不断向前递推，未被筛去的就是素数

bool prime[N];
void init(){
    for(int i=2;i<N;i++) prime[i]=true;
    for(int i=2;i<N;i++){
        if(prime[i]){
            for(int j=2*i;j<N;j+=i){
                prime[j]=false;
            }
        }
    }
}

当然，min25筛其实并不用到埃氏筛，但是其借鉴了一点埃氏筛的思想

引理1

对于给定的正整数n和所有的正整数，$\lfloor \frac{n}{m}\rfloor$的取值只有$O(\sqrt{n})$种

证明. 当$m\le \sqrt{n}$时满足条件的m只有$O(\sqrt{n})$种。当$m>\sqrt{n}$时，$0\le \lfloor \frac{n}{m}\rfloor <\sqrt{n}$满足这一条件的$\lfloor \frac{n}{m}\rfloor$也只有$O(\sqrt{n})$种，故得证

min25筛的分析与推导

我们将问题泛化，令$f(x)={\sigma }_0(x^k)$，由于${\sigma }_0(x)$函数为积性函数，易证*f(x)*也为积性函数。且因为$f(x^c)=c\ast k+1$，因此积性函数$f(x)$满足上面所提到的min25筛求解积性函数前缀和的条件。现要求积性函数$f(x)$的前缀和。

为了方便公式的书写，我们需要先引入一些定义
$$\begin{gathered} P:素数集合 \\ lpf(n):n的最小质因数 \\ {P}_i:埃氏筛法筛了前i个质数之后剩下的自然数集合 \\ {p}_i:第i个质数 \end{gathered}$$

将原前缀和的求和方式进行转化
$$\sum _{i=1}^{n}f(i)=1+\sum _{2\le p\le np\in P}\sum _{2\le x\le nlpf(x)=p}f(x)$$
由于*f(x)*为积性函数，故
$$1+\sum _{2\le p\le np\in P}\sum _{2\le x\le nlpf(x)=p}f(x)=1+\sum _{2\le p^e\le n,e\ge 1,p\in P}f(p^e)\left(1+\sum _{2\le x\le \lfloor \frac{n}{p^e}\rfloor ,lpf(x)>p}f(x)\right)$$
显然，对于一个合数来说，其最小的质数定然不会超过$\sqrt{n}$,因此上式就可以进行进一步的转化
$$\begin{array}{rl}\sum _{i=1}^{n}f(i)& =1+\sum _{2\le p\le np\in P}\sum _{2\le x\le nlpf(x)=p}f(x)\\ & =1+\sum _{2\le p^e\le n,e\ge 1,p\in P}f(p^e)\left(1+\sum _{2\le x\le \lfloor \frac{n}{p^e}\rfloor ,lpf(x)>p}f(x)\right)\\ & =1+\sum _{2\le p\le \sqrt{n}2\le p^e\le n,e\ge 1,p\in P}f(p^e)\left(1+\sum _{2\le x\le \lfloor \frac{n}{p^e}\rfloor ,lpf(x)>p}f(x)\right)+\sum _{\sqrt{n}<p\le n,p\in P}f(p)\end{array}$$
但是直接求这个有一定的难度，为此引入两个新的函数
$$\begin{gathered} g(n,m)=\sum _{2\le x\le n,lpf(x)>m}f(x) \\ h(n)=\sum _{2\le p\le n,p\in P}f(p) \end{gathered}$$
原式即是$g(n,0)$

为此我们对g(n,m)进行推导
$$\begin{array}{rl}g(n,m)& =\sum _{2\le x\le n,lpf(x)>m}f(x)\\ & =\sum _{m<p\le \sqrt{n},p^e\le n,e\ge 1,p\in P}f(p^e)\left(1+\sum _{2\le x\le \lfloor \frac{n}{p^e}\rfloor ,lpf(x)>p}f(x)\right)+\sum _{\sqrt{n}<p\le n,p\in P}f(p)\\ & =\sum _{m<p\le \sqrt{n},p^e\le n,e\ge 1,p\in P}f(p^e)\left(1+g(\lfloor \frac{n}{p^e},p\rfloor )\right)+h(n)-h(\sqrt{n})\end{array}$$
但是这样不便于式子之间的转化，因此将原式中的素数部分的函数值的递推转化一下
$$\begin{array}{rl}g(n,m)& =\sum _{m<p\le \sqrt{n},p^e\le n,e\ge 1,p\in P}f(p^e)\left([e>1]+g(\lfloor \frac{n}{p^e},p\rfloor )\right)+\sum _{m<p\le \sqrt{n},p\in P}f(x)+h(n)-h(\sqrt{n})\\ & =\sum _{m<p\le \sqrt{n},p^e\le n,e\ge 1,p\in P}f(p^e)\left([e>1]+g(\lfloor \frac{n}{p^e},p\rfloor )\right)+h(n)-h(m)\end{array}$$
对于我们想要得到的$g(n,0)$我们需要枚举所有的小于等于$\sqrt{n}$的素数作为m并将其对应的$g(n,m)$加上才可以最后得到答案。需要注意的是当$n\le m^2$时，$g(n,m)=0$这是显然可以得知的，当$n\le m^2$时，1到n之间，不存在最小质因数还大于m的数字。

由于前面的条件$f(p^e)$是可以快速计算的，因此接下来只要预处理出h就可以进行递推了。

接下来考虑如何求$h$

首先可以发现所有需要用到的h的标号均为0到$\sqrt{n}$之间的质数，或为$\lfloor \frac{n}{p^e}\rfloor$

现在证明：对于所有需要用到的$h(i)$其一定满足，存在m使得$\lfloor \frac{n}{m}\rfloor =i$。

$$\begin{array}{rl}证：& 令\lfloor \frac{n}{p}\rfloor =k，p为任意的小于等于n的数\\ & 则有n=p\ast k+b(b<p)\\ & n÷k=p\dots \dots b\\ & \because p\le \sqrt{n}\\ & \therefore k\ge \sqrt{n}\\ & \therefore k>b\\ & \therefore \lfloor \frac{n}{k}\rfloor =p\end{array}$$
由此得出，所有小于等于$\sqrt{n}$的数字都是满足上述的可以通过$\lfloor \frac{n}{m}\rfloor$得到的，而所有0到$\sqrt{n}$之间的数字自然也就满足这个条件，而$\lfloor \frac{n}{p^e}\rfloor$则是也显然满足这个条件，因此得证

根据引理1，满足这种条件的数字共有$O(\sqrt{n})$个.由于h函数是对f函数求素数前缀和，且上面说到了，所有满足可以用min25筛求的积性函数都需要满足$f(p)$可以通过多项式表示，也即f函数可以表示为$f(p)=\sum a_t{p}^{b_t}$，不妨令$f_t(p)=a_t{p}^{b_t}$,$f_t$函数的前缀和为$h_t(i)$因此h函数也就可以表示为$h(i)={\sum }_t{h}_t(i)$。进而，对于每个h(i)我们只需要求出所有的$h_t(i)$即可。

到了这里我们的问题就变成了：给定n和t，对所有满足$i=\lfloor \frac{n}{m}\rfloor$，求：
$$h_t(i)=\sum _{p\le i,p\in P}{p}^{b_t}$$
为了求这个函数我们需要引入一个新的函数$h_{t,i}^{\prime }(j)$,其含义为埃氏筛法筛出前i个数之后，前j个数字种剩余的数字的$b_t$次方和
$$\begin{gathered} h_{t,i}^{\prime }(j)=\sum _{1\le p\le j \\ p=1orp\in Porlpf(p)>p_i}{p}^{b_t} \end{gathered}$$
那么函数$h_t(j)=h_{t,i}(j)(j<p_{i+1}\ast p_{i+1})$

对于$h_{t,i}^{\prime }(j)$有转移式
$$h_{t,i}(j)=h_{t,0}^{\prime }(j)-\sum _{k=0}^{k<i}{p}_{k+1}^{b_t}(h_{t,k}^{\prime }(\lfloor \frac{j}{p_{k+1}}\rfloor )-h_{t,k}^{\prime }(p_k))$$
由此，h函数的求值就结束了。

由此就可以得到积性函数f的前缀和了。上述算法的复杂度为$O(\frac{n^{\frac{3}{4}}}{log(n)})$,现在还有一种新版的min25筛可以优化到$O(n^{\frac{2}{3}})$但是实际应用上，由于常数巨大，实际上的运行时间并不优化多少，且代码冗长，故此处不予介绍。

模板及解释

//问题描述定义：
//定义sigma(n)为n的正因数的数量，求f(n,k)=sum_{i=1}^n (sigma(i^k))
#include<stdio.h>
#include<math.h>

using namespace std;
typedef unsigned long long ull;
typedef long long ll;

const int N=200005;

int T,S,pr[N],pc;
ll n,num[N],m,K;
ull g[N];
bool fl[N];
// 给定一个数字X求出其为第几个可以得到有效的g的数字
inline int ID(ll x){return x<=S?x:m-n/x+1;}

ull f(ll n,int i){
    if(n<1||pr[i]>n)return 0;
    ull ret=g[ID(n)]-(i-1)*(K+1);
    while((ll)pr[i]*pr[i]<=n){
        int p=pr[i];
        ull e=K+1,t=n/p;
        while(t>=p)ret+=f(t,i+1)*e+e+K,t/=p,e+=K;
        // ret= sum{sigma(p^es)([n>1]+f(n/p^es,p)+g[n]-g[num[i-1]]}
        // 因为对于函数g(n,m)当n<=m^2,g(n,m)为0
        // 且根据sigma函数的性质，对于质数p，有sigma(p^es)=es*k+1
        // 所以 ret+=(es*k+1)*f(n/p^es,p)+((es+1)k+1)(1<=es&& n/p^es>p)即当前项的sigma(p^e)*f(n/p^e,p)加上下一项的sigma(p^e)
        // 这样的做法避免了e<1也就不必进行特判
        i++;
    }
    return ret;
}
//g[i]即小于i的所有质数的sigma函数求和
//根绝前面的推导，只有当i可以通过[n/m]得到时，其才有用
ull solve(ll n){
    int i,p,x;ull y;
    S=sqrt(n);
    while((ll)S*S<=n)S++;
    while((ll)S*S>n)S--;
    while(m)num[m--]=0;
    for(i=1;i<=S;i++)num[++m]=i;
    for(i=S;i>=1;i--)if(n/i>S)num[++m]=n/i;
    for(i=1;i<=m;i++)g[i]=num[i]-1;
    //此处g[i]为小于等于第i个满足可以通过[n/m]得到的数的素因子的数量
    //故先减去“1”因为1一定不为素因子且无法被后续操作筛去
    x=1;y=0;
    for(p=2;p<=S;p++)if(g[p]!=g[p-1]){
        while(num[x]<(ll)p*p)x++;
        //令g'(i,j)为埃氏筛法筛出前j个素数后，1到j之间剩余数的数量
        //则有g'(i,j)=g'(0,j)-sum_{k=0,k<i}(g'(k,[j/p[k+1]])-g'(k,p[k]))
        //其中p[0]=1,p[i](i>0)为第i个质数
        //则g(j)即为g'(i,j)使得有p[i]<j&&p[i+1]>j
        //其中g'(k,p[k])即为k+1
        //由于以上的递推式仅仅作用到num[x]>=p^2为止，故仅仅将更新到x
        //由于每次都要减去尚未去除当前素数的g以充当g'因此从大往小更新
        for(i=m;i>=x;i--)g[i]-=g[ID(num[i]/p)]-y;
        y++;
    }
    for(i=1;i<=m;i++)g[i]*=K+1;
    return f(n,1)+1;
}

int main(){
    int i,j;
    //素数筛
    for(i=2;i<N;i++)if(!fl[i]){
        pr[++pc]=i;
        for(j=i+i;j<N;j+=i)fl[j]=1;
    }
    for(scanf("%d",&T);T--;){
        scanf("%lld%lld",&n,&K);
        printf("%llu\n",solve(n));
    }
    return 0;
}

参考资料

1. 朱震霆的国家集训队论文： 《一些特殊的数论函数求和问题》, 国家集训队2018论文集
2. 新版min25筛（O(n^(2/3))）详解 https://zhuanlan.zhihu.com/p/60378354
3. yyb的博客 https://www.cnblogs.com/cjyyb/p/9185093.html