Min_25 筛 & Min_26 筛 & zzt 求和法 给（人能看懂的）代码的 学习笔记

看不懂别人博客里写的代码，所以只好自己实现常数超大的版本了，，，

记号：

$$\begin{array}{rl}\mathbf{P}& :\text{质数集}\\ p_k& :\text{第 k 个质数}\\ \text{lpf}(n)& :\text{n 的最小质因子}\\ x/y& :\lfloor {\displaystyle \frac{x}{y}}\rfloor \end{array}$$

前置知识：

$n$ 的块筛，指 $n/i$ 的 $O(\sqrt{n})$ 种取值。

用整除分块可以轻松求出 $n$ 的块筛，并给它们编号：

L = sqrt(n);
for (int l = 1, r; l <= n; l = r + 1)
{
    r = n / (n / l), v[++_] = n / l;
    if (v[_] <= L)
        o1[v[_]] = _;
    else
        o2[n / v[_]] = _;
}
//此时 v 中从大到小存储 n 的块筛
int O(int x) { return x <= L ? o1[x] : o2[n / x]; } //返回 x 是从大到小第几个块筛

给定数论函数 $f$。$\mathrm{\forall }p\in \mathbf{P}$，$f(p)$ 是关于 $p$ 的少项式，$f(p^k)$ 可以快速求值。给定 $n$，求 $\sum _{i=1}^{n}f(i)$。

Part I 求前缀质数处的和

即求 $\sum _{p\in \mathbf{P},p\le n}f(p)$。

Min_25 筛为啥叫 ex埃筛 呢？

上面提到，$f(p)$ 是关于 $p$ 的少项式（假设有 $o$ 项），即 $f(p)=\sum _{i=1}^o{a}_i{p}^{s_i}$，则 $\sum _{p\in \mathbf{P},p\le n}f(p)=\sum _{i=1}^o{a}_i\sum _{p\in \mathbf{P},p\le n}{p}^{s_i}$。

于是对于少项式的每一项，我们只需要求 $\sum _{p\in \mathbf{P},p\le n}{p}^{s_i}$，设 $g(n)=n^{s_i}$（$g$ 是完全积性函数！）。

模拟埃筛的过程，设 $G_{k,n}$ 表示 $\sum _{i=1}^n[\text{lpf}(i)>p_k\vee i\in \mathbf{P}]g(i)$，即第 $k$ 轮埃筛（筛去 $p_k$ 的倍数）后剩下的位置的 $g$ 值之和。

显然有初始值 $G_{0,n}=\sum _{i=2}^{n}g(i)$，我们钦定 $1$ 初始时就被筛掉了。

埃筛只用筛 $\sqrt{n}$ 轮，所以设 $p_c$ 为 $\sqrt{n}$ 以内最大的质数，则答案为 $G_{c,n}$。

对于 $k\le c$，预处理 $s_k=\sum _{i=1}^{k}g(p_i)$。

Min_25

考虑算 $G_{k,n}$。

筛完 $k-1$ 轮，剩下的位置的 $g$ 值之和为 $G_{k-1,n}$。

第 $k$ 轮需要筛掉 $p_k$ 的倍数，也就是说 $G_{k,n}=G_{k-1,n}-\sum _{i=1}^n[\text{lpf}(i)=p_k]g(i)$。

考虑后面的 $\sum _{i=1}^n[\text{lpf}(i)=p_k]g(i)$，$\text{lpf}(i)=p_k$ 当且仅当 $\text{lpf}\left({\displaystyle \frac{i}{p_k}}\right)\ge p_k$，

枚举 ${\displaystyle \frac{i}{p_k}}$，则 $\sum _{i=1}^n[\text{lpf}(i)=p_k]g(i)=\sum _{i=1}^{n/p_k}[\text{lpf}(i)\ge p_k]g(i\times p_k)=g(p_k)\sum _{i=1}^{n/p_k}[\text{lpf}(i)\ge p_k]g(i)$。

考虑后面的 $\sum _{i=1}^{n/p_k}[\text{lpf}(i)\ge p_k]g(i)$，它是 $G_{k-1,n/p_k}$ 吗？并不是！

可以发现，$G_{k-1,n/p_k}$ 比 $\sum _{i=1}^{n/p_k}[\text{lpf}(i)\ge p_k]g(i)$ 多统计了小于 $p_k$ 的素数处的值，即 $\sum _{i=1}^{n/p_k}[\text{lpf}(i)\ge p_k]g(i)=G_{k-1,n/p_k}-s_{k-1}$。

依次回代，可得 $G_{k,n}=G_{k-1,n}-g(p_k)(G_{k-1,n/p_k}-s_{k-1})$。

这里的 $n$ 只需取到块筛位置，而 $k$ 只需取到 $\pi (\sqrt{n})$，则总复杂度为 $O(\sum _{i=1}^{\sqrt{n}}{\displaystyle \frac{\sqrt{i}}{\log i}}+\sum _{i=1}^{\sqrt{n}}{\displaystyle \frac{\sqrt{n/i}}{\log n/i}})=O\left({\displaystyle \frac{n^{3/4}}{\log n}}\right)$。

实现时 $G$ 可以滚，整除太慢可以对 $i\le \sqrt{n}$ 预处理 $d_i={\displaystyle \frac{1}{i}}$（没错，$d_i$ 是 double），然后要算 $a/b$ 时算 $\lfloor a\times d_b+{10}^{-6}\rfloor$（加 eps 防止掉精度）即可。

给出洛谷 P5493 的代码：

