Min_25 筛 & Min_26 筛 & zzt 求和法 给（人能看懂的）代码的 学习笔记
看不懂别人博客里写的代码，所以只好自己实现常数超大的版本了，，，
记号：
前置知识：
\(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\)。\(\forall p\in\mathbf P\)，\(f(p)\) 是关于 \(p\) 的少项式，\(f(p^k)\) 可以快速求值。给定 \(n\)，求 \(\sum\limits_{i=1}^nf(i)\)。
Part I 求前缀质数处的和
即求 \(\sum\limits_{p\in\mathbf P,p\le n}f(p)\)。
Min_25 筛为啥叫 ex埃筛 呢？
上面提到，\(f(p)\) 是关于 \(p\) 的少项式（假设有 \(o\) 项），即 \(f(p)=\sum\limits_{i=1}^oa_ip^{s_i}\)，则 \(\sum\limits_{p\in\mathbf P,p\le n}f(p)=\sum\limits_{i=1}^oa_i\sum\limits_{p\in\mathbf P,p\le n}p^{s_i}\)。
于是对于少项式的每一项，我们只需要求 \(\sum\limits_{p\in\mathbf P,p\le n}p^{s_i}\)，设 \(g(n)=n^{s_i}\)（\(g\) 是完全积性函数！）。
模拟埃筛的过程，设 \(G_{k,n}\) 表示 \(\sum\limits_{i=1}^n[\text{lpf}(i)>p_k\vee i\in\mathbf P]g(i)\)，即第 \(k\) 轮埃筛（筛去 \(p_k\) 的倍数）后剩下的位置的 \(g\) 值之和。
显然有初始值 \(G_{0,n}=\sum\limits_{i=2}^ng(i)\)，我们钦定 \(1\) 初始时就被筛掉了。
埃筛只用筛 \(\sqrt n\) 轮，所以设 \(p_c\) 为 \(\sqrt n\) 以内最大的质数，则答案为 \(G_{c,n}\)。
对于 \(k\le c\)，预处理 \(s_k=\sum\limits_{i=1}^kg(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\limits_{i=1}^n[\text{lpf}(i)=p_k]g(i)\)。
考虑后面的 \(\sum\limits_{i=1}^n[\text{lpf}(i)=p_k]g(i)\)，\(\text{lpf}(i)=p_k\) 当且仅当 \(\text{lpf}\left(\dfrac i{p_k}\right)\ge p_k\)，
枚举 \(\dfrac i{p_k}\)，则 \(\sum\limits_{i=1}^n[\text{lpf}(i)=p_k]g(i)=\sum\limits_{i=1}^{n/p_k}[\text{lpf}(i)\ge p_k]g(i\times p_k)=g(p_k)\sum\limits_{i=1}^{n/p_k}[\text{lpf}(i)\ge p_k]g(i)\)。
考虑后面的 \(\sum\limits_{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\limits_{i=1}^{n/p_k}[\text{lpf}(i)\ge p_k]g(i)\) 多统计了小于 \(p_k\) 的素数处的值，即 \(\sum\limits_{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\left(\sum\limits_{i=1}^{\sqrt n}\dfrac{\sqrt i}{\log i}+\sum\limits_{i=1}^{\sqrt n}\dfrac{\sqrt{n/i}}{\log{n/i}}\right)=O\left(\dfrac{n^{3/4}}{\log n}\right)\)。
实现时 \(G\) 可以滚，整除太慢可以对 \(i\le\sqrt n\) 预处理 \(d_i=\dfrac 1i\)（没错，\(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(\dfrac{n^{3/4}}{\log n}\right)\) 也太菜了，连杜都不如，能不能 给 力 点 啊？
考虑 根 号 分 治。
设置阈值 \(T\)。求 \(G_{k,n}\) 时，对 \(n\ge T\) 用上面的转移式转移，
对 \(n<T\)，根据 \(G_{k,n}=G_{k-1,n}-\sum\limits_{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\limits_{i=1}^{n/T}\dfrac{\sqrt{n/i}}{\log{n/i}}\right)=O\left(\dfrac n{\sqrt T\log n}\right)\) 次，复杂度 \(O\left(\dfrac n{\sqrt T}\right)\)。
总复杂度 \(O(\left(\dfrac n{\sqrt T}+T\log T\right)\)，取 \(T=O\left(\left(\dfrac n{\log n}\right)^{\frac 23}\right)\) 得到最优复杂度 \(O(n^{\frac 23}\log^{\frac 13}n)\)……
……这是不是比原来还菜了？
继续根号分治，求 \(G_{k,n}\) 时对 \(p_k<T_2\) 只用转移式而不用树状数组，复杂度 \(O(\pi(T_2)\sqrt n)=O\left(\dfrac{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 16}\) 爆搜到的 \(i\) 就只有大概 \(O\left(\dfrac T{\log n}\right)\) 个了，
于是树状数组上的修改贡献的复杂度为 \(O(T)\)，而查询贡献的复杂度仍为 \(O\left(\dfrac n{\sqrt T}\right)\)，
总复杂度 \(O\left(\dfrac n{\sqrt T}+T\right)\)，取 \(T=O(n^{\frac 23})\) 得到最优复杂度 \(O(n^{\frac 23})\)，和杜一样了！
给出洛谷 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\Leftrightarrow p_k^3\le n\) 时才需要在树状数组上查询 \(G_{k-1,n/p_k}\)，
则树状数组上的查询共有 \(O\left(\sum\limits_{i=1}^{n/T}\dfrac{\sqrt[3]{n/i}}{\log n/i}\right)=O\left(\dfrac n{T^\frac 23\log n}\right)\) 次， 复杂度 \(O\left(\dfrac n{T^\frac 23}\right)\)，
而 \(n\ge T\) 的转移就有 \(O\left(\dfrac n{\sqrt T\log n}\right)>O\left(\dfrac n{T^\frac 23}\right)\) 次，所以瓶颈在枚举所有状态进行转移，而不在树状数组。
树状数组上的修改贡献的复杂度仍为 \(O(T)\)，所以总复杂度 \(O\left(\dfrac n{\sqrt T\log n}+T\right)\)，
取 \(T=O\left(\left(\dfrac n{\log n}\right)^{\frac 23}\right)\) 得到最优复杂度 \(O\left(\left(\dfrac n{\log n}\right)^{\frac 23}\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\limits_{i=1}^nf(i)\)。
仿照上面的 \(G\)，我们设 \(F_{k,n}=\sum\limits_{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\limits_{p_k\le p_i\le n}\sum\limits_{c\ge 1,p_i^c\le n}f(p_i^c)F_{i+1,n/p_i^c}\)，
最小值因子 \(p_i\) 要枚举到 \(n\)，复杂度炸了。考虑把其中的质数拿出来单独算，则剩下的合数的最小质因子只需枚举到 \(\sqrt n\)，则
（最后一个等号成立的原因：考虑 \(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\limits_{i=1}^kg(p_i)\)，可以发现 \(p_k\le\sqrt n\) 所以 \(s\) 够用。
至于 \(\sum\limits_{p\in\mathbf P,p\le n}f(p)\)……这不是 Part I 吗？
直接按照转移式开搜，边界条件 \(p_k>n\) 时 \(F_{k,n}=0\)，记忆化都不用，复杂度可以证明是 \(O\left(\dfrac{n^{\frac 34}}{\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\limits_{p_{t_3}\le p_i\le\sqrt n}\left(f(p_i^2)+f(p_i)\left(\sum\limits_{j\in\mathbf P,j\le n/p_i}f(j)-s_i\right)\right)+\sum\limits_{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\limits_{i=1}^n[\text{lpf(i)}=p_k]f(i)=F_{k+1,n}+\sum\limits_{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 23}\) 时会枚举第一个 \(\sum\)，复杂度 \(O\left(\sum\limits_{i=1}^{\sqrt[3]n}\dfrac{\sqrt{n/i}}{\log n/i}\right)=O\left(\dfrac{n^{\frac 23}}{\log n}\right)\)。
接着，由 \(F_{t_3,i}\) 递推到 \(F_{t_6,i}\) 的过程与 Part I 中由 \(G_{t_6,i}\) 递推到 \(G_{t_3,i}\) 的过程分析方法一样，复杂度 \(O(n^{\frac 23})\)。
最后由 \(F_{t_6,i}\) 递推到 \(F_{1,i}\) 时，需要枚举 \(\sqrt[6]n\) 以内的素数在 \(n\) 以内的幂，以及所有块筛，
jijidawang 证明了前者有 \(O\left(\dfrac{\sqrt[6]n}{\log n}\right)\) 个，而块筛有 \(O(\sqrt n)\) 个，复杂度就是二者之积 \(O\left(\dfrac{n^{\frac 23}}{\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\) 满足 \(\forall p\in\mathbf P,g(p)=f(p)\)，已知 \(g\) 在 \(n\) 的所有块筛处的前缀和，求 \(\sum\limits_{i=1}^n F(i)\)。
做法：
称质因子次数最低二次的数为 Powerful Number（PN）。显然所有 PN 都能表示成 \(a^2b^3\) 的形式。
考虑 \(n\) 以内 PN 的个数。有 \(O(\sum\limits_{i=1}^{\sqrt n}\sqrt[3]{n/i^2})=O(\sqrt n)\) 个。
设积性函数 \(h\) 满足 \(g*h=f\)，则 \(\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{aligned} \sum\limits_{i=1}^nf(i)&=\sum\limits_{d|i}h(d)g\left(\dfrac id\right)\\ &=\sum\limits_{d=1}^nh(d)\sum\limits_{i=1}^{n/d}g(i)\\ &=\sum\limits_{d\in\mathbf{PN},d\le n}h(d)\sum\limits_{i=1}^{n/d}g(i) \end{aligned} \]爆搜 \(d\in\mathbf{PN},d\le n\)，统计答案即可。
注意跑 PN 筛的前提是已知 \(g\) 在 \(n\) 的所有块筛处的前缀和，
所以在 Min_25 筛板子题构造 \(g(p)=p\varphi(p)\) 并使用 PN 筛的复杂度不是 \(O(\sqrt n)\)，
因为求 \(g\) 在 \(n\) 的所有块筛处的前缀和需要跑一遍 \(O(n^{\frac 23})\) 的杜教筛。
