#!/usr/bin/env python3
"""
verify_P194.py — Numerical verification for Addendum 194.

Topic: Identification-rule alternatives for Δ_{S³} eigenvalues → lepton masses.

Five rules are checked:
  R1. Density-weighted eigenvalues: E_ell * <sin chi>_ell
  R2. Heat-kernel identification: t* = ln(R*) / (E2-E1)
  R3. Spectral-zeta identification: s* = ln(R*) / ln(E2/E1)
  R4. GON-lapse identification: two prescriptions (avg, WKB turning point)
  R5. Moment-hierarchy identification: exp(mu1/mu0 * (E2-E1))

Expected output (exit 0):
  All checks PASS.
  R1 density-weighted ratio ≈ 2.571  (factor ~80 short)
  R2 t* ≈ 1.0663,  pi/3 miss = 1.83%
  R3 s* ≈ 5.4358,  7pi/4 miss = 1.14%
  R4a avg lapse ratio ≈ 1.037  (excluded)
  R4b WKB ratio ≈ 1.225  (excluded)
  R5 exp(MU_RATIO*5) ≈ 52.81  (factor 3.9 short)

Copyright: Léon Fernando Vlegels. License: MIT. May 2026.
"""
import sys
import math

# =========================================================================
# High-precision via mpmath (55 digits)
# =========================================================================
try:
    from mpmath import mp, mpf, pi, exp, log, sqrt, sin, cos, fabs, nstr
    mp.dps = 55
    USE_MP = True
except ImportError:
    # Fallback to standard math
    import math as _math
    USE_MP = False
    mpf = float
    pi = _math.pi
    exp = _math.exp
    log = _math.log
    sqrt = _math.sqrt
    sin = _math.sin
    cos = _math.cos
    fabs = abs
    def nstr(x, n): return str(round(float(x), n))

# =========================================================================
# TOE constants
# =========================================================================
ALPHA_INV   = 4*pi**3 + pi**2 + pi        # ≈ 137.036
BREATH      = pi * ALPHA_INV              # ≈ 430.5
FRAC_EDGE   = pi / ALPHA_INV              # ≈ 0.02293
FRAC_BOUND  = pi**2 / ALPHA_INV           # ≈ 0.07202
FRAC_BULK   = 4*pi**3 / ALPHA_INV         # ≈ 0.9053
OMEGA_0     = pi**3 / 4                   # ≈ 7.7516

def mu_n(n):
    """Closed-form moment: mu_n = 16pi^3/(n+4) + 3pi^2/(n+3) + 2pi/(n+2)."""
    return 16*pi**3/(n+4) + 3*pi**2/(n+3) + 2*pi/(n+2)

MU_0    = mu_n(0)       # = alpha^{-1}
MU_1    = mu_n(1)
MU_2    = mu_n(2)
MU_RATIO = MU_1 / MU_0  # ≈ 0.79334

# S3 eigenvalues E_ell = ell*(ell+2)
def E(l): return l*(l+2)
E1, E2, E3 = E(1), E(2), E(3)   # 3, 8, 15

# Experimental target
TARGET = mpf("206.768")          # m_mu / m_e

PASS_COUNT = 0
FAIL_COUNT = 0

def check(name, condition, msg=""):
    global PASS_COUNT, FAIL_COUNT
    if condition:
        PASS_COUNT += 1
    else:
        FAIL_COUNT += 1
    n = PASS_COUNT + FAIL_COUNT
    print(f"  [{'PASS' if condition else 'FAIL'}] {n:>2}. {name}"
          + (f"  [{msg}]" if msg else ""))

# =========================================================================
# Preamble: TOE constant verification
# =========================================================================
print("=" * 72)
print("verify_P194.py  —  Identification-rule alternatives for S³ eigenvalues")
print(f"Precision: {'mpmath mp.dps=55' if USE_MP else 'standard float'}")
print("=" * 72)

print("\n--- TOE constants ---")
check("alpha_inv = 4pi^3+pi^2+pi ≈ 137.036",
      fabs(ALPHA_INV - mpf("137.036")) < mpf("0.001"),
      f"alpha_inv = {float(ALPHA_INV):.8f}")
check("MU_0 = alpha_inv  (monad closure)",
      fabs(MU_0 - ALPHA_INV) < mpf("1e-12"),
      f"MU_0 = {float(MU_0):.8f}")
check("MU_1 ≈ 108.717",
      fabs(MU_1 - mpf("108.717")) < mpf("0.001"),
      f"MU_1 = {float(MU_1):.8f}")
