专栏文章

Min_25 筛

ccyh_toby
算法·理论参与者 35已保存评论 39

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
39 条
当前快照
1 份
快照标识符
@mhz5rzz4
此快照首次捕获于
2025/11/15 01:55
4 个月前
此快照最后确认于
2025/11/29 05:24
3 个月前
查看原文

简介

对于积性函数 f(x)f(x)，如果能快速知道 f(pk)f(p^k) 的值（其中 pp 为质数，下同），且 f(p)f(p) 可以表示为项数较少的多项式（或能表示成若干项使得每一项都是完全积性函数），那么 min_25\rm min\_25 筛可以在 O(n34log⁡n)O(\frac{n^{\frac{3}{4}}}{\log n}) 时间复杂度内求出 ∑i=1nf(i)\sum_{i=1}^nf(i)。
由于其思想与埃氏筛类似，所以也被称作扩展埃氏筛。

内容

第一步

推导

目标：求 g(n)=∑p≤nf(p)g(n)=\sum_{p\le n} f(p) 。
不妨设 f(p)f(p) 是完全积性函数，如果不是可以尝试拆成若干项完全积性函数，分别求然后相加。
首先要线性筛求出 n\sqrt n 以内的质数。
g(n)g(n) 很难直接求解，考虑用 DP 计算。设 g(n,j)=∑i=1nf(i)[i是质数或其最小质因子>pj]g(n,j)=\sum_{i=1}^nf(i)[i\text{是质数或其最小质因子}>p_j]，其中 pjp_j 表示第 jj 个质数，那么我们要的就是 g(n,k),k为最小的满足pk≥ng(n,k),k\text{为最小的满足}p_k\ge \sqrt n。考虑从 j−1j-1 变到 jj ，那么最小质因子为 pjp_j 的合数会被筛掉，那么它们的贡献要减去。则有转移
g(n,j)=g(n,j−1)−f(pj)(g(⌊npj⌋,j−1)−g(pj−1,j−1))g(n,j)=g(n,j-1)-f(p_j)\left(g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right)-g(p_{j-1},j-1)\right)
系数 f(pj)f(p_j) 表示由于 f(p)f(p) 是完全积性函数，所以可以把它从后面提出来。 g(⌊npj⌋,j−1)g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right) 表示考虑所有 pjp_j 的倍数，它们除以 pjp_j 之后，最小质因子 >pj−1>p_{j-1} 的合数以及所有质数的贡献，应当减去。但是，这些质数中可能有 ≤pj−1\le p_{j-1} 的，它们在之前就被筛掉过了，所以要加回来，也就是 g(pj−1,j−1)g(p_{j-1},j-1)。
由于有公式 ⌊⌊ab⌋c⌋=⌊abc⌋\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\rfloor=\lfloor\frac{a}{bc}\rfloor ，因此容易发现上述式子只会用到形如 ⌊nx⌋,x≤n\lfloor\frac{n}{x}\rfloor,x\le n 的点处的 DP 值，即第一项的状态数是 O(n)O(\sqrt n)（实际实现的时候注意状态数是 2n2\sqrt n）。我们预处理出这 O(n)O(\sqrt n) 个数，把他们离散化，顺带求出 g(x,0)g(x,0)，然后 DP 即可。
注意到上述转移式中还有一项 g(pj−1,j−1)g(p_{j-1},j-1)，我们只处理了所有形如 ⌊nx⌋\lfloor\frac{n}{x}\rfloor 的数，是否会漏掉某些质数 pp 呢？其实不会漏。注意到，能用来转移的 pp 一定满足 p≤np \le \sqrt n。我们只需要证明 ∀k≤n,∃x,k=⌊nx⌋\forall k \le \sqrt n, \exists x,k=\lfloor\frac{n}{x}\rfloor。
证明：
  1. 若 k=⌊n⌋k=\lfloor\sqrt n\rfloor，设 n=k2+dn=k^2+d，则由于 k2≤n<(k+1)2k^2\le n < (k+1)^2，故 d∈[0,2k]d\in[0,2k]。
    1. d∈[0,k)d\in [0,k)，那么 ⌊nk⌋=k+⌊dk⌋=k\lfloor\frac{n}{k}\rfloor=k+\lfloor\frac{d}{k}\rfloor=k。
    2. d∈[k,2k]d\in [k,2k]，那么 ⌊nk+1⌋=⌊k2+k+d−kk+1⌋=k+⌊d−kk+1⌋=k\lfloor\frac{n}{k+1}\rfloor=\lfloor\frac{k^2+k+d-k}{k+1}\rfloor=k+\lfloor\frac{d-k}{k+1}\rfloor=k。
  2. 若 k<⌊n⌋k < \lfloor\sqrt n\rfloor 即 k≤n−1k \le \sqrt n - 1，假设存在 i,⌊ni+1⌋<k<⌊ni⌋i,\lfloor\frac{n}{i+1}\rfloor<k<\lfloor\frac{n}{i}\rfloor，此时 kk 恰好夹在两个连续的 ⌊nx⌋\lfloor\frac{n}{x}\rfloor 之间，即不可被表出。则 ni+1<k\frac{n}{i+1}<k，故 n<k(i+1)n<k(i+1)，从而 k<⌊ni⌋<⌊k(i+1)i⌋=k+⌊ki⌋k<\lfloor\frac{n}{i}\rfloor<\lfloor\frac{k(i+1)}{i}\rfloor=k+\lfloor\frac{k}{i}\rfloor。另一方面，ni+1<k<n\frac{n}{i+1}<k<\sqrt n，所以 i+1>ni+1>\sqrt n，于是 i>n−1≥ki>\sqrt n - 1 \ge k，因此 ⌊ki⌋=0\lfloor\frac{k}{i}\rfloor=0，于是得到 k<kk<k，故假设不成立，原命题成立。
