前置数学符号 在展开讨论之前，先约定本文所使用的符号体系：

• 记 lpf(x) 为整数 x 的最小质因子
• 记 p_k 为第 k 个质数
• 记 π(n) 为不超过 n 的质数总数

算法核心思想

Min_25筛法是一种用于计算数论函数前缀和的巧妙算法。其基本思路是将求和过程拆解为质数部分与合数部分分别处理。设我们要求解的是：

∑i=1n [i ∈ ℙ] f(i)

这里的 ℙ 表示质数集合。由于 n 的规模通常很大，直接枚举所有质数显然不切实际。Min_25筛采用了一种"拟合"策略：构造一个完全积性函数 F，使得 F(p) = f(p) 对所有质数 p 成立。

有了完全积性函数 F 后，原问题转化为：

∑i=1n [i ∈ ℙ] F(i)

容斥原理的应用

直接计算质数处的函数值很困难，我们转而采用容斥思想：用 ∑<sub>i=1</sub><sup>n</sup> F(i) 减去所有合数处的 F 值。注意到任意合数 x 都满足 lpf(x) ≤ ⌊√x⌋，这意味着我们可以采用类似埃拉托斯特尼筛法的思路，遍历所有不超过 ⌊√n⌋ 的质数，逐步扣除它们的贡献。

动态规划 formulation

定义状态：

g(n, x) = ∑i=1n [i ∈ ℙ 或 lpf(i) > px] F(i)

即统计在 n 范围内，所有质数以及最小质因子大于 p_x 的数对应的 F 值之和。

考虑从 g(n, x-1) 转移到 g(n, x)。当限制条件从 lpf > p<sub>x-1</sub> 变为 lpf > p_x 时，被剔除的恰好是那些最小质因子恰好等于 p_x 的数。因此：

Δ = g(n, x-1) - g(n, x) = ∑i=1n [lpf(i) = p_x] F(i)

接下来分析 lpf(i) = p_x 的结构。由于 F 是完全积性函数，我们可以提取一个因子 F(p_x)，剩下的部分需要满足 lpf ≥ p_x。考察 g(⌊n/p_x⌋, x-1)，发现它比目标多计算了那些质因子小于 p_x 的质数处的 F 值。预先处理 sp_x = ∑<sub>i=1</sub><sup>x</sup> F(p_i)，即可得到递推公式：

g(n, x) = g(n, x-1) - F(p_x) · (g(⌊n/p_x⌋, x-1) - spx-1)

根据杜教筛的理论，我们仅需要 g(⌊n/i⌋, x) 中满足 p_x² ≤ ⌊n/i⌋ 的那些状态。可以证明，对这些关键点进行动态规划，时间复杂度为 O(n^{3/4} / log n)。

离散化处理

实现时需要对 ⌊n/i⌋ 的不同取值进行编号。采用两个大小为 O(√n) 的数组 idx_small 和 idx_large：当 ⌊n/i⌋ ≤ √n 时，将编号存入 idx_small[⌊n/i⌋]；否则令 j = ⌊n / ⌊n/i⌋⌋，将编号存入 idx_large[j] 中。

最终我们只需要 g(n, π(⌊√n⌋))，即所有质数处的 F 值之和。由于第二维只与其前一个状态有关，可以采用滚动数组优化空间。

处理原始函数 f

定义另一状态：

S(n, x) = ∑i=1n [lpf(i) > p_x] f(i)

将求和分为质数和合数两部分考虑。

质数部分的贡献

对于质数部分，贡献为：

g(n, π(⌊√n⌋)) - sp_x

合数部分的贡献

由于 f 不一定是完全积性函数，无法单独提取质因子。这里采用枚举质数和指数的策略：先枚举 p_i (i > x)，再枚举其指数 c，这样可以将 f(p_i^c) 完整提取出来。剩余部分应满足 lpf > p_i，即 S(⌊n/p_i^c⌋, i)。需要注意的是，当 c > 1 时，需要补回 f(p_i^c)；当 c = 1 时，f(p_i) 已在质数部分计算过，无需重复。

综合得到递推式：

S(n, x) = g(n, π(⌊√n⌋)) - sp_x + ∑i=x+1p_i²≤n ∑c≥1, p_i^c≤n f(p_i^c) · (S(⌊n/p_i^c⌋, i) + [c > 1])

该部分的时间复杂度被证明是 O(n^{1-ε})。

最终答案为 S(n, 0) + f(1)。

算法要求

Min_25筛法对目标函数 f 有以下要求：

• 能够快速计算 F(p) = f(p)，且 F 的前缀和可快速求得
• 能够快速计算 f(p^c)（c ≥ 1）

这些条件在实际应用中通常容易满足。若 f(p) 可以分解为若干完全积性函数的线性组合，则对应地组合 g 和 sp 的值即可。

应用实例

示例一：Luogu P5325 本题要求 f(p) = p² - p。构造 F₁(p) = p 和 F₂(p) = p² 即可求解。

