spb/prime-mystery-engine

1#!/usr/bin/env python32# =============================================================================3# core.py — Prime Mystery Engine: core primitives4# Author: Simon-Pierre Boucher — contact@spboucher.ai5# =============================================================================6# Building blocks (Phase 0 of the protocol):7#   - primes_upto(n)          : numpy Sieve of Eratosthenes (odd-only)8#   - segmented_primes(lo,hi) : segmented sieve, constant memory per segment9#   - is_prime(n)             : deterministic Miller-Rabin for n < 3.317e2410#   - bpsw(n)                 : Baillie-PSW (MR base 2 + strong Lucas)11#   - prime_count(n)          : pi(n) via sieve (validation helper)12# All routines are pure-Python/numpy, no hidden prime tables beyond the13# 12 deterministic MR bases (which are witnesses, not a prime list).14# =============================================================================1516from __future__ import annotations  # keep annotations lazy (Python 3.9 nodes)1718import numpy as np1920# Deterministic Miller-Rabin witness set: correct for all n < 3.317e2421# (Sorenson & Webster 2015). Documented per the Golden Rule.22_MR_BASES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)23_MR_LIMIT = 3_317_044_064_679_887_385_961_981242526def primes_upto(n: int) -> np.ndarray:27    """All primes <= n, via odd-only Sieve of Eratosthenes (numpy)."""28    if n < 2:29        return np.array([], dtype=np.int64)30    if n < 3:31        return np.array([2], dtype=np.int64)32    # index i represents the odd number 2i+1; index 0 (=1) is not prime33    size = (n + 1) // 234    sieve = np.ones(size, dtype=bool)35    sieve[0] = False36    for i in range(1, (int(n ** 0.5) + 1) // 2 + 1):37        if sieve[i]:38            p = 2 * i + 139            start = (p * p) // 240            sieve[start::p] = False41    odds = 2 * np.nonzero(sieve)[0].astype(np.int64) + 142    return np.concatenate(([np.int64(2)], odds))434445def sieve_segment(lo: int, hi: int, base_primes: np.ndarray | None = None) -> np.ndarray:46    """Primes in [lo, hi) via segmented sieve. base_primes must cover sqrt(hi)."""47    if hi <= 2:48        return np.array([], dtype=np.int64)49    lo = max(lo, 2)50    if base_primes is None:51        base_primes = primes_upto(int(hi ** 0.5) + 1)52    seg = np.ones(hi - lo, dtype=bool)53    for p in base_primes:54        p = int(p)55        if p * p >= hi:56            break57        start = max(p * p, ((lo + p - 1) // p) * p)58        seg[start - lo::p] = False59    if lo <= 1:60        seg[: 2 - lo] = False61    return lo + np.nonzero(seg)[0].astype(np.int64)626364def segmented_primes(lo: int, hi: int, segment_size: int = 10_000_000):65    """Yield numpy arrays of primes covering [lo, hi) segment by segment."""66    base = primes_upto(int(hi ** 0.5) + 1)67    for start in range(lo, hi, segment_size):68        end = min(start + segment_size, hi)69        yield sieve_segment(start, end, base)707172def prime_count(n: int) -> int:73    """pi(n), computed by segmented sieve (validation helper)."""74    total = 075    for chunk in segmented_primes(2, n + 1):76        total += len(chunk)77    return total787980def _mr_witness(n: int, a: int, d: int, r: int) -> bool:81    """True if a is a Miller-Rabin witness for compositeness of n."""82    x = pow(a, d, n)83    if x == 1 or x == n - 1:84        return False85    for _ in range(r - 1):86        x = x * x % n87        if x == n - 1:88            return False89    return True909192def is_prime(n: int) -> bool:93    """Deterministic primality for n < 3.317e24 (Miller-Rabin, 12 bases).94    For larger n, falls back to BPSW (documented: no known counterexample,95    verified exhaustively below 2^64)."""96    if n < 2:97        return False98    for p in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):99        if n % p == 0:100            return n == p101    d, r = n - 1, 0102    while d % 2 == 0:103        d //= 2104        r += 1105    if n < _MR_LIMIT:106        return not any(_mr_witness(n, a, d, r) for a in _MR_BASES if a % n)107    return bpsw(n)108109110def _jacobi(a: int, n: int) -> int:111    """Jacobi symbol (a/n), n odd positive."""112    a %= n113    result = 1114    while a:115        while a % 2 == 0:116            a //= 2117            if n % 8 in (3, 5):118                result = -result119        a, n = n, a120        if a % 4 == 3 and n % 4 == 3:121            result = -result122        a %= n123    return result if n == 1 else 0124125126def _strong_lucas(n: int) -> bool:127    """Strong Lucas probable-prime test (Selfridge parameters). n odd, >2,128    not a perfect square, no small factors."""129    # Selfridge: first D in 5,-7,9,-11,... with Jacobi(D/n) = -1130    D = 5131    while True:132        j = _jacobi(D % n, n)133        if j == -1:134            break135        if j == 0 and abs(D) != n:136            return False137        D = -D - 2 if D > 0 else -D + 2138        if abs(D) > 1_000_000:  # would indicate a square slipped through139            raise ArithmeticError("no suitable D found; n likely a square")140    Q = (1 - D) // 4141    # factor n+1 = d * 2^s142    d, s = n + 1, 0143    while d % 2 == 0:144        d //= 2145        s += 1146    # Lucas sequences U_d, V_d by binary ladder147    U, V, k = 0, 2, 0148    Qk = 1149    inv2 = (n + 1) // 2  # inverse of 2 mod n (n odd)150    for bit in bin(d)[2:]:151        # double: U_{2k}=U_k V_k ; V_{2k}=V_k^2 - 2 Q^k152        U, V = U * V % n, (V * V - 2 * Qk) % n153        Qk = Qk * Qk % n154        if bit == "1":155            # increment: U_{k+1}=(P U + V)/2 ; V_{k+1}=(D U + P V)/2 ; P=1156            U, V = (U + V) * inv2 % n, (D * U + V) * inv2 % n157            Qk = Qk * Q % n158    if U == 0 or V == 0:159        return True160    for _ in range(s - 1):161        V = (V * V - 2 * Qk) % n162        if V == 0:163            return True164        Qk = Qk * Qk % n165    return False166167168def _is_square(n: int) -> bool:169    r = math_isqrt(n)170    return r * r == n171172173from math import isqrt as math_isqrt  # noqa: E402174175176def bpsw(n: int) -> bool:177    """Baillie-PSW: strong MR base 2 + strong Lucas (Selfridge).178    No known counterexample; deterministic below 2^64."""179    if n < 2:180        return False181    for p in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):182        if n % p == 0:183            return n == p184    d, r = n - 1, 0185    while d % 2 == 0:186        d //= 2187        r += 1188    if _mr_witness(n, 2, d, r):189        return False190    if _is_square(n):191        return False192    return _strong_lucas(n)193194195def gaps_in_range(lo: int, hi: int, segment_size: int = 50_000_000):196    """Iterate consecutive prime gaps over [lo, hi). Yields (p, gap) where197    gap = next_prime(p) - p, for all consecutive pairs with p in [lo, hi)."""198    prev = None199    for chunk in segmented_primes(max(lo, 2), hi, segment_size):200        if len(chunk) == 0:201            continue202        if prev is not None:203            yield int(prev), int(chunk[0] - prev)204        diffs = np.diff(chunk)205        for p, g in zip(chunk[:-1], diffs):206            yield int(p), int(g)207        prev = chunk[-1]208