check("MU_RATIO ≈ 0.79334",
      fabs(MU_RATIO - mpf("0.79334")) < mpf("0.00001"),
      f"MU_RATIO = {float(MU_RATIO):.8f}")
check("E1=3, E2=8, E3=15",
      E1==3 and E2==8 and E3==15,
      f"E1={E1}, E2={E2}, E3={E3}")
check("FRAC_EDGE + FRAC_BOUND + FRAC_BULK = 1",
      fabs(FRAC_EDGE + FRAC_BOUND + FRAC_BULK - 1) < mpf("1e-14"),
      f"sum = {float(FRAC_EDGE+FRAC_BOUND+FRAC_BULK):.15f}")

# =========================================================================
# RULE 0: Canonical identification
# =========================================================================
print("\n--- R0: Canonical identification ---")
r0_ratio = mpf(E2) / mpf(E1)
check("R0: E2/E1 = 8/3",
      fabs(r0_ratio - mpf(8)/3) < mpf("1e-14"),
      f"ratio = {float(r0_ratio):.8f}")
check("R0: ratio is factor ~78 short of 207",
      fabs(TARGET / r0_ratio - mpf("77.5")) < mpf("0.1"),
      f"factor = {float(TARGET/r0_ratio):.2f}")

r0b_ratio = exp(MU_RATIO * (E2 - E1))
check("R0b: P18 exp(MU_RATIO*5) ≈ 52.81",
      fabs(r0b_ratio - mpf("52.81")) < mpf("0.01"),
      f"ratio = {float(r0b_ratio):.6f}")
check("R0b: factor ~3.9 short of 207",
      fabs(TARGET / r0b_ratio - mpf("3.9")) < mpf("0.05"),
      f"factor = {float(TARGET/r0b_ratio):.3f}")

# =========================================================================
# RULE 1: Density-weighted eigenvalues
# Fiber density rho_ell = <sin chi>_ell = integral numerically
# =========================================================================
print("\n--- R1: Density-weighted eigenvalues ---")

import numpy as np

def gegenbauer_C1(l, x):
    """Gegenbauer polynomial C^{alpha=1}_l(x). Recurrence: C^1_l = 2x*C^1_{l-1} - C^1_{l-2}."""
    x = np.asarray(x, dtype=float)
    if l == 0: return np.ones_like(x)
    if l == 1: return 2*x
    c_prev = np.ones_like(x)
    c_curr = 2*x
    for _ in range(2, l+1):
        c_next = 2*x*c_curr - c_prev
        c_prev, c_curr = c_curr, c_next
    return c_curr

def fiber_density(l, n_pts=100000):
    """Compute <sin chi>_l via numerical quadrature."""
    chi = np.linspace(1e-6, np.pi - 1e-6, n_pts)
    C = gegenbauer_C1(l, np.cos(chi))
    psi_sq = C**2 * np.sin(chi)**2        # eigenfunction^2 * measure on S3
    norm = np.trapezoid(psi_sq, chi)
    fd = np.trapezoid(psi_sq * np.sin(chi), chi) / norm
    return fd

fd1 = fiber_density(1)
fd2 = fiber_density(2)
fd3 = fiber_density(3)

check("R1: <sin chi>_1 ≈ 0.6791 (fiber density ell=1)",
      abs(fd1 - 0.6791) < 0.0002,
      f"fd1 = {fd1:.8f}")
check("R1: <sin chi>_2 ≈ 0.6548 (fiber density ell=2)",
      abs(fd2 - 0.6548) < 0.0002,
      f"fd2 = {fd2:.8f}")
check("R1: fd1 > fd2  (densities converge from above)",
      fd1 > fd2,
      f"fd1-fd2 = {fd1-fd2:.6f}")
check("R1: fd2 < fd3  (monotone convergence)",
      fd2 > fd3,
      f"fd2-fd3 = {fd2-fd3:.6f}")

r1_ratio = (E2 * fd2) / (E1 * fd1)
check("R1: density-weighted ratio ≈ 2.571",
      abs(r1_ratio - 2.571) < 0.002,
      f"ratio = {r1_ratio:.6f}")
check("R1: ratio < canonical E2/E1 = 8/3",
      r1_ratio < float(r0_ratio),
      f"{r1_ratio:.4f} < {float(r0_ratio):.4f}")
