Min_25筛

目录

• -1.参考资料
• 0.定义
• 1.算法
• 2.复杂度分析
• 3.实现

被 蒟蒻tyy 按在地上摩擦/dk

-1.参考资料

oi wiki

0.定义

定义 (x/y=lfloorfrac{x}{y} floor)。

定义

[operatorname{isprime}(n)=egin{cases}1&n prime\0&otherwiseend{cases} ]

定义 (p_k) 为第 (k) 个质数，且 (p_0=1)。

定义 (operatorname{lpf}(n)=min{p|[p|n]})。

1.算法

min_25 筛可以在 (O(frac{n^frac{3}{4}}{log n})) 或者 (O(n^{1-epsilon}))时间内计算出积性函数前缀和问题。

要求：(f(p^c)) 是关于 (p) 的项数较少的多项式，可以快速求值。

考虑 dp。设

[F_{prime}(n)=sum_{2le ple n}f(p) ]

[F_k(n)=sum_{i=2}^{n}f(i)[p_kle operatorname{lpf}(i)] ]

枚举 (i)最小质因子，有：

[egin{aligned}F_k(n)&=sum_{i=2}^{n}f(i)[p_kle operatorname{lpf}(i)]\ &=sum_{ige k,p_i^2le n}sum_{cge 1,p_i^cle n}f(p_i^c)([c>1]+F_{i+1}(n/p_i^c))+sum_{ige k,p_ile n}f(p_i)\ &=sum_{ige k,p_i^2le n}sum_{cge 1,p_i^cle n}f(p_i^c)([c>1]+F_{i+1}(n/p_i^c))+F_{prime}(n)-F_{prime}(p_{k-1})\ &=sum_{ige k,p_i^2le n}sum_{cge 1,p_i^{c+1}le n}(f(p_i^c)F_{i+1}(n/p_i^c)+f(p_i^{c+1}))+F_{prime}(n)-F_{prime}(p_{k-1})\ end{aligned} ]

最后一步是因为

[sum_{cge 1,p_i^cle n}f(p_i^c)([c>1]+F_{i+1}(n/p_i^c)) ]

[=sum_{cge 1,p_i^cle n}f(p_i^c)[c>1]+sum_{cge 1,p_i^cle n}f(p_i^c)F_{i+1}(n/p_i^c) ]

而因为对于满足 (p_i^cle n<p_i^{c+1}) 的 (c)，有 (1le n/p_i^c<p_i<p_{i+1})，所以 (F_{i+1}(n/p_i^c)=0)，所以

[=sum_{cge 1,p_i^{c+1}le n}f(p_i^{c+1})+sum_{cge 1,p_i^{c+1}le n}f(p_i^c)F_{i+1}(n/p_i^c) ]

[=sum_{cge 1,p_i^{c+1}le n}(f(p_i^c)F_{i+1}(n/p_i^c)+f(p_i^{c+1})) ]

下面有两种处理方法：

1. 直接按照定义式递推。
2. 从大到小枚举 (p) ，注意到仅当 (p^2<n) 时的转移增加值不为 (0)，所以按照递推式后缀和优化即可

下面考虑如何计算 (F_{prime}(n))。

观察求 (F_k(n)) 的过程，我们会发现只有 (1,2,dots,sqrt n,n/sqrt n,n/(sqrt n-1),dots,n) 共 (2sqrt n) 个数有用。

因为 (f(p)) 是关于 (p) 的多项式，按照次数拆分后，我们只要考虑 (f(p)=p^s) 的情况，于是问题就转化为了：

给定 (n,s,g(p)=p^s)，对所有 (m=n/i)，求 (sum_{ple m}g(p))。

设

[G_k(n)=sum_{i=1}^{n}[p_k<operatorname{lpf}(i)lor operatorname{isprime(i)}]g(i) ]

对任意合数 (x)，都有 (operatorname{lpf}le sqrt x)，所以所求 (F_{prime}(n)=G_{lfloorsqrt n floor}(n))。

显然有 (G_0(n)=sum_{i=2}^{n}g(i))，考虑如何转移。

