Sàng Min-25 là một kỹ thuật mạnh mẽ được sử dụng để tính tổng tiền tố của các hàm nhân tính \(f(i)\) trong thời gian cận tuyến tính (sub-linear). Cụ thể, thuật toán này có thể giải quyết bài toán tìm \(\sum_{i=1}^n f(i)\) với độ phức tạp thời gian khoảng \(\mathcal{O}(\frac{n^{3/4}}{\log n})\) và không gian \(\mathcal{O}(\sqrt{n})\).
Để áp dụng sàng Min-25, hàm nhân tính \(f(i)\) cần thỏa mãn hai điều kiện cơ bản:
- Tại các số nguyên tố \(p\), \(f(p)\) có thể được biểu diễn dưới dạng một đa thức theo \(p\).
- Tại các lũy thừa của số nguyên tố \(p^c\), giá trị \(f(p^c)\) có thể được tính toán nhanh chóng (thường là \(\mathcal{O}(1)\)).
1. Nguyên lý hoạt động
Ý tưởng cốt lõi của thuật toán là chia các số trong khoảng \([1, n]\) thành hai nhóm: số nguyên tố và hợp số. Các hợp số sẽ được phân loại dựa trên nhân tử nguyên tố nhỏ nhất của chúng.
Gọi \(p_j\) là số nguyên tố thứ \(j\). Ta định nghĩa hàm \(S(n, k)\) là tổng của \(f(i)\) với mọi \(i \in [1, n]\) sao cho nhân tử nguyên tố nhỏ nhất của \(i\) lớn hơn \(p_k\).
Công thức truy hồi của \(S(n, k)\) như sau:
\[ S(n, k) = Q(n) - \sum_{i=1}^k f(p_i) + \sum_{i=k+1}^{p_i^2 \le n} \sum_{c=1}^{p_i^{c+1} \le n} \left( f(p_i^c) S(\lfloor n/p_i^c \rfloor, i) + f(p_i^{c+1}) \right) \]
Trong đó \(Q(n) = \sum_{p \le n} f(p)\) là tổng giá trị của hàm tại các số nguyên tố không vượt quá \(n\). Kết quả cuối cùng của bài toán sẽ là \(S(n, 0) + f(1)\).
Để tính \(Q(n)\), vì \(f(p)\) là một đa thức, ta có thể tính tổng từng đơn thức của nó. Giả sử ta cần tính tổng các số nguyên tố \(p^d\). Ta định nghĩa \(G(n, k)\) là tổng các \(i^d\) với \(i \in [1, n]\) sao cho \(i\) là số nguyên tố hoặc nhân tử nguyên tố nhỏ nhất của \(i\) lớn hơn \(p_k\).
Công thức chuyển trạng thái của \(G(n, k)\):
\[ G(n, k) = \begin{cases} G(n, k-1) & \text{nếu } p_k^2 > n \\ G(n, k-1) - p_k^d \left( G(\lfloor n/p_k \rfloor, k-1) - \sum_{i=1}^{k-1} p_i^d \right) & \text{nếu } p_k^2 \le n \end{cases} \]
Điều kiện biên là \(G(n, 0) = \sum_{i=2}^n i^d\), có thể tính bằng các công thức tổng lũy thừa thông thường.
2. Cài đặt thuật toán
Dưới đây là mã nguồn minh họa cách tính tổng hàm nhân tính bằng phương pháp đệ quy. Trong ví dụ này, ta xét một hàm nhân tính giả định nơi \(f(p) = p^2 - p\).
#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
typedef long long ll;
const int MOD = 1e9 + 7;
ll N;
int sqrtN;
vector<int> primes;
bool is_prime[200005];
ll sp1[200005], sp2[200005];
ll values[200005];
int id1[200005], id2[200005];
ll g1[200005], g2[200005];
int m = 0;
void sieve(int n) {
fill(is_prime + 2, is_prime + n + 1, true);
for (int i = 2; i <= n; i++) {
if (is_prime[i]) {
primes.push_back(i);
sp1[primes.size()] = (sp1[primes.size() - 1] + i) % MOD;
sp2[primes.size()] = (sp2[primes.size() - 1] + (ll)i * i) % MOD;
}
for (int p : primes) {
if (i * p > n) break;
is_prime[i * p] = false;
if (i % p == 0) break;
}
}
}
inline int get_id(ll x) {
return x <= sqrtN ? id1[x] : id2[N / x];
}
ll S(ll x, int y) {
if (x <= 1 || (y <= (int)primes.size() && primes[y - 1] > x)) return 0;
int k = get_id(x);
ll res = (g2[k] - g1[k] + MOD) % MOD;
res = (res - (sp2[y - 1] - sp1[y - 1] + MOD) % MOD + MOD) % MOD;
for (int i = y; i <= (int)primes.size() && (ll)primes[i - 1] * primes[i - 1] ≤ x; i++) {
ll pe = primes[i - 1];
for (int e = 1; pe * primes[i - 1] ≤ x; e++, pe *= primes[i - 1]) {
ll f_pe = (pe % MOD * (pe % MOD - 1 + MOD)) % MOD;
ll f_pe_next = ((pe * primes[i - 1]) % MOD * ((pe * primes[i - 1]) % MOD - 1 + MOD)) % MOD;
res = (res + f_pe * S(x / pe, i + 1) % MOD + f_pe_next) % MOD;
}
}
return res;
}
int main() {
cin >> N;
sqrtN = sqrt(N);
sieve(sqrtN);
for (ll i = 1, j; i <= N; i = j + 1) {
j = N / (N / i);
values[++m] = N / i;
ll val = values[m] % MOD;
g1[m] = val * (val + 1) % MOD * 500000004 % MOD - 1; // sum of i
g2[m] = val * (val + 1) % MOD * (2 * val + 1) % MOD * 166666668 % MOD - 1; // sum of i^2
if (values[m] <= sqrtN) id1[values[m]] = m;
else id2[N / values[m]] = m;
}
for (int j = 1; j <= (int)primes.size(); j++) {
for (int i = 1; i <= m && (ll)primes[j - 1] * primes[j - 1] <= values[i]; i++) {
int k = get_id(values[i] / primes[j - 1]);
g1[i] = (g1[i] - (ll)primes[j - 1] * (g1[k] - sp1[j - 1] + MOD) % MOD + MOD) % MOD;
g2[i] = (g2[i] - (ll)primes[j - 1] * primes[j - 1] % MOD * (g2[k] - sp2[j - 1] + MOD) % MOD + MOD) % MOD;
}
}
cout << (S(N, 1) + 1) % MOD << endl;
return 0;
}
3. Phân tích độ phức tạp
Độ phức tạp không gian của thuật toán là \(\mathcal{O}(\sqrt{n})\) do chúng ta chỉ cần lưu trữ các giá trị tại các điểm \(\lfloor n/i \rfloor\).
Về thời gian, quá trình tính \(G(n, k)\) và \(S(n, k)\) đều có độ phức tạp được chứng minh là xấp xỉ \(\mathcal{O}(n^{3/4}/\log n)\). Trong thực tế, với \(n = 10^{10}\), thuật toán thường chạy trong khoảng dưới 1 giây nếu được tối ưu tốt, nhanh hơn đáng kể so với các phương pháp như Sàng Du (Du's Sieve) trong nhiều trường hợp.
4. Ví dụ ứng dụng: Tổng hàm số lượng ước số
Giả sử cần tính \(\sum_{i=1}^n \sigma_0(i^k)\), trong đó \(\sigma_0(x)\) là số lượng ước của \(x\).
Vì \(\sigma_0(i^k)\) là hàm nhân tính, ta có:
- Với \(p\) là số nguyên tố, \(f(p) = \sigma_0(p^k) = k + 1\). Đây là một đa thức bậc 0.
- Với \(p^c\), \(f(p^c) = \sigma_0(p^{ck}) = ck + 1\).
Khi đó, bước tính \(G(n, k)\) đơn giản là đếm số lượng số nguyên tố (hàm \(f'(p) = 1\)). Sau đó áp dụng công thức \(S(n, k)\) để thu được kết quả.
// Đoạn mã tính S cho hàm sigma_0(i^k)
ll calculate_S(ll x, int y, int K) {
if (x <= 1 || primes[y - 1] > x) return 0;
int k = get_id(x);
ll res = (ll)(g[k] - (y - 1)) * (K + 1) % MOD;
for (int i = y; i <= (int)primes.size() && (ll)primes[i - 1] * primes[i - 1] <= x; i++) {
ll pe = primes[i - 1];
for (int e = 1; pe * primes[i - 1] ≤ x; e++, pe *= primes[i - 1]) {
res = (res + (ll)(e * K + 1) * calculate_S(x / pe, i + 1, K) + ( (e + 1) * K + 1 )) % MOD;
}
}
return res;
}
Phương pháp này cho phép xử lý các giá trị \(n\) lên tới \(10^{11}\) hoặc thậm chí \(10^{12}\) tùy thuộc vào giới hạn thời gian và tính chất cụ thể của hàm \(f\).