min_25 筛
已知一积性函数 \(f(x)\),\(f(p)\) 为关于 \(p\) 的多项式,且 \(f(p^k)\) 可以快速算出来,求 \(\sum_{i=1}^nf(i)\)。
例:求 \(\sum_{i=1}^n \sigma_0(i^k)\),其中 \(n = 10^{10}\)
step1 : 求质数的积性函数前缀和
求 \(1 ... n\) 以内的质数个数。
我们设 \(f(n,m)\) 表示 \(1...n\) 中的所有质数或者不含最小质因子为 \(p_1...p_m\) (\(p_i\) 表示第 \(i\) 个质数)的数的个数。
那么有:
\[f(n,m) = f(n,m-1) + (f(\lfloor \frac{n}{p_m} \rfloor,m-1) - (m-1)) \]意义为:考虑从 \(f(n,m-1)\) 处转移过来,那么我们只需要筛掉最小质因子恰好是 \(p_m\) 的数。那么我们把所有要筛的数提出来个 \(p_m\),剩下的数的最小质因子为 \(p_1...p_{m-1}\)。注意删掉“质数"部分。
更普遍的情况。
类似地,我们设 \(f(n,m)\) 表示 \(1...n\) 中的所有质数或者不含最小质因子为 \(p_1...p_m\) (\(p_i\) 表示第 \(i\) 个质数)的数的函数和。这时候我们需要认为合数也有一个假的和质数一样的函数和。那么有:
\[f(n,m) = f(n,m-1) + F(p_m)(f(\lfloor \frac{n}{p_m} \rfloor,m-1) - psum(m-1)) \](函数重名了,用 \(F(n)\) 代替了积性函数;\(psum\) 为预处理出来的质数函数前缀和)
由于是多项式,我们可以拆成 \(F(p) = p^k\) 的形式,最终拼起来。这样,假函数成为完全积性函数,直接提出来 \(F(p_m)\) 就是对的了。
复杂度证明
见那篇博客。
发现 \(p_m^2 \le n\)
可以积分证明,复杂度为 \(O(n^{0.75} / \log n)\)
step2 : 求所有数的积性函数和
我们设 \(g(n,m)\) 表示 \(1...n\) 中最小质因子大于 \(p_m\) 的数的(真)积性函数和。
那么有:
\[g(n,m) = f(n,m) - psum(m-1) + \sum_{p_k^e} F(p_k^e)(g(\lfloor \frac{n}{p^e} \rfloor,k) + [e \not= 1]) \]意义为:分质数,合数考虑贡献。质数的贡献我们已经搞好了,为 \(f(n,m) - psum(m-1)\),合数的话我们需要枚举最小质因子及其次数,即枚举 \(k > m\) 的 \(p_k\),再枚举次数 \(e\)。由于是积性函数,我们可以直接将 \(F(p_k^e)\) 提出来。值得注意的是,形如 \(p^k,k > 1\) 的数我们既没有在质数处枚举到,也没有在合数处提出来 \(p_k^e\) 后枚举到,所以需要加上个 \([e \not= 1]\)。
复杂度:不会证。状态数和上面一样是 \(O(n^{0.75} / \log n)\) 的,但是转移复杂度比较玄学,据说总复杂度是 \(O(n^{0.75} / \log n)\) 的,也有的说是 \(O(n^{1 - \epsilon})\) 的。
代码:
//divcntk
int block, itot, pos[N];
ull id[N], n, K, f[N];
inline int get_id(ull x) {
return pos[x <= block ? x : block + n / x];
}
inline void get_f() {
for (ull l = 1, r; l <= n; l = r + 1) {
ull zhi = n / l;
r = n / zhi;
id[++itot] = zhi;
pos[zhi <= block ? zhi : block + n / zhi] = itot;
f[itot] = zhi - 1;
}
for (int j = 1; 1ull * pri[j] * pri[j] <= n; ++j) {
for (int i = 1; id[i] >= 1ull * pri[j] * pri[j]; ++i) {
f[i] = f[i] - f[get_id(id[i] / pri[j])] + j - 1;
}
}
}
ull get_g(ull n, int j) {
if (pri[j] >= n) return 0;
ull res = (f[get_id(n)] - j) * (K + 1);//bug
for (int k = j + 1; 1ull * pri[k] * pri[k] <= n; ++k) {
ull tmp = pri[k];
for (int e = 1; tmp <= n; ++e, tmp = tmp * pri[k]) {
res += (K * e + 1) * (get_g(n / tmp, k) + (e > 1));
}
}
return res;
}
inline void work() {
read(n), read(K); block = sqrt(n);//bug
get_f();
ull res = get_g(n, 0) + 1;
printf("%llu\n", res);
}