Min_25 筛

怎么都开Min_25筛了那我也开

一些定义：

• \(n\)：前缀和上界
• \(p_i\)：第 \(i\) 个质数（定义 \(p_0 = 1\)）
• \(p_{\max}\)：\(< \sqrt n\) 的最大质数
• \(\text{mindiv}(n)\)：\(n\) 的最小质因子
• \(\mathbb{P}\)：质数集合

如无特殊注明，以下使用的 \(p\) 都属于 \(\mathbb{P}\)。

是更快且泛用性更广的筛子（相较于杜教筛而言）。它的复杂度在 \(\le 10^{13}\) 范围内被证明是常数很小的 \(O\left(\dfrac{n^{3/4}}{\log n}\right)\)，在实际应用中跑得比 \(O(n ^ {\frac 23})\) 都快。具体复杂度的证明可以看zzt的2018集训队论文。这里是我的阅读随笔
他（指Min_25）解决的是这么一类积性函数 \(f\) 的前缀和问题：

1. 积性
2. 在 \(p^k\) 处取值为一个关于 \(p\) 的多项式
3. \(f(p^k)\) 可以快速算得

目标：\(\sum\limits_{i = 1}^n f(i)\)

首先进行式子的拆解：

\[\sum_{i = 1}^n f(i) = \sum_{\text{p为质数}}^nf(p) + \sum_{\text{k为合数}}^n f(k) + f(1) \]

分别处理两个部分。你问 f(1)去哪了？最后加1不就行了

质数部分

在质数处的快速计算需要设辅助函数 \(g(x)\)，要求 \(g\) 完全积性，且在质数处的取值和 \(f\) 相同。容易发现由于原函数在质数处取值为多项式，我们总可以设一个或多个 \(g\) 表示原函数。当 \(g\) 不止一个时分别计算贡献后相加即可。
于是我们可以单独讨论一个完全积性函数 \(f(x)\) 对质数部分答案的贡献。

我们设一个二元函数 \(g(n, k)\)，定义如下：

\[g(n,k) = \sum_{i \in \mathbb{P}\ \lor \ \text{mindiv}(i) > p_k } ^n f(i) \]

可以发现，质数部分的答案为 \(g(n,p_{\max})\)，即小于等于 \(n\) 的所有质数。我们发现这样 \(\lor\) 后面的东西没用，但是带着更好求一些。具体见下。
考虑通过容斥一下求出 \(g(n,k-1)\) 转移到 \(g(n,k)\) 的转移方程。

我们发现，如果 \(k\) 增大，相应会删除一些数，这些数就是满足最小质因子为 \(p_k\) 的合数。考虑如何表示这些被删掉的数字，从共性下手。
我们提出一个 \(p_k\) 作为他们共有的最小质因子，即，应该被删掉的数字在除掉 \(p_k\) 后不会有 \(\le p_{k-1}\) 的最小质因子，即满足 \(\text{mindiv}(i) > p_{k-1}\) 且小于等于 \(\lfloor \frac n {p_k} \rfloor\) 的合数。
我们发现，这部分数字对答案的贡献可以被形式化地表示为下式：

\[f(p_k) \times (g(\lfloor \frac n {p_k} \rfloor , k - 1) - g(p_{k-1}, k-1)) \]

其中减号前面为质数与合数的共同贡献，减号后面为质数单独的贡献。
因此我们有状态转移方程：

\[g(n, k) = g(n, k-1) - f(p_k) \times (g(\lfloor \frac n {p_k} \rfloor , k - 1) - g(p_{k-1}, k-1)) \]

此处表明了 \(f\) 为完全积性函数的意义。质因子 \(p_k\) 可以被直接提出，当且仅当 \(f\) 完全积性。

考虑怎么转移。由数论分块可知我们所需的所有 \(\lfloor \frac n {p_k} \rfloor\) 状态是 \(O(\sqrt n)\) 级别的，且 \(p_{k} \le p_{\max}\)，因此总状态是 \(O(\sqrt n)\) 的。因此我们可以直接通过数论分块处理 \(i \le \sqrt{n}\) 的所有 \(g(i, 0)\) 与 \(g(\lfloor \frac n i \rfloor, 0)\) 取值。随后朴素地转移状态，复杂度正确。这里不需要写二维数组，只需要从大到小更新数值就行。
对于值的存储，我们可以单独开两个 \(\text{id}\) 数组，分别存储 \(x \le \sqrt n\) 时 \(x\) 对应值数组下标与 \(x > \sqrt n\) 时 \(\lfloor \frac n x \rfloor\) 对应值数组下标。容易发现两个 \(\text{id}\) 数组长度都是 \(O(\sqrt n)\)，这样实现了快速离散化。

