Min_25 筛
由于其由 Min_25 发明并最早开始使用,故称「Min_25 筛」。
从此种筛法的思想方法来说,其又被称为「Extended Eratosthenes Sieve」。
其可以在 O\left(\frac{n^{\frac{3}{4}}}{\log{n}}\right) 或 \Theta\left(n^{1 - \epsilon}\right) 的时间复杂度下解决一类 积性函数 的前缀和问题。
要求: f(p) 是一个关于 p 的项数较少的多项式或可以快速求值; f(p^{c}) 可以快速求值。
记号¶
- 如无特别说明,本节中所有记为 p 的变量的取值集合均为全体质数。
- x / y := \left\lfloor\frac{x}{y}\right\rfloor
- \operatorname{isprime}(n) := [ |\{d : d \mid n\}| = 2 ] ,即 n 为质数时其值为 1 ,否则为 0 。
- p_{k} :全体质数中第 k 小的质数(如: p_{1} = 2, p_{2} = 3 )。特别地,令 p_{0} = 1 。
- \operatorname{lpf}(n) := [1 < n] \min\{p : p \mid n\} + [1 = n] ,即 n 的最小质因数。特别地, n=1 时,其值为 1 。
- F_{\mathrm{prime}}(n) := \sum_{2 \le p \le n} f(p)
- F_{k}(n) := \sum_{i = 2}^{n} [p_{k} \le \operatorname{lpf}(i)] f(i)
具体方法¶
观察 F_{k}(n) 的定义,可以发现答案即为 F_{1}(n) + f(1) = F_{1}(n) + 1 。
考虑如何求出 F_{k}(n) 。通过枚举每个 i 的最小质因子及其次数可以得到递推式:
最后一步推导基于这样一个事实:对于满足 p_{i}^{c} \le n < p_{i}^{c + 1} 的 c ,有 p_{i}^{c + 1} > n \iff n / p_{i}^{c} < p_{i} < p_{i + 1} ,故 F_{i + 1}\left(n / p_{i}^{c}\right) = 0 。
其边界值即为 F_{k}(n) = 0 (p_{k} > n) 。
假设现在已经求出了所有的 F_{\mathrm{prime}}(n) ,那么有两种方式可以求出所有的 F_{k}(n) :
- 直接按照递推式计算。
- 从大到小枚举 p 转移,仅当 p^{2} < n 时转移增加值不为零,故按照递推式后缀和优化即可。
现在考虑如何计算 F_{\mathrm{prime}}{(n)} 。
观察求 F_{k}(n) 的过程,容易发现 F_{\mathrm{prime}} 有且仅有 1, 2, \dots, \left\lfloor\sqrt{n}\right\rfloor, n / \sqrt{n}, \dots, n / 2, n 这 O(\sqrt{n}) 处的点值是有用的。
一般情况下, f(p) 是一个关于 p 的低次多项式,可以表示为 f(p) = \sum a_{i} p^{c_{i}} 。
那么对于每个 p^{c_{i}} ,其对 F_{\mathrm{prime}}(n) 的贡献即为 a_{i} \sum_{2 \le p \le n} p^{c_{i}} 。
分开考虑每个 p^{c_{i}} 的贡献,问题就转变为了:给定 n, s, g(p) = p^{s} ,对所有的 m = n / i ,求 \sum_{p \le m} g(p) 。
Notice: g(p) = p^{s} 是完全积性函数!
于是设 G_{k}(n) := \sum_{i = 1}^{n} \left[p_{k} < \operatorname{lpf}(i) \lor \operatorname{isprime}(i)\right] g(i) ,即埃筛第 k 轮筛完后剩下的数的 g 值之和。
对于一个合数 x ,必定有 \operatorname{lpf}(x) \le \sqrt{x} ,则 F_{\mathrm{prime}} = G_{\left\lfloor\sqrt{n}\right\rfloor} ,故只需筛到 G_{\left\lfloor\sqrt{n}\right\rfloor} 即可。
考虑 G 的边界值,显然为 G_{0}(n) = \sum_{i = 2}^{n} g(i) 。(还记得吗?特别约定了 p_{0} = 1 )
对于转移,考虑埃筛的过程,分开讨论每部分的贡献,有:
- 对于 n < p_{k}^{2} 的部分, G 值不变,即 G_{k}(n) = G_{k - 1}(n) 。
- 对于 p_{k}^{2} \le n 的部分,被筛掉的数必有质因子 p_{k} ,即 -g(p_{k}) G_{k - 1}(n / p_{k}) 。
- 对于第二部分,由于 p_{k}^{2} \le n \iff p_{k} \le n / p_{k} ,故会有 \operatorname{lpf}(i) < p_{k} 的 i 被减去。这部分应当加回来,即 g(p_{k}) G_{k - 1}(p_{k - 1}) 。
则有:
复杂度分析¶
对于 F_{k}(n) 的计算,其第一种方法的时间复杂度被证明为 O\left(n^{1 - \epsilon}\right) (见 zzt 集训队论文 2.3);
对于第二种方法,其本质即为洲阁筛的第二部分,在洲阁论文中也有提及(6.5.4),其时间复杂度被证明为 O\left(\frac{n^{\frac{3}{4}}}{\log{n}}\right) 。
对于 F_{\mathrm{prime}}(n) 的计算,事实上,其实现与洲阁筛第一部分是相同的。
考虑对于每个 m = n / i ,只有在枚举满足 p_{k}^{2} \le m 的 p_{k} 转移时会对时间复杂度产生贡献,则时间复杂度可估计为:
对于空间复杂度,可以发现不论是 F_{k} 还是 F_{\mathrm{prime}} ,其均只在 n / i 处取有效点值,共 O(\sqrt{n}) 个。
则可以使用 杜教筛一节中介绍的 trick 来将空间复杂度优化至 O(\sqrt{n}) 。
有关代码实现¶
对于 F_{k}(n) 的计算,我们实现时一般选择实现难度较低的第一种方法,其在数据规模较小时往往比第二种方法的表现要好;
对于 F_{\mathrm{prime}}(n) 的计算,直接按递推式实现即可。
对于 p_{k}^{2} \le n ,可以用线性筛预处理出 s_{k} := F_{\mathrm{prime}}(p_{k}) 来替代 F_{k} 递推式中的 F_{\mathrm{prime}}(p_{k - 1}) 。
相应地, G 递推式中的 G_{k - 1}(p_{k - 1}) = \sum_{i = 1}^{k - 1} g(p_{i}) 也可以用此方法预处理。
用 Extended Eratosthenes Sieve 求 积性函数 f 的前缀和时,应当明确以下几点:
- 如何快速(一般是线性时间复杂度)筛出前 \sqrt{n} 个 f 值;
- f(p) 的多项式表示;
- 如何快速求出 f(p^{c}) 。
明确上述几点之后按顺序实现以下几部分即可:
- 筛出 [1, \sqrt{n}] 内的质数与前 \sqrt{n} 个 f 值;
- 对 f(p) 多项式表示中的每一项筛出对应的 G ,合并得到 F_{\mathrm{prime}} 的所有 O(\sqrt{n}) 个有用点值;
- 按照 F_{k} 的递推式实现递归,求出 F_{1}(n) 。
例题¶
求莫比乌斯函数的前缀和¶
求 \displaystyle \sum_{i = 1}^{n} \mu(i) 。
易知 f(p) = -1 。则 g(p) = -1, G_{0}(n) = \sum_{i = 2}^{n} g(i) = -n + 1 。
直接筛即可得到 F_{\mathrm{prime}} 的所有 O(\sqrt{n}) 个所需点值。
求欧拉函数的前缀和¶
求 \displaystyle \sum_{i = 1}^{n} \varphi(i) 。
首先易知 f(p) = p - 1 。
对于 f(p) 的一次项 (p) ,有 g(p) = p, G_{0}(n) = \sum_{i = 2}^{n} g(i) = \frac{(n + 2) (n - 1)}{2} ;
对于 f(p) 的常数项 (-1) ,有 g(p) = -1, G_{0}(n) = \sum_{i = 2}^{n} g(i) = -n + 1 。
筛两次加起来即可得到 F_{\mathrm{prime}} 的所有 O(\sqrt{n}) 个所需点值。
「LOJ #6053」简单的函数¶
给定 f(n) :
易知 f(p) = p - 1 + 2[p = 2] 。则按照筛 \varphi 的方法筛,对 2 讨论一下即可。
此处给出一种 C++ 实现:
参考代码
/* 「LOJ #6053」简单的函数 */
#include <algorithm>
#include <cmath>
#include <cstdio>
using i64 = long long;
constexpr int maxs = 200000; // 2sqrt(n)
constexpr int mod = 1000000007;
template <typename x_t, typename y_t>
inline void inc(x_t &x, const y_t &y) {
x += y;
(mod <= x) && (x -= mod);
}
template <typename x_t, typename y_t>
inline void dec(x_t &x, const y_t &y) {
x -= y;
(x < 0) && (x += mod);
}
template <typename x_t, typename y_t>
inline int sum(const x_t &x, const y_t &y) {
return x + y < mod ? x + y : (x + y - mod);
}
template <typename x_t, typename y_t>
inline int sub(const x_t &x, const y_t &y) {
return x < y ? x - y + mod : (x - y);
}
template <typename _Tp>
inline int div2(const _Tp &x) {
return ((x & 1) ? x + mod : x) >> 1;
}
template <typename _Tp>
inline i64 sqrll(const _Tp &x) {
return (i64)x * x;
}
int pri[maxs / 7], lpf[maxs + 1], spri[maxs + 1], pcnt;
inline void sieve(const int &n) {
for (int i = 2; i <= n; ++i) {
if (lpf[i] == 0)
pri[lpf[i] = ++pcnt] = i, spri[pcnt] = sum(spri[pcnt - 1], i);
for (int j = 1, v; j <= lpf[i] && (v = i * pri[j]) <= n; ++j) lpf[v] = j;
}
}
i64 global_n;
int lim;
int le[maxs + 1], // x \le \sqrt{n}
ge[maxs + 1]; // x > \sqrt{n}
#define idx(v) (v <= lim ? le[v] : ge[global_n / v])
int G[maxs + 1][2], Fprime[maxs + 1];
i64 lis[maxs + 1];
int cnt;
inline void init(const i64 &n) {
for (i64 i = 1, j, v; i <= n; i = n / j + 1) {
j = n / i;
v = j % mod;
lis[++cnt] = j;
idx(j) = cnt;
G[cnt][0] = sub(v, 1ll);
G[cnt][1] = div2((i64)(v + 2ll) * (v - 1ll) % mod);
}
}
inline void calcFprime() {
for (int k = 1; k <= pcnt; ++k) {
const int p = pri[k];
const i64 sqrp = sqrll(p);
for (int i = 1; lis[i] >= sqrp; ++i) {
const i64 v = lis[i] / p;
const int id = idx(v);
dec(G[i][0], sub(G[id][0], k - 1));
dec(G[i][1], (i64)p * sub(G[id][1], spri[k - 1]) % mod);
}
}
/* F_prime = G_1 - G_0 */
for (int i = 1; i <= cnt; ++i) Fprime[i] = sub(G[i][1], G[i][0]);
}
inline int f_p(const int &p, const int &c) {
/* f(p^{c}) = p xor c */
return p xor c;
}
int F(const int &k, const i64 &n) {
if (n < pri[k] || n <= 1) return 0;
const int id = idx(n);
i64 ans = Fprime[id] - (spri[k - 1] - (k - 1));
if (k == 1) ans += 2;
for (int i = k; i <= pcnt && sqrll(pri[i]) <= n; ++i) {
i64 pw = pri[i], pw2 = sqrll(pw);
for (int c = 1; pw2 <= n; ++c, pw = pw2, pw2 *= pri[i])
ans +=
((i64)f_p(pri[i], c) * F(i + 1, n / pw) + f_p(pri[i], c + 1)) % mod;
}
return ans % mod;
}
int main() {
scanf("%lld", &global_n);
lim = sqrt(global_n);
sieve(lim + 1000);
init(global_n);
calcFprime();
printf("%lld\n", (F(1, global_n) + 1ll + mod) % mod);
return 0;
}
build本页面最近更新:,更新历史
edit发现错误?想一起完善? 在 GitHub 上编辑此页!
people本页面贡献者:TrisolarisHD, Xeonacid
copyright本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供,附加条款亦可能应用