简介

min25筛能在在亚线性时间复杂度求出一类积性函数的前缀和，前提是这个积性函数在质数和质数的幂等位置的函数值比较好求。它借助埃氏筛的思想，将原问题转化为在质因子处的相似的子问题，从而得到一个递推式，达到快速求解的目的。

算法原理

一些符号

为了方便描述，定义以下符号：
设我们的问题是求F(n)∑i=1nf(i)F(n)\sum\limits_{i=1}^nf(i)F(n)i=1∑n​f(i)，假设f(i)f(i)f(i)是一个完全积性函数。

• minp(i)minp(i)minp(i)，表示iii的最小质因子
• prkpr_kprk​，表示nnn以内第kkk个质数
• ∣pr(n)∣|pr(n)|∣pr(n)∣，表示nnn以内质数的个数

g函数

g(n,i)g(n,i)g(n,i)表示在埃氏筛中，在nnn以内前iii个质数筛完后，剩下的数的fff值之和。剩下的数包含所有质数，以及最小质因子大于pripr_ipri​的合数。

g(n,i)=∑(j∈prime)∧(minp(j)>pri)f(j)g(n,i)=\sum\limits_{(j\in prime)\wedge(minp(j)>pr_i)}f(j)g(n,i)=(j∈prime)∧(minp(j)>pri​)∑​f(j)

由此可知g(n,∣pr(n)∣)g(n,|pr(n)|)g(n,∣pr(n)∣)表示做完埃氏筛后未被删掉的数（即nnn以内所有质数）的fff值之和，g(prk,k)g(pr_k,k)g(prk​,k)表示前kkk个质数的fff值之和。

我们可以用如下递推式求g(n,i)g(n,i)g(n,i)。

