#!/usr/bin/env python3
"""
verify_P016.py -- Paper 16: shell completion and Mersenne-prime prediction.

This verifier checks the explicit shell arithmetic, prime-exponent tests, and
the statistical-validation table against the displayed model definition.
"""

from __future__ import annotations

import math
import sys
from pathlib import Path

PASS = FAIL = 0
_N = 0


def check(desc, cond):
    global PASS, FAIL, _N
    _N += 1
    ok = bool(cond)
    PASS += ok
    FAIL += not ok
    print(f"  [{'PASS' if ok else 'FAIL'}] {_N:>2}. {desc}")
    return ok


class Verifier:
    """Same check semantics as verify_common.Verifier; modern output style."""

    def __init__(self, name):
        print(name)

    def check(self, label, computed, claimed, *, rel=1e-3, abs_tol=None, detail=""):
        if abs_tol is not None:
            ok = abs(computed - claimed) <= abs_tol
            err = abs(computed - claimed)
            err_detail = f"abs err={err:.6g}, tol={abs_tol:.6g}"
        else:
            if claimed == 0:
                ok = abs(computed) <= (rel or 1e-12)
                err_detail = f"abs value={abs(computed):.6g}, tol={rel:.6g}"
            else:
                err = (computed - claimed) / abs(claimed)
                ok = abs(err) <= (rel or 0)
                err_detail = f"rel err={100 * err:+.6g}%, tol={100 * (rel or 0):.6g}%"
        return self.record(label, ok, computed, claimed, err_detail + (f"; {detail}" if detail else ""))

    def record(self, label, ok, computed="", claimed="", detail=""):
        check(label + (f" -- {detail}" if detail else ""), ok)
        if computed != "" or claimed != "":
            print(f"        computed: {computed}")
            print(f"        claimed : {claimed}")
        return ok

    def summary(self):
        print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
        return 1 if FAIL else 0


v = Verifier("P016 -- Shell Completion")

ROOT = Path(__file__).resolve().parents[2]
TEX = (ROOT / "toe" / "16_Paper_ShellCompletion.tex").read_text()

PI = math.pi
ALPHA_INV = 4.0 * PI**3 + PI**2 + PI
ALPHA = 1.0 / ALPHA_INV
KAPPA = ALPHA ** 1.25
KAPPA_PRINTED = 0.00213306

# Known Mersenne-prime exponents through the pre-2024 M51.  The entry
# 42,643,801 is easy to miss because it was discovered after larger ones.
M51_EXPONENTS = [
    2,
    3,
    5,
    7,
    13,
    17,
    19,
    31,
    61,
    89,
    107,
    127,
    521,
    607,
    1279,
    2203,
    2281,
    3217,
    4253,
    4423,
    9689,
    9941,
    11213,
    19937,
    21701,
    23209,
    44497,
    86243,
    110503,
    132049,
    216091,
    756839,
    859433,
    1257787,
    1398269,
    2976221,
    3021377,
    6972593,
    13466917,
    20996011,
    24036583,
    25964951,
    30402457,
    32582657,
    37156667,
    42643801,
    43112609,
    57885161,
    74207281,
    77232917,
    82589933,
]

PREDICTED_TRIPLET = [160_964_569, 260_446_093, 421_410_673]
ACTUAL_2024_M52 = 136_279_841
PHI = (1.0 + math.sqrt(5.0)) / 2.0


def shell_distance(value: float) -> float:
    frac = value - math.floor(value)
    return min(frac, 1.0 - frac)


