min_25 筛

一个亚线性筛，复杂度大概是 $O(n34log⁡n)O(\frac{n ^{\frac{3}{4}}}{ \log n})$ 。

使用 $min\_25$ 求前缀和，有两个基本特征：① 积性函数，② 满足质数点为多项式。

算法思路

给定 $\leq 10 ^ {11}$ ，设 $f (n)$ 为一个积性函数，质数点为多项式，求 $∑i=1nf(i)\sum\limits_{i = 1} ^{n} f(i)$ 。

先求出 $n\sqrt n$ 内的所有质数，以及 $f (i)$ 。
求出 $∑i=1n[i∈prime]f(i)\sum\limits_{i = 1} ^{n} [i \in prime]f(i)$ ， $p_i$ 为第 $i$ 个质因子。

我们设 $g(n,j)=∑i=1n[i∈primeormin_prime_of_i>pj]f(i)g(n, j) = \sum\limits_{i = 1} ^{n}[i \in prime\ or\ min\_prime\_of\_i > p_j]f(i)$ ，质数或者最小质因子大于 $p_j$ 的 $f$ 值之和。

考虑从 $g (n, i - 1)$ 推到 $g (n, i)$ ，我们也就是要减去最小质因子等于 $p_i$ 的贡献，
$\left\{\begin{matrix}\\ g(n, i - 1)& p_i \times p_i > n\\ g(n, i - 1) - f(p_i) \times \left(g(\frac{n}{p_i}, i - 1) - \sum\limits_{j = 1} ^{i - 1}f(p_j) \right) & p_i \times p_i \leq n\\ \end{matrix}\right.$
考虑如何递归计算积性函数前缀和 $∑i=1nf(i)\sum\limits_{i = 1} ^{n} f(i)$ 。

设 $S(n,j)=∑i=1n[min_prime_of_i≥pj]f(i)S(n, j) = \sum\limits_{i = 1} ^{n} [min\_prime\_of\_i \geq p_j] f(i)$ ，最小质因子大于等于 $p_j$ 的 $f$ 值之和，所以这个积性函数的前缀和即为 $S (n, 1) + f (1)$ 。

对于 $S (n, j)$ ，中质数的贡献我们是已知的 $\sum\limits_{i = 1} ^{n}[i \in prime]f(i) - \sum\limits_{i = 1} ^{j - 1} f(p_i)$ ，接下来考虑计算合数对答案的贡献。

我们枚举合数的最小质因子以及最小质因子的次幂，考虑从 $S (n, j + 1)$ 推到 $S (n, j)$ ，也就是在原来的基础上加上，
$∑i=jprimei×primei≤n∑k=1primeik+1≤n(S(nprimeik,j+1)f(primeik)+f(primeik+1))\sum_{i = j} ^{prime_i \times prime_i \leq n} \sum_{k = 1} ^{prime_i ^ {k + 1} \leq n} \left(S(\frac{n}{prime_i ^ k}, j + 1)f(prime_i ^ k) + f(prime_i ^ {k + 1})\right)\\$
所以有
$\sum_{i = 1} ^{j - 1} f(p_i) +\sum_{i = j} ^{prime_i \times prime_i \leq n} \sum_{k = 1} ^{prime_i ^ {k + 1} \leq n} \left(S(\frac{n}{prime_i ^ k}, j + 1)f(prime_i ^ k) + f(prime_i ^ {k + 1})\right)\\$
这里的 $g (n)$ 是最后筛得的，也就是里面存的是质数 $f$ 的前缀和。

这里解释一下那后面一长串的东西，为什么要加上 $f(prime_i ^{k + 1})$ ，而且上界为 $primeik+1≤nprime_i ^{k + 1} \leq n$ ，

由于我们枚举的是最小质因子及幂次，那么一定有 $n > prime_i ^ k$ ，所以有上界为 $primeik+1≤nprime_i ^{k + 1} \leq n$ ，

至于为什么要加上 $f(prime_i ^ {k + 1})$ ，由于我们定义的函数中 $f (1)$ 在执行 $min\_25$ 的过程中始终为 $0$ ，

当我们枚举的幂次为 $k + 1$ 时，我们会漏了 $f(prime_i ^{k + 1})$ 这一项，所以得在前一步给他加上去。
优化一点点？如何递推计算积性函数前缀和 $∑i=1nf(i)\sum\limits_{i = 1} ^{n} f(i)$ 。

不就是把上面的 $S (n, j)$ ，从 $S (n, j + 1)$ 一步一步往下推吗？

对于 $nprimejk≥primej,nprimejk<primej+12\frac{n}{prime_j ^ k} \geq prime_j, \frac{n}{prime_j ^ k} < prime_{j + 1} ^ 2$ ，这一步，中我们递归会返回的是 $\sum\limits_{i = 1} ^{j}f(p_i)$ ，

也就是大于 $prime_j$ 的质数的部分，按照 $S (n, j)$ 的定义，我们在递推的时候，得到的却是 $0$ 。

所以我们修改一下递推时候的 $S (n, j)$ 的定义：设 $S(n,j)=∑i=1n[i∈primeormin_prime_of_i≥pj]f(i)S(n, j) = \sum\limits_{i = 1} ^{n}[i \in prime\ or\ min\_prime\_of\_i \geq p_j]f(i)$ ，

我们还是同样的枚举最小质因子，及其次幂，考虑从 $S (n, j + 1)$ 推到 $S (n, j)$ ，也就是要加上
$∑i=jprimei×primei≤n∑k=1primeik+1≤n((S(nprimeik,j+1)−∑k=1jf(pi))f(primeik)+f(primeik+1))\sum_{i = j} ^{prime_i \times prime_i \leq n} \sum_{k = 1} ^{prime_i ^{k + 1} \leq n} \left(\left(S(\frac{n}{prime_i ^ k}, j + 1) - \sum_{k = 1} ^{j} f(p_i) \right)f(prime_i ^ k) + f(prime_i ^{k + 1})\right)\\$

P5325 【模板】Min_25筛

有积性函数 $f (x)$ ，且 $\in primes$ ，求 $∑i=1nf(i)\sum\limits_{i = 1} ^{n} f(i)$ 。

#include <bits/stdc++.h>using namespace std;namespace min_25 {const int N = 1e6 + 10, mod = 1e9 + 7, inv2 = 500000004, inv6 = 166666668;int prime[N], id1[N], id2[N], m, T, cnt;long long a[N], sum1[N], g1[N], sum2[N], g2[N], s[N], n;bool st[N];inline int ID(long long x) {return x <= T ? id1[x] : id2[n / x];}inline long long f(long long x) {x %= mod;return x * (x - 1) % mod;}inline long long calc1(long long x) {x %= mod;return (x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod - 1 + mod) % mod;}inline long long calc2(long long x) {x %= mod;return (x * (x + 1) % mod * inv2 % mod - 1 + mod) % mod;}void init() {T = sqrt(n + 0.5);for(int i = 2; i <= T; i++) {if(!st[i]) {prime[++cnt] = i;sum1[cnt] = (sum1[cnt - 1] + 1ll * i * i) % mod;sum2[cnt] = (sum2[cnt - 1] + i) % mod;}for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {break;}}}for (long long l = 1, r; l <= n; l = r + 1) {r = n / (n / l);a[++m] = n / l;a[m] <= T ? id1[a[m]] = m : id2[n / a[m]] = m;g1[m] = calc1(a[m]);g2[m] = calc2(a[m]);}for (int j = 1; j <= cnt && 1ll * prime[j] * prime[j] <= n; j++) {for (int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {g1[i] = ((g1[i] - 1ll * prime[j] * prime[j] % mod * (g1[ID(a[i] / prime[j])] - sum1[j - 1]) % mod) % mod + mod) % mod;g2[i] = ((g2[i] - 1ll * prime[j] * (g2[ID(a[i] / prime[j])] - sum2[j - 1]) % mod) % mod + mod) % mod;}}for (int i = 1; i <= m; i++) {s[i] = (g1[i] - g2[i] + mod) % mod;}for (int j = cnt; j >= 1; j--) {for (int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {for (long long cur = prime[j]; cur * prime[j] <= a[i]; cur *= prime[j]) {s[i] = (s[i] + f(cur % mod) * (s[ID(a[i] / cur)] - sum1[j] + sum2[j]) % mod + f(cur * prime[j])) % mod;s[i] = (s[i] + mod) % mod;}}}printf("%lld\n", s[ID(n)] + 1);}}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);scanf("%lld", &min_25::n);min_25::init();return 0;
}

