Đếm số nguyên tố đến N lớn (Lucy/Meissel)

Bài toán: cho $N$ rất lớn (tới ${10}^{11}$, ${10}^{12}$), tính $\pi (N)$ là số số nguyên tố không vượt quá $N$. Sàng Eratosthenes cần $O(N)$ bộ nhớ nên không dùng được khi $N>{10}^9$ (một mảng bool cỡ ${10}^{11}$ là 100 GB). Thuật toán Lucy_Hedgehog (đặt theo tên người đăng nó trên diễn đàn Project Euler) tính $\pi (N)$ trong $O(N^{3/4})$ thời gian và chỉ $O(\sqrt{N})$ bộ nhớ: với $N={10}^{11}$ chạy dưới $0,1$ giây, với $N={10}^{12}$ khoảng $0,4$ giây. Cùng khung đó tính được tổng các số nguyên tố, đếm số nguyên tố theo lớp đồng dư, và đếm số lượng số có ước nguyên tố nhỏ nhất bằng $p$ cho trước.

Ý tưởng / Trực giác

Nhận xét 1: chỉ có $O(\sqrt{N})$ giá trị $\lfloor N/i\rfloor$ khác nhau

Xét tập $V=\{\lfloor N/i\rfloor :1\le i\le N\}$. Với $i\le \sqrt{N}$ có không quá $\sqrt{N}$ giá trị; với $i>\sqrt{N}$ thì $\lfloor N/i\rfloor <\sqrt{N}$, lại không quá $\sqrt{N}$ giá trị. Vậy $|V|\le 2\sqrt{N}$.

Tập $V$ có một tính chất "đóng" rất quý: nếu $v\in V$ thì $\lfloor v/p\rfloor \in V$ với mọi $p$, vì

$$\lfloor \frac{\lfloor N/i\rfloor }{p}\rfloor =\lfloor \frac{N}{ip}\rfloor .$$

Nghĩa là nếu ta tính một đại lượng $S(v)$ cho mọi $v\in V$, mà công thức truy hồi của $S(v)$ chỉ cần $S(\lfloor v/p\rfloor )$ và $S$ tại giá trị nhỏ, thì mọi thứ cần dùng đều nằm trong $V$. Ta sẽ tính $\pi (v)$ cho cả $2\sqrt{N}$ giá trị này cùng lúc, và đáp án là $\pi (N)$ (ứng với $i=1$).

Nhận xét 2: sàng Eratosthenes "theo số lượng"

Nhớ lại sàng Eratosthenes: duyệt số nguyên tố $p$ tăng dần, gạch bỏ các bội $p^2,p^2+p,\dots$. Sau khi xử lý xong mọi $p\le \sqrt{N}$, số còn lại trong $[2,N]$ đúng là các số nguyên tố.

Lucy không gạch từng số, mà chỉ theo dõi số lượng còn lại. Định nghĩa

$$S(v,p)=\#\{2\le x\le v:x\text{ nguyên tố, hoặc ước nguyên tố nhỏ nhất của }x\text{ lớn hơn }p\}.$$

Tức là $S(v,p)$ là số phần tử còn sống trong đoạn $[2,v]$ sau khi sàng bằng mọi số nguyên tố $\le p$. Ta có:

• Khởi tạo $S(v,1)=v-1$ (chưa gạch gì, còn đủ $2,3,\dots ,v$).
• Khi $p^2>v$: sàng thêm $p$ không gạch gì trong $[2,v]$, nên $S(v,p)=S(v,p-1)$. Đặc biệt $S(v,\lfloor \sqrt{v}\rfloor )=\pi (v)$.
• Nếu $p$ không phải số nguyên tố: không có gì thay đổi, $S(v,p)=S(v,p-1)$.

Công thức truy hồi và vì sao đúng

Với $p$ nguyên tố và $p^2\le v$:

$$\boxed{S(v,p)=S(v,p-1)-[S(\lfloor v/p\rfloor ,p-1)-S(p-1,p-1)]}$$

Lượng bị trừ phải bằng số phần tử bị gạch ở bước $p$, tức số hợp số $x\le v$ có ước nguyên tố nhỏ nhất đúng bằng $p$. Mỗi số như vậy viết duy nhất được thành $x=p·m$ với:

