Meissel–Lehmer 算法

前言

推荐先行阅读我的blog文章————Min_25 筛

什么是Meissel–Lehmer 算法

Meissel-Lehmer 算法是一种基于 \(ϕ\) 函数的的快速计算前缀质数个数（当然也可以推广到前缀和质数幂次）的算法。

【模板例题】Meissel–Lehmer 算法

给定整数 \(n\)，求出 \(\pi(n)\) 的值。
\(\pi(n)\) 表示 \(1 \sim n\) 的整数中质数的个数。
对于 \(100\%\) 的数据，\(1 \leq n \leq 10^{13}\)。

解法

前置知识：Eratosthenes 筛法。Eratosthenes 筛法是最简单的筛法之一

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;
}

我们只需要知道在标星号的那一行代码执行时有多少个 is_prime 从 true 变为 false 。如果我们可以维护这个值，我们就可以知道从 1 到 n 的素数个数了。据说，这个东西在 Project Euler 上曾经有一个所谓 Lucy DP 的东西可以解决，但是这个不是现在我们想要的。

定义： 记 is_prime 为 \(a\)，\(S(v,p)\) 为 \(a_j ( 2\le j\le v)\)在 \(i=p\)时处理完毕后为 true 的个数。

然后上标处理，虽然时空复杂度达不到线性，但很快了。

之后修为Min_25筛发
在 min_25 的代码中的一些变量解释：

\(\texttt{smalls}[i] = S(i, p)\)
\(\texttt{roughs}[i]\)：没有被筛掉的第 \(i\) 个整数。
其大小由一个变量 \(\texttt s\) 维护。
把 \(1\) 留下，过筛时使用的素数会消失。
\(\texttt{larges}[i] = S(\lfloor n/\texttt{roughs}[i]\rfloor, p)\)
\(\texttt{pc} = \pi(p-1)\)
\(\texttt{skip}[i]\)：\(i\) 如果被筛出来的话就是真的。
但是不更新偶数。

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;
}

仔细看上面的代码应该能大致了解算法，要是不懂可以看这里我的一些解释：

roughs 的更新：roughs[..s]是未过筛的整数。ns 是下一个循环中的 s。ns = 0，在 roughs[ns] 中有值的话 ns 就自增。
larges[smalls[d]-pc] 的含义：smalls[d] - pc 是 “未过筛剩下的整数中 d 是第几个数字?” 的意思。
smalls[i] 的更新：保持 \(j = \frac i p\)被更新。当 \(i \ge pj\) 时，smalls[i] 被更新，它的值被减去 smalls[j] - pc，也就是被减去 smalls[i / p] - pc。

完整代码

#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;
}