ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

min_25筛瞎写

2020-05-14 21:54:42  阅读:348  来源: 互联网

标签:lfloor 25 min int ll 筛瞎写 include mo define


%%%gmh差不多一年前学会min_25筛
%%%某初一大佬似乎已经会了min_25筛
菜哭了

约定

以下记\(P\)为素数集合,\(P(n)\)为所有小于等于\(n\)的素数的集合。
\(minp(x)\)表示\(x\)的最小质因子

问题

这种什么筛之类的,多是求积性函数的前缀和的算法。
min_25筛能做的积性函数要满足这两个条件:

  1. 对于所有\(p\in P\),\(f(p)\)能够用多项式表示出来。
  2. 对于所有\(p \in P,k\in N\),\(f(p^k)\),\(f(p^k)\)能够快速计算。

介绍

min_25筛具体分为两个部分:求素数的答案,和求所有数的答案(不包括\(1\)的答案,\(1\)的答案最后自己加上)。

算\(g\)

先说素数部分怎么求。
分别枚举每一项的系数,不同的系数分开来算。
以下记素数\(x\)的函数为\(f(x)\)(指某一特定的系数下)

设\(g(n,j)=\sum_{i \in P(n) 或 minp(i)>p_j}f(i)\)
考虑从\(g(n,j-1)\)得到\(g(n,j)\)
想想\(g(n,j-1)\)得到\(g(n,j)\)要去掉些什么,根据概念可知要去掉所有满足\(minp(x)=p_j\)的合数\(x\)。

这个合数最小是\(p_j^2\),于是当\(n<p_j^2\)时,有\(g(n,j)=g(n,j-1)\)

古爷:这就是min_25筛的精髓!
有了这个东西,意味着大于\(\sqrt n\)的质数都没必要算,
\(j\)最多取到\(|P(\sqrt n)|\)就可以了。

当\(n\ge p_j^2\)时,考虑将\(p_j\)的倍数筛去。
由于是积性函数,所以可以将\(f(p_j)\)提出来,剩下的就是\(g(\lfloor\frac{n}{p_j} \rfloor,j-1)\)。
真的是吗?其实这样还算多了小于\(p_j\)的所有质数的贡献,即\(g(p_j-1,j-1)\),减去即可。
于是这个时候\(g(n,j)=g(n,j-1)+f(p_j)(g(\lfloor\frac{n}{p_j} \rfloor,j-1)-g(p_j-1,j-1))\)

注意一点,上面说“由于是积性函数”,这里应该是广义的积性函数(即便不是题目想要求的),也就是\(f(x)\),形式为\(f(x)=x^k\)
这个步骤下我们最终只是需要素数的贡献和,至于合数的贡献,我们把它当做广义积性函数来算。
合数的贡献只是中间值,计算到最后是会消掉的。

临界:当\(j=0\)的时候,直接用自然数幂和来计算。注意减去\(1\)的贡献。

显然,\(g(n,|P(n)|)\)就是所有小于等于\(n\)的质数的答案。

联想?

形象化地理解这个东西:
当从\(P(n,j-1)\)转移到\(P(n,j)\)的时候,将\(minp(x)=p_j\)的合数\(x\)去掉。
这个东西看起来很像埃氏筛法。


求\(S\)

接下来是计算所有数的答案。
这里的\(F(x)\)直接表示\(x\)的贡献(不需要分具体是哪一位)。
设\(S(n,j)=\sum_{minp(i)\ge p_j} f(i)\)
先列出式子:
\(S(n,j)=g(n,|P(n)|)-\sum_{k<j} F(k)+\sum_{k\ge j,p_k^2\le n}\sum_{e>0,p_k^{e+1}\le n}(F(p_k^e)S(\lfloor \frac{n}{p_k^e}\rfloor,k+1)+F(p_k^{e+1}))\)
前面的\(g(n,|P(n)|)-\sum_{k<j} f(k)\)是所有大于等于\(p_j\)的素数的答案。
后面就是枚举最小质因子\(p_k\)及其系数\(e\),提出来之后计算。
\(S(\lfloor \frac{n}{p_k^e}\rfloor,k+1)\)如果不为\(0\),则至少要满足\(p_k\le\lfloor\frac{n}{p_k^e}\rfloor\),移项一下就得到了\(p_k^{e+1}\le n\)
最后面之所以要加上那个东西,是因为\(S(\lfloor \frac{n}{p_k^e}\rfloor,k+1)\)中并不包括\(p_k\),所以要单独计算\(p_k^{e+1}\)的贡献。

