Min_25

初步理解。

我们要求的是积性函数 \(f(p)\) 的前缀和。\(f(p)\) 要是一个能快速求值的东西（一般为多项式），\(f(p^c)\) 能快速求值。

\(p_k\)：第 k 小的质数。

\(minp(n)\)：\(n\) 的最小质因数。

\(F_p(n) = \sum_{2 \leq p \leq n} f(p)\)，p 为质数。\(p_0 = 1\)。

\(F(k, n) = \sum_{i = 2}^{n} [p_k \leq minp(i)] f(i)\)。

我们先把问题转化为求 \(F(1, n) + f(1) = F(1, n) + 1\)。

如何求 \(F(k, n)\)？枚举每个数的最小质因数：

\[\begin{aligned} F(k, n) &= \sum_{2 \leq i \leq n} [p_k \leq minp(i)] f(i) \\ &= \sum_{\substack{k \leq i \\ p_i^2 \leq n}} \sum_{\substack{c \geq 1 \\ p_i^c \leq n}} f(p_i^c)([c > 1] + F(i+1, n / p_i^c)) + \sum_{\substack{k \leq i \\ p_i \leq n}} f(p_i) \\ &= \sum_{\substack{k \leq i \\ p_i^2 \leq n}} \sum_{\substack{c \geq 1 \\ p_i^c \leq n}} f(p_i^c)([c > 1] + F(i+1, n / p_i^c)) + F_p(n) - F_p(p_{k - 1}) \end{aligned} \]

由于满足 \(p_i^c \leq n < p_i^{c + 1}\)，有 \(n / p_i^c < p_i < p_{i + 1}\)，所以 \(F(i + 1, n / p_i^c) = 0\)。可以进一步化简成 \(\sum_{\substack{k \leq i \\ p_i^2 \leq n}} \sum_{\substack{c \geq 1 \\ p_i^{c + 1} \leq n}} f(p_i^c)F(i+1, n / p_i^c) + f(p_i^{c + 1})\)。边界值为 \(F(k, n) = 0\)，其中 \(p_k > n\)。

问题转化为求 \(F_p(n)\)。

发现递归求 \(F(k, n)\) 的过程中最多只会用到 1, 2, ... , \(\sqrt{n}\)，\(n / \sqrt{n}\), ... , \(n\) 这 \(2 \sqrt{n}\) 个地方的值。考虑只求出这些地方的 \(F_p\)。由于 \(f(p)\) 一般是多项式，我们把多项式每一项拆开单独考虑其对整体的贡献。

假设当前考虑 \(g(p) = p^s\)，要求 \(m\) 属于 \(2 \sqrt{n}\) 个值中，所有 \(\sum_{2 \leq p \leq m} g(p)\)。

设 \(G(k, n)\) 为埃式筛法 \(k\) 轮后剩下的数的 \(g(i)\) 的值之和。即 \(G(k, n) = \sum_{2 \leq i \leq n} [i \in \mathbb{P} \ \vee \ p_k < minp(i) ] g(i)\)。

• 对于 \(p_k > n\) 的部分，没有贡献。

• 对于 \(p_k \leq n\) 的部分， \(p_k\) 会被筛掉，所以减少 \(g(p_k) \times G(k - 1, n / p_k)\) （ \(g(p)\) 是完全积性函数）。

• 但是这样会多减掉小于 \(p_k\) 的质数的部分，所以要加回去 \(g(p_k) \times G(k - 1, p_{k - 1})\)。

所以有 \(G(k, n) = G(k - 1, n) - [p_k^2 \leq n]g(p_k)(G(k - 1, n / p_k) - G(k - 1, p_{k - 1}))\)。

总结一下，我们需要预处理出 \(\sqrt{n}\) 级别的所有素数以及素数的若干次方的前缀和。然后递推计算 \(G(k, n)\)，合并成 \(F_p(n)\)。最后递归计算 \(F(1, n)\)。注意到有效值只有 \(2 \sqrt{n}\) 个，开两个数组记录这些位置就好了。

LG5325

i64 n; int lim;
int primes[S], tot = 0, minp[S];
i64 sp1[S], sp2[S];
i64 val[S];
int cnt, id1[S], id2[S];
i64 G[S][2], Fp[S];

void sieve(int n) {
    tot = 0;
    for (int i = 2; i <= n; ++i) {
        if (minp[i] == 0) {
            minp[i] = i;
            primes[++tot] = i;
        }
        for (int j = 1; j <= tot && 1ll * i * primes[j] <= n; ++j) {
            int p = primes[j];
            minp[i * p] = p;
            if (p == minp[i]) {
                continue;
            }
        }
    }
    for (int i = 1; i <= tot; ++i) {
        sp1[i] = (sp1[i - 1] + primes[i]) % mod;
        sp2[i] = (sp2[i - 1] + 1ll * primes[i] * primes[i] % mod) % mod;
    }
}

int idx(i64 x) {
    return x <= lim ? id1[x] : id2[n / x];
}

void init(i64 n) {
    cnt = 0;
    for (i64 l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        i64 v = n / l;
        val[++cnt] = v;
        if (v <= lim) {
            id1[v] = cnt;
        } else {
            id2[n / v] = cnt;
        }
        G[cnt][0] = 1ll * ((v + 2) % mod) * inv2 % mod * ((v - 1) % mod) % mod;
        G[cnt][1] = 1ll * v % mod * ((v + 1) % mod) % mod * ((v * 2 + 1) % mod) % mod * inv6 % mod;
        G[cnt][1] = (G[cnt][1] - 1ll + mod) % mod;
    }
    for (int i = 1; i <= tot; ++i) {
        int p = primes[i];
        for (int j = 1; j <= cnt && 1ll * p * p <= val[j]; ++j) {
            i64 v = val[j];
            G[j][0] = (G[j][0] - 1ll * p * ((G[idx(v / p)][0] - sp1[i - 1] + mod) % mod) % mod + mod) % mod;
            G[j][1] = (G[j][1] - 1ll * p * p % mod * ((G[idx(v / p)][1] - sp2[i - 1] + mod) % mod) % mod + mod) % mod;
        }
    }
    for (int i = 1; i <= cnt; ++i) {
        Fp[i] = (G[i][1] - G[i][0] + mod) % mod;
    }
}

i64 calc(i64 k, i64 n) {
    if (n < primes[k] || n <= 1) {
        return 0;
    }
    i64 res = (Fp[idx(n)] - (sp2[k - 1] - sp1[k - 1] + mod) % mod + mod) % mod;
    for (int i = k; i <= tot && 1ll * primes[i] * primes[i] <= n; ++i) {
        int p = primes[i];
        i64 pw = p;
        for (i64 c = 1; 1ll * p * pw <= n; ++c, pw = 1ll * p * pw) {
            i64 v = pw % mod, v2 = 1ll * v * v % mod;
            v2 = (v2 - v + mod) % mod;
            res = (res + 1ll * v2 * calc(i + 1, n / pw) % mod) % mod;
            v = 1ll * pw * p % mod, v2 = 1ll * v * v % mod;
            v2 = (v2 - v + mod) % mod;
            res = (res + v2) % mod;
        }
    }
    return res % mod;
}

void solve() {
    std::cin >> n;
    lim = sqrt(n);
    init(n);
    i64 ans = (calc(1, n) + 1ll) % mod;
    std::cout << ans % mod << "\n";
}