check("R1: ratio factor ~80 short of 207",
      abs(float(TARGET)/r1_ratio - 80) < 2,
      f"factor = {float(TARGET)/r1_ratio:.1f}")

# =========================================================================
# RULE 2: Heat-kernel identification
# =========================================================================
print("\n--- R2: Heat-kernel identification ---")

dE = E2 - E1   # = 5
t_star = log(TARGET) / dE
check("R2: t* = ln(206.768)/5 ≈ 1.06632",
      fabs(t_star - mpf("1.0663")) < mpf("0.0001"),
      f"t* = {float(t_star):.10f}")

# Verify: exp(5 * t*) = TARGET exactly
ratio_at_tstar = exp(dE * t_star)
check("R2: exp(5*t*) = TARGET  (by construction)",
      fabs(ratio_at_tstar - TARGET) < mpf("1e-10"),
      f"exp(5*t*) = {float(ratio_at_tstar):.6f}")

# Compare t* with TOE constants
pct_pi3 = float((t_star / (pi/3) - 1) * 100)
check("R2: t* vs pi/3 — within 2%",
      abs(pct_pi3) < 2.0,
      f"pi/3 = {float(pi/3):.8f},  miss = {pct_pi3:+.2f}%")

pct_frac_b = float((t_star / FRAC_BOUND - 1) * 100)
check("R2: t* / FRAC_BOUNDARY ≈ 14.8  (not a TOE integer)",
      abs(float(t_star / FRAC_BOUND) - 14.8) < 0.2,
      f"ratio = {float(t_star/FRAC_BOUND):.4f}")

# Ratio at t = pi/3
ratio_at_pi3 = exp(5 * pi / 3)
check("R2: exp(5*pi/3) ≈ 187.5  (1.1*pi/3 gives 9% short)",
      fabs(ratio_at_pi3 - mpf("187.5")) < mpf("0.5"),
      f"exp(5*pi/3) = {float(ratio_at_pi3):.4f}")

# t* in units of BREATH_PERIOD
t_star_breath = float(t_star / BREATH)
check("R2: t*/BREATH_PERIOD ≈ 2.5e-3  (tiny fraction — not natural)",
      abs(t_star_breath - 2.48e-3) < 0.05e-3,
      f"t*/BREATH = {t_star_breath:.6e}")

# =========================================================================
# RULE 3: Spectral-zeta identification
# =========================================================================
print("\n--- R3: Spectral-zeta identification ---")

s_star = log(TARGET) / log(mpf(E2) / mpf(E1))
check("R3: s* = ln(206.768)/ln(8/3) ≈ 5.4358",
      fabs(s_star - mpf("5.4358")) < mpf("0.0001"),
      f"s* = {float(s_star):.10f}")

# Verify: (8/3)^s* = TARGET  (bare)
ratio_at_sstar = (mpf(E2)/mpf(E1))**s_star
check("R3: (8/3)^s* = TARGET  (by construction)",
      fabs(ratio_at_sstar - TARGET) < mpf("1e-8"),
      f"ratio = {float(ratio_at_sstar):.6f}")

# Compare s* with TOE candidates
pct_7pi4 = float((s_star / (7*pi/4) - 1) * 100)
check("R3: s* vs 7*pi/4  — within 2%",
      abs(pct_7pi4) < 2.0,
      f"7*pi/4 = {float(7*pi/4):.8f},  miss = {pct_7pi4:+.2f}%")

pct_2pi1 = float((s_star / (2*pi - 1) - 1) * 100)
check("R3: s* vs 2*pi-1  — miss > 2%",
      abs(pct_2pi1) > 2.0,
      f"2*pi-1 = {float(2*pi-1):.6f},  miss = {pct_2pi1:+.2f}%")

# Ratio at s = 7*pi/4
ratio_7pi4 = (mpf(E2)/mpf(E1))**(7*pi/4)
check("R3: (8/3)^(7pi/4) ≈ 219.7  (6.3% over target)",
      fabs(ratio_7pi4 - mpf("219.7")) < mpf("0.2"),
      f"(8/3)^(7pi/4) = {float(ratio_7pi4):.4f}")

# Ratio at s = 5 (integer)
ratio_s5 = (mpf(E2)/mpf(E1))**5
check("R3: (8/3)^5 ≈ 134.8  (35% short, = 8^5/3^5 = 32768/243)",
      fabs(ratio_s5 - mpf("134.8")) < mpf("0.1"),
      f"(8/3)^5 = {float(ratio_s5):.4f}")

