简述
时隔一年之后再学一个新的筛法是不是没救了啊。
主要原理
min25 筛的主要思路来源于埃式筛法，通过神奇的 DP 数组进行预处理之后，再用很短的时间就可以求出整个前缀和了。最终的主要思路是快速计算 \(f(p^k)\) 来表示要算的积性函数 \(h(x)\)，所以 min25 筛求解 \(h(x)\) 可能需要计算多个不同的 \(f(p^k)\) 前缀和来进行对 \(h(x)\) 的组合表示。
首先需要用线性筛算出足够多的 \(pre[n] = \sum_{i = 1}^n prime_i^k\)。假设我们现在要求的一个 \(f(p^k)\) 为 \(f(x) = \sum_{i = 1} [i \in prime] i^k\)，那么我们可以来设一个 DP （别问我为什么要这么干，真的不知道咋想出来的）。
定义 \(g(n, i)\) 为前 \(n\) 个「质数或者是最小值因子大于 \(prime_i\) 的数」的 \(k\) 次方之和，可以理解为埃式筛法筛完 \(i\) 轮后剩下的局部解。最后 \(\sum_{i = 1}^n [i \in primes] i^k = g(n, \sqrt{n})\)。
乍一看如果直接 DP 的话我们的时间复杂度是 \(\Theta(n \sqrt n)\) 的，但是我们发现其实这个 DP 的转移不需要依赖于 \(n – 1\)，而是从 \(\lfloor \frac{n}{primes_i} \rfloor\) 转移。所以最后的层数最多 \(\sqrt n\)层，然而时间复杂度是 \( \Theta(\frac{n^{\frac{3}{4}}}{\log_2 n}) \)，我不会算（zzt 的论文写了证明，神仙可以去看看）。
如何进行转移呢？如果假设这个 \(primes_i^2 > n\) 了，那么显然这一轮埃筛只筛了空气出来；否则，我们需要认真思考下怎样去筛掉那些最小值因子恰好为 \(primes_i\) 的数的 \(k\) 次方。
首先，我们已经知道要筛掉的数是最小值因子恰好为 \(primes_i\) 的数的 \(k\) 次方。假设为 \(x^k\)，那么可以被表示为 \(primes_i^k x’^k\)。前面那个 \(primes_i^k\) 可以直接提出来，关键是怎么去算 \(x’^k\)。可以考虑直接拿 \(n\) 去暴力除掉一个 \(primes_i\)，拿 \(g(\lfloor \frac{n}{primes_i} \rfloor, i – 1)\) 即可。
不过要注意一点，最小值因子恰好为 \(primes_i\) 的数的 \(k\) 次方中，虽然除去了 \(primes_i)\) 得到 \(g(\lfloor \frac{n}{primes_i} \rfloor, i – 1)\)，但是也要注意保持条件「最小值因子恰好为 \(primes_i\) 的数」，我们需要从 \(g(\lfloor \frac{n}{primes_i} \rfloor, i – 1)\) 去掉前 \(i – 1\) 个质数的 \(k\) 次方来保持答案的正确性。
整理得到：
\[ g(n, i) = \begin{cases} g(n, i – 1), primes_i^2 > n \\ g(n, i – 1) – primes_i^k (g(\lfloor \frac{n}{primes_i} \rfloor, i – 1) – pre[i – 1]) \end{cases} \]
所以你会发现，那个 \(n\) 其实就被优化成了 \(\sqrt n\)。我们只需要存 \(\sqrt n\) 个即可。
算完这个之后，考虑算最后的答案。设计任务 \(solve(n, nid)\)，考虑当前要算的数的最小值因子编号 \(nid\)，然后先考虑把质数的答案算干净，为 \( g(n, maxPrimeId) – pre[nid – 1] \)。再考虑算合数的贡献，枚举下一个质因子 \(primes_j\)，因为是积性函数，所以从 \(solve(\frac{n}{primes_j}, j + 1)\) 拿到之后直接乘上一个 \(f(primes_j^k)\)。
总结为：
\[ solve(n, nid) = g(n, maxPrimeId) – pre[nid – 1] + \sum_{j = i}^{primes_j^2 \leq n} \sum_{k = 1}^{primes_j^{k + 1} \leq n} (solve(\lfloor \frac{n}{primes_j^k} \rfloor, nid + 1) * f(primes_j^k) + f(primes_j^{k + 1}) ) \]
代码
// min25.cpp
// This code file is a template of min25 sieve with rich comments for easier understanding;
// This code is based on the solution of problem https://www.luogu.com.cn/problem/P5325;
// Mathematical background:
// The goal is to get {\sum_{i = 1}^n f(i), i \in [1, 10^10]}, which f(p^k) = p^k(p^k - 1).
// Of course, the function f is multiplicative function.
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAX_N = 2e5 + 200, mod = 1e9 + 7;
// primes[] : used by linear sieve;
// primePre[] : \sum_{i = 0} [i \in prime] i;
// primePre[] : \sum_{i = 0} [i \in prime] i^2;
// idx[2][] : locations;
// facts[] & ftot : factor stored;
// g[2][] : the original g[][] with space saved;
int primes[MAX_N], tot, primePre[MAX_N], primePre2[MAX_N], idx[2][MAX_N], ftot, g[2][MAX_N];
ll facts[MAX_N];
bool vis[MAX_N];
ll n;
// helper function to locate the factor address;
struct locator
{
int &operator[](const ll &x) { return x < MAX_N ? idx[0][x] : idx[1][n / x]; }
} loc;
int fpow(int bas, int tim)
{
int ret = 1;
while (tim)
{
if (tim & 1)
ret = 1LL * ret * bas % mod;
bas = 1LL * bas * bas % mod;
tim >>= 1;
}
return ret;
}
const int inv2 = fpow(2, mod - 2), inv3 = fpow(3, mod - 2);
// preprocess;
void sieve()
{
for (int i = 2; i < MAX_N; i++)
{
if (!vis[i])
primes[++tot] = i, primePre[tot] = (0LL + primePre[tot - 1] + i) % mod,
primePre2[tot] = (0LL + primePre2[tot - 1] + 1LL * i * i % mod) % mod;
for (int j = 1; j <= tot && 1LL * i * primes[j] < MAX_N; j++)
{
vis[i * primes[j]] = true;
if (i % primes[j] == 0)
break;
}
}
}
int sieve(ll nbound, int nid)
{
if (primes[nid] >= nbound)
return 0;
// answer would be splited into 2 parts;
// the \sum_{i = 1} f(i) [i \in primes] part, and the rest part;
int id = loc[nbound], nppart = 0, ppart = (0LL + g[1][id] + mod - g[0][id] + mod - (0LL + primePre2[nid] + mod - primePre[nid]) % mod) % mod;
for (int i = nid + 1; i <= tot && 1LL * primes[i] * primes[i] <= nbound; i++)
for (ll j = 1, acc = primes[i]; acc <= nbound; j++, acc *= primes[i])
nppart = (0LL + nppart + 1LL * acc % mod * ((acc - 1) % mod) % mod * ((0LL + sieve(nbound / acc, i) + (j != 1)) % mod) % mod) % mod;
return (0LL + nppart + ppart) % mod;
}
int main()
{
scanf("%lld", &n), sieve();
for (ll i = 1, gx; i <= n; i = gx + 1)
{
gx = n / (n / i);
int id = ++ftot;
facts[id] = n / i;
// x^1;
g[0][ftot] = (1LL * facts[id] % mod * ((facts[id] + 1) % mod) % mod * inv2 % mod + mod - 1) % mod;
// x^2;
g[1][ftot] = (1LL * facts[id] % mod * ((facts[id] + 1) % mod) % mod * inv2 % mod * (2LL * facts[id] % mod + 1) % mod * inv3 % mod + mod - 1) % mod;
loc[facts[id]] = id;
}
// do dp;
for (int i = 1; i <= tot; i++)
for (int j = 1; j <= ftot && 1LL * primes[i] * primes[i] <= facts[j]; j++)
{
int pos = loc[facts[j] / primes[i]];
g[0][j] = (0LL + g[0][j] + mod - 1LL * primes[i] * ((0LL + g[0][pos] + mod - primePre[i - 1]) % mod) % mod) % mod;
g[1][j] = (0LL + g[1][j] + mod - 1LL * primes[i] * primes[i] % mod * ((0LL + g[1][pos] + mod - primePre2[i - 1]) % mod) % mod) % mod;
}
printf("%lld\n", (1LL + sieve(n, 0)) % mod);
return 0;
}
例题 A：经典积性函数筛 – BZOJ 3944 Sum
先考虑欧拉函数 \(\varphi(x)\) 的前缀和。发现 \(\varphi(p^k) = p^k – p^{k – 1} = p^{k – 1}(p – 1)\)，所以我们只需要维护质数积和即可。
再考虑莫比乌斯函数 \(\mu(x)\)，其中发现合数肯定都是 \(0\)，唯一贡献的只有 \(x = \prod p_i\)，所以只需要算区间内的质数个数即可，因为只有 \(\mu(p^k), k = 1\) 有贡献。
// BZ3944.cpp
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAX_N = 2e5 + 200;
int T, primes[MAX_N], tot, idx[2][MAX_N];
ll g[2][MAX_N], pre[2][MAX_N], n, facts[MAX_N], ftot;
bool vis[MAX_N];
struct locator
{
int &operator[](const ll &x) { return x < MAX_N ? idx[0][x] : idx[1][n / x]; }
} loc;
void sieve()
{
for (int i = 2; i < MAX_N; i++)
{
if (!vis[i])
primes[++tot] = i, pre[0][tot] = pre[0][tot - 1] + i, pre[1][tot] = -tot;
for (int j = 1; j <= tot && 1LL * i * primes[j] < MAX_N; j++)
{
vis[i * primes[j]] = true;
if (i % primes[j] == 0)
break;
}
}
}
ll sieve_phi(ll nbound, int nid)
{
if (primes[nid] >= nbound || nbound <= 1)
return 0;
int id = loc[nbound];
ll ret = g[0][id] - g[1][id] - pre[0][nid] + nid;
for (int i = nid + 1; i <= tot && 1LL * primes[i] * primes[i] <= nbound; i++)
for (ll acc = primes[i] * primes[i], phi_acc = primes[i] - 1; acc <= nbound;)
{
ret += 1LL * phi_acc * sieve_phi(nbound / (acc / primes[i]), i) + phi_acc * primes[i];
acc *= primes[i], phi_acc *= primes[i];
}
return ret;
}
ll sieve_mu(ll nbound, int nid)
{
if (primes[nid] >= nbound || nbound <= 1)
return 0;
int id = loc[nbound];
ll ret = -(g[1][id] - nid);
for (int i = nid + 1; i <= tot && 1LL * primes[i] * primes[i] <= nbound; i++)
ret -= sieve_mu(nbound / primes[i], i);
return ret;
}
void solve()
{
ftot = 0;
for (ll i = 1, gx; i <= n; i = gx + 1)
{
gx = n / (n / i);
int id = ++ftot;
facts[id] = n / i;
// prime sum & number;
g[0][ftot] = (facts[id] + 2) * (facts[id] - 1) / 2, g[1][ftot] = facts[id] - 1;
// prime num;
loc[facts[id]] = id;
}
for (int i = 1; i <= tot; i++)
for (int j = 1; j <= ftot && 1LL * primes[i] * primes[i] <= facts[j]; j++)
{
int pos = loc[facts[j] / primes[i]];
g[0][j] -= 1LL * primes[i] * (0LL + g[0][pos] - pre[0][i - 1]);
g[1][j] -= g[1][pos] - i + 1;
}
if (n == 0)
puts("0 0");
else
printf("%lld %lld\n", sieve_phi(n, 0) + 1, sieve_mu(n, 0) + 1);
}
int main()
{
scanf("%d", &T), sieve();
while (T--)
scanf("%lld", &n), solve();
return 0;
}