#include <cmath>
#include <cstdio>
#define int long long
bool b[10000050];
double d[10000050];
int n, k, M, L, c, _, Z, _f[20], _o[20], _v[20], v[10000050], o1[10000050], o2[10000050], g[10000050], p[10000050], s[10000050];
int P(int x, int y)
{
    int q = 1;
    for (x %= M; y; y >>= 1, x = x * x % M)
        if (y & 1)
            q = q * x % M;
    return q;
}
int Q(int n) //拉插求 2^k+...+n^k
{
    if (n <= k + 2)
        return _f[n];
    else
    {
        static int _l[20], _r[20];
        int q = 0, _z = 1;
        _l[0] = _r[k + 3] = 1;
        for (int i = 1; i <= k + 2; ++i)
            _l[i] = _l[i - 1] * ((n - i) % M) % M;
        for (int i = k + 2; i >= 1; --i)
            _r[i] = _r[i + 1] * ((n - i) % M) % M;
        for (int i = 1; i <= k + 2; ++i)
            q = (q + _f[i] * _l[i - 1] % M * _r[i + 1] % M * _v[i - 1] % M * _v[k + 2 - i] % M * (k + 2 - i & 1 ? M - 1 : 1)) % M;
        return q;
    }
}
int O(int x) { return x <= L ? o1[x] : o2[n / x]; } //返回 x 是从大到小第几个块筛
signed main()
{
    scanf("%lld%lld%lld", &n, &k, &M);
    for (int i = _o[0] = 1; i <= 15; ++i)
        _o[i] = _o[i - 1] * i % M;
    _v[15] = P(_o[15], M - 2);
    for (int i = 14; i >= 0; --i)
        _v[i] = _v[i + 1] * (i + 1) % M;
    for (int i = 2; i <= k + 2; ++i)
        _f[i] = (_f[i - 1] + P(i, k)) % M;
    L = sqrt(n), d[1] = 1;
    for (int i = 2; i <= L; ++i)
    {
        if (!b[i])
            p[++c] = i, s[c] = (s[c - 1] + P(i, k)) % M; //记得预处理 s
        for (int j = 1; i * p[j] <= L; ++j)
        {
            b[i * p[j]] = 1;
            if (!(i % p[j]))
                break;
        }
        d[i] = 1.0 / i; //还有预处理 d
    }
    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = n / (n / l), v[++_] = n / l;
        if (v[_] <= L)
            o1[v[_]] = _;
        else
            o2[n / v[_]] = _;
        g[_] = Q(n / l); //g[0][n/l]=2^k+...+(n/l)^k, n/l 的编号为 _ 所以存在 g[_] 处
    }
    //此时 v 中从大到小存储 n 的块筛
    for (int j = 1; j <= c; ++j)
        for (int i = 1; p[j] * p[j] <= v[i]; ++i)
            g[i] = (g[i] + M - (s[j] + M - s[j - 1]) * (g[O(v[i] * d[p[j]] + 1e-6)] + M - s[j - 1]) % M) % M;
            // g[i] 即为博客中的 g[j][v[i]], v[i] 的编号为 i 所以存在 g[i] 处
            // 根据博客，g[j][v[i]]=g[j+1][v[i]]-p[j]^k*(g[j+1][v[i]/p[j]]-s[j-1]), 可以在上面的代码中一一对应
    for (int i = 1; i <= L; ++i)
        Z = (Z + g[O(n * d[i] + 1e-6)] * i % M * i) % M;
    printf("%lld", Z);
    return 0;
}

经测试，这份代码 10s 内能跑到 $n=6\times {10}^{12}$ 左右。

Min_26

$O\left({\displaystyle \frac{n^{3/4}}{\log n}}\right)$ 也太菜了，连杜都不如，能不能 给 力 点 啊？

考虑 根 号 分 治。

设置阈值 $T$。求 $G_{k,n}$ 时，对 $n\ge T$ 用上面的转移式转移，

对 $n<T$，根据 $G_{k,n}=G_{k-1,n}-\sum _{i=1}^n[\text{lpf}(i)=p_k]g(i)$，爆搜出所有 $i\le n,\text{lpf}(i)=p_k$ 的 $i$ 来转移……

诶等等啊，本来每次转移是 $O(1)$ 的，现在 $n<T$ 的每次转移还要爆搜一遍，那不是更菜了吗？

所以我们不能对每个 $n$ 爆搜所有 $i\le n,\text{lpf}(i)=p_k$ 的 $i$，

但我们可以爆搜所有 $i<T,\text{lpf}(i)=p_k$ 的 $i$，然后对每个 $i$ 更新所有 $i\le n$ 的 $n$，一次爆搜完成所有转移……

？那你复杂度不是没变吗

所以更新所有 $i\le n$ 的 $n$ 时不能枚举所有 $i\le n$ 的 $n$ 来更新，注意到这是一个后缀减，所以可以用 树 状 数 组 来优化。

分析一下复杂度，首先每个数只会被爆搜到一次，所以树状数组上的修改只有 $O(T)$ 次，复杂度 $O(T\log T)$，

对于 $n>T$，每次转移需要一个树状数组上的查询，转移有 $O\left(\sum _{i=1}^{n/T}{\displaystyle \frac{\sqrt{n/i}}{\log n/i}}\right)=O\left({\displaystyle \frac{n}{\sqrt{T}\log n}}\right)$ 次，复杂度 $O\left({\displaystyle \frac{n}{\sqrt{T}}}\right)$。