# =========================================================================
# RULE 4: GON-lapse identification
# =========================================================================
print("\n--- R4: GON-lapse identification ---")

# Prescription A: eigenfunction-averaged lapse = fiber density
r4a_ratio = fd1 / fd2
check("R4a: avg lapse ratio fd1/fd2 ≈ 1.037",
      abs(r4a_ratio - 1.037) < 0.002,
      f"ratio = {r4a_ratio:.6f}")
check("R4a: ratio << 207  (excluded by magnitude)",
      r4a_ratio < 2.0,
      f"ratio = {r4a_ratio:.4f}")

# Prescription B: WKB turning point sin(chi_tp) = sqrt(ell/(ell+2))
import math as _math
wkb1 = _math.sqrt(1/3)      # = 1/sqrt(3)
wkb2 = _math.sqrt(2/4)      # = 1/sqrt(2)
check("R4b: WKB lapse ell=1: sqrt(1/3) ≈ 0.5774",
      abs(wkb1 - 0.5774) < 0.0001,
      f"m1_wkb = {wkb1:.8f}")
check("R4b: WKB lapse ell=2: sqrt(1/2) ≈ 0.7071",
      abs(wkb2 - 0.7071) < 0.0001,
      f"m2_wkb = {wkb2:.8f}")
check("R4b: m2_wkb > m1_wkb  (wrong direction: heavier has larger lapse)",
      wkb2 > wkb1,
      f"m2/m1 = {wkb2/wkb1:.6f}")
r4b_ratio = wkb2 / wkb1
check("R4b: WKB ratio m2/m1 = sqrt(3/2) ≈ 1.2247",
      abs(r4b_ratio - _math.sqrt(1.5)) < 1e-6,
      f"ratio = {r4b_ratio:.8f}")
check("R4b: ratio << 207  (excluded)",
      r4b_ratio < 2.0,
      f"ratio = {r4b_ratio:.4f}")

# To achieve ratio=207 with m1≤1: m2 ≤ 1/207 → β_2 ≈ 0.0024 rad  (near Hopf pole)
m2_needed = 1.0 / float(TARGET)
check("R4: producing ratio=207 requires m2 ≤ 1/207 ≈ 0.0048",
      abs(m2_needed - 0.00483) < 0.0001,
      f"m2_needed = {m2_needed:.6f}")
beta2_needed = _math.asin(m2_needed) / 2
check("R4: β_2 ≈ 0.0024 rad  (unphysically near Hopf pole β=0)",
      abs(beta2_needed - 0.00242) < 0.0001,
      f"β_2 = {beta2_needed:.6f} rad")
check("R4: β_2 << FRAC_BOUNDARY  (no natural origin in TOE layers)",
      beta2_needed < float(FRAC_BOUND) / 10,
      f"β_2/FRAC_BOUNDARY = {beta2_needed/float(FRAC_BOUND):.4f}")

# =========================================================================
# RULE 5: Moment-hierarchy identification
# =========================================================================
print("\n--- R5: Moment-hierarchy identification ---")

check("R5: moment closed form: mu_0 = alpha_inv",
      fabs(mu_n(0) - ALPHA_INV) < mpf("1e-12"),
      f"mu_0 = {float(mu_n(0)):.8f}")
check("R5: mu_0 > mu_1 > mu_2  (decreasing)",
      mu_n(0) > mu_n(1) > mu_n(2),
      f"mu: {float(mu_n(0)):.4f} > {float(mu_n(1)):.4f} > {float(mu_n(2)):.4f}")

for n in range(5):
    check(f"R5: mu_{n} > 0",
          mu_n(n) > 0,
          f"mu_{n} = {float(mu_n(n)):.6f}")

# Consecutive ratios
mu01 = float(mu_n(0) / mu_n(1))
mu12 = float(mu_n(1) / mu_n(2))
check("R5: mu0/mu1 ≈ 1.260  (far from 207)",
      abs(mu01 - 1.260) < 0.001,
      f"mu0/mu1 = {mu01:.6f}")
check("R5: mu1/mu2 ≈ 1.206",
      abs(mu12 - 1.206) < 0.001,
      f"mu1/mu2 = {mu12:.6f}")
check("R5: all moment ratios in [1.1, 1.3]  (nowhere near 207)",
      all(1.1 < float(mu_n(i)/mu_n(i+1)) < 1.3 for i in range(5)),
      "ratios bounded far from target")