所以，上述担心就是多虑了。
我看到很多博客、题解都额外预处理了 n\sqrt n 内所有素数处的函数值的前缀和，但因为我“粗心”，代码中直接调用 gg 数组，却顺利通过了一堆题目，写笔记时意识到这一点，于是琢磨出了奥妙所在。其实，在看 zzt 集训队论文时注意到他说只需要用到 ⌊nx⌋\lfloor\frac{n}{x}\rfloor 处的值，就心生疑惑，现在终于明白了。因此，在我看来，预处理根号内素数处的函数前缀和的值是不必要的。

时间复杂度

∑i≤nO(π(i))+∑i≤nO(π(ni))=∑i≤nO(π(ni))=∑i≤nO(nilog⁡ni)=O(∫1nnxlog⁡nxdx)=O(n34log⁡n)\begin{aligned}&\sum_{i\le \sqrt n}O(\pi(\sqrt i))+\sum_{i\le \sqrt n}O(\pi(\sqrt\frac{n}{i}))\\ =&\sum_{i\le \sqrt n}O(\pi(\sqrt\frac{n}{i}))\\ =&\sum_{i\le \sqrt n}O\left(\frac{\sqrt\frac{n}{i}}{\log\sqrt\frac{n}{i}}\right)\\ =&O\left(\int_1^{\sqrt n} \frac{\sqrt\frac{n}{x}}{\log\sqrt\frac{n}{x}}{\rm d}x\right)\\ =&O\left(\frac{n^\frac{3}{4}}{\log n}\right) \end{aligned}

第二步

推导

目标：求 S(n)=∑i≤nf(i)S(n)=\sum_{i\le n} f(i)。与第一步类似，设 S(n,j)=∑i=1nf(i)[i的最小质因子>pj]S(n,j)=\sum_{i=1}^nf(i)[i\text{的最小质因子}>p_j]。但此处 ff 不需要再拆分成单项式，直接是原函数即可（因为不需要依赖于完全积性，只需要积性即可）（但要能快速计算 f(pk)f(p^k) 的值）。

方法一