#include <bits/stdc++.h>
using namespace std;
using int64 = long long;
const int MAXN = 2000005;
const int MOD = 1e9 + 7;

int64 n, sqrt_n;
int64 values[MAXN << 1];
int prime_cnt;
int primes[MAXN / 10];
int64 sum_p1[MAXN / 10];
int64 sum_p2[MAXN / 10];
int total_vals;
int idx_small[MAXN];
int idx_large[MAXN];
int64 dp1[MAXN << 1];
int64 dp2[MAXN << 1];
bool is_composite[MAXN];

int64 addmod(int64 a, int64 b) {
    a += b;
    if (a >= MOD) a -= MOD;
    return a;
}
int64 submod(int64 a, int64 b) {
    a -= b;
    if (a < 0) a += MOD;
    return a;
}
int64 mulmod(int64 a, int64 b) {
    return (a * b) % MOD;
}
int64 mod_pow(int64 a, int64 e) {
    int64 r = 1;
    while (e) {
        if (e & 1) r = mulmod(r, a);
        a = mulmod(a, a);
        e >>= 1;
    }
    return r;
}

void init_sieve(int limit) {
    for (int i = 2; i <= limit; ++i) {
        if (!is_composite[i]) {
            primes[++prime_cnt] = i;
            sum_p1[prime_cnt] = addmod(sum_p1[prime_cnt - 1], i);
            sum_p2[prime_cnt] = addmod(sum_p2[prime_cnt - 1], mulmod(i, i));
        }
        for (int j = 1; j <= prime_cnt && i * primes[j] <= limit; ++j) {
            is_composite[i * primes[j]] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

int64 compute_prefix_sum(int64 x) {
    return mulmod(mulmod(x % MOD, (x + 1) % MOD), 500000004);
}

int64 compute_square_sum(int64 x) {
    x %= MOD;
    return mulmod(mulmod(x, (x + 1) % MOD), mulmod((2 * x + 1) % MOD, 166666668));
}

int64 solve_recursive(int64 m, int pos) {
    if (primes[pos] >= m) return 0;
    int index = (m <= sqrt_n) ? idx_small[m] : idx_large[n / m];
    int64 result = submod(submod(dp2[index], dp1[index]), 
                          submod(sum_p2[pos], sum_p1[pos]));
    
    for (int i = pos + 1; i <= prime_cnt && 1LL * primes[i] * primes[i] <= m; ++i) {
        int64 p = primes[i];
        for (int64 power = p, exp = 1; power <= m; power *= p, ++exp) {
            int64 p_mod = power % MOD;
            int64 term = mulmod(mulmod(p_mod, submod(p_mod, 1)), 
                               addmod(solve_recursive(m / power, i), exp > 1));
            result = addmod(result, term);
        }
    }
    return result;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    
    cin >> n;
    sqrt_n = sqrt(n);
    init_sieve(sqrt_n);
    
    for (int64 left = 1, right, val; left <= n; left = right + 1) {
        val = n / left;
        right = n / val;
        values[++total_vals] = val;
        
        int64 v_mod = val % MOD;
        dp1[total_vals] = submod(compute_prefix_sum(v_mod), 1);
        dp2[total_vals] = submod(compute_square_sum(v_mod), 1);
        
        if (val <= sqrt_n) idx_small[val] = total_vals;
        else idx_large[n / val] = total_vals;
    }
    
    for (int i = 1; i <= prime_cnt; ++i) {
        int p = primes[i];
        int64 p_mod = p % MOD;
        int64 p2_mod = mulmod(p_mod, p_mod);
        
        for (int j = 1; j <= total_vals && 1LL * p * p <= values[j]; ++j) {
            int64 div_result = values[j] / p;
            int target_idx = (div_result <= sqrt_n) ? idx_small[div_result] 
                                                    : idx_large[n / div_result];
            
            dp1[j] = submod(dp1[j], mulmod(p_mod, submod(dp1[target_idx], sum_p1[i - 1])));
            dp2[j] = submod(dp2[j], mulmod(p2_mod, submod(dp2[target_idx], sum_p2[i - 1])));
        }
    }
    
    cout << addmod(solve_recursive(n, 0), 1) << '\n';
    return 0;
}

示例二：SPOJ DIVCNTK 本题要求计算 f(n) = d(n^k)，其中 d(x) 表示约数个数。可以证明 f 是积性函数，且满足 f(p) = k + 1、f(p^c) = c·k + 1。使用 Min_25 筛，取 F = 1，最后乘以 k + 1 即可。

示例三：Luogu P5493 质数前缀统计 仅需使用 Min_25 筛的前半部分进行动态规划。边界条件 g(n, 0) 是自然数幂和，可通过 O(k) 的拉格朗日插值求得。总时间复杂度 O(n^{3/4} / log n + √n · k)。

示例四：LOJ 6235 区间素数个数 令 F(p) = 1，直接使用 Min_25 筛的前半部分进行 DP 即可。时间复杂度 O(n^{3/4} / log n + √n)。

示例五：SPOJ DIVFACT4 要计算 d(n!)，即 n! 的约数个数。经典结论表明，n! 中质数 p 的指数为 ∑_{c≥1} ⌊n/p^c⌋。因此：

d(n!) = ∏i=1π(n) (1 + ∑c≥1 ⌊n/p_i^c⌋)

对于大于 √n 的质数，它们仅贡献 1 + ⌊n/p⌋，可以通过整除分块处理：对每个块 [l, r]，先用 Min_25 筛求出该区间内的质数个数 cnt，然后给答案乘上 (1 + ⌊n/l⌋)^{cnt}。时间复杂度 O(n^{3/4} / log n + √n · log n)。

对于不超过 √n 的质数，直接暴力计算即可，共有 O(√n / log n) 个。

由于模数可能很大，需要使用快速乘法避免溢出。

#include <bits/stdc++.h>
using namespace std;
using int64 = long long;
using i128 = __int128_t;

const int MAXN = 5000005;
int64 mod;

int prime_cnt;
int primes[MAXN];
int prime_idx[MAXN];
bool is_composite[MAXN];
int sqrt_n;
int64 n;
int t;

int total_vals;
int64 seq[MAXN];
int idx_small[MAXN];
int idx_large[MAXN];
int64 g[MAXN];

int64 mulmod(int64 a, int64 b) {
    return (i128)a * b % mod;
}

int64 powmod(int64 a, int64 e) {
    int64 r = 1 % mod;
    a %= mod;
    while (e) {
        if (e & 1) r = mulmod(r, a);
        a = mulmod(a, a);
        e >>= 1;
    }
    return r;
}

void build_sieve(int limit) {
    is_composite[1] = true;
    for (int i = 2; i <= limit; ++i) {
        if (!is_composite[i]) {
            prime_idx[i] = ++prime_cnt;
            primes[prime_cnt] = i;
        }
        for (int j = 1; j <= prime_cnt && i * primes[j] <= limit; ++j) {
            is_composite[i * primes[j]] = true;
            prime_idx[i * primes[j]] = j;
            if (i % primes[j] == 0) break;
        }
    }
}

inline int get_index(int64 x) {
    return (x <= sqrt_n) ? idx_small[x] : idx_large[n / x];
}

int64 count_factor(int64 base, int64 p) {
    int64 res = 0;
    while (base >= p) {
        base /= p;
        res += base % mod;
    }
    return res % mod;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    
    build_sieve(MAXN - 5);
    cin >> t;
    
    while (t--) {
        cin >> n >> mod;
        sqrt_n = sqrt(n);
        total_vals = 0;
        int64 answer = 1 % mod;
        
        for (int64 left = 1, right, val; left <= n; left = right + 1) {
            val = n / left;
            right = n / val;
            seq[++total_vals] = val;
            g[total_vals] = (val - 1) % mod;
            
            if (val <= sqrt_n) idx_small[val] = total_vals;
            else idx_large[right] = total_vals;
        }
        
        for (int i = 1; i <= prime_cnt && primes[i] <= sqrt_n; ++i) {
            for (int j = 1; j <= total_vals && 1LL * primes[i] * primes[i] <= seq[j]; ++j) {
                int pos = get_index(seq[j] / primes[i]);
                g[j] -= (g[pos] - (i - 1));
                if (g[j] < 0) g[j] += mod;
            }
        }
        
        for (int i = 1; i <= prime_cnt && primes[i] <= sqrt_n; ++i) {
            answer = mulmod(answer, addmod(count_factor(n, primes[i]), 1));
        }
        
        for (int64 left = sqrt_n + 1, right, val; left <= n; left = right + 1) {
            right = n / (val = n / left);
            int64 cnt = g[get_index(right)] - g[get_index(left - 1)];
            if (cnt < 0) cnt += mod;
            answer = mulmod(answer, powmod((val + 1) % mod, cnt));
        }
        
        cout << answer << '\n';
    }
    return 0;
}

示例六：SPOJ APS2 这是 Min_25 筛的一种非平凡应用。需要分别处理质数和合数。

对于质数部分，要求 ∑i=1n [i ∈ ℙ] i，直接令 F(p) = p 并使用 Min_25 筛的前半部分即可。

对于合数部分，考察 Min_25 筛前半部分的 DP 公式：

g(n, x) = g(n, x-1) - F(p_x) · (g(⌊n/p_x⌋, x-1) - spx-1)

注意到 g(⌊n/p_x⌋, x-1) - sp<sub>x-1</sub> 正好统计了 n 以内最小质因子等于 p_x 的数的个数。在 DP 过程中累加这一结果，即可得到合数部分的答案。

总时间复杂度 O(n^{3/4} / log n + √n)。

算法小结 Min_25 筛法是一种在数论函数求和中极为强大的工具，尤其适用于处理与质数相关的复杂求和问题。其核心在于通过构造完全积性函数进行拟合，结合容斥原理与动态规划，在 O(n^{3/4} / log n) 的时间复杂度内完成计算。该算法在竞赛编程和数论研究中都有广泛的应用价值。