质数前缀统计
求出 \(n\) 以内所有质数的 \(c\) 次方之和。
考虑埃氏筛，每次用一个质数枚举其的所有倍数，筛去所有不合法的数。
一些定义以及解释
- \(p_{\min}(x)\)：\(x\) 的最小质因子，\(p_{\min}(1) = +\infty\)，如果 \(x\) 非质数，那么 \(p_{\min}(x) \le \sqrt{x}\)。
- \(m：\left\lfloor \sqrt{n}\right\rfloor\) 以内质数的个数。
- \(p_k\)：第 \(k\) 个质数。
- \(P_k\)：前 \(k\) 个质数构成的集合。
- \(S(n, k)\)：\(n\) 以内的数用前 \(k\) 个质数删除了不合法的数后，剩下的数的集合，也就是 \(S(n, k) = \{x \in \mathbb{N}^* | x \le n \wedge p_{\min}(x) \ge p_k\}\)。同时，我们注意到，对于 \(S(n, m)\) 就是所有还没有遍历的质数。
- \(D(n, k)\)：\(n\) 以内的数中用第 \(k\) 个质数删去的数的集合，即 \(D(n, k) = S(n, k - 1) - S(n, k)\)。
- \(f(n, k)\)：\(\displaystyle \sum_{x \in S(n, k)} x^c\)，于是题目要求的答案就是 \(f(n, m) - 1 + \displaystyle \sum_{x \in P_m} x^c\)。
推导过程
\(S(n, k) = S(n, k - 1) - p_k S(\left\lfloor \frac{n}{p_k}\right\rfloor , k - 1)\)
也就是考虑被 \(p_k\) 删除的数，这个可以根据 \(p_k S(\left\lfloor \frac{n}{p_k}\right\rfloor , k - 1)\) 与 \(D(n, k)\) 是个双射来证明。
于是，我们可以得到 \(f\) 的计算式子：
\(f(n, k) = f(n, k - 1) - p_k^c f(\left\lfloor\frac{n}{p_k}\right\rfloor, k - 1)\)
如果我们对于每个质数都直接用上面的式子计算一遍，那么总共有 \(\frac{\sqrt{n}}{\log n}\) 个质数，而根据整除分块，我们实际上要计算的值有 \(\sqrt{n}\) 个，于是复杂度就是 \(\mathcal O(\frac{n}{\log n})\) 了！
我们可以考虑优化，我们会注意到一个这样的事情，对于 \(k > m\)，\(S(n, k) = S(n, k - 1) - \{p_k\}\)，因为 \(S\) 还在改变，于是我们还要计算那些 \(k > m\) 的值，如果我们让 \(S\) 不改变，那么是否会有一些优化呢？
我们可以定义一些新的东西：
- \(T(n, k) = S(n, k) +P_k\)
- \(g(n, k) = f(n, k) + \displaystyle \sum_{x \in P_m} x^c\)，最后答案就是 \(g(n, m) - 1\)
于是我们可以利用之前 \(S\) 的式子得到：
\(T(n, k) = T(n, k - 1) - p_k \left(T(\left\lfloor \frac{n}{p_k}\right\rfloor , k - 1) - P_{k - 1}\right) + \{p_k\}\)
又我们注意到，\(T(p_k - 1, k - 1)\) 就是 \(P_{k - 1} + {1}\)，因此我们也可以得到 \(g\) 的式子：
\(g(n, k) = g(n, k - 1) - \left(g(\left\lfloor \frac{n}{p_k}\right\rfloor , k - 1) - g(p_k - 1, k - 1)\right)\)
至于 \(g(n, 0) = \sum_{i \le n} i^c\)，可以在 \(\mathcal O(c)\) 的时间里用拉格朗日插值求出。
然后复杂度就是 \(\mathcal O(\dfrac{n^{3 / 4}}{\log n})\)，不太会证。
模板 tea：Luogu5493 质数前缀统计
#define in read<int>()
#define lin read<ll>()
#define rep(i, x, y) for(int i = (x); i <= (y); i++)
#define per(i, x, y) for(int i = (x); i >= (y); i--)
using namespace std;
using ll = long long;
template < typename T > void chkmax(T &x, const T &y) { x = x > y ? x : y; }
template < typename T > void chkmin(T &x, const T &y) { x = x < y ? x : y; }
const int N = 1e6 + 10;
const int inf = 0x7fffffff;
int mod, C;
ll lim, n;
ll g0[N], g1[N];
// g0[i] 是 g(i), g1[i] 是 g(n / i)
ll qp(ll x, int t) { x %= mod; ll res = 1; for(; t; t = t >> 1, x = x * x % mod) if(t & 1) res = res * x % mod; return res; }
namespace calcg {
const int N = 20;
int y[N], fac[N], ifac[N], pre[N], suf[N], m;
void init() {
m = C + 2;
fac[0] = 1; rep(i, 1, m) fac[i] = 1ll * fac[i - 1] * i % mod;
ifac[m] = qp(fac[m], mod - 2);
per(i, m - 1, 0) ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
rep(i, 1, m) y[i] = qp(i, C), y[i] = (y[i] + y[i - 1]) % mod;
}
int calc(ll v) {
v %= mod;
pre[0] = 1; rep(i, 1, m) pre[i] = 1ll * pre[i - 1] * (v - i + mod) % mod;
suf[m + 1] = 1; per(i, m, 1) suf[i] = 1ll * suf[i + 1] * (v - i + mod) % mod;
ll res = 0;
rep(i, 1, m) {
res += 1ll * pre[i - 1] * suf[i + 1] % mod * ifac[i - 1] % mod
* ifac[m - i] % mod * y[i] % mod * (m - i & 1 ? mod - 1 : 1) % mod;
res %= mod;
} return res;
}
}
using calcg :: calc;
int main() {
n = lin, C = in, mod = in, lim = sqrtl(n);
calcg :: init();
rep(i, 1, lim) g0[i] = (calc(i) - 1) % mod, g1[i] = (calc(n / i) - 1) % mod;
rep(i, 2, lim) {
g0[i] = (g0[i] + mod) % mod; if(g0[i] == g0[i - 1]) continue;
ll pc = qp(i, C), pg = g0[i - 1];
int l1 = lim / i, l2 = min(n / i / i, lim); chkmin(l1, l2);
rep(j, 1, l1) g1[j] = (g1[j] - pc * (g1[j * i] - pg)) % mod;
ll tt = n / i;
if(tt < inf) { // 卡常
int t = tt;
rep(j, l1 + 1, l2) g1[j] = (g1[j] - pc * (g0[t / j] - pg)) % mod;
} else {
ll t = tt;
rep(j, l1 + 1, l2) g1[j] = (g1[j] - pc * (g0[t / j] - pg)) % mod;
}
ll st = 1ll * i * i;
per(j, lim, st) g0[j] = (g0[j] - pc * (g0[j / i] - pg)) % mod;
}
ll ans = 0;
rep(i, 1, lim) {
g1[i] = (g1[i] + mod) % mod;
ans += 1ll * i * i % mod * g1[i] % mod; ans %= mod;
} printf("%lld\n", ans);
return 0;
}
当然，这个算法可以进行一波魔改，然后就可以做 LOJ6028「from CommonAnts」质数计数 II 了，代码。
Min_25 筛
给出一个积性函数 \(F\)：
- \(F(p)\) 为关于 \(p\) 的简单多项式。
- \(F(p^k)\) 能够较快速求出。
- 求出 \(\sum_{i = 1}^n F(i)\) 。
按照 \(k\) 从大到小的顺序，考虑 \(p_{\min} = k\) 的数，那么记 \(S(n, k) = \{x \in \mathbb{N}^* | x \le n \wedge p_{\min}(x) \ge p_k\}\)，记 \(f(n, k) = \sum_{x \in S(n, k)} F(x)\)，最后答案就是 \(f(n, 1)\)。
于是有 \(\displaystyle f(n, k) = \sum_{i = k}^{p_i \le n} \sum_{v} F(p_i ^ v)f\left(\left\lfloor\frac{n}{p_i ^ v}\right\rfloor, k+ 1\right)\)，当 \(p_i \ge \sqrt{n}\) 的时候，后面这个式子的值就是 \(F(p_i)\)，于是可以用之前预处理的质数前缀和解决。
复杂度为 \(\mathcal O(\dfrac{n^{3 / 4}}{\log n})\)，不太会证。
比上面的更好的写法 \(\downarrow\)
#include <bits/stdc++.h>
#define in read<int>()
#define lin read<ll>()
#define rep(i, x, y) for(int i = (x); i <= (y); i++)
#define per(i, x, y) for(int i = (x); i >= (y); i--)
using namespace std;
using ll = long long;
const int N = 1e6 + 10;
const int mod = 1e9 + 7;
const int inv2 = (mod + 1) >> 1;
const int inv3 = (mod + 1) / 3;
ll n, lim, sp1[N], sp2[N], g1[N], g2[N], w[N];
int ind1[N], ind2[N], tot;
int pri[N / 10], pnum;
bool v[N];
void shai(int l) {
rep(i, 2, l) {
if(!v[i]) {
pri[++pnum] = i;
sp1[pnum] = (sp1[pnum - 1] + i) % mod;
sp2[pnum] = (sp2[pnum - 1] + 1ll * i * i) % mod;
}
rep(j, 1, pnum) {
if(pri[j] * i > l) break;
v[pri[j] * i] = true; if(i % pri[j] == 0) break;
}
}
}
int S(ll x, int y) {
if(x < pri[y]) return 0;
ll res = 0;
res += g2[x <= lim ? ind1[x] : ind2[n / x]] - sp2[y - 1];
res -= g1[x <= lim ? ind1[x] : ind2[n / x]] - sp1[y - 1];
res = (res % mod + mod) % mod;
rep(k, y, pnum) {
if(1ll * pri[k] * pri[k] > x) break;
ll cur = pri[k];
for(int c = 1; cur <= x; cur *= pri[k], c++) {
res += cur % mod * ((cur - 1) % mod) % mod * (S(x / cur, k + 1) + (c != 1)) % mod;
res %= mod;
} // 枚举计算的次数，当 c!=1 时，应该包括 1，c = 1 的时候，我们已经计算了。
} return res % mod;
}
int main() {
n = lin, lim = sqrtl(n);
shai(lim * 2); // 稍微多筛一点
for(ll i = 1; i <= n; i++) {
++tot, w[tot] = (n / i); ll v = w[tot] % mod;
g1[tot] = (v * (v + 1) / 2 - 1) % mod;
g2[tot] = (v * (v + 1) % mod * (v * 2 + 1) % mod * inv2 % mod * inv3 % mod - 1) % mod;
if(w[tot] <= lim) ind1[w[tot]] = tot; else ind2[n / w[tot]] = tot;
i = n / (n / i);
}
rep(i, 1, pnum) {
int p = pri[i]; ll x1 = sp1[i - 1], x2 = sp2[i - 1];
rep(x, 1, tot) {
if(1ll * p * p > w[x]) break;
int y = w[x] / p <= lim ? ind1[w[x] / p] : ind2[n / (w[x] / p)];
g1[x] = (g1[x] - 1ll * p * (g1[y] - x1) % mod + mod) % mod;
g2[x] = (g2[x] - 1ll * p * p % mod * (g2[y] - x2) % mod + mod) % mod;
}
}
printf("%d\n", (S(n, 1) + 1) % mod);
return 0;
}