总复杂度 $O(({\displaystyle \frac{n}{\sqrt{T}}}+T\log T)$，取 $T=O\left({\left({\displaystyle \frac{n}{\log n}}\right)}^{\frac{2}{3}}\right)$ 得到最优复杂度 $O(n^{\frac{2}{3}}{\log }^{\frac{1}{3}}n)$……

……这是不是比原来还菜了？

继续根号分治，求 $G_{k,n}$ 时对 $p_k<T_2$ 只用转移式而不用树状数组，复杂度 $O(\pi (T_2)\sqrt{n})=O\left({\displaystyle \frac{T_2\sqrt{n}}{\log T_2}}\right)$。

此时我们只需对 $p_k\ge T_2$ 爆搜所有 $i<T,\text{lpf}(i)=p_k$ 的 $i$，取个 $T_2=n^{\frac{1}{6}}$ 爆搜到的 $i$ 就只有大概 $O\left({\displaystyle \frac{T}{\log n}}\right)$ 个了，

于是树状数组上的修改贡献的复杂度为 $O(T)$，而查询贡献的复杂度仍为 $O\left({\displaystyle \frac{n}{\sqrt{T}}}\right)$，

总复杂度 $O({\displaystyle \frac{n}{\sqrt{T}}}+T)$，取 $T=O(n^{\frac{2}{3}})$ 得到最优复杂度 $O(n^{\frac{2}{3}})$，和杜一样了！

给出洛谷 P5493 的代码：

#include <cmath>
#include <cstdio>
#define int long long
bool b[10000050];
double d[10000050];
int n, k, M, L, c, _, Z, n6, n3, n23, t6, t3, _f[20], _o[20], _v[20], v[10000050], o1[10000050], o2[10000050], g[10000050], p[10000050], s[10000050], C[10000050];
int P(int x, int y)
{
    int q = 1;
    for (x %= M; y; y >>= 1, x = x * x % M)
        if (y & 1)
            q = q * x % M;
    return q;
}
int Q(int n) //拉插求 2^k+...+n^k
{
    if (n <= k + 2)
        return _f[n];
    else
    {
        static int _l[20], _r[20];
        int q = 0, _z = 1;
        _l[0] = _r[k + 3] = 1;
        for (int i = 1; i <= k + 2; ++i)
            _l[i] = _l[i - 1] * ((n - i) % M) % M;
        for (int i = k + 2; i >= 1; --i)
            _r[i] = _r[i + 1] * ((n - i) % M) % M;
        for (int i = 1; i <= k + 2; ++i)
            q = (q + _f[i] * _l[i - 1] % M * _r[i + 1] % M * _v[i - 1] % M * _v[k + 2 - i] % M * (k + 2 - i & 1 ? M - 1 : 1)) % M;
        return q;
    }
}
int O(int x) { return x <= L ? o1[x] : o2[n / x]; } //返回 x 是从大到小第几个块筛
//块筛是从大到小存的，所以对 >=i 的块筛减去一个数是前缀减
//后缀差分一下，前缀减单点求值变为单点减后缀求和
void G(int x, int k)
{
    for (; x; x &= x - 1)
        C[x] = (C[x] + k) % M;
}
int W(int x)
{
    int q = 0;
    for (; x <= _; x += x & -x)
        q = (q + C[x]) % M;
    return q;
}
void D(int u, int fu, int l)
{
    G(O(n / (n / u)), (M - fu) % M); //找 >=u 的第一个块筛，设这个块筛是 n/o，则 o<=n/u，这个块筛为 n/(n/u)，
    //>=u 的所有块筛，即这个块筛及其之前的块筛，g 值都要减去 fu
    for (int i = l; i <= c && u * p[i] < n23; ++i)
        for (int j = p[i], fj = 1; u * j < n23; j *= p[i])
            fj = fj * (s[i] + M - s[i - 1]) % M, D(u * j, fu * fj % M, i + 1);
}
void U(int u) //爆搜 lpf=p[u] 的数，注意 p[u] 本身不应该被筛掉，所以要做一次单点加把之后的单点减抵消掉
{
    G(O(n / (n / p[u])), (s[u] + M - s[u - 1]) % M);
    for (int i = p[u], fi = 1; i < n23; i *= p[u])
        fi = fi * (s[u] + M - s[u - 1]) % M, D(i, fi, u + 1);
}
signed main()
{
    scanf("%lld%lld%lld", &n, &k, &M);
    L = sqrt(n), n6 = pow(n, 1.0 / 6), n23 = pow(n, 2.0 / 3), n3 = sqrt(n23);
    //注意这里需要保证用 n3 内的质数可以筛完 n23，即要保证 n3=sqrt(n23)，所以写 n3=pow(n,1.0/3) 会有问题
    for (int i = _o[0] = 1; i <= 15; ++i)
        _o[i] = _o[i - 1] * i % M;
    _v[15] = P(_o[15], M - 2);
    for (int i = 14; i >= 0; --i)
        _v[i] = _v[i + 1] * (i + 1) % M;
    for (int i = 2; i <= k + 2; ++i)
        _f[i] = (_f[i - 1] + P(i, k)) % M;
    d[1] = t6 = t3 = 1; //t6,t3 表示 n6,n3 内最大质数的编号+1（我下面的循环是左闭右开的，所以要+1）
    for (int i = 2; i <= L; ++i)
    {
        if (!b[i])
            p[++c] = i, s[c] = (s[c - 1] + P(i, k)) % M; //记得预处理 s
        for (int j = 1; i * p[j] <= L; ++j)
        {
            b[i * p[j]] = 1;
            if (!(i % p[j]))
                break;
        }
        if (i == n6)
            t6 += c;
        if (i == n3)
            t3 += c;
        d[i] = 1.0 / i; //还有预处理 d
    }
    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = n / (n / l), v[++_] = n / l;
        if (v[_] <= L)
            o1[v[_]] = _;
        else
            o2[n / v[_]] = _;
        g[_] = Q(n / l); //g[0][n/l]=2^k+...+(n/l)^k, n/l 的编号为 _ 所以存在 g[_] 处
    }
    //此时 v 中从大到小存储 n 的块筛
    for (int j = 1; j < t6; ++j) //对于 p[j]<T2 只用转移式而不用树状数组，这里是左闭右开所以 t6 要加一
        for (int i = 1; p[j] * p[j] <= v[i]; ++i)
            g[i] = (g[i] + M - (s[j] + M - s[j - 1]) * (g[O(v[i] * d[p[j]] + 1e-6)] + M - s[j - 1]) % M) % M;
    for (int i = _; v[i] < n23; --i)
        C[i] = (g[i] + M - g[i + (i & -i)]) % M; //O(n) 建树状数组
    for (int j = t6; j < t3; ++j)                //这里是左闭右开所以 t3 要加一
    {
        for (int i = 1; n23 <= v[i]; ++i) //v[i]>=T 用转移式求 g[i]
        {
            int o = v[i] * d[p[j]] + 1e-6;
            if (o < n23) //v[i]/p[j]<T 时，g[v[i]/p[j]] 在树状数组上
                g[i] = (g[i] + M - (s[j] + M - s[j - 1]) * (W(O(o)) + M - s[j - 1]) % M) % M;
            else //v[i]/p[j]>=T 时，g[v[i]/p[j]] 是用转移式得到的，直接在 g 上查
                g[i] = (g[i] + M - (s[j] + M - s[j - 1]) * (g[O(o)] + M - s[j - 1]) % M) % M;
        }
        U(j); //爆搜 lpf=p[j] 的数，更新 v[i]<T 的 g[i]
    }
    for (int i = _; v[i] < n23; --i)
        g[i] = W(i);              //用 sqrt(T) 以内的质数筛完之后，所有 v[i]<T 的 g[i] 不会再变，从树状数组上取出 g 值
    for (int j = t3; j <= c; ++j) //对于 p[j]>=sqrt(T) 只有 v[i]>=T 的 g[i] 会受影响，直接用转移式求 g[i]
        for (int i = 1; p[j] * p[j] <= v[i]; ++i)
            g[i] = (g[i] + M - (s[j] + M - s[j - 1]) * (g[O(v[i] * d[p[j]] + 1e-6)] + M - s[j - 1]) % M) % M;
    for (int i = 1; i <= L; ++i)
        Z = (Z + g[O(n * d[i] + 1e-6)] * i % M * i) % M;
    printf("%lld", Z);
    return 0;
}