考虑把贡献拆成质数的和合数的，合数枚举最小质因子以及次数，于是有转移：
S(n,j)=g(n)−g(pj)+∑j<k,pk≤n,1≤e,pke≤nf(pke)(S(⌊npke⌋,k)+[e≠1])S(n,j)=g(n)-g(p_j)+\sum_{j<k,p_k\le \sqrt n,1\le e,p_k^e\le n} f(p_k^e)\left(S\left(\left\lfloor\frac{n}{p_k^e}\right\rfloor,k\right)+[e\neq 1]\right)
最后一项 [e≠1][e\neq 1] 的意思是，对于 e=1e=1 的情况，SS 没有计算 11 贡献，刚好，因为此时 pk×1p_k\times 1 是质数，其贡献在之前计算过；对于 e>1e > 1 的情况，pke×1p_k^e \times 1 是合数，贡献算漏了，要补上。直接暴力递归计算（并且不需要记忆化）。

方法二

也是把贡献拆成质数和合数，只是采用类似于第一步的递推方式：
S(n,j)=f(pj+1)+S(n,j+1)+∑pj+1≤n,1≤e,pj+1e≤nf(pj+1e)(S(⌊npj+1e⌋,j+1)+[e≠1])S(n,j)=f(p_{j+1})+S(n,j+1)+\sum_{p_{j+1}\le \sqrt n,1\le e, p_{j+1}^e\le n}f(p_{j+1}^e)\left(S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)+[e\neq 1]\right)
但直接这样转移复杂度不对，我们要严格控制只有 ≤n\le \sqrt n 的质数才发生转移。即，对于 pj+1>np_{j+1}>\sqrt n 的状态 S(n,j)S(n,j)，不显式地计算出来，用到的时候特判为 g(n)−g(pj)g(n)-g(p_j)。所以，我们要更新状态 S(n,j)S(n,j) 时，一定有 pj+1≤np_{j+1}\le \sqrt n，此时需要用到 S(n,j+1)S(n,j+1)，即使要特判为 g(n)−g(pj+1)g(n)-g(p_{j+1})，由步骤一中的论述可知 g(pj+1)g(p_{j+1}) 已经计算出，所以不会有问题。需要注意 n=2,3n=2,3 的状态不会被任何 pp 更新到，所以也要特判。一种较为简洁的方法是，像方法一的代码一样写一个函数 S(x,y)S(x,y)，只不过不用来递归，而是用来实现各种特判。
但是，这样代码不太好些，且实现出来常数比较大。我实现了一份：code。耗时几乎是别人的两倍。

方法二改进

既然采用了类似第一步的递推方式，干脆直接套用第一步的状态设计。设 S(n,j)=∑i=1nf(i)[i是质数或其最小质因子>pj]S(n,j)=\sum_{i=1}^nf(i)[i\text{是质数或其最小质因子}>p_j] ，那么参考第一步的方程，有转移：
S(n,j)=S(n,j+1)+∑pj+1≤n,1≤e,pj+1e≤nf(pj+1e)(S(⌊npj+1e⌋,j+1)−g(min⁡{⌊npj+1e⌋pj+1})+[e≠1])S(n,j)=S(n,j+1)+\sum_{p_{j+1}\le \sqrt n,1\le e, p_{j+1}^e\le n}f(p_{j+1}^e)\left(S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)-g(\min\{\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor p_{j+1}\})+[e\neq 1]\right)
最后面的 −g(… )-g(\dots) 是因为，S(⌊npj+1e⌋,j+1)S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right) 中包含了 ≤pj+1\le p_{j+1} 的质数的贡献，与转移式的意义不符，应当减去。但是 gg 里面有个 min⁡\min 非常令人不爽，并且为了卡常，我们可以进一步简化上式：
S(n,j)=S(n,j+1)+∑pj+1≤n,1≤e,pj+1e+1≤nf(pj+1e)(S(⌊npj+1e⌋,j+1)−g(pj+1))+f(pj+1e+1)S(n,j)=S(n,j+1)+\sum_{p_{j+1}\le \sqrt n,1\le e, p_{j+1}^{e+1}\le n}f(p_{j+1}^e)\left(S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)-g(p_{j+1})\right)+f(p_{j+1}^{e+1})
这是因为当 pj+1e+1>np_{j+1}^{e+1}>n 的时候 s(… )−g(pj+1)s(\dots)-g(p_{j+1}) 这一项不会产生贡献。（其实方法一也可以这样简化式子。）
所以只需要设置 S(n,+∞)S(n,+\infty) 初值为 g(n)g(n)，即可进行和第一步几乎一样的 DP 过程了。
此外，注意这种方法不仅计算出了 S(n)S(n)，也顺带计算出了所有 S(⌊n/x⌋)S(\lfloor n/x\rfloor)，这是方法一做不到的。

