Mở đầu

Nên đọc bài viết blog của tôi trước đó — Min_25 sieve

Thuật toán Meissel-Lehmer là gì

Thuật toán Meissel-Lehmer là phương pháp tính nhanh số lượng số nguyên tố trong dãy số từ 1 đến n (có thể mở rộng cho tổng lũy thừa số nguyên tố).

[Bài toán mẫu] Thuật toán Meissel-Lehmer

Cho số nguyên n, hãy tính giá trị π(n). π(n) là số lượng số nguyên tố trong khoảng 1 đến n. Với 100% dữ liệu, 1 ≤ n ≤ 10^13.

Giải pháp

Kiến thức nền tảng: Sàng Eratosthenes. Sàng Eratosthenes là phương pháp sàng đơn giản nhất.

int Eratosthenes(int n) {
  int p = 0;
  for (int i = 0; i <= n; ++i) is_prime[i] = 1;
  is_prime[0] = is_prime[1] = 0;
  for (int i = 2; i <= n; ++i) {
    if (is_prime[i]) {
      prime[p++] = i;
      if ((long long)i * i <= n) for (int j = i * i; j <= n; j += i)
          is_prime[j] = 0; //(*)
    }
  }
  return p;
}

Chúng ta chỉ cần biết có bao nhiêu giá trị is_prime chuyển từ true sang false tại dòng đánh dấu. Nếu có thể theo dõi giá trị này, ta có thể xác định số lượng số nguyên tố từ 1 đến n. Có thông tin cho rằng phương pháp Lucy DP trên Project Euler có thể giải quyết vấn đề này, nhưng đó không phải là điều chúng ta cần.

Định nghĩa: Gọi is_prime là a, S(v,p) là số lượng a_j (2 ≤ j ≤ v) còn true sau khi xử lý đến i=p.

Tiếp tục xử lý, dù độ phức tạp thời gian và không gian không đạt tuyến tính nhưng tốc độ vẫn rất nhanh.

Sau đó là phương pháp Min_25 sieve Giải thích một số biến trong mã Min_25:

smalls[i] = S(i, p) roughs[i]: số nguyên thứ i chưa bị loại bỏ. Kích thước được duy trì bởi biến s. Giữ lại số 1, các số nguyên tố được sử dụng để sàng sẽ bị xóa. larges[i] = S(⌊n/roughs[i]⌋, p) pc = π(p-1) skip[i]: nếu i bị loại bỏ thì giá trị này là true. Không cập nhật số chẵn.

using i64 = long long;
int isqrt(i64 n) {
  return sqrtl(n);
}
__attribute__((target("avx"), optimize("O3", "unroll-loops")))
i64 prime_pi(const i64 N) {
  if (N <= 1) return 0;
  if (N == 2) return 1;
  const int v = isqrt(N);
  int s = (v + 1) / 2;
  vector<int> smalls(s); for (int i = 1; i < s; ++i) smalls[i] = i;
  vector<int> roughs(s); for (int i = 0; i < s; ++i) roughs[i] = 2 * i + 1;
  vector<i64> larges(s); for (int i = 0; i < s; ++i) larges[i] = (N / (2 * i + 1) - 1) / 2;
  vector<bool> skip(v + 1);
  const auto divide = [] (i64 n, i64 d) -> int { return double(n) / d; };
  const auto half = [] (int n) -> int { return (n - 1) >> 1; };
  int pc = 0;
  for (int p = 3; p <= v; p += 2) if (!skip[p]) {
    int q = p * p;
    if (i64(q) * q > N) break;
    skip[p] = true;
    for (int i = q; i <= v; i += 2 * p) skip[i] = true;
    int ns = 0;
    for (int k = 0; k < s; ++k) {
      int i = roughs[k];
      if (skip[i]) continue;
      i64 d = i64(i) * p;
      larges[ns] = larges[k] - (d <= v ? larges[smalls[d >> 1] - pc] : smalls[half(divide(N, d))]) + pc;
      roughs[ns++] = i;
    }
    s = ns;
    for (int i = half(v), j = ((v / p) - 1) | 1; j >= p; j -= 2) {
      int c = smalls[j >> 1] - pc;
      for (int e = (j * p) >> 1; i >= e; --i) smalls[i] -= c;
    }
    ++pc;
  }
  larges[0] += i64(s + 2 * (pc - 1)) * (s - 1) / 2;
  for (int k = 1; k < s; ++k) larges[0] -= larges[k];
  for (int l = 1; l < s; ++l) {
    int q = roughs[l];
    i64 M = N / q;
    int e = smalls[half(M / q)] - pc;
    if (e < l + 1) break;
    i64 t = 0;
    for (int k = l + 1; k <= e; ++k) t += smalls[half(divide(M, roughs[k]))];
    larges[0] += t - i64(e - l) * (pc + l - 1);
  }
  return larges[0] + 1;
}

