Min_25 筛

# 简介 对于**积性函数** $f(x)$，如果能快速知道 $f(p^k)$ 的值（其中 $p$ 为质数，下同），且 $f(p)$ 可以表示为项数较少的**多项式**（或能表示成若干项使得每一项都是**完全积性函数**），那么 $\rm min\_25$ 筛可以在 $O(\frac{n^{\frac{3}{4}}}{\log n})$ 时间复杂度内求出 $\sum_{i=1}^nf(i)$。 由于其思想与埃氏筛类似，所以也被称作扩展埃氏筛。 # 内容 ## 第一步 ### 推导 目标：求 $g(n)=\sum_{p\le n} f(p)$ 。 不妨设 $f(p)$ 是完全积性函数，如果不是可以尝试拆成若干项完全积性函数，分别求然后相加。 首先要线性筛求出 $\sqrt n$ 以内的质数。 $g(n)$ 很难直接求解，考虑用 DP 计算。设 $g(n,j)=\sum_{i=1}^nf(i)[i\text{是质数或其最小质因子}>p_j]$，其中 $p_j$ 表示第 $j$ 个质数，那么我们要的就是 $g(n,k),k\text{为最小的满足}p_k\ge \sqrt n$。考虑从 $j-1$ 变到 $j$ ，那么最小质因子为 $p_j$ 的合数会被筛掉，那么它们的贡献要减去。则有转移 $$ 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(p_j)$ 表示由于 $f(p)$ 是完全积性函数，所以可以把它从后面提出来。 $g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right)$ 表示考虑所有 $p_j$ 的倍数，它们除以 $p_j$ 之后，最小质因子 $>p_{j-1}$ 的合数**以及所有质数**的贡献，应当减去。但是，这些**质数**中可能有 $\le p_{j-1}$ 的，它们在之前就被筛掉过了，所以要加回来，也就是 $g(p_{j-1},j-1)$。 由于有公式 $\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\rfloor=\lfloor\frac{a}{bc}\rfloor$ ，因此容易发现上述式子只会用到形如 $\lfloor\frac{n}{x}\rfloor,x\le n$ 的点处的 DP 值，即第一项的状态数是 $O(\sqrt n)$（实际实现的时候注意状态数是 $2\sqrt n$）。我们预处理出这 $O(\sqrt n)$ 个数，把他们离散化，顺带求出 $g(x,0)$，然后 DP 即可。 注意到上述转移式中还有一项 $g(p_{j-1},j-1)$，我们只处理了所有形如 $\lfloor\frac{n}{x}\rfloor$ 的数，是否会漏掉某些质数 $p$ 呢？其实不会漏。注意到，能用来转移的 $p$ 一定满足 $p \le \sqrt n$。我们只需要证明 $\forall k \le \sqrt n, \exists x,k=\lfloor\frac{n}{x}\rfloor$。 > 证明： > > 1. 若 $k=\lfloor\sqrt n\rfloor$，设 $n=k^2+d$，则由于 $k^2\le n < (k+1)^2$，故 $d\in[0,2k]$。 > 1. $d\in [0,k)$，那么 $\lfloor\frac{n}{k}\rfloor=k+\lfloor\frac{d}{k}\rfloor=k$。 > 2. $d\in [k,2k]$，那么 $\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 < \lfloor\sqrt n\rfloor$ 即 $k \le \sqrt n - 1$，假设存在 $i,\lfloor\frac{n}{i+1}\rfloor<k<\lfloor\frac{n}{i}\rfloor$，此时 $k$ 恰好夹在两个连续的 $\lfloor\frac{n}{x}\rfloor$ 之间，即不可被表出。则 $\frac{n}{i+1}<k$，故 $n<k(i+1)$，从而 $k<\lfloor\frac{n}{i}\rfloor<\lfloor\frac{k(i+1)}{i}\rfloor=k+\lfloor\frac{k}{i}\rfloor$。另一方面，$\frac{n}{i+1}<k<\sqrt n$，所以 $i+1>\sqrt n$，于是 $i>\sqrt n - 1 \ge k$，因此 $\lfloor\frac{k}{i}\rfloor=0$，于是得到 $k<k$，故假设不成立，原命题成立。 所以，上述担心就是多虑了。 > 我看到很多博客、题解都额外预处理了 $\sqrt n$ 内所有素数处的函数值的前缀和，但因为我“粗心”，代码中直接调用 $g$ 数组，却顺利通过了一堆题目，写笔记时意识到这一点，于是琢磨出了奥妙所在。其实，在看 zzt 集训队论文时注意到他说只需要用到 $\lfloor\frac{n}{x}\rfloor$ 处的值，就心生疑惑，现在终于明白了。因此，在我看来，**预处理根号内素数处的函数前缀和的值是不必要的**。 ### 时间复杂度 $$ \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)=\sum_{i\le n} f(i)$。与第一步类似，设 $S(n,j)=\sum_{i=1}^nf(i)[i\text{的最小质因子}>p_j]$。但此处 $f$ 不需要再拆分成单项式，直接是原函数即可（因为不需要依赖于**完全积性**，只需要**积性**即可）（但要能快速计算 $f(p^k)$ 的值）。 #### 方法一 考虑把贡献拆成质数的和合数的，合数枚举最小质因子以及次数，于是有转移： $$ 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\neq 1]$ 的意思是，对于 $e=1$ 的情况，$S$ 没有计算 $1$ 贡献，刚好，因为此时 $p_k\times 1$ 是质数，其贡献在之前计算过；对于 $e > 1$ 的情况，$p_k^e \times 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) $$ 但直接这样转移复杂度不对，我们要严格控制只有 $\le \sqrt n$ 的质数才发生转移。即，对于 $p_{j+1}>\sqrt n$ 的状态 $S(n,j)$，不显式地计算出来，用到的时候特判为 $g(n)-g(p_j)$。所以，我们要更新状态 $S(n,j)$ 时，一定有 $p_{j+1}\le \sqrt n$，此时需要用到 $S(n,j+1)$，即使要特判为 $g(n)-g(p_{j+1})$，由步骤一中的论述可知 $g(p_{j+1})$ 已经计算出，所以不会有问题。需要注意 $n=2,3$ 的状态不会被任何 $p$ 更新到，所以也要特判。一种较为简洁的方法是，像方法一的代码一样写一个函数 $S(x,y)$，只不过不用来递归，而是用来实现各种特判。 **但是，这样代码不太好些，且实现出来常数比较大**。我实现了一份：[code](https://loj.ac/s/1344247)。耗时几乎是别人的两倍。 #### 方法二改进 既然采用了类似第一步的递推方式，干脆直接套用第一步的状态设计。设 $S(n,j)=\sum_{i=1}^nf(i)[i\text{是质数或其最小质因子}>p_j]$ ，那么参考第一步的方程，有转移： $$ 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(\dots)$ 是因为，$S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)$ 中包含了 $\le p_{j+1}$ 的质数的贡献，与转移式的意义不符，应当减去。但是 $g$ 里面有个 $\min$ 非常令人不爽，~~并且为了卡常~~，我们可以进一步简化上式： $$ 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}) $$ 这是因为当 $p_{j+1}^{e+1}>n$ 的时候 $s(\dots)-g(p_{j+1})$ 这一项不会产生贡献。（其实方法一也可以这样简化式子。） 所以只需要设置 $S(n,+\infty)$ 初值为 $g(n)$，即可进行和第一步几乎一样的 DP 过程了。 **此外，注意这种方法不仅计算出了 $S(n)$，也顺带计算出了所有 $S(\lfloor n/x\rfloor)$，这是方法一做不到的。** ### 时间复杂度 对于第一种方法，zzt 大佬在集训队论文中证明，一般情况下其复杂度将会是 $O(n^{1-\epsilon})$，其中 $\epsilon$ 代表一个无穷小量。但是，他也证明了，**在 $n \le 10^{13}$ 时，复杂度是 $O\left(\frac{n^\frac{3}{4}}{\log n}\right)$！** 对于第二种方法，其复杂度和第一步类似；虽然要枚举质数的指数，但随 $p$ 的增大，指数的枚举次数下降速率急剧趋缓，因而不会对复杂度产生太大影响，~~视为常数（？）~~ 听说其密度确实仍然是 $O(\frac{1}{\log n})$ 级别。复杂度 $O\left(\frac{n^\frac{3}{4}}{\log n}\right)$。不过第一种方法常数更小，实际运行效可能会高那么一点点。 ## 实现细节和代码 第一步中要使空间为 $O(\sqrt n)$，可以考虑根号分治表示下标。即，对于 $x< \sqrt n$，用 $x$ 映射到下标；对于 $x>\sqrt n$，用 $n/x$ 映射到下标。 [P5325【模板】Min\_25 筛](https://www.luogu.com.cn/problem/P5325) 第一步把原函数拆成 $g_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; } ``` [评测记录](https://www.luogu.com.cn/record/66603549) ### 方法二 使用改进后的递推方式。 ```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; } ``` [评测记录](https://www.luogu.com.cn/record/66603825) 可以看出方法二实际效率确实劣于方法一。 # 例题 - [loj6235区间素数个数](https://loj.ac/p/6235) 求质数个数，联想到 min\_25 筛的第一步。相当于求 $f(p)=1$ 这个函数在质数处取值的前缀和。 [评测记录](https://loj.ac/s/1343058) - [loj6027质数计数 I](https://loj.ac/p/6027) / [loj6028质数计数 II](https://loj.ac/p/6028) 容易发现两题本质相同，以第二道为例。 求合法的质数个数，容易联想到 min\_25 筛的第一步。可以设 $f(n,j,k)$ 表示前 $n$ 个数中，最小质因子 $>p_j$，且 $\bmod m=k$ 的数的个数。转移的时候，设 $p_j\bmod m=c$，那么枚举 $k\in[0,m-1]$，把 $k$ 的状态贡献到 $kc\bmod m$ 的状态即可。 [评测记录](https://loj.ac/s/1343604) - [loj6053简单的函数](https://loj.ac/p/6053) 考虑函数质数处的取值。容易发现，$f(p)=\begin{cases}p+1 &p=2\\p-1 &p\neq 2\end{cases}$。 第一步中，应当直接把 $f(p)$ 当成 $p-1$，然后拆成 $g_1(p)=p,g_0(p)=1$ 计算，最后整合起来，并特判含 $2$ 的前缀和。**注意，我们应当保证拆成的若干个函数是完全积性函数，因此 $g_0(p)=-1$ 这样的拆分是不可以的。** 第二步没有特别之处，套用模板即可。 [评测记录](https://loj.ac/s/1344737) - [loj6181某个套路求和题](https://loj.ac/p/6181) 观察发现，当 $\mu(n)=0$ 时，$f(n)=0$。当 $n=1$ 时 $f(n)=1$。接下来考虑 $\mu(n)\neq 0$ 且 $n>1$ 的情况。 根据 $\mu$ 的定义知，此时 $n$ 可以表示成 $p_1p_2\dots p_m$。而 $f(n)$ 取 $+1$ 还是 $-1$ 取决于 $\sum_{d|n}[\mu(d)=-1]$ 的奇偶性。我们发现，$\mu(d)=-1$ 当且仅当 $d$ 中包含奇数个质因子，因此上式等价于 $\sum_{i=2k+1,1\le i \le m} {m\choose i}$。由简单组合知识可知这个式子等于 $2^{m-1}$。于是，当 $m=1$ 即 $n$ 为质数时，$f(n)=-1$，否则 $f(n)=1$。 我们惊奇地发现 $f(n)=\begin{cases}\mu^2(n) & n\text{为合数}\\ -1 &n\text{为质数}\end{cases}$。所以要求 $\sum_{i=1}^nf(i)$，就是求 $\sum_{i=1}^n\mu^2(i)-2\cdot\sum_{i=1}^n[i\text{是质数}]$。 前一项是积性函数前缀和，可以 min\_25 筛。后一项显然套用 min\_25 筛的第一步即可。 [评测记录](https://loj.ac/s/1343555) - [loj6682梦中的数论](https://loj.ac/p/6682) 容易知道就是求 $\sum_{i=1}^n{d(i)\choose 2}$，其中 $d(n)$ 表示 $n$ 的因数个数。展开变成 $\frac{1}{2}\left(\sum_{i=1}^nd^2(i)-\sum_{i=1}^nd(i)\right)$。 对于第一项，$d^2(n)$ 是积性函数，且 $d^2(p)=4,d^2(p^k)=(k+1)^2$，考虑 min\_25 筛。第一步要求函数是完全积性函数，所以提一个 $4$ 出来，$d'(p)=1$，最后再乘回去。 对于第二项，枚举因数：$\sum_{d=1}^n\lfloor\frac{n}{d}\rfloor$。数论分块即可。 [评测记录](https://loj.ac/s/1343773) - [loj6783/6784/6785简单的函数 的一堆加强版](https://loj.ac/p/6784) 就是说要知道所有 $S(\lfloor\frac{n}{x}\rfloor)$ 的值。刚好 min\_25 筛的第二种写法可以求出这些东西。 [评测记录](https://loj.ac/s/1344627) 不过 $10^{13}$ 的版本似乎不是为 min\_25 筛准备的（或者我的写法常数太大？），反正我觉得很卡长，尽管有大佬的代码跑得飞快。最后卡了一晚上+一早上，终于离 TL $500ms$ 左右通过了，但代码几乎惨不忍睹（。[评测记录](https://loj.ac/s/1344689) 有没有大佬教教我正解的啊（悲） # 小结 好像重要的点在简介里都提到了。 # 参考资料 - <https://www.luogu.com.cn/blog/wucstdio/solution-p5325> - <https://oi-wiki.org/math/number-theory/min-25/> - $2018$ 集训队论文：朱震霆 《一些特殊的数论函数求和问题 》 - $2016$ 集训队论文：任之洲《积性函数求和的几种方法》