• $p\le m\le \lfloor v/p\rfloor$;
• $m$ không có ước nguyên tố nào nhỏ hơn $p$ (nếu không thì ước nhỏ nhất của $x$ nhỏ hơn $p$).

Đếm các $m$ như vậy:

• $S(\lfloor v/p\rfloor ,p-1)$ đếm các $m\in [2,\lfloor v/p\rfloor ]$ là số nguyên tố hoặc không có ước nguyên tố $<p$. Tập này gồm hai phần: các số nguyên tố nhỏ hơn $p$ (bị đếm vì là số nguyên tố, nhưng không hợp lệ) và mọi $m\ge p$ không có ước nguyên tố $<p$ (hợp lệ).
• Số nguyên tố nhỏ hơn $p$ có đúng $\pi (p-1)$ số, và $\pi (p-1)=S(p-1,p-1)$: mọi số nguyên tố $\le \sqrt{p-1}$ đều $\le p-1$, nên sau khi sàng bằng mọi số nguyên tố $\le p-1$ thì đoạn $[2,p-1]$ đã được sàng xong, chỉ còn các số nguyên tố.

Trừ đi là ra đúng số $m$ hợp lệ, nên công thức đúng. Chú ý $\lfloor v/p\rfloor \ge p>p-1$ nên hiệu luôn không âm.

Cài đặt thực tế chỉ giữ một mảng $S[v]$ cho $v\in V$ và cập nhật tại chỗ theo từng $p$. Để $S[\lfloor v/p\rfloor ]$ vẫn là giá trị của bước trước ($p-1$) khi dùng, ta duyệt $v$ giảm dần: $\lfloor v/p\rfloor <v$ nên nó chưa bị sửa ở bước $p$. Còn $S[p-1]$ thì không bị sửa ở bước $p$ vì $p-1<p^2$.

Nhận biết $p$ nguyên tố cũng không cần sàng riêng: $p$ nguyên tố khi và chỉ khi $S[p]\ne S[p-1]$ (lúc xét tới $p$, hai giá trị này đã là $\pi (p)$ và $\pi (p-1)$).

Ví dụ chạy tay

$N=30$, $\lfloor \sqrt{30}\rfloor =5$. Tập $V$ gồm các giá trị $\lfloor 30/i\rfloor$ với $i=1.\mathrm{.5}$ (phần "lớn") và $1.\mathrm{.5}$ (phần "nhỏ", các $v\le \lfloor \sqrt{30}\rfloor =5$):

V (giảm dần):  30  15  10   7   6 |  5   4   3   2   1
               --- large (30/i) -- | ------ small ------

Khởi tạo $S[v]=v-1$:

v      : 30  15  10   7   6   5   4   3   2   1
S[v]   : 29  14   9   6   5   4   3   2   1   0

Bước $p=2$ ($S[2]=1\ne S[1]=0$ nên 2 nguyên tố; below $=S[1]=0$). Cập nhật mọi $v\ge 4$, giảm dần:

v=30: S[30] -= S[15] - 0 = 14   -> 15    (gạch 4,6,...,30: 14 số chẵn)
v=15: S[15] -= S[7]  - 0 = 6    -> 8
v=10: S[10] -= S[5]  - 0 = 4    -> 5
v=7 : S[7]  -= S[3]  - 0 = 2    -> 4
v=6 : S[6]  -= S[3]  - 0 = 2    -> 3
v=5 : S[5]  -= S[2]  - 0 = 1    -> 3
v=4 : S[4]  -= S[2]  - 0 = 1    -> 2
         ^ S[15], S[7], ... vẫn là giá trị cũ vì duyệt v giảm dần

v      : 30  15  10   7   6   5   4   3   2   1
S[v]   : 15   8   5   4   3   3   2   2   1   0

Bước $p=3$ ($S[3]=2\ne S[2]=1$; below $=S[2]=1$). Cập nhật $v\ge 9$:

v=30: S[30] -= S[10] - 1 = 4    -> 11    (gạch 9, 15, 21, 27)
v=15: S[15] -= S[5]  - 1 = 2    -> 6     (gạch 9, 15)
v=10: S[10] -= S[3]  - 1 = 1    -> 4     (gạch 9)