g(n,i)={g(n,i−1)pri2>ng(n,i−1)−f(pri)×(g(⌊npri⌋,i−1)−g(pri−1,i−1))pri2≤ng(n,i)= \begin{aligned} \left\{\begin{matrix} g(n,i-1)\qquad &pr_i^2>n \\ \\ g(n,i-1)-f(pr_i)\times(g(\lfloor\dfrac{n}{pr_i}\rfloor,i-1)-g(pr_{i-1},i-1)) \qquad &pr_i^2\leq n \end{matrix}\right. \end{aligned}g(n,i)=⎩⎨⎧​g(n,i−1)g(n,i−1)−f(pri​)×(g(⌊pri​n​⌋,i−1)−g(pri−1​,i−1))​pri2​>npri2​≤n​​

对于第一部分，因为pri2>npr_i^2>npri2​>n，所以在埃氏筛中pripr_ipri​无法筛掉任何一个数，所以等于上一个状态。在埃氏筛中，一个合数只会被它的最小质因数筛掉，而质数不会被筛。

对于第二部分，g(n,i)g(n,i)g(n,i)所包含的fff值，等于g(n,i−1)g(n,i-1)g(n,i−1)所包含的fff值减去pripr_ipri​筛掉的数的fff值，即减去最小质因子为pripr_ipri​的合数的fff值。

我们发现这个式子也可以写成含ggg的表达式，即f(pri)×(g(⌊npri⌋,i−1)−g(pri−1,i−1))f(pr_i)\times(g(\lfloor\dfrac{n}{pr_i}\rfloor,i-1)-g(pr_{i-1},i-1))f(pri​)×(g(⌊pri​n​⌋,i−1)−g(pri−1​,i−1))。其中g(⌊npri⌋,i−1)g(\lfloor\dfrac{n}{pr_i}\rfloor,i-1)g(⌊pri​n​⌋,i−1)是⌊npri⌋\lfloor\dfrac{n}{pr_i}\rfloor⌊pri​n​⌋以内用前i−1i-1i−1个质数筛过后的fff值之和，它包含两个部分：

• 最小质因子大于pri−1pr_{i-1}pri−1​的数的fff值
• 前i−1i-1i−1个质数的fff值，即g(pri−1,i−1)g(pr_{i-1},i-1)g(pri−1​,i−1)

那么减掉第二部分后，得到的是⌊npri⌋\lfloor\dfrac{n}{pr_i}\rfloor⌊pri​n​⌋以内最小质因子大于pri−1pr_{i-1}pri−1​的数的fff值。如果再乘上pripr_ipri​，得到的是nnn以内最小质因子等于pripr_ipri​的合数的fff值之和，这部分就是pripr_ipri​筛去的数的fff值之和。

但是只有当fff为完全积性函数时，f(pri)f(pr_i)f(pri​)才能提取出来。而这里fff是积性函数，却不一定是完全积性函数，我们怎么能够这样推出式子呢？而且，这个递推式的初始条件是g(n,0)g(n,0)g(n,0)，而g(n,0)g(n,0)g(n,0)正是题目要求的，怎么可以用来做初始值呢？

其实在这里的g(n,0)g(n,0)g(n,0)所求的并不是真正的fff值之和，它求的是另一个完全积性函数的和，假设这个函数为f′f'f′。在质数位置和质数幂的位置，f′f'f′值与fff值相等；而在其他位置，f′f'f′值与fff值不一定相等。设置f′f'f′的目的，是因为这里需要完全积性函数。用f′f'f′代替fff，即可使递推式成立。

下文要求的ggg值中只包含fff在质数和质数幂处的值，其他地方的值都被筛掉了，所以求出的答案是正确的。

递推式已经推出来了，接下来看看如何实现。

ggg函数有两个参数直接做的话，时间复杂度和空间复杂度都很大。我们发现第二维是可以滚动的，所以空间上可以省掉第二维。
令si=g(pri,i)s_i=g(pr_i,i)si​=g(pri​,i)，则sis_isi​可以在求质数时求出。又因为⌊⌊na⌋b⌋=⌊nab⌋\lfloor\dfrac{\lfloor\frac na\rfloor}{b}\rfloor=\lfloor\dfrac{n}{ab}\rfloor⌊b⌊an​⌋​⌋=⌊abn​⌋，所以在递推中所有ggg值的参数都在数论分块的2n2\sqrt n2n​个值之中。所以我们一开始将这2n2\sqrt n2n​个数先存起来，然后在上面进行滚动求值。

code

void init(){x=sqrt(n)+1;dd();//埃氏筛for(long long l=1,r;l<=n;l=r+1){r=n/(n/l);v[++vt]=n/l;//记下数论分块中可能用到的值if(n/l<=x) f[0][n/l]=vt;else f[1][l]=vt;long long t=n/l;g[vt]=...;//求g(v[vt],0)处的值}long long w;for(int i=1;i<=p1;i++){w=1ll*pr[i]*pr[i];long long t;for(int j=1;j<=vt&&w<=v[j];j++){t=v[j]/pr[i];if(t<=x) t=f[0][t];else t=f[1][n/t];g[j]=(...);//求g(j,i)的值}}//滚动求g(n,i)
}

S函数

设S(n,i)S(n,i)S(n,i)表示nnn以内的最小质因子大于pripr_ipri​的数的fff值之和，则有

S(n,i)=g(n,∣pr(n)∣)−g(pri,∣pr(n)∣)+∑j>i∑prjk≤nf(prjk)×(S(⌊nprjk⌋,j)+[k>1])S(n,i)=g(n,|pr(n)|)-g(pr_i,|pr(n)|)+\sum\limits_{j>i}\sum\limits_{pr_j^k\leq n}f(pr_j^k)\times (S(\lfloor\dfrac{n}{pr_j^k}\rfloor,j)+[k>1])S(n,i)=g(n,∣pr(n)∣)−g(pri​,∣pr(n)∣)+j>i∑​prjk​≤n∑​f(prjk​)×(S(⌊prjk​n​⌋,j)+[k>1])

其中g(n,∣pr(n)∣)−g(pri,∣pr(n)∣)g(n,|pr(n)|)-g(pr_i,|pr(n)|)g(n,∣pr(n)∣)−g(pri​,∣pr(n)∣)表示nnn以内大于pripr_ipri​的质数的fff值之和。后面相当于求nnn以内不是质数的最小质因子大于pripr_ipri​的数的fff值之和。当k=1k=1k=1时，f(prj)f(pr_j)f(prj​)在前面已经被计算过了，所以不用加1；当k>1k>1k>1时，prjkpr_j^kprjk​不是质数，所以要加上。

SSS的初始值为S(n,∣pr(n)∣)=0S(n,|pr(n)|)=0S(n,∣pr(n)∣)=0。

我们要求的是S(n,0)+f(1)S(n,0)+f(1)S(n,0)+f(1)。

观察上述的递推式，前两项已经求出。最后一项有两个求和号，分别枚举质因数和它的幂数。

当prj2>npr_j^2>nprj2​>n时，在k=1k=1k=1时S(⌊nprjk⌋,j)S(\lfloor\dfrac{n}{pr_j^k}\rfloor,j)S(⌊prjk​n​⌋,j)显然为000（⌊nprjk⌋​以内的质因子。而S(⌊nprjk⌋,j)S(\lfloor\dfrac{n}{pr_j^k}\rfloor,j)S(⌊prjk​n​⌋,j)中第一个参数也只能取数论分块的2n2\sqrt n2n​个值。

code

long long S(long long p,int q){if(pr[q]>=p) return 0;long long td=(p<=x?f[0][p]:f[1][n/p]);long long re=((g[td]-s[q])%mod+mod)%mod;//前两项for(int j=q+1;j<=p1;j++){if(1ll*pr[j]*pr[j]>p) break;for(long long o=1,pq=pr[j];pq<=p;++o,pq=pq*pr[j]){re=(re+1ll*(...)%mod*(S(p/pq,j)%mod+(o>1))%mod)%mod;//最后一项}}return re;
}

思考

观察S(n,i)S(n,i)S(n,i)的递推式，将其改成如下式子

S(n,i)=∑j>i∑prjk≤nf(prjk)×(S(⌊nprjk⌋,j)+1)S(n,i)=\sum\limits_{j>i}\sum\limits_{pr_j^k\leq n}f(pr_j^k)\times (S(\lfloor\dfrac{n}{pr_j^k}\rfloor,j)+1)S(n,i)=j>i∑​prjk​≤n∑​f(prjk​)×(S(⌊prjk​n​⌋,j)+1)

递推式意义不变，而且不用求ggg函数了。那为什么我们不能用这个式子来求SSS呢？

注意当+[k>1]+[k>1]+[k>1]变为+1+1+1时，若k=1k=1k=1，prj2>npr_j^2>nprj2​>n，则prjpr_jprj​也有贡献，这就要求我们求出nnn以内的所有质数。一般的nnn都是10910^9109以上的级别的，这显然是不可行的。

总结

min25筛的时间复杂度为O(n34ln⁡n)O(\dfrac{n^{\frac 34}}{\ln n})O(lnnn43​​)，相对于杜教筛来说，min25筛要快一些，而且不用找g,hg,hg,h以及求它们的前缀和。

有时候题目中的fff是一个多项式，我们需要把它拆开，对每一项都用ggg的递推式求一次，最后再求和。

例题

洛谷P5325 【模板】Min_25筛

f(pk)=pk(pk−1)=p2k−pkf(p^k)=p^k(p^k-1)=p^{2k}-p^kf(pk)=pk(pk−1)=p2k−pk，令f1(x)=x2f_1(x)=x^2f1​(x)=x2，f2(x)=xf_2(x)=xf2​(x)=x，则f(pk)=f1(pk)+f2(pk)f(p^k)=f_1(p^k)+f_2(p^k)f(pk)=f1​(pk)+f2​(pk)。因为f1f_1f1​和f2f_2f2​都是完全积性函数，所以我们可以分别将f1,f2f_1,f_2f1​,f2​带入ggg的递推式中，做min25筛。

code

#include
using namespace std;
const int N=200000;
const long long mod=1000000007;
int p1,vt,z[N+5],pr[N+5];
long long n,x,ny,v[N+5],f[2][N+5],s1[N+5],s2[N+5],g1[N+5],g2[N+5];
long long mi(long long t,long long v){if(!v) return 1;long long re=mi(t,v/2);re=re*re%mod;if(v&1) re=re*t%mod;return re; 
}
void dd(){for(int i=2;i<=N;i++){if(!z[i]){pr[++p1]=i;s1[p1]=(s1[p1-1]+1ll*i*i%mod)%mod;s2[p1]=(s2[p1-1]+i)%mod;}for(int j=1;j<=p1&&i*pr[j]<=N;j++){z[i*pr[j]]=1;if(i%pr[j]==0) break;}}
}
void init(){x=sqrt(n)+1;ny=mi(6ll,mod-2);dd();for(long long l=1,r;l<=n;l=r+1){r=n/(n/l);v[++vt]=n/l;if(n/l<=x) f[0][n/l]=vt;else f[1][l]=vt;long long t=(n/l)%mod;g1[vt]=(t*(t+1)%mod*(2*t+1)%mod*ny%mod-1+mod)%mod;g2[vt]=(t*(t+1)/2%mod-1+mod)%mod;}long long w;for(int i=1;i<=p1;i++){w=1ll*pr[i]*pr[i];long long t;for(int j=1;j<=vt&&w<=v[j];j++){t=v[j]/pr[i];if(t<=x) t=f[0][t];else t=f[1][n/t];g1[j]=(g1[j]-1ll*pr[i]*pr[i]%mod*(g1[t]-s1[i-1])%mod+mod)%mod;g2[j]=(g2[j]-1ll*pr[i]*(g2[t]-s2[i-1])%mod+mod)%mod;}}
}
long long S(long long p,int q){if(pr[q]>=p) return 0;long long td=(p<=x?f[0][p]:f[1][n/p]);long long re=((g1[td]-s1[q]-g2[td]+s2[q])%mod+mod)%mod;for(int j=q+1;j<=p1;j++){if(1ll*pr[j]*pr[j]>p) break;for(long long o=1,pq=pr[j],t;pq<=p;++o,pq=pq*pr[j]){t=pq%mod;re=(re+t*(t-1)%mod*(S(p/pq,j)%mod+(o>1))%mod)%mod;}}return re;
}
int main()
{scanf("%lld",&n);init();printf("%lld",(S(n,0)+1)%mod);return 0;
}