时间复杂度

对于第一种方法，zzt 大佬在集训队论文中证明，一般情况下其复杂度将会是 O(n1−ϵ)O(n^{1-\epsilon})，其中 ϵ\epsilon 代表一个无穷小量。但是，他也证明了，在 n≤1013n \le 10^{13} 时，复杂度是 O(n34log⁡n)O\left(\frac{n^\frac{3}{4}}{\log n}\right)！
对于第二种方法，其复杂度和第一步类似；虽然要枚举质数的指数，但随 pp 的增大，指数的枚举次数下降速率急剧趋缓，因而不会对复杂度产生太大影响，视为常数（？） 听说其密度确实仍然是 O(1log⁡n)O(\frac{1}{\log n}) 级别。复杂度 O(n34log⁡n)O\left(\frac{n^\frac{3}{4}}{\log n}\right)。不过第一种方法常数更小，实际运行效可能会高那么一点点。

实现细节和代码

第一步中要使空间为 O(n)O(\sqrt n)，可以考虑根号分治表示下标。即，对于 x<nx< \sqrt n，用 xx 映射到下标；对于 x>nx>\sqrt n，用 n/xn/x 映射到下标。
P5325【模板】Min_25 筛
第一步把原函数拆成 g1(p)=p,g2(p)=p2g_1(p) = p,g_2(p)=p^2 两个完全积性函数计算。

方法一

第二步中直接递归计算，但是用了化简过的递推式。
CPP
#include <bits/stdc++.h>

using namespace std;

const int N = 2e5 + 5, mod = 1e9 + 7, i2 = 5e8 + 4, i6 = 166666668;
typedef long long ll;

ll n;
int id1[N], id2[N], m;
int p[N], vis[N], cnt;
ll g1[N], g2[N], v[N];

inline int get(ll x) { return x < N ? id1[x] : id2[n/x]; }
inline ll S1(ll x) { return x %= mod, x * (x + 1) % mod * i2 % mod; }
inline ll S2(ll x) { return x %= mod, x * (x + 1) % mod * (2 * x + 1) % mod * i6 % mod; }
inline ll sq(ll x) { return x %= mod, x * x % mod; }
inline ll F(ll x) { return x %= mod, (sq(x) - x + mod) % mod; }

inline void init(int n) {
	for (int i = 2; i <= n; i++) {
		if (!vis[i]) p[++cnt] = i;
		for (int j = 1; j <= cnt && p[j] <= n / i; j++) {
			vis[i*p[j]] = 1;
			if (i % p[j] == 0) break;
		}
	}
} 

ll S(ll x, int y) {
	if (p[y] >= x) return 0;
	ll res = (g2[get(x)] - g1[get(x)] - g2[get(p[y])] + g1[get(p[y])] + 2 * mod) % mod;
	for (int i = y + 1; i <= cnt && p[i] <= x / p[i]; i++) {
		ll w = p[i];
		for (int j = 1; w <= x / p[i]; j++, w = w * p[i])//注意后面有 x / w, 所以此处不能取模 
			res = (res + F(w) * S(x / w, i) % mod + F(w * p[i])) % mod;
	}
	return res;
}