经测试，这份代码在 10s 内能跑到 $n=1.5\times {10}^{13}$ 左右。

zzt 求和法

能 不 能 再 给 力 点 啊 ？

考虑减少一些树状数组上的查询。

注意到对于 $n<T$，$p_k^2>n$ 时 $G_{k,n}=G_{k-1,n}$，也就是说树状数组上 $n$ 处的值不会变了，

此时不妨把 $n$ 处的值存下来，以后就不用在树状数组上查询了。

这样的话，对于 $n\ge T$，算 $G_{k,n}$ 时只有在 $p_k^2\le n/p_k\iff p_k^3\le n$ 时才需要在树状数组上查询 $G_{k-1,n/p_k}$，

则树状数组上的查询共有 $O\left(\sum _{i=1}^{n/T}{\displaystyle \frac{\sqrt[3]{n/i}}{\log n/i}}\right)=O\left({\displaystyle \frac{n}{T^{\frac{2}{3}}\log n}}\right)$ 次， 复杂度 $O\left({\displaystyle \frac{n}{T^{\frac{2}{3}}}}\right)$，

而 $n\ge T$ 的转移就有 $O\left({\displaystyle \frac{n}{\sqrt{T}\log n}}\right)>O\left({\displaystyle \frac{n}{T^{\frac{2}{3}}}}\right)$ 次，所以瓶颈在枚举所有状态进行转移，而不在树状数组。

树状数组上的修改贡献的复杂度仍为 $O(T)$，所以总复杂度 $O({\displaystyle \frac{n}{\sqrt{T}\log n}}+T)$，

取 $T=O\left({\left({\displaystyle \frac{n}{\log n}}\right)}^{\frac{2}{3}}\right)$ 得到最优复杂度 $O\left({\left({\displaystyle \frac{n}{\log n}}\right)}^{\frac{2}{3}}\right)$。

给出洛谷 P5493 的代码：

