Min_25筛

蒟蒻tyy 按在地上摩擦/dk

-1.参考资料

oi wiki

0.定义

定义 \(x/y=\lfloor\frac{x}{y}\rfloor\)

定义

\[\operatorname{isprime}(n)=\begin{cases}1&n\ prime\\0&otherwise\end{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_{2\le p\le n}f(p) \]

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

枚举 \(i\)最小质因子,有:

\[\begin{aligned}F_k(n)&=\sum_{i=2}^{n}f(i)[p_k\le \operatorname{lpf}(i)]\\ &=\sum_{i\ge k,p_i^2\le n}\sum_{c\ge 1,p_i^c\le n}f(p_i^c)([c>1]+F_{i+1}(n/p_i^c))+\sum_{i\ge k,p_i\le n}f(p_i)\\ &=\sum_{i\ge k,p_i^2\le n}\sum_{c\ge 1,p_i^c\le n}f(p_i^c)([c>1]+F_{i+1}(n/p_i^c))+F_{prime}(n)-F_{prime}(p_{k-1})\\ &=\sum_{i\ge k,p_i^2\le n}\sum_{c\ge 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_{c\ge 1,p_i^c\le n}f(p_i^c)([c>1]+F_{i+1}(n/p_i^c)) \]

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

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

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

\[=\sum_{c\ge 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\)\(2\sqrt n\) 个数有用。

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

给定 \(n,s,g(p)=p^s\),对所有 \(m=n/i\),求 \(\sum_{p\le 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_{\lfloor\sqrt n\rfloor}(n)\)

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

考虑每部分的贡献:

  1. 对于 \(n<p_k^2\) 的部分,\(G\) 值不变,贡献 \(G_{k-1}(n)\)
  2. 对于 \(p_k^2\le n\) 的部分,被筛掉的数必有质因子 \(p_k\),贡献 \(-g(p_k)G_{k-1}(n/p_k)\)
  3. 对于第二部分,由于 \(p_k^2\le n\Leftrightarrow p_k\le 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^2\le n]g(p_k)(G_{k-1}(n/p_k)-G_{k-1}(p_{k-1})) \]

2.复杂度分析

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

对于 \(F_{prime}(n)\) 的计算,可以证明复杂度为 \(O\left(\frac{n^\frac{3}{4}}{\log n}\right)\)

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\n",(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)\)
posted @ 2020-10-07 23:08  happydef  阅读(217)  评论(3编辑  收藏  举报