Min_25 筛 学习笔记
zhiyangfan
·
·
个人记录
Min_25 筛
Min_25 筛是一种亚线性筛法，可以在 \mathcal{O}(\frac{n^{\frac{3}{4}}}{\log n}) 的时间复杂度下快速算出形如：
\sum_{i=1}^n f(i)
的值，不过一般比较好实现的方法被证明复杂度是 \mathcal{O}(n^{1-\epsilon})。
本文约定 \mathbb{P} 表示质数集，\mathbb{P}' 表示所有 \le \sqrt{n} 的质数集。
对于函数的限制
在了解一个筛法的具体内容前，先了解其适用范围是有助于推导过程时的理解和应用时不至于用错的。对于 Min_25 筛，其要求求和的 f(i) 满足如下性质：
- 积性函数，即对于任意的 \gcd(x,y)=1，都有 f(x)f(y)=f(xy)。
-
-
接下来，我们就找一个满足以上条件的函数并以它为例子进行介绍吧。
P5325 【模板】Min_25筛
定义积性函数 f(x)，其中 f(p^k)=p^k(p^k-1),p\in\mathbb{P}，求：
\sum_{i=1}^nf(i)\bmod{10^9+7}
(1\le n\le 10^{10})
对于和式的恰当拆分
首先我们注意到许多性质是针对质数来的，所以我们考虑把质数从求和式中单独拉出来：
\sum_{i=1}^nf(i)=\sum_{p=1}^n[p\in\mathbb{P}]f(p)+\sum_{i=1}^nf(i)[p\notin\mathbb{P}\operatorname{or}p=1]
但对于后面的和式怎么办呢？我们想到因为是积性函数，所以合数可以提取出来质因数与质数扯上关系，枚举合数的 \rm LPF（\rm lowest\ prime\ factor，最小质因数）及其次数有：
\sum_{p=1}^n[p\in\mathbb{P}]f(p)+\sum_{p=1}^{\sqrt{n}}[1\le p^e\le n\operatorname{and}p\in\mathbb{P}]f(p^e)\left(\sum_{i=1}^{\lfloor\frac{n}{p^e}\rfloor}[\mathrm{LPF}_i>p]f(i)\right)
提取出来 p^e 后如果原数满足 \mathrm{LPF}=p 则剩下的数应该满足 \mathrm{LPF}_i>p。
引入递推思想
（不知道 Min_25 神仙咋想到的啊Orz）考虑递推：设 g(n,j) 表示：
\sum_{i=1}^n i^k[i\in\mathbb{P}\operatorname{or}\mathrm{LPF}_i>p_j]
这里 i^k 是把 f(p) 作为多项式拆分成单项式的结果，所以如果 f(p) 有 l 项，那我们就要求 l 个 g(n,j)，并用它们的线性组合组合出 f(p),p\in\mathbb{P}。
接下来的问题就是为 $g(n,j)$ 找到一个递归式，也即如何实现 $g(n,j-1)\rightarrow g(n,j)$ 的转移。考虑这个转移的实质其实是把限制拉大了，所以一定是去掉一些数，而去掉的数应该满足 $\mathrm{LPF}=p_j$。既然满足 $\mathrm{LPF}=p_j$，我们就试着提取一个 $p_j$ 出来：
$$g(n,j)=g(n,j-1)-p_j^k\left(g\left(\dfrac{n}{p_j},j-1\right)-g(p_{j-1},j-1)\right)$$
解释一下这个式子。首先要发现，$i^k$ 是完全积性函数，即 $\forall x,y$ 都有 $f(x)f(y)=f(xy)$，所以即使 $\dfrac{n}{p_j}$ 不一定与 $p_j$ 互质也能提取（而上文 $f$ 仅为积性函数，所以必须一口气提取完所有的 $p$ 才能保证互质）。然后分别看一下括号里面的两项。$g\left(\left\lfloor\dfrac{n}{p_j}\right\rfloor,j-1\right)$ 关于的是提取完 $p_j$ 后的数，而一个数如果想 $\mathrm{LPF}=p_j$，提取完 $p_j$ 后其 $\rm LPF$ 应该 $> p_{j-1}$。但这有个问题，$\le p_{j-1}$ 的质数会被统计到这里面，但是它们并不满足 $\mathrm{LPF}>p_{j-1}$ 所以要减去。$g(p_{j-1},j-1)$ 就是干这个的，显然 $\le p_{j-1}$ 的数不可能满足 $\mathrm{LPF}>p_{j-1}$，所以这个统计的就是 $\le p_{j-1}$ 的质数。
### 求递归式
我们刚刚求出了关于 $g(n,j)$ 的递归式，但发现求 $g(n,j)$ 还得用到 $g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right)$，求 $g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right)$ 又得用到 $g\left(\left\lfloor\frac{\left\lfloor\frac{n}{p_j}\right\rfloor}{p_k}\right\rfloor,j-2\right)$ 等等。以此类推，难道我们对于所有 $1\sim n$ 都求出 $g(n,j)$ 吗？显然复杂度不能接受。但我们发现上述过程中出现了一个很熟悉的结构，有结论：
$$\left\lfloor\dfrac{\left\lfloor\dfrac{n}{a}\right\rfloor}{b}\right\rfloor=\left\lfloor\dfrac{n}{ab}\right\rfloor$$
那容易发现，不管递推多少次，我们需要的值总能被表示为 $g\left(\left\lfloor\dfrac{n}{x}\right\rfloor,j\right),x\in\mathbb{N}_+$ 的形式，我们只需要对 $\left\lfloor\dfrac{n}{x}\right\rfloor$ 求 $g$ 就好了，而它仅有 $\mathcal{O}(\sqrt{n})$ 种，可以接受。
有一个实现上的小问题，我们在套上述递归式计算的时候，肯定要给每个 $\left\lfloor\dfrac{n}{x}\right\rfloor$ 一个编号，但在递归的时候我们需要实现从 $\left\lfloor\dfrac{n}{x}\right\rfloor$ 到其编号的转化，而前者的值域是 $[0,n]$，如果用 `std::map` 的话复杂度会多一个 $\log$。不过也不难实现 $\mathcal{O}(1)$ 的转化，考虑对于所有 $\le \sqrt{n}$ 的 $\left\lfloor\dfrac{n}{x}\right\rfloor$ 开一个桶记录 $\left\lfloor\dfrac{n}{x}\right\rfloor\rightarrow \mathrm{index}_{\lfloor\frac{n}{x}\rfloor}$ 的关系，对于 $>\sqrt{n}$ 的 $\left\lfloor\dfrac{n}{x}\right\rfloor$ 用另一个桶记录 $\left\lfloor\dfrac{n}{\left\lfloor\dfrac{n}{x}\right\rfloor}\right\rfloor\rightarrow \mathrm{index}_{\lfloor\frac{n}{x}\rfloor}$ 的关系。这样键值的值域就保持在 $[0,\sqrt{n}]$ 了，可以数组存。
### 另一个递归式
考虑设 $S(n,x)$ 表示：
$$\sum_{i=1}^n f(i)[\mathrm{LPF}_i>p_x]$$
说人话就是所有 $\mathrm{LPF}>p_x$ 的数求和，则显然答案就是 $S(n,0)+1$（因为一开始求 $g$ 的时候把 $1$ 排除在外了）。接下来考虑怎么计算它。注意到符合条件的数有两种：
- $> p_x$ 的质数。
- $\mathrm{LPF}_i>p_x$ 的合数。
首先是 $> p_x$ 的质数，这可以通过我们刚刚算出的 $g$ 求，即：
$$g(n,|\mathbb{P}'|)-sp_x$$
其中 $sp_x$ 表示 $\sum_{p=1}^{p_x}p^k[p\in\mathbb{P}]$，可以通过线性筛在 $\mathcal{O}(\sqrt{n})$ 的时间前缀和出来。
然后是 $\mathrm{LPF}_i>p_x$ 的合数，这时候我们考虑枚举 $\rm LPF$ 及其次数，得到：
$$\sum_{p_k^e\le n\operatorname{and}k>x}f(p_k^e)\left(S\left(\dfrac{n}{p^e_k},k\right)+[e\ne 1]\right)$$
这里的提取跟最开始拆分 $g$ 的思路类似，只不过这次因为不是完全积性函数是全部提取完的，$\rm LPF$ 应该要 $>p_k$。然后是 $[e\ne 1]$，显然 $f(p_k^e)$ 应该算上，但这样提取是不会被 $S$ 包含的，应该 $+1$ 把 $f(p_k^e)$ 包含进去。但当 $e=1$ 的时候，$f(p_k^e)$ 是质数，已经被算过了，不能加上。
这样把上面俩式子加一起就可以算了，这里采用递归的形式实现最方便，由于某个我并不知道的定理，可以不用记忆化也能保证复杂度。最后记得 $+f(1)$ 把之前没算上的 $1$ 给补上。
### 总结
Min_25 筛的基本步骤：
1. 将待求函数在质数情况下的表达式转化成 $\sum i^k$ 的形式。
1. 线性筛出 $1\sim \sqrt{n}$ 之间的质数，并算出其 $k$ 次方前缀和。
1. 对于形如 $\lfloor\frac{n}{x}\rfloor$ 的数算出 $g(\lfloor\frac{n}{x}\rfloor,0)$，注意去掉 $1^k$。
1. 套递归式算 $g(n,j)$，这里 $g$ 可以用滚动数组.
1. 递归计算 $S(n,x)$，无需记忆化。答案即为 $S(n,0)+f(1)$。
板子题代码：
```cpp
#include <cmath>
#include <cstdio>
const int N = 1e6 + 10, mod = 1e9 + 7, inv3 = 333333336, inv2 = 500000004;
typedef long long ll; ll n, w[N]; int p[N], vis[N], sp1[N], sp2[N], tp, sqrn;
// 线性筛出 <= \sqrt{n} 的质数，并求出 p^k 的前缀和
inline void getP(int m)
{
vis[1] = 1;
for (int i = 2; i <= m; ++i)
{
if (!vis[i]) p[++tp] = i;
for (int j = 1; j <= tp && (ll)p[j] * i <= m; ++j)
{
vis[i * p[j]] = 1;
if (i % p[j] == 0) break;
}
}
for (int i = 1; i <= tp; ++i) sp1[i] = (sp1[i - 1] + p[i]) % mod;
for (int i = 1; i <= tp; ++i) sp2[i] = (sp2[i - 1] + (ll)p[i] * p[i] % mod) % mod;
}
int ind1[N], ind2[N], tot, g1[N], g2[N];
int S(ll x, int y)
{
if (p[y] >= x) return 0;
int k = x <= sqrn ? ind1[x] : ind2[n / x];
int ans = ((ll)g2[k] - g1[k] - sp2[y] + sp1[y] + mod + mod) % mod;
for (int i = y + 1; i <= tp && (ll)p[i] * p[i] <= x; ++i)
{
ll pe = p[i];
for (int e = 1; pe <= x; ++e, pe *= p[i])
{
int xx = pe % mod;
ans = (ans + (ll)xx * (xx - 1) % mod * (S(x / pe, i) + (e != 1)) % mod) % mod;
}
}
return ans;
}
int main()
{
scanf("%lld", &n); sqrn = sqrt(n); getP(sqrn);
// 求出 g(\lfloor\frac{n}{x}\rfloor,0) 的值，并存一下 \lfloor\frac{n}{x}\rfloor 对应的下标
for (ll l = 1, r, tn; l <= n; l = r + 1)
{
r = n / (n / l); w[++tot] = n / l; tn = w[tot] % mod;
g1[tot] = tn * (tn + 1) % mod * inv2 % mod - 1;
g2[tot] = tn * (tn + 1) % mod * (2 * tn + 1) % mod * inv2 % mod * inv3 % mod - 1;
if (w[tot] <= sqrn) ind1[w[tot]] = tot;
else ind2[n / w[tot]] = tot;
}
// 根据递归式算出 g(n,j) 的值
for (int i = 1; i <= tp; ++i)
for (int j = 1; j <= tot && (ll)p[i] * p[i] <= w[j]; ++j)
{
int k = (w[j] / p[i]) <= sqrn ? ind1[w[j] / p[i]] : ind2[n / (w[j] / p[i])];
(g1[j] -= (ll)p[i] * (g1[k] - sp1[i - 1] + mod) % mod) %= mod;
(g2[j] -= (ll)p[i] * p[i] % mod * (g2[k] - sp2[i - 1] + mod) % mod) %= mod;
(g1[j] += mod) %= mod; (g2[j] += mod) %= mod;
}
// 递归求 S(n,0)，根据某个定理，不需要记忆化也能保证复杂度
printf("%d\n", (S(n, 0) + 1) % mod); return 0;
}
```