# P18 exponential formula
r5_ratio = exp(MU_RATIO * (E2 - E1))
check("R5: exp(MU_RATIO*5) ≈ 52.81",
      fabs(r5_ratio - mpf("52.81")) < mpf("0.02"),
      f"exp(mu1/mu0 * 5) = {float(r5_ratio):.6f}")
check("R5: factor ~3.9 short of 207",
      fabs(TARGET / r5_ratio - mpf("3.9")) < mpf("0.05"),
      f"factor = {float(TARGET/r5_ratio):.4f}")

# Gap needed
delta_E_needed = log(TARGET) / MU_RATIO
check("R5: to reach 207, need ΔE = ln(207)/MU_RATIO ≈ 6.720",
      fabs(delta_E_needed - mpf("6.720")) < mpf("0.002"),
      f"ΔE* = {float(delta_E_needed):.6f}")
check("R5: actual S3 gap ΔE = E2-E1 = 5  (factor 1.34 short)",
      abs(float(delta_E_needed) / (E2-E1) - 1.34) < 0.01,
      f"ΔE*/ΔE = {float(delta_E_needed)/(E2-E1):.4f}")

# Equivalence: t* = MU_RATIO needed to close gap = same as R2's t*
MU_RATIO_needed = float(log(TARGET)) / (E2-E1)
check("R5/R2 equivalence: MU_RATIO_needed = t* (both = ln(207)/5)",
      abs(MU_RATIO_needed - float(t_star)) < 1e-10,
      f"MU_RATIO_needed = {MU_RATIO_needed:.8f} = t* = {float(t_star):.8f}")

# =========================================================================
# Cross-checks and near-miss characterisation
# =========================================================================
print("\n--- Cross-checks and near-miss characterisation ---")

# R2/R3 near-misses
check("Near-miss R2: |t* - pi/3| < 0.02",
      fabs(t_star - pi/3) < mpf("0.02"),
      f"|t* - pi/3| = {float(fabs(t_star-pi/3)):.6f}")
check("Near-miss R3: |s* - 7pi/4| < 0.07",
      fabs(s_star - 7*pi/4) < mpf("0.07"),
      f"|s* - 7pi/4| = {float(fabs(s_star-7*pi/4)):.6f}")

# s* > 5 (spectral gap dimension interpretation)
check("R3: s* > 5  (not the integer 5 = E2-E1 = dim(B4)+1)",
      s_star > 5,
      f"s* = {float(s_star):.4f} > 5")

# R4 is categorically excluded: max ratio bounded
check("R4 excluded: both lapse ratios << 207",
      r4a_ratio < 5 and r4b_ratio < 5,
      f"R4a={r4a_ratio:.4f}, R4b={r4b_ratio:.4f}")

# R1 is categorically bounded: fd1, fd2 both in (0.6, 0.7)
check("R1 bounded: fiber densities in (0.6, 0.7) for ell=1,2",
      0.6 < fd1 < 0.7 and 0.6 < fd2 < 0.7,
      f"fd1={fd1:.4f}, fd2={fd2:.4f}")

# Summary: only R2/R3 can in principle reach 207, but require unpinned parameter
print()
print("  Classification summary:")
print(f"    R0:  {float(r0_ratio):.4f}  (canonical; factor {float(TARGET/r0_ratio):.1f} short)")
print(f"    R0b: {float(r0b_ratio):.4f}  (P18 exp; factor {float(TARGET/r0b_ratio):.2f} short)")
print(f"    R1:  {r1_ratio:.4f}  (density-weighted; factor {float(TARGET)/r1_ratio:.1f} short)")
print(f"    R2:  t*={float(t_star):.6f}  (near pi/3={float(pi/3):.6f}, miss {pct_pi3:+.2f}%)")
print(f"    R3:  s*={float(s_star):.6f}  (near 7pi/4={float(7*pi/4):.6f}, miss {pct_7pi4:+.2f}%)")
print(f"    R4a: {r4a_ratio:.4f}  (lapse avg; EXCLUDED)")
print(f"    R4b: {r4b_ratio:.4f}  (WKB lapse; EXCLUDED)")
print(f"    R5:  {float(r5_ratio):.4f}  (P18 moment-exp; factor {float(TARGET/r5_ratio):.2f} short)")

# =========================================================================
# Final tally
# =========================================================================
print(f"\n{'='*60}\nRESULT: {PASS_COUNT} PASS / {FAIL_COUNT} FAIL")
sys.exit(0 if FAIL_COUNT == 0 else 1)