int main()
{
	scanf("%lld", &n), init(sqrt(n) + 1);
	for (ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l), v[++m] = n / l;
		if (v[m] < N) id1[v[m]] = m;
		else id2[n/v[m]] = m;
		g1[m] = (S1(v[m]) - 1 + mod) % mod, g2[m] = (S2(v[m]) - 1 + mod) % mod;
	}
	for (int j = 1; j <= cnt; j++) {
		for (int i = 1; i <= m && p[j] <= v[i] / p[j]; i++) {
			g1[i] = (g1[i] - p[j] * (g1[get(v[i]/p[j])] - g1[get(p[j-1])]) % mod + mod) % mod;
			g2[i] = (g2[i] - sq(p[j]) * (g2[get(v[i]/p[j])] - g2[get(p[j-1])]) % mod + mod) % mod;
		}
	}
	printf("%lld\n", (S(n, 0) + 1) % mod);
	return 0;
}
评测记录

方法二

使用改进后的递推方式。
CPP
#include <bits/stdc++.h>

using namespace std;

const int N = 2e5 + 5, mod = 1e9 + 7, i2 = 5e8 + 4, i6 = 166666668;
typedef long long ll;

ll n;
int id1[N], id2[N], m;
int p[N], vis[N], cnt;
ll s[N], g1[N], g2[N], v[N];

inline int get(ll x) { return x < N ? id1[x] : id2[n/x]; }
inline ll S1(ll x) { return x %= mod, x * (x + 1) % mod * i2 % mod; }
inline ll S2(ll x) { return x %= mod, x * (x + 1) % mod * (2 * x + 1) % mod * i6 % mod; }
inline ll sq(ll x) { return x %= mod, x * x % mod; }
inline ll F(ll x) { return x %= mod, (sq(x) - x + mod) % mod; }

inline void init(int n) {
	for (int i = 2; i <= n; i++) {
		if (!vis[i]) p[++cnt] = i;
		for (int j = 1; j <= cnt && p[j] <= n / i; j++) {
			vis[i*p[j]] = 1;
			if (i % p[j] == 0) break;
		}
	}
} 

int main()
{
	scanf("%lld", &n), init(sqrt(n) + 1);
	for (ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l), v[++m] = n / l;
		if (v[m] < N) id1[v[m]] = m;
		else id2[n/v[m]] = m;
		g1[m] = (S1(v[m]) - 1 + mod) % mod, g2[m] = (S2(v[m]) - 1 + mod) % mod;
	}
	for (int j = 1; j <= cnt; j++) {
		for (int i = 1; i <= m && p[j] <= v[i] / p[j]; i++) {
			g1[i] = (g1[i] - p[j] * (g1[get(v[i]/p[j])] - g1[get(p[j-1])]) % mod + mod) % mod;
			g2[i] = (g2[i] - sq(p[j]) * (g2[get(v[i]/p[j])] - g2[get(p[j-1])]) % mod + mod) % mod;
		}
	}
	for (int i = 1; i <= m; i++) s[i] = g1[i] = (g2[i] - g1[i] + mod) % mod;
	for (int j = cnt; j >= 1; j--) {
		for (int i = 1; i <= m && p[j] <= v[i] / p[j]; i++) {
			ll w = p[j];
			for (int k = 1; w <= v[i] / p[j]; k++, w *= p[j])
				s[i] = (s[i] + F(w) * (s[get(v[i]/w)] - g1[get(p[j])] + mod) % mod + F(w * p[j])) % mod;
		}
	}
	printf("%lld\n", (s[get(n)] + 1) % mod);
	return 0;
}
评测记录
可以看出方法二实际效率确实劣于方法一。

例题