#6235. 区间素数个数

直接是 $min\_25$ 的 $g$ 函数求解。

#include <bits/stdc++.h>using namespace std;typedef long long ll;const int N = 1e6 + 10;int prime[N], id1[N], id2[N], m, cnt, T;bool st[N];ll a[N], sum[N], g[N], n;inline int ID(ll x) {return x <= T ? id1[x] : id2[n / x];
}void init() {T = sqrt(n + 0.5);for (int i = 2; i <= T; i++) {if (!st[i]) {prime[++cnt] = i;sum[cnt] = sum[cnt - 1] + 1;}for (int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {st[i * prime[j]] = 1;if (i % prime[j] == 0) {break;}}}for (ll l = 1, r; l <= n; l = r + 1) {r = n / (n / l);a[++m] = n / l;if (a[m] <= T) {id1[a[m]] = m;}else {id2[n / a[m]] = m;}g[m] = a[m] - 1;}for (int j = 1; j <= cnt && 1ll * prime[j] * prime[j] <= n; j++) {for (int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {g[i] -= g[ID(a[i] / prime[j])] - sum[j - 1];}}
}ll solve(ll n) {return g[ID(n)];
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);scanf("%lld", &n);init();printf("%lld\n", solve(n));return 0;
}

#572. 「LibreOJ Round #11」Misaka Network 与求和

$f (n)$ 为次大质因子， $f (6) = 3, f (4) = 2, f (1) = 0, f (p) = 1$ 。
$∑i=1n∑j=1nf(gcd⁡(i,j))k∑d=1nf(d)k∑i=1n∑j=1n[gcd(i,j)=d]∑d=1nf(d)k∑i=1nd∑j=1nd[gcd(i,j)=1]∑d=1nf(d)k(2×∑i=1nd∑j=1i[gcd(i,j)=1]−1)∑d=1nf(d)k(2∑i=1ndϕ(i)−1)\sum_{i = 1} ^{n} \sum_{j = 1} ^{n} f(\gcd(i, j)) ^ k\\ \sum_{d = 1} ^{n} f(d) ^ k \sum_{i = 1} ^{n} \sum_{j = 1} ^{n}[gcd(i, j) = d]\\ \sum_{d = 1} ^{n} f(d) ^ k \sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{\frac{n}{d}}[gcd(i, j) = 1]\\ \sum_{d = 1} ^{n} f(d) ^ k \left( 2 \times\sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{i}[gcd(i, j) = 1] - 1 \right)\\ \sum_{d = 1} ^{n} f(d) ^ k \left(2 \sum_{i = 1} ^{\frac{n}{d}}\phi(i) - 1\right)\\$
设 $\sum\limits_{i = 1} ^{n} \phi(i), F(n) = f(n) ^ k$ ，而有：
$∑d=1nF(d)(2×S(nd)−1)\sum_{d = 1} ^{n} F(d) \left(2 \times S(\frac{n}{d}) - 1 \right)\\$
对于 $\times S(\frac{n}{d}) - 1$ ，可以通过杜教筛求得，考虑如何求前面的 $∑i=1nF(i)\sum\limits_{i = 1} ^{n} F(i)$ ，对前项考虑 $min\_25$ 求解。

我们定义 $S(n,j)=∑i=1n[i∈primeormin_primeofi≥primej]f(i)S(n, j) = \sum\limits_{i = 1} ^{n} [i \in prime\ or\ min\_prime\ of\ i \geq prime_j] f(i)$ ，也就是质数或者，最小质因子大于等于 $prime_j$ 的答案。

每一轮我们枚举最小质因子及其幂次，考虑这一部分合数的贡献，分两种：

$prime_j$ 是次小质因子。
$prime_j$ 不是次小质因子。

#include <bits/stdc++.h>using namespace std;typedef unsigned int uint;uint quick_pow(uint a, int n) {uint ans = 1;while (n) {if (n & 1) {ans *= a;}a *= a;n >>= 1;}return ans;
}int n, k;namespace min_25 {const int N = 1e6 + 10;int prime[N], a[N], id1[N], id2[N], m, cnt, T;uint g[N], sum[N], s[N], f[N];bool st[N];inline int ID(int x) {return x <= T ? id1[x] : id2[n / x];}void init() {T = sqrt(n + 0.5);for(int i = 2; i <= T; i++) {if(!st[i]) {prime[++cnt] = i;f[cnt] = quick_pow(i, k);sum[cnt] = sum[cnt - 1] + 1;}for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {break;}}}for(int l = 1, r; l <= n; l = r + 1) {r = n / (n / l);a[++m] = n / l;a[m] <= T ? id1[a[m]] = m : id2[n / a[m]] = m;g[m] = a[m] - 1;}for(int j = 1; j <= cnt; j++) {for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {g[i] -= g[ID(a[i] / prime[j])] - sum[j - 1];}}for (int i = 1; i <= m; i++) {s[i] = g[i];}for (int j = cnt; j >= 1; j--) {for (int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {for (int cur = prime[j]; 1ll * cur * prime[j] <= a[i]; cur *= prime[j]) {s[i] += (s[ID(a[i] / cur)] - g[ID(a[i] / cur)]) + (g[ID(a[i] / cur)] - sum[j - 1]) * f[j];/*合数答案贡献分为两类：一、当前枚举的质数不是次小的，二、当前枚举的质数是次小的，也就是形如prime[j] ^ k * b, b是一个质数且这个质数大于 prime[j], 最后再加上 f(prime[j] ^ {k + 1})。*/}}}}uint solve(int x) {return x <= 1 ? 0 : s[ID(x)];}
}namespace Djs {const int N = 2e6 + 10;int prime[N], cnt;uint phi[N];bool st[N];void init() {phi[1] = 1;for (int i = 2; i < N; i++) {if (!st[i]) {prime[++cnt] = i;phi[i] = i - 1;}for (int j = 1; j <= cnt && 1ll * i * prime[j] < N; j++) {st[i * prime[j]] = 1;if (i % prime[j] == 0) {phi[i * prime[j]] = phi[i] * prime[j];break;}phi[i * prime[j]] = phi[i] * (prime[j] - 1);}}for (int i = 1; i < N; i++) {phi[i] += phi[i - 1];}}unordered_map<int, uint> mp;uint Djs_ans(int n) {if (n < N) {return phi[n];}if (mp.count(n)) {return mp[n];}uint ans = 1ll * n * (n + 1) / 2;for (uint l = 2, r; l <= n; l = r + 1) {r = n / (n / l);ans -= (r - l + 1) * Djs_ans(n / l);}return mp[n] = ans;}
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);scanf("%d %d", &n, &k);min_25::init();Djs::init();uint ans = 0;for (uint l = 1, r; l <= n; l = r + 1) {r = n / (n / l);ans += (min_25::solve(r) - min_25::solve(l - 1)) * (2 * Djs::Djs_ans(n / l) - 1);}cout << ans << "\n";return 0;
}

