min_25筛瞎写

Posted 2021-02-19 jz-597

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

约定

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

问题

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

1. 对于所有(pin P)，(f(p))能够用多项式表示出来。
2. 对于所有(p in P,kin 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)|)就可以了。

当(nge p_j^2)时，考虑将(p_j)的倍数筛去。
由于是积性函数，所以可以将(f(p_j))提出来，剩下的就是(g(lfloorfrac{n}{p_j} floor,j-1))。
真的是吗？其实这样还算多了小于(p_j)的所有质数的贡献，即(g(p_j-1,j-1))，减去即可。
于是这个时候(g(n,j)=g(n,j-1)+f(p_j)(g(lfloorfrac{n}{p_j} floor,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_{kge j,p_k^2le n}sum_{e>0,p_k^{e+1}le n}(F(p_k^e)S(lfloor frac{n}{p_k^e} floor,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} floor,k+1))如果不为(0)，则至少要满足(p_klelfloorfrac{n}{p_k^e} floor)，移项一下就得到了(p_k^{e+1}le n)
最后面之所以要加上那个东西，是因为(S(lfloor frac{n}{p_k^e} floor,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(lfloorfrac{n}{p_j^e} floor)+F(p_j^{e+1})))
哇哈哈看起来似乎连(g)都不用求
然后得到了古爷的指导：这个东西，(O(n))起步。
后面的那一段求和只有(p_j^2le n)的时候是要做的，但是前面的那个要做(|P(n)|)遍。

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

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

待填坑。

---

实现？

观察一下式子，可以发现，要计算到的(n)都满足(nin{lfloorfrac{N}{i} floor|iin[1,N],i in Z})
(N)是题目给出的(N)。
自己证
可以发现这个集合的大小不超过(2sqrt N)
考虑预处理出每个(g(n,|P(n)|))。
简写为(g_n)
先考虑怎么存储。对于每个(n)，可以按小于等于或者大于(sqrt N)讨论，后者借助(lfloorfrac{N}{n} floor)来存储。可以证明(lfloorfrac{N}{n} floor)对于不同的(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
",(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
",(1+S(_n,1))%mo);
	return 0;
}