v      : 30  15  10   7   6   5   4   3   2   1
S[v]   : 11   6   4   4   3   3   2   2   1   0

Bước $p=4$: $S[4]=2=S[3]$ nên 4 là hợp số, bỏ qua.

Bước $p=5$ ($S[5]=3\ne S[4]=2$; below $=S[4]=2$). Chỉ $v=30\ge 25$:

v=30: S[30] -= S[6] - 2 = 1     -> 10    (gạch 25; còn 35 > 30)

v      : 30  15  10   7   6   5   4   3   2   1
S[v]   : 10   6   4   4   3   3   2   2   1   0
         ^^

Kết quả $\pi (30)=10$ (các số $2,3,5,7,11,13,17,19,23,29$). Tiện thể ta có luôn $\pi (15)=6$, $\pi (10)=4$, $\pi (7)=4$: cả bảng là $\pi (v)$ cho mọi $v\in V$.

Cài đặt

Lưu $S$ bằng hai mảng small / large

Không dùng map (chậm $O(\log )$ mỗi lần truy cập). Với $r=\lfloor \sqrt{N}\rfloor$:

• small[v] lưu $S(v)$ với $v\le r$;
• large[i] lưu $S(\lfloor N/i\rfloor )$ với $i\le r$.

Khi cần $S(\lfloor v/p\rfloor )$ với $v=\lfloor N/i\rfloor$: vì $\lfloor v/p\rfloor =\lfloor N/(ip)\rfloor$, nếu $ip\le r$ thì đọc large[i*p], ngược lại $\lfloor N/(ip)\rfloor <\sqrt{N}$ nên đọc small[N/(i*p)]. Mọi truy cập đều là đọc mảng trực tiếp theo chỉ số; phép chia chỉ còn ở N / (i*p) (khi $ip>r$), ở v / p trong vòng small, và một lần n / p2 cho mỗi số nguyên tố $p$.

Thứ tự "$v$ giảm dần" tương ứng với: duyệt large[i] theo $i$ tăng dần, rồi small[v] theo $v$ giảm dần.

C++: đếm số nguyên tố

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

// pi(n) by Lucy_Hedgehog's method: O(n^{3/4}) time, O(sqrt n) memory.
ll countPrimes(ll n) {
    if (n < 2) return 0;
    ll r = sqrtl((long double)n);
    while (r * r > n) r--;                 // fix floating-point error
    while ((r + 1) * (r + 1) <= n) r++;
    // small[v] = S(v) for v <= r;  large[i] = S(n / i) for i <= r
    vector<ll> small(r + 1), large(r + 1);
    for (ll v = 1; v <= r; v++) small[v] = v - 1;
    for (ll i = 1; i <= r; i++) large[i] = n / i - 1;

    for (ll p = 2; p <= r; p++) {
        if (small[p] == small[p - 1]) continue;   // p is composite
        ll below = small[p - 1];                  // pi(p - 1)
        ll p2 = p * p;
        // large values: v = n / i >= p^2  <=>  i <= n / p^2
        ll iEnd = min(r, n / p2);
        for (ll i = 1; i <= iEnd; i++) {
            ll ip = i * p;                        // S(v / p) = S(n / (i * p))
            ll sub = (ip <= r) ? large[ip] : small[n / ip];
            large[i] -= sub - below;
        }
        // small values, decreasing, only v >= p^2
        for (ll v = r; v >= p2; v--)
            small[v] -= small[v / p] - below;
    }
    return large[1];                              // S(n / 1) = pi(n)
}

int main() {
    ll n;
    cin >> n;
    cout << countPrimes(n) << "\n";
}

Kiểm chứng: $\pi ({10}^{10})=455052511$, $\pi ({10}^{11})=4118054813$, $\pi ({10}^{12})=37607912018$.

C++: tổng các số nguyên tố (theo modulo)

Thay "đếm" bằng "cộng": khởi tạo $G(v)=2+3+\cdots +v$ và mỗi bước trừ $p·[G(\lfloor v/p\rfloor )-G(p-1)]$. Hệ số $p$ xuất hiện vì mỗi số bị gạch có dạng $x=p·m$ nên đóng góp $p·m$. Công thức chỉ đúng khi hàm được cộng (ở đây $f(x)=x$; tương tự với $f(x)=x^k$ hoặc $f(x)=1$) là nhân tính hoàn toàn: $f(pm)=f(p)f(m)$ với mọi $p,m$.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MOD = 1e9 + 7;