1847 奇怪的数学题

给定 $n, k$ ，要我们求 $∑i=1n∑j=1nsgcd(i,j)k\sum\limits_{i = 1} ^{n} \sum\limits_{j = 1} ^{n}sgcd(i, j) ^ k$ ， $s g c d (i, j)$ 是 $g c d (i, j)$ 的次大公约数， $sgcd(i,j)=gcd⁡(i,j)min_primeofgcd⁡(i,j)sgcd(i, j) = \frac{\gcd(i, j)}{min\_prime\ of\ \gcd(i, j)}$ ，另 $s g c d (i, j) = f (g c d (i, j))$ 。
$∑i=1n∑j=1nf(gcd(i,j))k∑d=1nf(d)k∑i=1n∑j=1n[gcd(i,j)=d]∑d=1nf(d)k∑i=1nd∑j=1nd[gcd(i,j)=1]∑d=1nf(d)k(2×∑i=1nd∑j=1i[gcd(i,j)=1]−1)∑d=1nf(d)k(2×∑i=1ndϕ(i)−1)\sum_{i = 1} ^{n} \sum_{j = 1} ^{n} f(gcd(i, j)) ^ k\\ \sum_{d = 1} ^{n} f(d) ^ k \sum_{i = 1} ^{n} \sum_{j = 1} ^{n} [gcd(i, j) = d]\\ \sum_{d = 1} ^{n} f(d) ^ k \sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{\frac{n}{d}}[gcd(i,j) = 1]\\ \sum_{d = 1} ^{n} f(d) ^ k \left(2 \times \sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{i} [gcd(i, j) = 1] - 1\right)\\ \sum_{d = 1} ^{n} f(d) ^ k \left(2 \times \sum_{i = 1} ^{\frac{n}{d}} \phi(i) -1 \right)\\$
后项可以用杜教筛求解，考虑前项的 $∑d=1nf(d)k\sum\limits_{d = 1} ^{n} f(d) ^ k$ ， $min\_25$ 中有枚举最小质因子的一步，下面定义 $f(n) = n ^ k$ 。

