前置知识
powerful number
定义：\(\text{powerful number}\) 是没有 \(1\) 次质因子的数，也简称为 \(\text{PN}\)。
结论：\(n\) 以内的 \(\text{PN}\) 个数为 \(O(\sqrt n)\)。
证明：
所有 \(\text{PN}\) 必然可表示成 \(a^2b^3\) 的形式。
那么：
\[\sum_{a=1}^{\sqrt n}(\frac{n}{a^2})^{\frac{1}{3}}=\sqrt n \]构造
构造 \(G(p^k)=F(p)^k,G\times H=F\)，可以注意到 \(G(P)\) 为完全积性函数。
那么：
\[\begin{aligned} F(p)&=G(1)H(p)+H(1)G(p)\\ &=G(p)+H(p) \end{aligned}\]所以，\(H(p)=0\)。
\[\begin{aligned} \sum_{i=1}^n F(i)&=\sum_{i=1}^n\sum_{j|i}G(j)H(\frac{i}{j})\\ &=\sum_{i=1}^nH(i)\sum_{j=1}^{\frac{n}{i}}G(j)\\ &=\sum_{i=1}^n[i\in PN]H(i)SG_{\frac{n}{i}}\\ \end{aligned}\]其中:
\[\begin{aligned} H(p)&=0\\ H(p^2)&=F(p^2)-F(p)^2\\ H(p^3)&=F(p^3)-F(p)F(p^2)\\ H(p^k)&=F(p^k)-F(p)F(p^{k-1}) \end{aligned}\]所以，我们只需要求出 \(G\) 的块筛即可。
块筛
考虑两种方式。
法一
若 \(G\) 可杜教筛。
可以直接使用杜教筛。
此时，这个筛法就是常称的 \(\text{PN}\) 筛。
法二
若 \(F(p)\) 为一个低次多项式。
考虑 \(S(n,a)\)。
\[\begin{aligned} S(n,a)&=\sum_{i=1}^n[i\in P ~or~ min_i > P_a] G(i)\\ &=S(n,a-1)-g(p_a)(S(\frac{n}{P_a},a-1)-S(P_a-1,a-1))\\ S(n,a-1)&=S(n,a)+g(p_a)(S(\frac{n}{P_a},a-1)-S(P_a-1,a-1))\\ \end{aligned}\]要求：\(G(x)\) 为完全积性函数。
若 \(\sum_{i=1}^n [i\in P]G(x)\) 很好求出，那么直接用第二个递推式即可。
否则，将多项式 \(G(x)\) 的每一项分别提出来，跑第一个递推式，合并后再用第二个递推式即可。
本质是 \(\text{min-25}\) 筛的前半部分。
完美利用了 \(G(x)\) 的性质。
个人更加推荐法二，这也可以看出这个筛法极强的拓展性。
速度
考虑 \(\text{PN}\) 的过程复杂度为 \(O(\sqrt n)\)。
块筛的过程要么是 \(O(n^{\frac{2}{3}})\) 或者 \(O(\frac{n^{\frac{3}{4}}}{\log n})\)（很可能不准）。
实际效率上比普通 \(\text{min-25}\) 筛快一点（板子题）。
\(10^{10}\)：\(0.2s\)。
\(10^{11}\)：\(0.6s\)。
\(10^{12}\)：\(3.0s\)。
代码难度个人认为比较简单，与杜教筛难度近似，感觉可以作为 \(\text{min-25}\) 筛的上位替代。
例题
P4213 【模板】杜教筛
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
const int N = 1e5 + 10;
i64 n, ans1, ans2, w[N];
i64 s1[N], s2[N], p1[N], p2[N];
int m, ct, tt, p[N], v[N], d1[N], d2[N];
inline void init_all(int m = N - 10) {
for (int i = 2; i <= m; i++) {
if (!v[i]) p[++ct] = i, p1[ct] = 1 + p1[ct - 1], p2[ct] = i + p2[ct - 1];
for (int j = 1; j <= ct && p[j] * i <= m; j++) {
v[p[j] * i] = 1; if (i % p[j] == 0) break;
}
}
}
inline int get(i64 x) { return (x <= m ? d1[x] : d2[n / x]); }
inline void init() {
m = sqrt(n) + 100, ct = tt = ans1 = ans2 = 0;
for (i64 l = 1, r; l <= n; l = r + 1) {
r = (n / (n / l)), w[++tt] = n / l;
w[tt] <= m ? d1[w[tt]] = tt : d2[r] = tt;
s1[tt] = w[tt], s2[tt] = (w[tt] + 1) * w[tt] / 2;
s1[tt]--, s2[tt]--;
}
while (1ll * p[ct + 1] * p[ct + 1] <= n) ct++;
for (int i = 1; i <= ct; i++) {
const i64 k1 = p[i];
const i64 k2 = k1 * k1;
int lim = 0;
while (lim < tt && k2 <= w[lim + 1]) lim++;
for (int j = 1; j <= lim; j++) {
if (k2 > w[j]) break;
int id = get(w[j] / k1);
s1[j] = (s1[j] - (s1[id] - p1[i - 1]));
s2[j] = (s2[j] - k1 * (s2[id] - p2[i - 1]));
}
}
for (int i = 1; i <= tt; i++) s1[i] = -s1[i], s2[i] = s2[i] + s1[i];
for (int i = 1; i <= tt; i++) p1[i] = -p1[i], p2[i] = p2[i] + p1[i];
for (int i = ct; i >= 1; i--) {
const i64 k1 = p[i];
const i64 k2 = k1 * k1;
int lim = 0;
while (lim < tt && k2 <= w[lim + 1]) lim++;
for (int j = lim; j >= 1; j--) {
int id = get(w[j] / k1);
s1[j] = (s1[j] + -1 * (s1[id] - p1[i - 1]));
s2[j] = (s2[j] + (k1 - 1) * (s2[id] - p2[i - 1]));
}
}
for (int i = 1; i <= tt; i++) s1[i]++, s2[i]++;
for (int i = 1; i <= tt; i++) p2[i] = p2[i] - p1[i], p1[i] = -p1[i];
}
inline void dfs(int now, i64 cur, i64 h1, i64 h2, i64&ans1, i64&ans2) {
ans1 = (ans1 + h1 * s1[get(n / cur)]);
ans2 = (ans2 + h2 * s2[get(n / cur)]);
for (int i = now; i <= ct; i++) {
const i64 k = p[i];
if (k * k > n / cur) break;
for (i64 x = k * k, y = 1, z = 2; x * cur <= n; x = x * k, y = y * k, z++) {
dfs(i + 1, cur * x, h1 * (z == 2 ? -1 : 0), h2 * (k - 1) * y, ans1, ans2);
}
}
}
int main() {
int t;
cin >> t, init_all();
while (t--) {
cin >> n, init();
dfs(1, 1, 1, 1, ans1, ans2);
cout << ans2 << " ";
cout << ans1 << "\n";
for (int i = 1; i <= tt; i++) s1[i] = s2[i] = d1[i] = d2[i] = w[i] = 0;
}
return 0;
}
P5325 【模板】Min_25 筛
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
const int N = 1e6 + 10;
const int mod = 1e9 + 7;
i64 n, ans, w[N];
i64 s1[N], s2[N], p1[N], p2[N];
int m, ct, tt, p[N], v[N], d1[N], d2[N];
inline int get(i64 x) { return (x <= m ? d1[x] : d2[n / x]); }
inline void init() {
m = sqrt(n) + 100;
for (i64 l = 1, r; l <= n; l = r + 1) {
r = (n / (n / l)), w[++tt] = n / l;
w[tt] <= m ? d1[w[tt]] = tt : d2[r] = tt;
s1[tt] = (__int128)(w[tt] + 1) * w[tt] * 500000004 % mod;
s2[tt] = (__int128)(2 * w[tt] + 1) * s1[tt] * 333333336 % mod;
s1[tt]--, s2[tt]--;
}
for (int i = 2; i <= m; i++) {
if (!v[i]) p[++ct] = i, p1[ct] = (i + p1[ct - 1]) % mod, p2[ct] = (1ll * i * i + p2[ct - 1]) % mod;
for (int j = 1; j <= ct && p[j] * i <= m; j++) {
v[p[j] * i] = 1;
if (i % p[j] == 0) break;
}
}
for (int i = 1; i <= ct; i++) {
const i64 k1 = p[i];
const i64 k2 = k1 * k1;
int lim = 0;
while (lim < tt && k2 <= w[lim + 1]) lim++;
for (int j = 1; j <= lim; j++) {
if (k2 > w[j]) break;
int id = get(w[j] / k1);
s1[j] = (s1[j] - k1 % mod * (s1[id] - p1[i - 1])) % mod;
s2[j] = (s2[j] - k2 % mod * (s2[id] - p2[i - 1])) % mod;
}
}
for (int i = 1; i <= tt; i++) s1[i] = (s2[i] - s1[i]) % mod;
for (int i = 1; i <= ct; i++) p1[i] = (p2[i] - p1[i]) % mod;
for (int i = ct; i >= 1; i--) {
const i64 k1 = p[i];
const i64 k2 = k1 * k1;
int lim = 0;
while (lim < tt && k2 <= w[lim + 1]) lim++;
for (int j = lim; j >= 1; j--) {
int id = get(w[j] / k1);
s1[j] = (s1[j] + k1 * (k1 - 1) % mod * (s1[id] - p1[i - 1])) % mod;
}
}
for (int i = 1; i <= tt; i++) s1[i]++;
}
inline void dfs(int now, i64 cur, i64 h, i64&ans) {
ans = (ans + h * s1[get(n / cur)]) % mod;
for (int i = now; i <= ct; i++) {
const i64 k = p[i];
if (k * k > n / cur) break;
for (i64 x = k * k, y = k; x * cur <= n; x = x * k, y = y * k) {
dfs(i + 1, cur * x, (__int128)h * x * (y + k - 2) % mod, ans);
}
}
}
int main() {
cin >> n, init();
dfs(1, 1, 1, ans);
cout << (ans + mod) % mod << "\n";
return 0;
}