Meissel–Lehmer 素数计数

What · 是什么

Meissel–Lehmer 算法在亚线性时间内计算素数计数函数 π(x)=#{p≤x}。思想：定义 φ(x,a)=所有≤x 且最小质因子>p_a 的数的个数，P_k(x,a)=恰有 k 个可重质因子且都>p_a 的数的个数。取 a=π(x^{1/3})，有 π(x)=φ(x,a)+a−1−P₂(x,a)。其中 φ 用递归/容斥（φ(x,b)=φ(x,b−1)−φ(⌊x/p_b⌋,b−1)）配合『普通叶子/特殊叶子』终止条件优化；P₂ 枚举质数对 (p,q) 用 π 前缀表求和。整体 O(x^{2/3}/log²x) 时间、O(x^{1/3}log³x) 空间（取 y=x^{1/3}）。更进一步的 Lehmer 公式与 Lagarias-Miller-Odlyzko (LMO) 变体达到 O(x^{2/3}/log²x) 级别。竞赛/科研中用于 n 极大（如 1e18）的素数计数。

Why · 为什么

How · 实现

伪代码

# 经典 Lehmer 版 (可查表工程实现):
phi(x, s):   # 最小质因子>p_s 的数个数
  if s == 0: return x
  if s == 1: return x - x//2
  if s == 2: return x - x//2 - x//3 + x//6
  if x < primes[s-1]**2: return pi(x) - s + 1
  return phi(x, s-1) - phi(x//primes[s-1], s-1)

pi(x):
  if x < N: return pi_arr[x]          # 预筛表
  a = pi(x**(1/3)); b = pi(isqrt(x))
  ans = phi(x, a) + a - 1 - P2(x, a)  # P2 枚举质数对
  return ans
# P2(x,a) = sum_{p in (p_a, sqrt(x)]} (pi(x//p) - pi(p) + 1)

Python 实现

"""Meissel-Lehmer 亚线性求 pi(x): 素数个数.
实现经典递归版 (用于中小 x 教学; 工程级加速见 OI Wiki 分块优化):
  pi(x) = phi(x, a) + a - 1 - P2(x, a),  a = pi(x^{1/3}), y = x^{1/3}
配合 phi(x, a) 递归 + P2 双素数对计数."""
import math


def sieve_primes(limit):
    is_p = [True] * (limit + 1)
    if limit >= 0: is_p[0] = False
    if limit >= 1: is_p[1] = False
    for i in range(2, int(limit ** 0.5) + 1):
        if is_p[i]:
            for j in range(i * i, limit + 1, i):
                is_p[j] = False
    return [i for i in range(2, limit + 1) if is_p[i]]


class MeisselLehmer:
    def __init__(self, maxn=10**8):
        # 预筛到 maxn, 供 pi 查询
        self.primes = sieve_primes(maxn)
        self.pi_arr = [0] * (maxn + 1)
        j = 0
        for i in range(maxn + 1):
            while j < len(self.primes) and self.primes[j] <= i:
                j += 1
            self.pi_arr[i] = j

    def phi(self, x, a):
        """<=x 且最小质因子 > p_a 的数的个数."""
        if a == 0:
            return x
        if a == 1:
            return x - x // 2
        if a == 2:
            return x - x // 2 - x // 3 + x // 6
        if x < self.primes[a - 1] ** 2:
            # 小 x: 直接数
            return self.pi_small(x) - a + 1
        # 递归: phi(x,a) = phi(x,a-1) - phi(x/p_a, a-1)
        return self.phi(x, a - 1) - self.phi(x // self.primes[a - 1], a - 1)

    def pi_small(self, x):
        if x < len(self.pi_arr):
            return self.pi_arr[x]
        raise ValueError("x too large")

    def pi(self, x):
        if x < len(self.pi_arr):
            return self.pi_arr[x]
        return self.pi_sublinear(x)

    def pi_sublinear(self, x):
        # 主定理: 取 y=x^{1/3}, a=pi(y)
        y = int(x ** (1.0 / 3.0))
        while (y + 1) ** 3 <= x: y += 1
        while y ** 3 > x: y -= 1
        a = self.pi_small(y)
        # P2(x,a) = sum_{p in (y, sqrt x], p<=q, pq<=x} 1
        p2 = 0
        lo = y + 1
        hi = int(math.isqrt(x))
        if lo <= hi:
            primes_in = [p for p in self.primes if lo <= p <= hi]
            for p in primes_in:
                p2 += self.pi_small(x // p) - self.pi_small(p) + 1
        phi_xa = self.phi(x, a)
        return phi_xa + a - 1 - p2


if __name__ == "__main__":
    ml = MeisselLehmer(10**6)
    # 已知 pi 值
    known = {0: 0, 1: 0, 2: 1, 3: 2, 10: 4, 100: 25, 1000: 168,
             10000: 1229, 100000: 9592, 10**6: 78498}
    for x, v in known.items():
        assert ml.pi_small(x) == v, (x, ml.pi_small(x), v)
        assert ml.pi(x) == v, (x, ml.pi(x), v)
    # 跨阈值的亚线性分支 (x > 10^6 走 phi/P2 路径), 用独立筛验证
    for x in [2 * 10**6, 5 * 10**6]:
        cnt = len(sieve_primes(x))
        got = ml.pi(x)
        assert got == cnt, (x, got, cnt)
    print("meissel-lehmer: all tests passed")

易错点

• φ 递归终止条件用错导致指数级展开
• x^{1/3}、√x 用浮点算导致边界差 1
• P2 枚举质数对时对 p 的上界/下界（p>p_a）弄错
• pi 查询超表未递归/未记忆化导致递归爆炸
• Möbius 前缀 S₀ 的容斥符号 μ(n) 写错
• 预筛表 N 太小导致频繁触发亚线性递归（常数差）

例题

信源

• OI Wiki✓ https://oi-wiki.org/math/number-theory/meissel-lehmer/
• 论文✓ https://dl.acm.org/doi/abs/10.1090/s0025-5718-96-00674-6
• 书籍✓ kactl / count_primes 模板（Meissel-Lehmer 工程实现）