#include <cmath>
#include <cstdio>
#define int long long
bool b[10000050];
double d[10000050];
int n, k, M, L, c, _, Z, n6, n3, n23, t6, t3, _f[20], _o[20], _v[20], v[10000050], o1[10000050], o2[10000050], g[10000050], p[10000050], s[10000050], C[10000050];
int P(int x, int y)
{
    int q = 1;
    for (x %= M; y; y >>= 1, x = x * x % M)
        if (y & 1)
            q = q * x % M;
    return q;
}
int Q(int n) //拉插求 2^k+...+n^k
{
    if (n <= k + 2)
        return _f[n];
    else
    {
        static int _l[20], _r[20];
        int q = 0, _z = 1;
        _l[0] = _r[k + 3] = 1;
        for (int i = 1; i <= k + 2; ++i)
            _l[i] = _l[i - 1] * ((n - i) % M) % M;
        for (int i = k + 2; i >= 1; --i)
            _r[i] = _r[i + 1] * ((n - i) % M) % M;
        for (int i = 1; i <= k + 2; ++i)
            q = (q + _f[i] * _l[i - 1] % M * _r[i + 1] % M * _v[i - 1] % M * _v[k + 2 - i] % M * (k + 2 - i & 1 ? M - 1 : 1)) % M;
        return q;
    }
}
int O(int x) { return x <= L ? o1[x] : o2[n / x]; } //返回 x 是从大到小第几个块筛
//块筛是从大到小存的，所以对 >=i 的块筛减去一个数是前缀减
//后缀差分一下，前缀减单点求值变为单点减后缀求和
void G(int x, int k)
{
    for (; x; x &= x - 1)
        C[x] = (C[x] + k) % M;
}
int W(int x)
{
    int q = 0;
    for (; x <= _; x += x & -x)
        q = (q + C[x]) % M;
    return q;
}
void D(int u, int fu, int l)
{
    G(O(n / (n / u)), (M - fu) % M); //找 >=u 的第一个块筛，设这个块筛是 n/o，则 o<=n/u，这个块筛为 n/(n/u)，
    //>=u 的所有块筛，即这个块筛及其之前的块筛，g 值都要减去 fu
    for (int i = l; i <= c && u * p[i] < n23; ++i)
        for (int j = p[i], fj = 1; u * j < n23; j *= p[i])
            fj = fj * (s[i] + M - s[i - 1]) % M, D(u * j, fu * fj % M, i + 1);
}
void U(int u) //爆搜 lpf=p[u] 的数，注意 p[u] 本身不应该被筛掉，所以要做一次单点加把之后的单点减抵消掉
{
    G(O(n / (n / p[u])), (s[u] + M - s[u - 1]) % M);
    for (int i = p[u], fi = 1; i < n23; i *= p[u])
        fi = fi * (s[u] + M - s[u - 1]) % M, D(i, fi, u + 1);
}
signed main()
{
    scanf("%lld%lld%lld", &n, &k, &M);
    L = sqrt(n), n6 = pow(n, 1.0 / 6), n23 = pow(n / log(n), 2.0 / 3), n3 = sqrt(n23);
    //注意这里需要保证用 n3 内的质数可以筛完 n23，即要保证 n3=sqrt(n23)，所以写 n3=pow(n,1.0/3) 会有问题
    for (int i = _o[0] = 1; i <= 15; ++i)
        _o[i] = _o[i - 1] * i % M;
    _v[15] = P(_o[15], M - 2);
    for (int i = 14; i >= 0; --i)
        _v[i] = _v[i + 1] * (i + 1) % M;
    for (int i = 2; i <= k + 2; ++i)
        _f[i] = (_f[i - 1] + P(i, k)) % M;
    d[1] = t6 = t3 = 1; //t6,t3 表示 n6,n3 内最大质数的编号+1（我下面的循环是左闭右开的，所以要+1）
    for (int i = 2; i <= L; ++i)
    {
        if (!b[i])
            p[++c] = i, s[c] = (s[c - 1] + P(i, k)) % M; //记得预处理 s
        for (int j = 1; i * p[j] <= L; ++j)
        {
            b[i * p[j]] = 1;
            if (!(i % p[j]))
                break;
        }
        if (i == n6)
            t6 += c;
        if (i == n3)
            t3 += c;
        d[i] = 1.0 / i; //还有预处理 d
    }
    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = n / (n / l), v[++_] = n / l;
        if (v[_] <= L)
            o1[v[_]] = _;
        else
            o2[n / v[_]] = _;
        g[_] = Q(n / l); //g[0][n/l]=2^k+...+(n/l)^k, n/l 的编号为 _ 所以存在 g[_] 处
    }
    //此时 v 中从大到小存储 n 的块筛
    for (int j = 1; j < t6; ++j) //对于 p[j]<T2 只用转移式而不用树状数组，这里是左闭右开所以 t6 要加一
        for (int i = 1; p[j] * p[j] <= v[i]; ++i)
            g[i] = (g[i] + M - (s[j] + M - s[j - 1]) * (g[O(v[i] * d[p[j]] + 1e-6)] + M - s[j - 1]) % M) % M;
    for (int i = _; v[i] < n23; --i)
        C[i] = (g[i] + M - g[i + (i & -i)]) % M; //O(n) 建树状数组
    int k = _;                                   //表示当前 i>k 的 g[i] 不会再变，已经存下来
    for (int j = t6; j < t3; ++j)                //这里是左闭右开所以 t3 要加一
    {
        for (; v[k] < p[j] * p[j]; --k)
            g[k] = W(k);                  //v[k]<p[j]^2 的 g[k] 不会再变，把它存下来
        for (int i = 1; n23 <= v[i]; ++i) //v[i]>=T 用转移式求 g[i]
        {
            int o = v[i] * d[p[j]] + 1e-6;
            if (o < n23 && O(o) <= k) //v[i]/p[j]<T 时，g[v[i]/p[j]] 在树状数组上
                                      //并且 O(o)>k 时 g[O(o)] 的值已经存下来，不用在树状数组上查询
                g[i] = (g[i] + M - (s[j] + M - s[j - 1]) * (W(O(o)) + M - s[j - 1]) % M) % M;
            else //v[i]/p[j]>=T 时，g[v[i]/p[j]] 是用转移式得到的，直接在 g 上查
                g[i] = (g[i] + M - (s[j] + M - s[j - 1]) * (g[O(o)] + M - s[j - 1]) % M) % M;
        }
        U(j); //爆搜 lpf=p[j] 的数，更新 v[i]<T 的 g[i]
    }
    for (; v[k] < n23; --k)
        g[k] = W(k);              //取出剩下的 v[k]<T 的 g[k] 值
    for (int j = t3; j <= c; ++j) //对于 p[j]>=sqrt(T) 只有 v[i]>=T 的 g[i] 会受影响，直接用转移式求 g[i]
        for (int i = 1; p[j] * p[j] <= v[i]; ++i)
            g[i] = (g[i] + M - (s[j] + M - s[j - 1]) * (g[O(v[i] * d[p[j]] + 1e-6)] + M - s[j - 1]) % M) % M;
    for (int i = 1; i <= L; ++i)
        Z = (Z + g[O(n * d[i] + 1e-6)] * i % M * i) % M;
    printf("%lld", Z);
    return 0;
}

经测试，这份代码在 10s 内能跑到 $n=2\times {10}^{13}$ 左右。

Part II 积性函数求和

……Min_25 筛好像是用来做积性函数求和的来着？

即求 $\sum _{i=1}^{n}f(i)$。