考虑每部分的贡献：

1. 对于 (n<p_k^2) 的部分，(G) 值不变，贡献 (G_{k-1}(n))。
2. 对于 (p_k^2le n) 的部分，被筛掉的数必有质因子 (p_k)，贡献 (-g(p_k)G_{k-1}(n/p_k))。
3. 对于第二部分，由于 (p_k^2le nLeftrightarrow p_kle n/p_k)，所以有一部分 (operatorname{lpf}(i)<p_k) 的 (i) 被减去，这部分要加回来，即 (g(p_k)G_{k-1}(p_{k-1}))。

所以

[G_k(n)=G_{k-1}(n)-[p_k^2le n]g(p_k)(G_{k-1}(n/p_k)-G_{k-1}(p_{k-1})) ]

2.复杂度分析

对于 (F_k(n)) 的计算，可以证明，第一种方法复杂度 (O(n^{1-epsilon}))，第二种方法复杂度 (Oleft(frac{n^frac{3}{4}}{log n} ight))。

对于 (F_{prime}(n)) 的计算，可以证明复杂度为 (Oleft(frac{n^frac{3}{4}}{log n} ight))。

3.实现

一般来说第一种方法代码难度低，而且数据范围较小时比第二种更快，所以有一般采用第一种方法实现。

这是洛谷模板的实现。

Code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll N=1000005,P=1e9+7,inv3=333333336;
ll n,num,cnt,tot,w[N],g1[N],g2[N],pr[N],s1[N],s2[N],ind1[N],ind2[N];
bool vst[N]; 
ll S(ll x,ll y){
	if(pr[y]>=x)return 0;
	ll k=(x<=num?ind1[x]:ind2[n/x]),ret=(g2[k]-g1[k]+P-(s2[y]-s1[y])+P)%P;
	for(ll i=y+1;i<=cnt&&pr[i]*pr[i]<=x;i++)
		for(ll e=1,pe=pr[i],t;pe<=x;e++,pe=pe*pr[i])
			ret=(ret+(t=pe%P)*(t-1)%P*(S(x/pe,i)+(e!=1)))%P;
	return ret%P;
} 
int main(){
	scanf("%lld",&n);
	num=sqrt(n);
	for(ll i=2;i<=num;i++){
		if(!vst[i])pr[++cnt]=i,s1[cnt]=(s1[cnt-1]+i)%P,s2[cnt]=(s2[cnt-1]+i*i)%P;
		for(ll j=1;j<=cnt&&i*pr[j]<=num;j++){
			vst[i*pr[j]]=1;
			if(i%pr[j]==0)break;
		}
	}
	for(ll l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		w[++tot]=n/l;
		g1[tot]=w[tot]%P;
		g2[tot]=g1[tot]*(g1[tot]+1)/2%P*(2*g1[tot]+1)%P*inv3%P-1;
		g1[tot]=g1[tot]*(g1[tot]+1)/2%P-1;
		if(n/l<=num)ind1[n/l]=tot;
		else ind2[r]=tot;
	}
	for(ll i=1;i<=cnt;i++)
		for(ll j=1,k;j<=tot&&pr[i]*pr[i]<=w[j];j++){
			k=(w[j]/pr[i]<=num?ind1[w[j]/pr[i]]:ind2[n/(w[j]/pr[i])]);
			g1[j]-=pr[i]*(g1[k]-s1[i-1]+P)%P;
			g2[j]-=pr[i]*pr[i]%P*(g2[k]-s2[i-1]+P)%P;
			g1[j]=(g1[j]%P+P)%P,g2[j]=(g2[j]%P+P)%P;
		}
	printf("%lld
",(S(n,0)+1)%P);
	return 0;
}

实现时依次做如下事情：

1. 筛出 ([1,sqrt n]) 内的质数和前 (sqrt n) 个 (f) 值；
2. 把 (f(p)) 关于 (p) 写成多项式表示，对每一项算出对应的 (G) 并合并算出 (F_{prime})。
3. 按照 (F_k(n)) 的递推式，用递归求出 (F_1(n))。