Min25筛学习笔记

一、简介

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

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

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

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

它的时间复杂度为：O(n34log⁡n)O(\frac{n^{\frac{3}{4}}}{\log{n}})，空间复杂度：O(n)O(\sqrt{n})

二、前置技能：埃氏筛

埃氏筛法每一轮找出第 ii 个质数，并把这个质数的倍数筛掉

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

第一次被 PiP_i 筛掉的 PiP_i 的倍数 xx 必然满足 x≥Pi2x\geq P_{i}^2

三、第一部分：处理质数点的前缀和

对于每个 ⌊nx⌋\left \lfloor \frac{n}{x} \right \rfloor，我们先来计算下面的式子：

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

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

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

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

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

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

设

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

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

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

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

那么我们要求的东西就是

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

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

  1. 若 Pj2>nP_{j}^2>n，那么第 jj 轮埃氏筛没有筛掉任何数

  2. 若 Pj2≤nP_{j}^2\leq n，那么第 jj 轮埃氏筛会筛掉

x≤⌊nPj⌋且x不含小于Pj的质因子{x\leq \left\lfloor\frac{n}{P_j}\right\rfloor}且x不含小于P_j的质因子

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

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

g(n,j)={g(n,j−1),Pj2>ng(n,j−1)−f(Pj)[g(⌊nPj⌋,j−1)−∑i=1j−1f(Pj)],Pj2≤ng(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−1j-1 轮埃氏筛后，区间 [1,Pj−1][1,P_{j-1}] 中就只剩下质数了，所以

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

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

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

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

四、第二部分：计算答案

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

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

质数点的答案为 g(n,+∞)g(n,+\infty) 减去小于 PjP_j 的质数的贡献，即

g(n,+∞)−∑i=2j−1f(Pi)=g(n,+∞)−g(Pj−1,+∞)g(n,+\infty)-\sum_{i=2}^{j-1}f(P_i)=g(n,+\infty)-g(P_{j-1},+\infty)

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

S(n,j)=g(n,+∞)−g(Pj−1,+∞)+∑k=jPk2≤n∑e=1Pke+1≤nS(⌊nPke⌋,k+1)f(Pke)+f(Pke+1)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)S(*,k+1) 保证了没有小于等于 PkP_k 的质因子，这样的话质数幂的情况需要单独考虑

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

五、代码实现

我们预处理所有 x=⌊ni⌋x=\left \lfloor \frac{n}{i} \right \rfloor 的值，从小到大记为 Pos1,Pos2,⋯ ,PosmPos_1,Pos_2,\cdots,Pos_{m}

  1. 对于 1≤x≤n1\leq x\leq\sqrt{n}，显然 x=Posxx=Pos_x

  2. 对于 x>nx>\sqrt{n}，观察下面的值：

{Posm=nPosm−1=⌊n2⌋⋯Posn=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.

我们发现 Posx=⌊nm−x+1⌋Pos_x=\left \lfloor \frac{n}{m-x+1} \right \rfloor，所以 x=m−⌊nPosx⌋+1x=m-\left \lfloor \frac{n}{Pos_x} \right \rfloor+1

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

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

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

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

六、练习题目

【LOJ6235】区间素数个数

第一部分中，取 f(p)=1f(p)=1 即可，求出的 g(n,+∞)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】简单的函数

对于质数PP的值

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

22 的情况单独处理，其余的就是套路了

#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)=σ0(Pk)=k+1f(P)=\sigma_0(P^k)=k+1

f(Pe)=σ0(Pek)=ek+1f(P^e)=\sigma_0(P^{ek})=ek+1

然后就无脑 Min25 筛了

至于对 2642^{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)=pf(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筛》——租酥雨

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