仿照上面的 $G$，我们设 $F_{k,n}=\sum _{i=1}^n[\text{lpf(i)}\ge p_k]f(i)$，则答案为 $F_{1,n}+1$（$f(1)$ 没有被统计到）。

Min_25

考虑 $F_{k,n}$ 统计到的数，枚举统计到的数的最小质因子及其次数，可得 $F_{k,n}=\sum _{p_k\le p_i\le n}\sum _{c\ge 1,p_i^c\le n}f(p_i^c)F_{i+1,n/p_i^c}$，

最小值因子 $p_i$ 要枚举到 $n$，复杂度炸了。考虑把其中的质数拿出来单独算，则剩下的合数的最小质因子只需枚举到 $\sqrt{n}$，则

$$\begin{array}{rl}{F}_{k,n}& =\sum _{p_k\le p_i\le \sqrt{n}}\sum _{c\ge 1,p_i^c\le n}f(p_i^c)([c>1]+F_{i+1,n/p_i^c})+\sum _{p_k\le p_i\le n}f(p_i)\\ & =\sum _{p_k\le p_i\le \sqrt{n}}\sum _{c\ge 1,p_i^{c+1}\le n}(f(p_i^c)F_{i+1,n/p_i^c}+f(p_i^{c+1}))+\sum _{p\in \mathbf{P},p\le n}f(p)-s_{k-1}\end{array}$$

（最后一个等号成立的原因：考虑 $p_i^c\le n<p_i^{c+1}$ 的 $c$，因为 $p_i^{c+1}>n$，所以 $n/p_i^c<p_i<p_{i+1}$，所以 $F_{i+1,n/p_i^c}=0$）

$s$ 就是之前预处理的 $s_k=\sum _{i=1}^{k}g(p_i)$，可以发现 $p_k\le \sqrt{n}$ 所以 $s$ 够用。

至于 $\sum _{p\in \mathbf{P},p\le n}f(p)$……这不是 Part I 吗？

直接按照转移式开搜，边界条件 $p_k>n$ 时 $F_{k,n}=0$，记忆化都不用，复杂度可以证明是 $O\left({\displaystyle \frac{n^{\frac{3}{4}}}{\log n}}\right)$。

给出洛谷 P5325 的代码：

#include <cmath>
#include <cstdio>
#define M 1000000007
#define _2 500000004
#define _6 166666668
#define int long long
bool b[10000050];
double d[10000050];
int n, L, c, _, v[10000050], o1[10000050], o2[10000050], p[10000050], g[10000050][2], s[10000050][2];
int O(int x) { return x <= L ? o1[x] : o2[n / x]; }
int F(int k, int n)
{
    if (p[k] > n)
        return 0;
    int q = (g[O(n)][1] + (M << 1) - g[O(n)][0] - s[k - 1][1] + s[k - 1][0]) % M;
    for (int i = k; i <= c && p[i] * p[i] <= n; ++i)        //枚举最小值因子 p[i]
        for (int j = p[i], o = n; j * p[i] <= n; j *= p[i]) //枚举 j=p[i]^c，保证 p[i]^(c+1)<=n
        {
            int u = j % M, v = j * p[i] % M;
            q = (q + u * (u + M - 1) % M * F(i + 1, o = o * d[p[i]] + 1e-6) + v * (v + M - 1)) % M;
            //答案加上 f(p[i]^c)*F(i+1,n/p[i]^c)+f(p[i]^(c+1))，可以在上面的代码中一一对应
        }
    return q;
}
signed main()
{
    scanf("%lld", &n);
    L = sqrt(n), d[1] = 1;
    for (int i = 2; i <= L; ++i)
    {
        if (!b[i])
            p[++c] = i, s[c][0] = (s[c - 1][0] + i) % M, s[c][1] = (s[c - 1][1] + i * i) % M;
        for (int j = 1; i * p[j] <= L; ++j)
        {
            b[i * p[j]] = 1;
            if (!(i % p[j]))
                break;
        }
        d[i] = 1.0 / i;
    }
    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = n / (n / l), v[++_] = n / l;
        if (v[_] <= L)
            o1[v[_]] = _;
        else
            o2[n / v[_]] = _;
        int u = n / l % M;
        g[_][0] = (u * (u + 1) % M * _2 + M - 1) % M;
        g[_][1] = (u * (u + 1) % M * (u << 1 | 1) % M * _6 + M - 1) % M;
    }
    for (int j = 1; j <= c; ++j)
        for (int i = 1; p[j] * p[j] <= v[i]; ++i)
        {
            int u = O(v[i] * d[p[j]] + 1e-6);
            g[i][0] = (g[i][0] + M - p[j] * (g[u][0] + M - s[j - 1][0]) % M) % M;
            g[i][1] = (g[i][1] + M - p[j] * p[j] % M * (g[u][1] + M - s[j - 1][1]) % M) % M;
        }
    //以上是 Part I，不再解释
    printf("%lld", (F(1, n) + 1) % M); //别忘了 +1
    return 0;
}

可见 Min_25 筛的码量并不大，具有一定的实战意义？

经测试，这份代码在 10s 内能跑到 $n=2.5\times {10}^{12}$ 左右。

Min_26

你应该可以猜到接下来会发生什么了（

考虑 根 号 分 治。

设 $t_3=\pi (\sqrt[3]{n})+1,t6=\pi (\sqrt[6]{n})+1$（与 Part I 给的代码中一样）。

考虑对 $n$ 的所有块筛 $i$，求出 $F_{t_3,i}$。

$F_{t_3,i}$ 统计到的数的 $\text{lpf}\ge p_{t_3}>\sqrt[3]{n}$，则其质因子（可重）最多只有两个，

把质数单独拿出来算，枚举有两个质因子的数的 $\text{lpf}$，则

