min_25筛瞎写

%%%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 ∈ { ⌊ N i ⌋ ∣ i ∈ [ 1 , N ] , i ∈ Z } n\in\{\lfloor\frac{N}{i}\rfloor|i\in[1,N],i \in Z\} n∈{⌊iN​⌋∣i∈[1,N],i∈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;
}