注意,这里的积性函数是“真”的积性函数
对于计算到的每个数,每个质因子都只会枚举一遍。
并不一定是广义积性函数。

答案就是\(S(n,1)+1\),后面那个是\(1\)的贡献。

自己瞎玩出来的求\(S\)做法

观察前面的那个式子,感觉好丑啊。
为什么还要枚举\(k\)?直接枚举\(p_j\)的系数就好了啊!
\(S(n,j)=S(n,j+1)+F(p_j)+\sum_{e>0,p_j^{e+1}\le n} (F(p_j^e)S(\lfloor\frac{n}{p_j^e}\rfloor)+F(p_j^{e+1}))\)
哇哈哈看起来似乎连\(g\)都不用求
然后得到了古爷的指导:这个东西,\(O(n)\)起步。
后面的那一段求和只有\(p_j^2\le n\)的时候是要做的,但是前面的那个要做\(|P(n)|\)遍。

于是就有个简单的优化:
当\(p_j^2>n\)的时候,直接计算\(\sum_{p_k^2>n} F(p_k)\),加上之后退出。
这个东西用\(g\)来求。所以还是要求\(g\)啊……

另一种高级的求\(S\)做法

待填坑。


实现?

观察一下式子,可以发现,要计算到的\(n\)都满足\(n\in\{\lfloor\frac{N}{i}\rfloor|i\in[1,N],i \in Z\}\)
\(N\)是题目给出的\(N\)。
自己证
可以发现这个集合的大小不超过\(2\sqrt N\)
考虑预处理出每个\(g(n,|P(n)|)\)。
简写为\(g_n\)
先考虑怎么存储。对于每个\(n\),可以按小于等于或者大于\(\sqrt N\)讨论,后者借助\(\lfloor\frac{N}{n}\rfloor\)来存储。可以证明\(\lfloor\frac{N}{n}\rfloor\)对于不同的\(n\)不重复。
接下来考虑模拟筛法来求:
枚举质数\(p_j\),然后计算当前这一层的\(g_n\)
具体见代码

求\(S\)的时候,不用记忆化,不用记忆化,不用记忆化!
另外,用前面那个方法就可以了,中间那个方法听说会慢,我自己脑补出来的那个方法和前面那个差不多(无法理解!为什么看起来优美那么多却没有地变快?)


实现

例题见洛谷板题。

大众做法

using namespace std;
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cassert>
#define M 100010
#define mo 1000000007
#define ll long long
ll qpow(ll x,ll y=mo-2){
	ll r=1;
	for (;y;y>>=1,x=x*x%mo)
		if (y&1)
			r=r*x%mo;
	return r;
}
ll _n,sq;
int p[M],np,sf[M],sg1[M],sg2[M];
bool inp[M];
int lsp[M];
#define f(p) ((ll)(p)*((p)-1))
ll w[M*2],cnt;
int id1[M],id2[M];
#define id(d) (d<=sq?id1[d]:id2[_n/d])
ll g1[M*2],g2[M*2];
void init(ll n){
	ll sq=sqrt(n);
	ll inv2=qpow(2),inv6=qpow(6);
	for (ll i=1;i<=n;++i){
		ll j=n/i;
		w[++cnt]=j;
		(j<=sq?id1[j]:id2[n/j])=cnt;
		i=n/j;
		j%=mo;
		g1[cnt]=(j*(j+1)%mo*inv2-1+mo)%mo;
		g2[cnt]=(j*(j+1)%mo*(2*j+1)%mo*inv6-1+mo)%mo;
	}
	for (int j=1;j<=np;++j)
		for (int i=1;i<=cnt && (ll)p[j]*p[j]<=w[i];++i){
			ll d=w[i]/p[j],k=id(d);
			g1[i]=((g1[i]-(ll)p[j]*(g1[k]-sg1[j-1]))%mo+mo)%mo;	
			g2[i]=((g2[i]-(ll)p[j]*p[j]%mo*(g2[k]-sg2[j-1]))%mo+mo)%mo;
		}
}
ll S(ll n,int j){
	ll tmp=id(n),r=(g2[tmp]-g1[tmp])-(j?sf[j-1]:0);
	r%=mo,r+=r<0?mo:0;
	if (j>sq)
		return r;
	for (int k=j;k<=np && (ll)p[k]*p[k]<=n;++k)
		for (ll w=p[k];w*p[k]<=n;w*=p[k])
			(r+=f(w%mo)%mo*S(n/w,k+1)+f(w*p[k]%mo))%=mo;
	return r;
}
int main(){
	scanf("%lld",&_n);
	sq=sqrt(_n);
	for (int i=2;i<=sq;++i){
		if (!inp[i]){
			p[++np]=i;
			sf[np]=(sf[np-1]+f(i))%mo;
			sg1[np]=(sg1[np-1]+i)%mo;
			sg2[np]=(sg2[np-1]+(ll)i*i)%mo;
		}
		lsp[i]=np;
		for (int j=1;j<=np && (ll)i*p[j]<=sq;++j){
			inp[i*p[j]]=1;
			if (i%p[j]==0)
				break;
		}
	}
	init(_n);
	printf("%lld\n",(1+S(_n,1))%mo);
	return 0;
}