def is_probable_prime(n: int) -> bool:
    if n < 2:
        return False
    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
    for p in small_primes:
        if n % p == 0:
            return n == p

    d = n - 1
    s = 0
    while d % 2 == 0:
        s += 1
        d //= 2

    # Deterministic for these < 2^32 exponents; deliberately stronger than
    # needed so the witness set is obvious in the audit output.
    for a in small_primes:
        if a >= n:
            continue
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(s - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True


def binomial_tail(n: int, observed: int, p: float) -> float:
    return sum(
        math.comb(n, k) * (p**k) * ((1.0 - p) ** (n - k))
        for k in range(observed, n + 1)
    )


s51 = sum(M51_EXPONENTS) * KAPPA
triplet_sum = sum(PREDICTED_TRIPLET)
s54 = s51 + triplet_sum * KAPPA
rounded_proof_s54 = 1_237_518.625 + triplet_sum * KAPPA_PRINTED

distances = []
running = 0.0
for exponent in M51_EXPONENTS:
    running += exponent * KAPPA
    distances.append(shell_distance(running))


v.check("alpha inverse", ALPHA_INV, 137.0363037769, rel=8e-12)
v.check("printed kappa", KAPPA, KAPPA_PRINTED, rel=2e-4)
v.record(
    "M51 exponent list has 51 entries",
    len(M51_EXPONENTS) == 51 and M51_EXPONENTS[-1] == 82_589_933,
    computed=f"count={len(M51_EXPONENTS)}, last={M51_EXPONENTS[-1]}",
    claimed="51 known exponents through p51=82,589,933",
)
v.check("S51 from exact kappa", s51, 1_237_518.625, abs_tol=2.0e-5)
v.check("triplet sum", triplet_sum, 842_821_335, rel=0)
v.check("Z3 coherence", triplet_sum % 3, 0, rel=0)
v.check("target shell from exact kappa", round(s54), 3_035_110, rel=0)
v.check("S54 exact shell distance", shell_distance(s54), 4.8e-6, rel=1e-3)
v.check("first phi ratio", PREDICTED_TRIPLET[1] / PREDICTED_TRIPLET[0], PHI, rel=3e-7)
v.check("second phi ratio", PREDICTED_TRIPLET[2] / PREDICTED_TRIPLET[1], PHI, rel=3e-7)
v.record(
    "predicted exponents pass Miller-Rabin witness set",
    all(is_probable_prime(p) for p in PREDICTED_TRIPLET),
    computed=[(p, is_probable_prime(p)) for p in PREDICTED_TRIPLET],
    claimed="all three exponents prime",
)
v.check("p54 decimal digits", math.floor(PREDICTED_TRIPLET[2] * math.log10(2.0)) + 1, 126_857_254, rel=0)
v.record(
    "prize threshold",
    PREDICTED_TRIPLET[2] * math.log10(2.0) > 1.0e8,
    computed=PREDICTED_TRIPLET[2] * math.log10(2.0),
    claimed="p54 log10(2) > 1e8",
)
v.check("target shell mod 5", 3_035_110 % 5, 0, rel=0)
v.check("target shell mod 13", 3_035_110 % 13, 0, rel=0)
v.record(
    "M52 discrepancy is acknowledged in the TeX",
    "differs from the announced" in TEX and "136{,}279{,}841" in TEX,
    computed="note-on-M52 paragraph present",
    claimed="paper flags predicted p52 differs from discovered M52",
)

v.check(
    "proof line works with printed kappa 0.00213306",
    rounded_proof_s54,
    3_035_109.99999,
    abs_tol=1.0e-4,
    detail="Expected fail: the headline shell closure uses the exact kappa, not the rounded number shown in the proof line.",
)
v.record(
    "theorem's p52 ordinal agrees with discovered M52",
    PREDICTED_TRIPLET[0] == ACTUAL_2024_M52,
    computed=f"predicted first triplet member={PREDICTED_TRIPLET[0]}, discovered M52={ACTUAL_2024_M52}",
    claimed="p52 = 160,964,569",
    detail="Expected fail: the paper acknowledges this later, but the theorem still labels the value p52.",
)
v.check(
    "observed d_n < 0.01 count",
    sum(d < 0.01 for d in distances),
    14,
    rel=0,
    detail="Expected fail under the cumulative S_n definition and the 51 listed exponents.",
)
v.check(
    "observed d_n < 0.05 count",
    sum(d < 0.05 for d in distances),
    19,
    rel=0,
    detail="Expected fail under the cumulative S_n definition and the 51 listed exponents.",
)
v.check(
    "mean shell distance",
    sum(distances) / len(distances),
    0.187,
    rel=2e-3,
    detail="Expected fail: the cumulative exact-kappa mean is about 0.221.",
)
v.check(
    "random expected count for d_n < 0.01 on [0,0.5]",
    51 * (0.01 / 0.5),
    10.2,
    rel=1e-12,
    detail="Expected fail: a uniform shell distance on [0,0.5] gives expectation 1.02, not 10.2.",
)
v.check(
    "random expected count for d_n < 0.05 on [0,0.5]",
    51 * (0.05 / 0.5),
    10.2,
    rel=1e-12,
    detail="Expected fail: this expectation is 5.1, not 10.2.",
)
v.check(
    "one-sided p-value for 14 distances below 0.01",
    binomial_tail(51, 14, 0.01 / 0.5),
    0.03,
    rel=0.1,
    detail="Expected fail: under the stated null, this tail is about 1e-12 and also does not match the observed count.",
)
v.record(
    "unique exhaustive search is reproducible from supplied data",
    False,
    computed="the TeX gives aggregate candidate counts but no search grid, tolerances beyond six bullet constraints, or candidate list",
    claimed="exhaustive search yields one unique triplet",
    detail="Expected reproducibility fail.",
)
v.record(
    "Mersenne primality of 2^p-1 is established by exponent primality",
    False,
    computed="Miller-Rabin verifies only that p is prime; it does not test whether 2^p-1 is a Mersenne prime",
    claimed="prediction is p54 = 421,410,673",
    detail="Expected logical-scope fail.",
)

sys.exit(v.summary())