还记得 $f(p_j) \left(g(\frac{n}{p_j}) - \sum\limits_{i = 1} ^{j - 1} p_j \right)$ ，的转移

后面减去的不就是最小质因子等于 $p_j$ 且为合数的函数值，那么 $g(npj−∑i=1j−1pj)g(\frac{n}{p_j} - \sum\limits_{i = 1} ^{j - 1}p_j)$ 就是除去 $f(p_j)$ 的，也就是我们需要的 $s g c d$ 了。

#include <bits/stdc++.h>using namespace std;typedef unsigned int uint;const int N = 1e6 + 10;uint quick_pow(uint a, int n) {uint ans = 1;while(n) {if(n & 1) ans = ans * a;a = a * a;n >>= 1;}return ans;
}namespace Sum {uint S[60][60], k;uint Sumk(int n) {uint ret = 0;for(int i = 1; i <= k; ++i) {uint s = 1;for(int j = 0; j <= i; ++j)if((n + 1 - j) % (i + 1)) s *= n + 1 - j;else s *= (n + 1 - j) / (i + 1);ret += S[k][i] *s;}return ret;}void init() {S[0][0] = 1;for(int i = 1; i <= k; i++)for(int j = 1; j <= i; j++)S[i][j] = S[i-1][j-1] + j * S[i-1][j];}     
}namespace Min_25 {uint prime[N], id1[N], id2[N], m, cnt, k, T;uint a[N], g1[N], sum1[N], g2[N], sum2[N], f[N], ans[N], n;bool st[N];int ID(int x) {return x <= T ? id1[x] : id2[n / x];}void init() {T = sqrt(n + 0.5);for(int i = 2; i <= T; i++) {if(!st[i]) {prime[++cnt] = i;sum1[cnt] = sum1[cnt - 1] + 1;f[cnt] = quick_pow(i, k);sum2[cnt] = sum2[cnt - 1] + f[cnt];}for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {break;}}}for(int l = 1, r; l <= n; l = r + 1) {r = n / (n / l);a[++m] = n / l;if(a[m] <= T) id1[a[m]] = m;else id2[n / a[m]] = m;g1[m] = a[m] - 1;g2[m] = Sum::Sumk(a[m]) - 1;}for(int j = 1; j <= cnt; j++) {for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {g1[i] -= (g1[ID(a[i] / prime[j])] - sum1[j - 1]);g2[i] -= f[j] * (g2[ID(a[i] / prime[j])] - sum2[j - 1]);ans[i] += g2[ID(a[i] / prime[j])] - sum2[j - 1];}}for(int i = 1; i <= m; i++) {ans[i] += g1[i];}}uint solve(int x) {if(x <= 1) return 0;return ans[ID(x)];}
}namespace Djs {uint prime[N], phi[N], cnt;bool st[N];void init() {phi[1] = 1;for(int i = 2; i < N; i++) {if(!st[i]) {prime[++cnt] = i;phi[i] = i - 1;}for(int j = 1; j <= cnt && 1ll * i * prime[j] < N; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {phi[i * prime[j]] = phi[i] * prime[j];break;}phi[i * prime[j]] = phi[i] * (prime[j] - 1);}}for(int i = 1; i < N; i++) {phi[i] += phi[i - 1];}}unordered_map<int, uint> ans_s;uint S(int n) {if(n < N) return phi[n];if(ans_s.count(n)) return ans_s[n];uint ans = 1ll * n * (n + 1) / 2;for(uint l = 2, r; l <= n; l = r + 1) {r = n / (n / l);ans -=(r - l + 1) * S(n / l);}return ans_s[n] = ans;}
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);int n, k;cin >> n >> k;Sum::k = k;Sum::init();Djs::init();Min_25::n = n, Min_25::k = k;Min_25::init();uint ans = 0;for(uint l = 1, r; l <= n; l = r + 1) {r = n / (n / l);ans += (Min_25::solve(r) - Min_25::solve(l - 1)) * (2 * Djs::S(n / l) - 1);}cout << ans << "\n";return 0;
}