Min_25 筛
前言
参考了贺了神 \(\color {balck} {\text z}\color{red}{\text{houkangyang}}\) 的blog。
Min_25 筛用于在亚线性时间复杂度（据说是 \(\mathcal O(\frac {n ^ {\frac 3 4}} {\log n})\)）下处理积性函数 \(f(x)\) 的前缀和，其中 \(f(x)\) 在质数处的值需满足是一个关于 \(x\) 低次多项式，并且在质数的幂次处可以快速算出 \(f(x)\) 的值，我们将其按 \(x\) 的系数拆分成多个单项式计算前缀和。
求 \(\sum_{i = 1} ^ N f(i)\)（注意下文 \(n\) 不是此处的 \(N\)）。
记号
定义 \(\mathcal P\) 为质数全集，
\(P_i\) 为第 \(i\) 个质数，
\(q(n)\) 为最小的 \(i\) 满足 \(P_i | n\)，
\(\frac n m\) 均为 \(\lfloor \frac n m \rfloor\)
思路
发现合数均可以变成多个质数的幂次的函数的积的形式，因此对于 \(i\) 分为质数和合数分类讨论，并以动态规划的思路求解。
先不管状态的空间和转移时间限制，以及不管 \(f(1)\) 的贡献，设 dp 状态 \(S(n, m) = \sum_{i = 2} ^ n f(i)\times[q(i) > m]\)。
我们的目标就是 \(S(N, 0) + f(1)\)。
质数
定义 \(g(n) = \sum_{i = 1} ^ n f(i) \times [i \in \mathcal P]\)。
质数在 \(S(n, m)\) 的贡献就是：\(g(n) - \sum_{i = 1} ^ m f(P_i)\)，在后面的计算过程中，因为我们会保证 \(m \le \sqrt n\)，所以这里 \(\sum_{i = 1} ^ m f(P_i)\) 可以前缀和简单预处理。
如何处理 \(g(n)\)，因为这是个积性函数，在合数处的取值不好处理，而我们只需要质数处的前缀和，所以若 \(f(x) = ax^k~(x\in \mathcal P)\) 我们设完全积性函数 \(f'(x)\) 满足 \(f'(x) = ax^k\)，只要保证质数处的取值满足 \(f(x) = f'(x)\)，我们再用动态规划建立状态 \(T(n, m) = \sum_{i = 1} ^ n f'(n) \times [i \in \mathcal P \or q(i) > m]\)。
观察如何从 \(T(n, m - 1)\) 转移到 \(T(n, m)\)，发现若 \(P_m > \sqrt n\)，两者是完全相等的，若 \(P_m \le \sqrt n\)，那么在 \(T(n, m - 1)\) 比 \(T(n, m)\) 会多出一部分贡献 \(f(x)~(x\notin \mathcal P\and q(x) = m)\)，于是就要减去这部分贡献，即 \(f(P_m) \times T(\frac n {P_m}, m - 1)\)，但是 \(T(\frac n {p_m}, m - 1)\) 会包含 \(\sum_{i = 1} ^ {m - 1} f'(P_i)\)，需要把这部分去除。
\(T(n, 0) = \sum_{i = 2} ^ n f(n)\)，最后的结果就是 \(g(n) = T(n, \text {maxp})\)。
合数
比较好处理，枚举合数的最小质因数和它的幂次直接转移，因为这个质因数上界不超过 \(\sqrt n\)，每次转移至少使 \(n\) 变成 \(\frac n 2\)，所以这里的复杂度不会超过 \(\mathcal O(\sqrt N \log N)\) （还能算到更低，但我不会qwq），有 \(S(n, m) = \sum_{i = m + 1} ^ {P_i \le \sqrt n}\sum_{j = 1} ^{P_i ^ j \le n} f(P_i ^ j)(S(\frac n {P_i ^ j}, i) + [j > 1])\)，注意这里特判 \(P_i ^ j\) 也是合数的情况。
合起来
实现
算 \(S(n, m)\) 的时候可以直接递归计算，但预处理 \(g(n)\) 即 \(T(n, m)\) 时无法算完所有状态。
但是我们发现递归计算 \(S(n, m)\) 的之后利用到的 \(g(n)\) 并不多，因为 \(\frac {\frac N m} k = \frac N {mk}\)，而对于 \(\frac N m\) 最多只有 \(2\sqrt N\) 个不同的数，我们只要计算出这 \(2\sqrt N\) 个不同的 \(g(n)\) 即可，在计算 \(T(n, m)\) 的时候也是同样 \(2\sqrt N\) 个不同的 \(g(n)\)，这样可以大大减少处理状态数，并且发现 \(T(n, m)\) 全由 \(T(k, m - 1)\) 转移而来，枚举 \(m\) 并将第二维压缩掉。
将这 \(2\sqrt N\) 个不同的数分类，其中 \(x \le \sqrt N\) 的存放在 id1[x]，\(x > \sqrt N\) 的存放在 id2[N / x]。
int a, pm[N + 5], top, g1[N + 5], g2[N + 5], id1[N + 5], id,
id2[N + 5], s[N + 5], sm1[N + 5], sm2[N + 5], B;
bool vis[N + 5];
void init() {
rep(i, 2, N) {
if(!vis[i]) pm[++top] = i;
for(int j = 1; j <= top && i * pm[j] <= N; ++j) {
vis[i * pm[j]] = 1;
if(i % pm[j] == 0) break;
}
}
}
const int inv2 = (mod + 1) / 2, inv6 = (mod + 1) / 6;
int s1(int n) {n %= mod; return (1 + n) * n % mod * inv2 % mod;}
int s2(int n) {n %= mod; return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;}
int get(int n) {return n <= B ? id1[n] : id2[a / n];}
int S(int n, int m) {
if(pm[m] > n) return 0;
int res = ((g2[get(n)] - g1[get(n)] + mod) - (sm2[m] - sm1[m])) % mod;
for(int i = m + 1; pm[i] * pm[i] <= n; ++i)
for(int k = pm[i]; k <= n; k *= pm[i])
(res += k % mod * (k % mod - 1) % mod * (S(n / k, i) + (k != pm[i])) % mod) %= mod;
return res;
}
signed main() {
// freopen("in1.in", "r", stdin);
init();
a = read();
B = sqrt(a);
for(int l = 1, r; l <= a; l = r + 1) {
r = a / (a / l);
s[++id] = a / l;
g1[id] = s1(s[id]) - 1;
g2[id] = s2(s[id]) - 1;
if(s[id] <= B) id1[s[id]] = id;
else id2[a / s[id]] = id;
}
rep(i, 1, top) {
sm1[i] = (sm1[i - 1] + pm[i]) % mod;
sm2[i] = (sm2[i - 1] + pm[i] * pm[i] % mod) % mod;
}
rep(i, 1, top) {
for(int j = 1; pm[i] * pm[i] <= s[j]; ++j) {
g1[j] = (g1[j] - pm[i] * (g1[get(s[j] / pm[i])] - sm1[i - 1] + mod) % mod + mod) % mod;
g2[j] = (g2[j] - pm[i] * pm[i] % mod * (g2[get(s[j] / pm[i])] - sm2[i - 1] + mod) % mod + mod) % mod;
}
}
cout << (S(a, 0) + 1) % mod;
return 0;
}