求质数个数，联想到 min_25 筛的第一步。相当于求 f(p)=1f(p)=1 这个函数在质数处取值的前缀和。
评测记录
容易发现两题本质相同，以第二道为例。
求合法的质数个数，容易联想到 min_25 筛的第一步。可以设 f(n,j,k)f(n,j,k) 表示前 nn 个数中，最小质因子 >pj>p_j，且  mod m=k\bmod m=k 的数的个数。转移的时候，设 pj mod m=cp_j\bmod m=c，那么枚举 k∈[0,m−1]k\in[0,m-1]，把 kk 的状态贡献到 kc mod mkc\bmod m 的状态即可。
评测记录
考虑函数质数处的取值。容易发现，f(p)={p+1p=2p−1p≠2f(p)=\begin{cases}p+1 &p=2\\p-1 &p\neq 2\end{cases}。
第一步中，应当直接把 f(p)f(p) 当成 p−1p-1，然后拆成 g1(p)=p,g0(p)=1g_1(p)=p,g_0(p)=1 计算，最后整合起来，并特判含 22 的前缀和。注意，我们应当保证拆成的若干个函数是完全积性函数，因此 g0(p)=−1g_0(p)=-1 这样的拆分是不可以的。
第二步没有特别之处，套用模板即可。
评测记录
观察发现，当 μ(n)=0\mu(n)=0 时，f(n)=0f(n)=0。当 n=1n=1 时 f(n)=1f(n)=1。接下来考虑 μ(n)≠0\mu(n)\neq 0 且 n>1n>1 的情况。
根据 μ\mu 的定义知，此时 nn 可以表示成 p1p2…pmp_1p_2\dots p_m。而 f(n)f(n) 取 +1+1 还是 −1-1 取决于 ∑d∣n[μ(d)=−1]\sum_{d|n}[\mu(d)=-1] 的奇偶性。我们发现，μ(d)=−1\mu(d)=-1 当且仅当 dd 中包含奇数个质因子，因此上式等价于 ∑i=2k+1,1≤i≤m(mi)\sum_{i=2k+1,1\le i \le m} {m\choose i}。由简单组合知识可知这个式子等于 2m−12^{m-1}。于是，当 m=1m=1 即 nn 为质数时，f(n)=−1f(n)=-1，否则 f(n)=1f(n)=1。
我们惊奇地发现 f(n)={μ2(n)n为合数−1n为质数f(n)=\begin{cases}\mu^2(n) & n\text{为合数}\\ -1 &n\text{为质数}\end{cases}。所以要求 ∑i=1nf(i)\sum_{i=1}^nf(i)，就是求 ∑i=1nμ2(i)−2⋅∑i=1n[i是质数]\sum_{i=1}^n\mu^2(i)-2\cdot\sum_{i=1}^n[i\text{是质数}]。
前一项是积性函数前缀和，可以 min_25 筛。后一项显然套用 min_25 筛的第一步即可。
评测记录
容易知道就是求 ∑i=1n(d(i)2)\sum_{i=1}^n{d(i)\choose 2}，其中 d(n)d(n) 表示 nn 的因数个数。展开变成 12(∑i=1nd2(i)−∑i=1nd(i))\frac{1}{2}\left(\sum_{i=1}^nd^2(i)-\sum_{i=1}^nd(i)\right)。
对于第一项，d2(n)d^2(n) 是积性函数，且 d2(p)=4,d2(pk)=(k+1)2d^2(p)=4,d^2(p^k)=(k+1)^2，考虑 min_25 筛。第一步要求函数是完全积性函数，所以提一个 44 出来，d′(p)=1d'(p)=1，最后再乘回去。
对于第二项，枚举因数：∑d=1n⌊nd⌋\sum_{d=1}^n\lfloor\frac{n}{d}\rfloor。数论分块即可。
评测记录
就是说要知道所有 S(⌊nx⌋)S(\lfloor\frac{n}{x}\rfloor) 的值。刚好 min_25 筛的第二种写法可以求出这些东西。
评测记录
不过 101310^{13} 的版本似乎不是为 min_25 筛准备的（或者我的写法常数太大？），反正我觉得很卡长，尽管有大佬的代码跑得飞快。最后卡了一晚上+一早上，终于离 TL 500ms500ms 左右通过了，但代码几乎惨不忍睹（。评测记录
有没有大佬教教我正解的啊（悲）

小结

好像重要的点在简介里都提到了。

参考资料

评论

共 39 条评论，欢迎与作者交流。

正在加载评论...