min25筛学习

简介

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

算法原理

一些符号

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

• $minp(i)$，表示$i$的最小质因子
• $p{r}_k$，表示$n$以内第$k$个质数
• $|pr(n)|$，表示$n$以内质数的个数

g函数

$g(n,i)$表示在埃氏筛中，在$n$以内前$i$个质数筛完后，剩下的数的$f$值之和。剩下的数包含所有质数，以及最小质因子大于$p{r}_i$的合数。

$$g(n,i)=\sum _{(j\in prime)\wedge (minp(j)>p{r}_i)}f(j)$$

由此可知$g(n,|pr(n)|)$表示做完埃氏筛后未被删掉的数（即$n$以内所有质数）的$f$值之和，$g(p{r}_k,k)$表示前$k$个质数的$f$值之和。

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

$$g(n,i)=\begin{array}{r}\{\begin{array}{cc}g(n,i-1)& p{r}_i^2>n\\ & \\ g(n,i-1)-f(p{r}_i)\times (g(\lfloor \frac{n}{p{r}_i}\rfloor ,i-1)-g(p{r}_{i-1},i-1))& p{r}_i^2\le n\end{array}\end{array}$$

对于第一部分，因为$p{r}_i^2>n$，所以在埃氏筛中$p{r}_i$无法筛掉任何一个数，所以等于上一个状态。在埃氏筛中，一个合数只会被它的最小质因数筛掉，而质数不会被筛。

对于第二部分，$g(n,i)$所包含的$f$值，等于$g(n,i-1)$所包含的$f$值减去$p{r}_i$筛掉的数的$f$值，即减去最小质因子为$p{r}_i$的合数的$f$值。

我们发现这个式子也可以写成含$g$的表达式，即$f(p{r}_i)\times (g(\lfloor \frac{n}{p{r}_i}\rfloor ,i-1)-g(p{r}_{i-1},i-1))$。其中$g(\lfloor \frac{n}{p{r}_i}\rfloor ,i-1)$是$\lfloor \frac{n}{p{r}_i}\rfloor$以内用前$i-1$个质数筛过后的$f$值之和，它包含两个部分：

• 最小质因子大于$p{r}_{i-1}$的数的$f$值
• 前$i-1$个质数的$f$值，即$g(p{r}_{i-1},i-1)$

那么减掉第二部分后，得到的是$\lfloor \frac{n}{p{r}_i}\rfloor$以内最小质因子大于$p{r}_{i-1}$的数的$f$值。如果再乘上$p{r}_i$，得到的是$n$以内最小质因子等于$p{r}_i$的合数的$f$值之和，这部分就是$p{r}_i$筛去的数的$f$值之和。

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

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

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

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

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

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)$表示$n$以内的最小质因子大于$p{r}_i$的数的$f$值之和，则有

$$S(n,i)=g(n,|pr(n)|)-g(p{r}_i,|pr(n)|)+\sum _{j>i}\sum _{p{r}_j^k\le n}f(p{r}_j^k)\times (S(\lfloor \frac{n}{p{r}_j^k}\rfloor ,j)+[k>1])$$

其中$g(n,|pr(n)|)-g(p{r}_i,|pr(n)|)$表示$n$以内大于$p{r}_i$的质数的$f$值之和。后面相当于求$n$以内不是质数的最小质因子大于$p{r}_i$的数的$f$值之和。当$k=1$时，$f(p{r}_j)$在前面已经被计算过了，所以不用加1；当$k>1$时，$p{r}_j^k$不是质数，所以要加上。

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

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

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

当$p{r}_j^2>n$时，在$k=1$时$S(\lfloor \frac{n}{p{r}_j^k}\rfloor ,j)$显然为$0$（$\lfloor \frac{n}{p{r}_j^k}\rfloor <p{r}_j$，不可能有数在$\lfloor \frac{n}{p{r}_j^k}\rfloor <p{r}_j$以内且最小质因子大于$p{r}_j$）。所以我们只要求$\sqrt{n}$以内的质因子。而$S(\lfloor \frac{n}{p{r}_j^k}\rfloor ,j)$中第一个参数也只能取数论分块的$2\sqrt{n}$个值。

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)=\sum _{j>i}\sum _{p{r}_j^k\le n}f(p{r}_j^k)\times (S(\lfloor \frac{n}{p{r}_j^k}\rfloor ,j)+1)$$

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

注意当$+[k>1]$变为$+1$时，若$k=1$，$p{r}_j^2>n$，则$p{r}_j$也有贡献，这就要求我们求出$n$以内的所有质数。一般的$n$都是$10^9$以上的级别的，这显然是不可行的。

总结

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

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

例题

洛谷P5325 【模板】Min_25筛

$f(p^k)=p^k(p^k-1)=p^{2k}-p^k$，令$f_1(x)=x^2$，$f_2(x)=x$，则$f(p^k)=f_1(p^k)+f_2(p^k)$。因为$f_1$和$f_2$都是完全积性函数，所以我们可以分别将$f_1,f_2$带入$g$的递推式中，做min25筛。

code

#include<bits/stdc++.h>
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;
}