Cirno_9's Blog

With a ruthless enemy

0%

3 岁小孩都能懂的 Min_25 筛速通教程

发表于 更新于 分类于 Valine：

目标：求 ∑i=1NF(x)\sum\limits_{i=1}^N F(x)，F(x)F(x) 积性并在素数处能拆成若干完全积性函数的和。

一些约定：fac(x)\mathrm{fac}(x) 表示 xx 最小质因子，pip_i 表示第 ii 个素数， ω(i)\omega(i) 表示 ii 的质因子个数，n/in/i 表示 ⌊ni⌋\left\lfloor\dfrac{n}{i}\right\rfloor，PP 表示素数集合。

再定义一个 H(n)=∑pi≤nF(pi)H(n) = \sum\limits_{p_i \leq n} F(p_i)，同时提醒读者注意文中 nn 与 NN 的区别。

Part 1 最小质因子分类

定义一个 S(n,k)=∑fac(x)≥pknF(x)S(n,k) = \sum\limits_{\mathrm{fac}(x) \ge p_k}^n F(x)，所求即为 S(N,1)S(N,1)。（Tip\rm Tip：fac(1)=+∞\mathrm{fac}(1) = +\infty）

需要做个递推式出来，考虑将其按照 fac\mathrm{fac} 分类来计算。具体而言枚举 pip_i 及其次数 cc，表示我要计算 fac(x)=pi\mathrm{fac}(x) = p_i 且其次数为 cc 的 ∑F(x)\sum F(x)：

S(n,k)=1+∑i=kω(n)∑c=1pic≤nF(pic)S(n/pic,i+1)S(n,k) = 1+\sum_{i=k}^{\omega(n)} \sum_{c=1}^{p_i^c \leq n} F(p_i^c)S(n/p_i^c,i+1)

观察到 pi≥np_i \ge \sqrt n 时，有贡献的只有 pip_i 一个数，那就分类计算：

S(n,k)=1+∑i=kω(n)∑c=1pic≤nF(pic)S(n/pic,i+1)+H(n)−H(n)S(n,k) = 1+\sum_{i=k}^{\omega(\sqrt n)} \sum_{c=1}^{p_i^c \leq n} F(p_i^c)S(n/p_i^c,i+1) + H(n) - H(\sqrt n)

这样的好处是只用计算 pk≤Np_k \leq \sqrt N 的状态。前面一半递归计算，H(x)H(\sqrt x) 可以暴力求，复杂度显然有一个 O(Nlog⁡N)O\left(\dfrac{N}{\log N}\right) 的上界但跑不满（听说在 N≤1013N \leq 10^{13} 时是 O(N0.75log⁡N)O\left( \dfrac{N^{0.75}}{\log N} \right) 的），问题转化为对每个 x=N/ix = N/i 求出 H(x)H(x)。

Part 2 质数点值和

把 FF 在素数处的取值拆成若干完全积性函数的和来算，假设拆出来的完全积性函数是 f1(x),f2(x)…f_1(x), f_2(x) \ldots。定义一个 hk(x)=∑pi≤nfk(pi)h_k(x) = \sum\limits_{p_i \leq n} f_k(p_i)，则 H(x)H(x) 可以用 ∑khk(x)\sum_k h_k(x) 来算，下面考虑 h(x)h(x) 与 f(x)f(x)。

令 g(n,k)=∑x=1n[x∈P∨fac(x)>pk]f(x)g(n,k) = \sum\limits_{x=1}^n [x\in P \lor \mathrm{fac}(x) > p_k]f(x)，显然 h(n)=g(n,ω(n))h(n) = g(n,\omega(n))。

gg 的第 kk 列可以由第 k−1k-1 列递推而来，具体而言 g(n,k)g(n,k) 是由 g(n,k−1)g(n,k-1) 去掉了 fac(x)=pk\mathrm{fac}(x) = p_k 的合数得到，利用 ff 的完全积性可以得到以下递推式：

g(n,k)={g(n,k−1)(pk>n)g(n,k−1)−f(pk)(g(n/pk,k−1)−h(pk−1))(pk≤n)g(n,k) = \begin{cases} g(n,k-1) &(p_k>\sqrt n)\\ g(n,k-1) - f(p_k)(g(n/p_k,k-1) - h(p_{k-1})) &(p_k \leq \sqrt n) \end{cases}

上面一半没卵用，只需要递推出 pk≤Np_k \leq \sqrt N 的点值即可，状态数一共 O(N0.75log⁡N)O\left( \dfrac{N^{0.75}}{\log N} \right) 个，线性筛出 N\sqrt N 以内的 ff 与 hh 即可求出所有 x=n/ix=n/i 的 h(x)h(x)。

代码实现

总而言之，认为复杂度是 O(run 1e10 is ok)O(\text{run 1e10 is ok}) 就行，大概。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include<bits/stdc++.h>
#define FOR(i,a,b) for(int i=a,i##i=b;i<=i##i;++i)
using namespace std;
typedef long long ll;
const int N=2e5+7,mod=1e9+7;
ll n,m,vl[N],g1[N],g2[N];
int isp[N],C,P[N],ct,o[N];
int id(ll x){return x<=m?x:C+1-n/x;}
ll S(ll n,int k){
	if(P[k-1]>=n)return 1;
	ll D=max((int)sqrt(n),P[k-1]),su=1+g2[id(n)]-g1[id(n)]-g2[D]+g1[D];
	FOR(i,k,o[D])for(ll p=P[i];p<=n;p*=P[i]){
		(su+=p%mod*(p%mod-1ll)%mod*S(n/p,i+1))%=mod;
	}
	return su;
}
int main(){
	cin>>n,m=sqrt(n),C=2*m;
	isp[1]=1;
	FOR(i,2,m){
		o[i]=o[i-1];
		if(!isp[i])++o[P[++ct]=i];
		for(int j=1;j<=ct&&i*P[j]<=m;++j){
			isp[i*P[j]]=1;
			if(i%P[j]==0)break;
		}
	}
	FOR(i,1,C){
		ll p=(vl[i]=i<=m?i:n/(C+1-i))%mod;
		g1[i]=p*(p+1ll)/2%mod,g2[i]=(p*2+1ll)*g1[i]%mod*(mod+1)/3%mod;
	}
	FOR(k,1,o[m])for(ll i=C,p=P[k];1ll*p*p<=vl[i];--i){
		ll E=id(vl[i]/p);
		(g1[i]+=p*(g1[p-1]-g1[E]))%=mod;
		(g2[i]+=(g2[p-1]-g2[E])*p%mod*p)%=mod;
	}
	cout<<(S(n,1)+mod)%mod;
	return 0;
}