Quan sát kỹ đoạn mã trên có thể hiểu được thuật toán, nếu chưa rõ có thể xem các giải thích của tôi:

Cập nhật roughs: roughs[..s] chứa các số chưa bị sàng. ns là giá trị s trong vòng lặp tiếp theo. ns = 0, khi roughs[ns] có giá trị thì ns tăng lên. larges[smalls[d]-pc] có nghĩa là: smalls[d] - pc là vị trí của d trong số các số chưa bị sàng. Cập nhật smalls: giữ nguyên j = i/p. Khi i ≥ pj, smalls[i] được cập nhật bằng cách trừ đi smalls[j] - pc, tức là số lượng số nguyên tố nhỏ hơn j.

Mã nguồn đầy đủ

#include <stdio.h>
#include <cmath>
#include <algorithm>
#include <vector>
#define qaq inline
typedef long long ll;
qaq int isqrt(ll n){
    return std::sqrt(n);
}
qaq ll half(ll n){
    return (n-1)>>1;
}
qaq ll divide(ll n,ll base){
    return double(n)/base;
}
ll piSieve(const ll n){
    if(n<=1) return 0LL;
    if(n==2) return 1LL;
    const int lim=isqrt(n);
    int vsz=(lim+1)>>1;
    std::vector<int> smalls(vsz);
    for(int cx=0;cx<vsz;++cx) smalls[cx]=cx;
    std::vector<int> roughs(vsz);
    for(int cx=0;cx<vsz;++cx) roughs[cx]=(cx<<1|1);
    std::vector<ll> larges(vsz);
    for(int cx=0;cx<vsz;++cx) larges[cx]=(n/(cx<<1|1)-1)>>1;
    std::vector<bool> skips(lim+1);
    int pCnt=0;
    for(int p=3;p<=lim;p+=2){
        if(skips[p]) continue;
        int p2=p*p;
        if(1LL*p2*p2>n) break;
        skips[p]=true;
        for(int cx=p2;cx<=lim;cx+=(p<<1))
            skips[cx]=true;
        int ns=0;
        for(int cx=0;cx<vsz;++cx){
            int cur=roughs[cx];
            if(skips[cur]) continue;
            ll d=1LL*cur*p;
            larges[ns]=larges[cx]-(d<=lim?larges[smalls[d>>1]-pCnt]
                                    :smalls[half(divide(n,d))])+pCnt;
            roughs[ns++]=cur;
        }
        vsz=ns;
        for(int cx=half(lim),cy=((lim/p)-1)|1;cy>=p;cy-=2){
            int cur=smalls[cy>>1]-pCnt;
            for(int cz=(cy*p)>>1;cz<=cx;--cx)
                smalls[cx]-=cur;
        }
        ++pCnt;
    }
    larges[0]+=1LL*(vsz+((pCnt-1)<<1))*(vsz-1)>>1;
    for(int cx=1;cx<vsz;++cx) larges[0]-=larges[cx];
    for(int cx=1;cx<vsz;++cx){
        int q=roughs[cx];
        ll m=n/q;
        int e=smalls[half(m/q)]-pCnt;
        if(e<cx+1) break;
        ll t=0;
        for(int cy=cx+1;cy<=e;++cy)
            t+=smalls[half(divide(m,roughs[cy]))];
        larges[0]+=t-1LL*(e-cx)*(pCnt+cx-1);
    }
    return larges[0]+1;
}
int main(){
    ll n;
    scanf("%lld",&n);
    printf("%lld\n",piSieve(n));
    return 0;
}