// Sum of primes <= n, modulo MOD.
ll sumPrimes(ll n) {
    if (n < 2) return 0;
    ll r = sqrtl((long double)n);
    while (r * r > n) r--;
    while ((r + 1) * (r + 1) <= n) r++;
    auto prefix = [](ll v) {             // (2 + 3 + ... + v) mod MOD
        ll a = v, b = v + 1;
        if (a % 2 == 0) a /= 2; else b /= 2;
        return ((a % MOD) * (b % MOD) % MOD - 1 + MOD) % MOD;
    };
    vector<ll> small(r + 1), large(r + 1), cnt(r + 1);
    for (ll v = 1; v <= r; v++) small[v] = prefix(v), cnt[v] = v - 1;
    for (ll i = 1; i <= r; i++) large[i] = prefix(n / i);

    for (ll p = 2; p <= r; p++) {
        if (cnt[p] == cnt[p - 1]) continue;       // exact counts detect primes
        ll below = small[p - 1], p2 = p * p, pm = p % MOD;
        ll iEnd = min(r, n / p2);
        for (ll i = 1; i <= iEnd; i++) {
            ll ip = i * p;
            ll sub = (ip <= r) ? large[ip] : small[n / ip];
            large[i] = ((large[i] - pm * ((sub - below + MOD) % MOD)) % MOD + MOD) % MOD;
        }
        for (ll v = r; v >= p2; v--) {
            small[v] = ((small[v] - pm * ((small[v / p] - below + MOD) % MOD)) % MOD + MOD) % MOD;
            cnt[v] -= cnt[v / p] - cnt[p - 1];    // keep exact pi(v) for small v
        }
    }
    return large[1];
}

int main() {
    ll n;
    cin >> n;
    cout << sumPrimes(n) << "\n";
}

Mảng cnt giữ $\pi (v)$ chính xác cho $v\le r$ chỉ để nhận biết số nguyên tố: sau khi lấy modulo, hai giá trị tổng có thể trùng nhau một cách tình cờ nên không dùng small[p] == small[p-1] được nữa. Kiểm chứng: tổng các số nguyên tố $\le {10}^{11}$ là $201467077743744681014$, chia ${10}^9+7$ dư $475146683$.

Python

Python đủ nhanh tới khoảng $N={10}^{10}$ (khoảng $1,5$ giây). Dùng dict cho gọn; muốn nhanh hơn thì đổi sang hai list small/large như bản C++.

from math import isqrt