我的辣鸡做法

using namespace std;
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cassert>
#define M 100010
#define mo 1000000007
#define ll long long
ll qpow(ll x,int y=mo-2){
	ll r=1;
	for (;y;y>>=1,x=x*x%mo)
		if (y&1)
			r=r*x%mo;
	return r;
}
ll _n,sq;
int p[M],np,sf[M],sg1[M],sg2[M];
bool inp[M];
int lsp[M];
#define f(p) ((ll)(p)*((p)-1))
ll w[M*2],cnt;
int id1[M],id2[M];
#define id(d) (d<=sq?id1[d]:id2[_n/d])
ll g1[M*2],g2[M*2];
void init(ll n){
	ll sq=sqrt(n);
	ll inv2=qpow(2),inv6=qpow(6);
	for (ll i=1;i<=n;++i){
		ll j=n/i;
		w[++cnt]=j;
		(j<=sq?id1[j]:id2[n/j])=cnt;
		i=n/j;
		j%=mo;
		g1[cnt]=(j*(j+1)%mo*inv2-1+mo)%mo;
		g2[cnt]=(j*(j+1)%mo*(2*j+1)%mo*inv6-1+mo)%mo;
	}
	for (int j=1;j<=np;++j)
		for (int i=1;i<=cnt && (ll)p[j]*p[j]<=w[i];++i){
			ll d=w[i]/p[j],k=id(d);
			g1[i]=((g1[i]-(ll)p[j]*(g1[k]-sg1[j-1]))%mo+mo)%mo;	
			g2[i]=((g2[i]-(ll)p[j]*p[j]%mo*(g2[k]-sg2[j-1]))%mo+mo)%mo;
		}
}
ll S(ll n,int j){
	ll r=0;
	if (j>np || (ll)p[j]*p[j]>n){
		int l=min(np,j-1),t=id(n);
		r+=g2[t]-g1[t]-sf[l];
		return (r%mo+mo)%mo;
	}
	r=(f(p[j]%mo)+S(n,j+1))%mo;
	ll w=p[j];
	for (int e=1;w*p[j]<=n;++e,w*=p[j])
		(r+=f(w%mo)%mo*S(n/w,j+1)+f(w*p[j]%mo))%=mo;
	return r;
}
int main(){
	freopen("in.txt","r",stdin);
	scanf("%lld",&_n);
	sq=sqrt(_n);
	for (int i=2;i<=sq;++i){
		if (!inp[i]){
			p[++np]=i;
			sf[np]=(sf[np-1]+f(i))%mo;
			sg1[np]=(sg1[np-1]+i)%mo;
			sg2[np]=(sg2[np-1]+(ll)i*i)%mo;
		}
		lsp[i]=np;
		for (int j=1;j<=np && (ll)i*p[j]<=sq;++j){
			inp[i*p[j]]=1;
			if (i%p[j]==0)
				break;
		}
	}
	init(_n);
	printf("%lld\n",(1+S(_n,1))%mo);
	return 0;
}

标签:lfloor,25,min,int,ll,筛瞎写,include,mo,define
来源: https://www.cnblogs.com/jz-597/p/12891605.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有