$F_{t_3,n}=\sum _{p_{t_3}\le p_i\le \sqrt{n}}(f(p_i^2)+f(p_i)(\sum _{j\in \mathbf{P},j\le n/p_i}f(j)-s_i))+\sum _{p\in \mathbf{P},p\le n}f(p)-s_{t_3-1}$。

接着，对 $n$ 的所有块筛 $i$ 求出 $F_{t_6,i}$。不妨设 $t_6\le k<t_3$。

可以发现， $F_{k,n}$ 比 $F_{k+1,n}$ 只多统计了 $\text{lpf}=p_k$ 的数，

枚举这些数中 $p_k$ 的次数，则 $F_{k,n}=F_{k+1,n}+\sum _{i=1}^n[\text{lpf(i)}=p_k]f(i)=F_{k+1,n}+\sum _{p_k^c\le n}f(p_k^c)(F_{k+1,n/p_k^c}+1)$。

直接转移复杂度显然不对，继续根号分治。

仿照 Part I 的做法，对 $n\ge T$ 用转移式转移，

对 $n<T$，爆搜所有 $i<T,\text{lpf(i)}=p_k$ 的 $i$，对每个 $i$ 更新所有 $i\le n$ 的 $n$，

然后这个还是后缀加，用树状数组维护。

最后，对 $n$ 的所有块筛 $i$ 求出 $F_{1,i}$。直接用刚刚的转移式转移即可。

分析一下复杂度，首先算 $F_{t_3,i}$ 时只有 $i\ge p_{t_3}^2>n^{\frac{2}{3}}$ 时会枚举第一个 $\sum$，复杂度 $O\left(\sum _{i=1}^{\sqrt[3]{n}}{\displaystyle \frac{\sqrt{n/i}}{\log n/i}}\right)=O\left({\displaystyle \frac{n^{\frac{2}{3}}}{\log n}}\right)$。

接着，由 $F_{t_3,i}$ 递推到 $F_{t_6,i}$ 的过程与 Part I 中由 $G_{t_6,i}$ 递推到 $G_{t_3,i}$ 的过程分析方法一样，复杂度 $O(n^{\frac{2}{3}})$。

最后由 $F_{t_6,i}$ 递推到 $F_{1,i}$ 时，需要枚举 $\sqrt[6]{n}$ 以内的素数在 $n$ 以内的幂，以及所有块筛，

jijidawang 证明了前者有 $O\left({\displaystyle \frac{\sqrt[6]{n}}{\log n}}\right)$ 个，而块筛有 $O(\sqrt{n})$ 个，复杂度就是二者之积 $O\left({\displaystyle \frac{n^{\frac{2}{3}}}{\log n}}\right)$。

给出洛谷 P5325 的代码：

