Min_25 筛学习笔记（伪）

2019-01-08

本文迁移自老博客，原始链接为 https://lookas2001.com/min_25-%e7%ad%9b%e5%ad%a6%e4%b9%a0%e7%ac%94%e8%ae%b0/

下面口胡，一些犄角旮旯地方的理解。

将所有数做质因数分解，考虑将数按照第一个质因子及其次数（即 p^e）分类，我们发现由于积性函数的性质，这一类的和可以共同提取出一个 f(p^e) ，即式子会变成 f(p^e)(f(a_1 / p^e) + f(a_2 / p^e) + ...) 。可以发现式子的后半部分就是一个很类似的问题。这大概就是 Min_25 筛的主要思想吧。

具体来讲是定义了两个求和函数，其中一个求和函数（即朱震霆论文中的 h ，网上大多数博客的 g）辅助另外一个求和函数（即朱震霆论文中的 g ，网上大多数博客的 S）求值。

那个辅助求值函数的求值过程非常类似埃氏筛的过程，故这个方法也称拓展埃氏筛法。

空间使用只需要 O(\sqrt n) ，因为每次递归形如（都是向下取整的除法） n -> n / i , n / i -> n / i / j ，而 n / i / j 可以等价为 n / (i * j) 。数论分块，下略。

复杂度啥的见朱震霆的论文吧。

下面是一些（应该没有锅的）数学公式，代码直接套这些公式就好了，模板在 Counting Divisors 。

积性函数： $f(x)$

要求有 $f(p^e) = \sum_k a_k(e)p^k$ 。

辅助求和函数： $h_k(n, j) = \sum_{2 \leq i \leq n 且 (i 的最小质因子 > p_j 或 i 是质数)} i^k$

若 $p_j > \sqrt n$ ，有 $h_k(n, j) = h_k(n, j - 1)$ ;

否则，有 $h_k(n, j) = h_k(n, j - 1) - p_j^k(h_k(\lfloor \frac n{p_j}\rfloor, j - 1) - h_k(p_j - 1, j - 1))$ 。

特别的，对于 $j = 0$ ，有 $h_k(n, j) = \sum_{2 \leq i \leq n} i^k$ 。

简记 $h_k(n) = h_k(n, j)$ 其中 $p_j > \sqrt n$ （即我们要的辅助求和函数最后没有合数项的）。

求和函数： $g(n, m) = \sum_{2 \leq i \leq n 且 (i 的最小质因子 > p_m)} f(i)$

可知 $g(n, m) = \sum_k a_k(1)(h(n) - h(p_m)) + \sum_{p_m < p_j\leq \sqrt n 且 e \geq 1 且 p_j^{e + 1} \leq n} (f(p_j^e) g(\lfloor \frac n{p_j^e}\rfloor, j) + f(p_j^{e + 1}))$

$p_0$ 的话当 $0$ 处理就好咯。

$g(n, 0)$ 即为所求。

#include <cstdio>
#include <cstring>

typedef unsigned long long u64;

const int MAX_SQRT_N = 1e5;

struct euler_sieve {
    static const int MAX_N = MAX_SQRT_N;
    bool f[MAX_N + 10];
    int p[MAX_N + 10];
    int p_n;

    inline void operator()(int n) {
        memset(f, 0, sizeof f);
        p_n = 0;
        for (int i = 2; i <= n; ++i) {
            if (!f[i]) {
                p[++p_n] = i;
            }
            for (int j = 1; p[j] * i <= n; ++j) {
                f[p[j] * i] = true;
                if (i % p[j] == 0) break;
            }
        }
    }
} es;

struct Min_25_sieve {
    u64 n;

    u64 b[2 * MAX_SQRT_N + 10]; // 所有可能的 n / i 值，值是从大到小的（废话）
    int b_n, sqrt_n;
    inline int locate(u64 i) { return i < sqrt_n ? b_n - i + 1 : n / i; } // 查找 b[j] == i 的位置 j

    int d0[2 * MAX_SQRT_N + 10]; // 存储 h 函数的值，数组命名个人癖好而已。 h 函数的取值只有 O(\sqrt n) 种。这里的第 i 个位置存储的不是 h(i) 而是 h(b[i])
    inline void pre_calc() { // 计算 h 函数
        for (int i = 1; i <= b_n; ++i) d0[i] = 1 * (b[i] - 1); // h(b[i], 0) ，j = 0 理解为没有“i 的最小质因子 > p_j”这条限制即可。这里以及下面的 1 * 代表的就是 i^k （本题中 k = 0（不是输入的那个 k 啊））
        for (int j = 1; (u64)es.p[j] * es.p[j] <= n; ++j) { // 枚举质数
            for (int i = 1; (u64)es.p[j] * es.p[j] <= b[i]; ++i) { // 枚举了有用的，即 b[i]
                d0[i] = d0[i] - 1 * (d0[locate(b[i] / es.p[j])] - d0[locate(es.p[j] - 1)]); // 一定有 p[j] - 1 < \sqrt n ，所以一定可以找到一个合法的位置（不会落在一个段内） // 由于 b[i] 的值是从大道小的，所以求值顺序是没有问题的，不用开另外一个缓存数组
            }
        }
    }

    u64 f_a0_k;
    inline u64 f_a0(int e) { return e * f_a0_k + 1; }
    inline u64 f(int p, int e) { return f_a0(e); }

    u64 calc(u64 n, int m) { // 计算 g 函数
        u64 rst = f_a0(1) * (d0[locate(n)] - d0[locate(es.p[m])]);
        for (int j = m + 1; (u64)es.p[j] * es.p[j] <= n; ++j) {
            u64 pe = es.p[j];
            for (int e = 1; (u64)pe * es.p[j] <= n; ++e, pe *= es.p[j]) {
                rst += f(es.p[j], e) * calc(n / pe, j) + f(es.p[j], e + 1);
            }
        }
        return rst;
    }

    inline u64 operator()(u64 n_, u64 f_a0_k_) {
        n = n_;
        f_a0_k = f_a0_k_;

        // 分块
        sqrt_n = 0; while ((u64)sqrt_n * sqrt_n < n) ++sqrt_n;
        b_n = 0; for (u64 i = 1; i <= n; i = n / (n / i) + 1) b[++b_n] = n / i;

        // 处理辅助求和函数
        pre_calc();

        // 求和
        es.p[0] = 0, d0[b_n + 1] = 0;
        return calc(n, 0) + 1; // + 1 是那个甩单的 1
    }
} mns;

int main() {
    // freopen("in.in", "r", stdin);

    int tc; scanf("%d", &tc);
    // int tc = 1;
    es(MAX_SQRT_N + 5); // 筛出质数
    while (tc--) {
        u64 n, f_a0_k; scanf("%llu%llu", &n, &f_a0_k);
        printf("%llu\n", mns(n, f_a0_k));
    }

    return 0;
}