luogu板子的代码

N = sqrt(n) ;
inline int get(register ll pos) { return pos < N ? id1[pos] : id2[n / pos]; } // 值的存储部分
inline ll S1(register ll x) { return x %= mod, x * (x + 1) / 2 % mod; }
inline ll S2(register ll x) { return x %= mod, x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod; }
// 两个完全积性函数

int prime[200001], cnt;
int vis[200001];
void init() { // 筛素数
    rep(i,2,N) {
        if (!vis[i]) prime[++cnt] = i;
        rep(j,1,cnt) { int tmp = i * prime[j]; if (tmp > N) break; vis[tmp] = 1; if(i % prime[j] == 0) break; }
    }
}

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;
	// 首先取得 g(k,0) 的值
}
rep(i,1,cnt) {
    rep(j,1,m) {
        if (1ll * prime[i] * prime[i] > v[j]) break;
        g1[j] = (g1[j] - 1ll * prime[i] * (g1[get(v[j] / prime[i])] - g1[get(prime[i - 1])]) % mod + mod) % mod;
        g2[j] = (g2[j] - 1ll * sq(prime[i]) * (g2[get(v[j] / prime[i])] - g2[get(prime[i - 1])]) % mod + mod) % mod;
		// 依照上式进行dp
    }
} 

统计答案

记 \(g(n, p_{\max})\) 为 \(g(n)\)。

我们再设一个函数 \(s(n, k)\)。这里的 \(f\) 函数指代原函数，不是自己构造的完全积性函数。

\[s(n, k) = \sum_{\text{mindiv}(i) > p_k}^n f(i) \]

容易发现 \(s(n,0)\) 即为最终答案。讨论 \(s(n, k)\) 如何求得。
仍然考虑将质数与合数的贡献分离。质数的贡献是所有大于 \(p_k\) 的质数，可以直接由 \(g\) 求得，即 \(g(n) - g(p_k)\)。合数的部分需要枚举质数的倍数求得，这也是 Min_25 筛被称为“拓展埃氏筛”的原因。我们枚举所有 \(k > j\) 的 \(p_k\) 以及其 \(e\) 次幂，每次贡献的数字为最小质因子为 \(p_k\) 且其次数为 \(e\) 的数，即除去 \(p_k^e\) 后满足 \(\text{mindiv}(j) > p_k\) 的所有 \(j \le n\)，因此最终统计答案时应考虑到 \(e\) 的边界为 \(p_k^{e+1} \le n\)。注意，我们在枚举次幂时总会计入 \(p_k\) 的贡献，而一个数字显然只能被计入一次，因此需要在计算时在 \(e>1\) 时减去1的贡献。形式化地，我们有状态转移方程：

\[s(n, k) = g(n) - g(p_k) + \sum_{p_k > p_j}^{p_{\max}} \sum_{e \ge 1 \land p_k ^{e+1} \le n} f(p_k^e) \times (s(\lfloor \frac n {p_k^e} \rfloor,k) + [e \neq 1]) \]

关于怎么转移……暴力
不进行记忆化的暴力转移的复杂度是对的。在 \(10^{11}\) 的数据下递归到被记忆化的部分的次数是 0。

第二部分luogu板子

inline ll S(ll x, int y) {
    if (prime[y] >= x) return 0;
    ll ret = (g2[get(x)] - g1[get(x)] - g2[get(prime[y])] + g1[get(prime[y])] + (mod<<1)) % mod;
    rep(i,  y+1, cnt) {
        if (1ll * prime[i] * prime[i] > x) break;
        for (register ll w = prime[i]; w * prime[i] <= x; w = w * prime[i]) {
            ret = (ret + 1ll * F1(w) * S(x / w, i) % mod + F1(w * prime[i])) % mod;
        }
    } return ret;
} 

printf("%lld\n",(S(n, 0) + 1) % mod);

未来闲话主题预定：Min_25 筛应用