#include <cmath>
#include <cstdio>
#define M 1000000007
#define _2 500000004
#define _6 166666668
#define int long long
bool b[10000050];
double d[10000050];
int n, L, c, _, n6, n3, n23, t6, t3, v[10000050], o1[10000050], o2[10000050], p[10000050], f[10000050], g[10000050][2], s[10000050][2], C[10000050];
int O(int x) { return x <= L ? o1[x] : o2[n / x]; }
void G(int x, int k)
{
    for (; x; x &= x - 1)
        C[x] = (C[x] + k) % M;
}
int W(int x)
{
    int q = 0;
    for (; x <= _; x += x & -x)
        q = (q + C[x]) % M;
    return q;
}
void D(int u, int fu, int l, int T) //与 Part I 不同的是，T=-1 时更新 f 而不是 g
{
    G(O(n / (n / u)), ~T ? (M - fu) % M : fu);
    for (int i = l; i <= c && u * p[i] < n23; ++i)
        for (int j = p[i], fj = 1; u * j < n23; j *= p[i])
            fj = fj = ~T ? fj * (s[i][T] + M - s[i - 1][T]) % M : j * (j - 1) % M, D(u * j, fu * fj % M, i + 1, T);
}
void U(int u, int T) //同上
{
    if (~T)
        G(O(n / (n / p[u])), (s[u][T] + M - s[u - 1][T]) % M);
    for (int i = p[u], fi = 1; i < n23; i *= p[u])
        fi = ~T ? fi * (s[u][T] + M - s[u - 1][T]) % M : i * (i - 1) % M, D(i, fi, u + 1, T);
}
signed main()
{
    scanf("%lld", &n);
    L = sqrt(n), n6 = pow(n, 1.0 / 6), n23 = pow(n, 2.0 / 3), n3 = sqrt(n23);
    d[1] = t6 = t3 = 1;
    for (int i = 2; i <= L + 1; ++i)
    {
        if (!b[i])
            p[++c] = i, s[c][0] = (s[c - 1][0] + i) % M, s[c][1] = (s[c - 1][1] + i * i) % M;
        for (int j = 1; i * p[j] <= L + 1; ++j)
        {
            b[i * p[j]] = 1;
            if (!(i % p[j]))
                break;
        }
        if (i == n6)
            t6 += c;
        if (i == n3)
            t3 += c;
        d[i] = 1.0 / i;
    }
    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = n / (n / l), v[++_] = n / l;
        if (v[_] <= L)
            o1[v[_]] = _;
        else
            o2[n / v[_]] = _;
        int u = n / l % M;
        g[_][0] = (u * (u + 1) % M * _2 + M - 1) % M;
        g[_][1] = (u * (u + 1) % M * (u << 1 | 1) % M * _6 + M - 1) % M;
    }
    for (int T = 0; T < 2; ++T)
    {
        for (int j = 1; j < t6; ++j)
            for (int i = 1; p[j] * p[j] <= v[i]; ++i)
                g[i][T] = (g[i][T] + M - (s[j][T] + M - s[j - 1][T]) * (g[O(v[i] * d[p[j]] + 1e-6)][T] + M - s[j - 1][T]) % M) % M;
        for (int i = _; v[i] < n23; --i)
            C[i] = (g[i][T] + M - g[i + (i & -i)][T]) % M;
        for (int j = t6; j < t3; ++j)
        {
            for (int i = 1; n23 <= v[i]; ++i)
            {
                int o = v[i] * d[p[j]] + 1e-6;
                if (o < n23)
                    g[i][T] = (g[i][T] + M - (s[j][T] + M - s[j - 1][T]) * (W(O(o)) + M - s[j - 1][T]) % M) % M;
                else
                    g[i][T] = (g[i][T] + M - (s[j][T] + M - s[j - 1][T]) * (g[O(o)][T] + M - s[j - 1][T]) % M) % M;
            }
            U(j, T);
        }
        for (int i = _; v[i] < n23; --i)
            g[i][T] = W(i);
        for (int j = t3; j <= c; ++j)
            for (int i = 1; p[j] * p[j] <= v[i]; ++i)
                g[i][T] = (g[i][T] + M - (s[j][T] + M - s[j - 1][T]) * (g[O(v[i] * d[p[j]] + 1e-6)][T] + M - s[j - 1][T]) % M) % M;
    }
    //以上是 Part I，不再解释
    for (int i = 1; v[i] >= p[t3]; ++i)
    {
        f[i] = (g[i][1] + (M << 1) - g[i][0] - s[t3 - 1][1] + s[t3 - 1][0]) % M; //初始化 f[t3][n]=g[n]-s[t3-1]
        for (int j = t3; j <= c && p[j] * p[j] <= v[i]; ++j)
        {
            int o = O(v[i] * d[p[j]] + 1e-6);
            f[i] = (f[i] + p[j] * p[j] % M * (p[j] * p[j] % M - 1) + p[j] * (p[j] - 1) % M * (g[o][1] + (M << 1) - g[o][0] - s[j][1] + s[j][0])) % M;
            //累加 f(p[j]^2)+f(p[j])*(g[n/p[j]]-s[j])，可以在上面的代码中一一对应
        }
    }
    for (int i = _; v[i] < n23; --i)
        C[i] = (f[i] + M - f[i + (i & -i)]) % M; //O(n) 建树状数组
    for (int j = t3 - 1; j >= t6; --j)
    {
        for (int i = 1; n23 <= v[i]; ++i)                      //v[i]>=T 用转移式求 f[i]
            for (int k = p[j], o = v[i]; k <= v[i]; k *= p[j]) //枚举 p[j]^c<=v[i]
            {
                o = o * d[p[j]] + 1e-6;
                if (o < n23) //v[i]/p[j]^c<T 时，f[v[i]/p[j]^c] 在树状数组上
                    f[i] = (f[i] + k % M * (k % M - 1) % M * (W(O(o)) + 1)) % M;
                else //v[i]/p[j]^c>=T 时，f[v[i]/p[j]^c] 是用转移式得到的，直接在 f 上查
                    f[i] = (f[i] + k % M * (k % M - 1) % M * (f[O(o)] + 1)) % M;
            }
        U(j, -1); //爆搜 lpf=p[j] 的数，更新 v[i]<T 的 f[i]
    }
    for (int i = _; v[i] < n23; --i)
        f[i] = W(i);                  //从树状数组上取出 v[i]<T 的 f[i]
    for (int j = t6 - 1; j >= 1; --j) //直接用转移式转移
        for (int i = 1; v[i] >= p[j]; ++i)
            for (int k = p[j], o = v[i]; k <= v[i]; k *= p[j])
                f[i] = (f[i] + k % M * (k % M - 1) % M * (f[O(o = o * d[p[j]] + 1e-6)] + 1)) % M;
    printf("%lld", (f[1] + 1) % M); //别忘了 +1
    return 0;
}

经测试，这份代码在 10s 内能跑到 $n=3.5\times {10}^{12}$ 左右。

值得一提的是，这份代码在模板题跑不过 Min_25，不知道为什么。

zzt 求和法

大 的 药 来 了

如果你会 PN 筛，你可以跳过下面的缩进部分。

PN 筛：有积性函数 $g$ 满足 $\mathrm{\forall }p\in \mathbf{P},g(p)=f(p)$，已知 $g$ 在 $n$ 的所有块筛处的前缀和，求 $\sum _{i=1}^{n}F(i)$。

做法：

称质因子次数最低二次的数为 Powerful Number（PN）。显然所有 PN 都能表示成 $a^2{b}^3$ 的形式。

考虑 $n$ 以内 PN 的个数。有 $O(\sum _{i=1}^{\sqrt{n}}\sqrt[3]{n/i^2})=O(\sqrt{n})$ 个。

设积性函数 $h$ 满足 $g\ast h=f$，则 $\mathrm{\forall }p\in \mathbf{P},f(p)=g(1)h(p)+g(p)h(1)=g(p)+h(p)$，

又因为 $g(p)=f(p)$，所以 $h(p)=0$，则 $h$ 只在 PN 处有值。则

$$\begin{array}{rl}\sum _{i=1}^{n}f(i)& =\sum _{d|i}h(d)g\left({\displaystyle \frac{i}{d}}\right)\\ & =\sum _{d=1}^{n}h(d)\sum _{i=1}^{n/d}g(i)\\ & =\sum _{d\in \mathbf{PN},d\le n}h(d)\sum _{i=1}^{n/d}g(i)\end{array}$$

爆搜 $d\in \mathbf{PN},d\le n$，统计答案即可。

注意跑 PN 筛的前提是已知 $g$ 在 $n$ 的所有块筛处的前缀和，

所以在 Min_25 筛板子题构造 $g(p)=p\phi (p)$ 并使用 PN 筛的复杂度不是 $O(\sqrt{n})$，

因为求 $g$ 在 $n$ 的所有块筛处的前缀和需要跑一遍 $O(n^{\frac{2}{3}})$ 的杜教筛。