def count_primes(n):
    """pi(n) by Lucy_Hedgehog's method."""
    if n < 2:
        return 0
    r = isqrt(n)
    vals = [n // i for i in range(1, r + 1)]        # large values, decreasing
    vals += list(range(vals[-1] - 1, 0, -1))        # small values, decreasing
    S = {v: v - 1 for v in vals}
    for p in range(2, r + 1):
        if S[p] == S[p - 1]:
            continue                                # p is composite
        below = S[p - 1]
        p2 = p * p
        for v in vals:                              # decreasing order
            if v < p2:
                break
            S[v] -= S[v // p] - below
    return S[n]

print(count_primes(int(input())))

Độ phức tạp

Bộ nhớ: $O(\sqrt{N})$, hai mảng độ dài $\lfloor \sqrt{N}\rfloor +1$. Với $N={10}^{12}$ chỉ khoảng $2·{10}^6$ số long long, tức 16 MB.

Thời gian: Ở bước $p$, ta cập nhật mọi $v\in V$ có $v\ge p^2$. Đổi thứ tự tính tổng: mỗi $v$ được cập nhật một lần cho mỗi số nguyên tố $p\le \sqrt{v}$, nên tổng công việc là

$$\sum _{v\in V}\pi (\sqrt{v})\le \sum _{v\in V}\sqrt{v}.$$

• Phần lớn, $v=N/i$ với $i\le \sqrt{N}$: ${\displaystyle {\sum }_{i\le \sqrt{N}}\sqrt{N/i}=\sqrt{N}{\sum }_{i\le \sqrt{N}}{i}^{-1/2}\approx \sqrt{N}·2N^{1/4}=2N^{3/4}}$.
• Phần nhỏ, $v\le \sqrt{N}$: ${\displaystyle {\sum }_{v\le \sqrt{N}}\sqrt{v}\approx \frac{2}{3}{N}^{3/4}}$.

Vậy thời gian $O(N^{3/4})$. Chính xác hơn, vì chỉ số nguyên tố mới gây cập nhật và $\pi (x)\approx x/\ln x$, tổng thực tế là $O(N^{3/4}/\log N)$. Đo thực tế số lần cập nhật: $3·{10}^6$ với $N={10}^9$, $8·{10}^7$ với $N={10}^{11}$, $4·{10}^8$ với $N={10}^{12}$; mỗi lần chỉ gồm vài lần đọc mảng, một phép trừ và nhiều nhất một phép chia 64 bit (n / ip ở vòng large khi $ip>r$, hoặc v / p ở vòng small), nên vẫn chạy dưới $0,5$ giây.

Vòng for p duyệt tới $\sqrt{N}$ cả số hợp, nhưng mỗi số hợp chỉ tốn $O(1)$ nên không ảnh hưởng.

⚠️ Lỗi thường gặp

• Tính $\lfloor \sqrt{N}\rfloor$ sai do số thực. double chỉ có 53 bit, nên (ll)sqrt(N) có thể lệch 1 khi $N$ lớn và sát một số chính phương (vd (ll)sqrt(1e18 - 1) ra ${10}^9$ thay vì $999999999$). Triệu chứng: thiếu hoặc thừa một giá trị trong $V$, đọc small[] ngoài mảng. Cách tránh: sau sqrtl luôn chỉnh bằng hai vòng while, hoặc dùng math.isqrt trong Python.
• Duyệt $v$ tăng dần. Khi đó $S[\lfloor v/p\rfloor ]$ đã bị sửa ở bước $p$ trước khi dùng, số bị gạch bị trừ thiếu. Triệu chứng: $\pi (30)$ ra $14$ thay vì $10$. Cách tránh: large[i] theo $i$ tăng (tức $v$ giảm), rồi small[v] theo $v$ giảm.
• Vòng for p dừng ở p < r thay vì p <= r. Khi $N=r^2$ (vd $N=25$), số $p=r$ vẫn phải gạch $p^2=N$. Triệu chứng: $\pi (25)$ ra $10$ thay vì $9$. Cách tránh: điều kiện p <= r.
• Tràn số khi tính i * p hoặc p * p. Với $N={10}^{12}$ thì $p\le {10}^6$, $p^2\le {10}^{12}$: vượt int. Cách tránh: mọi biến là long long. Với bản tính tổng, $p·G$ phải lấy modulo trước khi nhân, và tổng các số nguyên tố $\le {10}^{11}$ đã vượt $2^{64}$, nên hoặc lấy modulo hoặc dùng __int128.
• Nhận biết số nguyên tố bằng small[p] != small[p-1] trong bản có modulo. Hai tổng theo modulo có thể trùng tình cờ. Cách tránh: giữ song song một mảng đếm chính xác, hoặc sàng Eratosthenes riêng tới $\sqrt{N}$.
• Áp dụng cho hàm không nhân tính hoàn toàn. Công thức "trừ $f(p)·[\dots ]$" chỉ đúng nếu $f(pm)=f(p)f(m)$ với mọi $m$, như $f(x)=1$, $x$, $x^2$, hay ký tự ${\chi }_4$ trong phần mở rộng dưới. Với $\phi$, $\mu$ hay "số ước" thì không dùng trực tiếp được (xem phần Mở rộng).
• Dùng map/unordered_map trong C++. Chạy đúng nhưng chậm hơn 5 đến 20 lần, có thể TLE với $N={10}^{12}$. Cách tránh: hai mảng small/large.

Biến thể / Mở rộng

Đếm số theo ước nguyên tố nhỏ nhất

Lượng bị trừ ở bước $p$ khi $v=N$, tức $S(\lfloor N/p\rfloor ,p-1)-\pi (p-1)$, chính là số hợp số $\le N$ có ước nguyên tố nhỏ nhất bằng $p$. Cộng thêm $1$ (chính số $p$) ta được

$$\#\{2\le x\le N:spf(x)=p\}=S(\lfloor N/p\rfloor ,p-1)-\pi (p-1)+1,p\le \sqrt{N}.$$

Với số nguyên tố $p>\sqrt{N}$ thì đáp án là đúng $1$ (chỉ có $x=p$), vì mọi hợp số có ước nguyên tố nhỏ nhất $p$ đều $\ge p^2>N$: viết $x=kp$ với $k\ge 2$, nếu $k<p$ thì $k$ có ước nguyên tố nhỏ hơn $p$ (mâu thuẫn), còn nếu $k\ge p$ thì $x\ge p^2>N$. Ghi lại lượng trừ tại large[1] ở mỗi bước là có ngay bảng này, không tốn thêm độ phức tạp. Khi $N$ nhỏ (cỡ ${10}^7$), bảng tương tự lấy được trực tiếp từ mảng lp[] của sàng tuyến tính.

Đếm số nguyên tố theo lớp đồng dư mod 4

Hàm ${\chi }_4(x)=0,1,0,-1$ khi $x\equiv 0,1,2,3\left(\mathrm{mod}4\right)$ là nhân tính hoàn toàn. Cách đơn giản hơn là giữ hai mảng: $A(v)$ đếm số lẻ còn sống $\equiv 1\left(\mathrm{mod}4\right)$, $B(v)$ đếm số $\equiv 3\left(\mathrm{mod}4\right)$, bỏ qua $p=2$ (không đưa số chẵn vào ngay từ đầu):

• khởi tạo $A(v)=\lfloor (v-1)/4\rfloor$, $B(v)=\lfloor (v+1)/4\rfloor$;
• với $p\equiv 1\left(\mathrm{mod}4\right)$: $x=pm$ cùng lớp với $m$, nên $A$ trừ phần của $A$, $B$ trừ phần của $B$;
• với $p\equiv 3\left(\mathrm{mod}4\right)$: $x=pm$ đổi lớp so với $m$, nên $A$ trừ phần của $B$ và ngược lại.

Cùng ý tưởng mở rộng cho mọi modulo nhỏ $k$: giữ $\phi (k)$ mảng, mỗi số nguyên tố $p$ hoán vị các lớp theo phép nhân với $pmodk$.

Hàm nhân tính tổng quát: Min_25 và sàng Du

Muốn tính ${\sum }_{x\le N}f(x)$ cho $f$ nhân tính (không hoàn toàn) như $\phi$, $\mu$, số ước: bước 1 dùng đúng Lucy để có ${\sum }_{p\le v}f(p)$ cho mọi $v\in V$ (vì $f(p)$ thường là đa thức theo $p$), bước 2 cộng thêm các hợp số theo ước nguyên tố nhỏ nhất. Đó là sàng Min_25, chạy khoảng $O(N^{3/4}/\log N)$. Riêng tổng tiền tố của $\mu$ và $\phi$ còn có sàng Du (Du's sieve) dựa trên tích Dirichlet, $O(N^{2/3})$. Cả hai đều dựa trên cùng nhận xét "chỉ có $O(\sqrt{N})$ giá trị $\lfloor N/i\rfloor$".

Giới thiệu Meissel-Lehmer

Legendre định nghĩa $\varphi (x,a)$ là số số nguyên $1\le y\le x$ không chia hết cho $a$ số nguyên tố đầu tiên $p_1,\dots ,p_a$. Theo bao hàm, loại trừ, $\varphi$ có truy hồi

$$\varphi (x,a)=\varphi (x,a-1)-\varphi (\lfloor x/p_a\rfloor ,a-1),\varphi (x,0)=x.$$

Liên hệ chính xác với Lucy: với $p_a\le v$,

$$S(v,p_a)=\varphi (v,a)+a-1,$$

vì $\varphi (v,a)$ đếm số $1$ cùng mọi $x\ge 2$ có ước nguyên tố nhỏ nhất lớn hơn $p_a$, còn $S(v,p_a)$ bỏ số $1$ nhưng giữ thêm $a$ số nguyên tố $p_1,\dots ,p_a$. Meissel chọn $a=\pi (\sqrt[3]{x})$, $b=\pi (\sqrt{x})$. Các số $\le x$ không có ước nguyên tố nào trong $p_1..p_a$ gồm: số $1$, các số nguyên tố $p_{a+1},\dots ,p_{\pi (x)}$, và các tích đúng hai số nguyên tố $p_j{p}_k$ với $a<j\le k$ (ba thừa số thì vượt $x$). Với mỗi thừa số nhỏ $p_j$ (cần $p_j^2\le x$, tức $j\le b$), thừa số lớn $p_k$ chạy qua các số nguyên tố trong $[p_j,x/p_j]$, có $\pi (\lfloor x/p_j\rfloor )-(j-1)$ giá trị. Suy ra

$$\pi (x)=\varphi (x,a)+a-1-P_2(x,a),P_2(x,a)=\sum _{j=a+1}^b\left(\pi \right(\lfloor x/p_j\rfloor )-(j-1)).$$

$P_2$ cần $\pi$ tại các điểm $\le x^{2/3}$, lấy bằng một sàng đoạn; $\varphi (x,a)$ tính bằng đệ quy có cắt tỉa và bảng tra cho $a$ nhỏ. Lehmer, rồi Lagarias, Miller, Odlyzko cải tiến để đạt khoảng $O(x^{2/3})$ thời gian. Với $N\le {10}^{12}$ trong thi đấu, Lucy ngắn gọn hơn nhiều và đủ nhanh; Meissel-Lehmer chỉ đáng cài khi $N$ tới ${10}^{13}$ trở lên.

Khi nào dùng kỹ thuật nào

Bài tập luyện

Luyện kỹ thuật

• Số nguyên tố hay số đối xứng? (primepal) (Nâng cao) khởi động: cần $\pi (n)$ cho mọi $n$ tới khoảng $2·{10}^6$, đủ nhỏ để dùng sàng thường với mảng tiền tố; bài giúp phân biệt khi nào sàng là đủ và khi nào mới cần Lucy.
• Dãy số nguyên tố cùng nhau (kthcopr) (Kỳ cựu) đếm số $\le n$ không chia hết cho các ước nguyên tố của $p$ bằng bao hàm, loại trừ: đúng là hàm $\varphi (x,a)$ của Legendre (nền của Meissel-Lehmer) với tập số nguyên tố tùy ý, kết hợp chặt nhị phân.
• Số may mắn kép (dblhappy) (Kỳ cựu) định lý Fermat quy bài về đếm số nguyên tố $\equiv 1\left(\mathrm{mod}4\right)$ trong $[l,r]$, cộng thêm $1$ nếu số $2$ nằm trong đoạn (vì $2=1^2+1^2$); lời giải gốc sàng bitset tới $3·{10}^8$, nhưng Lucy với hai mảng theo lớp đồng dư mod 4 (phần Mở rộng) chạy tức thì và gần như không tốn bộ nhớ.

Trong ba bài trên, chỉ dblhappy áp dụng trực tiếp kỹ thuật Lucy; primepal (sàng thường) và kthcopr (bao hàm, loại trừ trên vài số nguyên tố) là bài liên quan để đối chiếu. Pool CTOJ hiện chưa có bài "thuần" đếm $\pi (N)$ hoặc tổng số nguyên tố với $N\ge {10}^{10}$; khi có sẽ bổ sung vào đây.

Đề VOI

• Xây dựng ma trận (voi22_matrix) (Chuyên gia) VOI 2022: điền tay vài hàng để thấy giá trị ô chỉ phụ thuộc vào thứ tự của ước nguyên tố nhỏ nhất của chỉ số hàng và cột; với $m,n\le {10}^9$, chạy Lucy riêng cho $N=m$ (chỉ số hàng) và $N=n$ (chỉ số cột): số chỉ số có cùng ước nguyên tố nhỏ nhất $p\le \sqrt{N}$ chính là lượng Lucy trừ đi ở bước $p$ (cộng 1), còn mỗi số nguyên tố $p>\sqrt{N}$ tạo một lớp có đúng một phần tử.