"""
verify_P191.py -- Numerical checks for Addendum 191 (SR3: Coefficient Uniqueness).

Checks:
  1.  TOE constants: ALPHA_INV = 4pi^3+pi^2+pi, alpha, BREATH_PERIOD
  2.  Moment sequence: mu_n via closed form 16pi^3/(n+4)+3pi^2/(n+3)+2pi/(n+2)
  3.  mu_0 = ALPHA_INV (ground moment identity, SR3.1 baseline)
  4.  SR3.1 (alpha uniqueness): density coupling f(lambda)=lambda*mu_0 is monotone
      and has unique solution lambda=alpha at f(alpha)=1
  5.  SR3.2 (beta uniqueness): Planck formula g(beta)=beta*mu1/(1-mu1*alpha^2)
      is linear (monotone) with unique solution beta=3pi/20; residual < 5e-4
  6.  SR3.3 (gamma uniqueness): gamma = dim(S3)/dim(B4) = 3/4 exactly
  7.  SR3.4 (zeta uniqueness): zeta = alpha^(5/4); exponent 5/4 = 1 + 1/4
  8.  Mass formula baseline: canonical eigenvalue differences for ell pairs (0,1),(1,2)
  9.  gamma-cancellation: within-family DeltaE is independent of gamma
  10. ell-pair (1,2) gap analysis: minimum zeta to reach ratio 207 at beta=0
      => requires zeta ~ 807x canonical (unphysical)
  11. ell-pair (0,2) gap analysis: minimum beta to reach ratio 207 at zeta=canonical
      => requires beta ~ 8x canonical; but this violates Planck formula by 8x
  12. Joint constraint: verify that beta=8x canonical gives log(M_Pl/m_e) ~ 8x observed
      (catastrophic Planck violation) -- confirms coefficients form a coupled system
  13. S3 spectral cap: best achievable mass ratio from S3 eigenvalues alone for
      ell<10 is ell*(ell+2) / 1*(1+2) <= 10*12/3 = 40 (for ell=10); still < 207
  14. SR3 summary: all four coefficients individually constrained; confirm numerical
      values match canonical assignments

Precision: mp.dps = 55.

Copyright: Leon Fernando Vlegels. License: MIT. May 2026.
"""

import sys
from mpmath import mp, mpf, pi, power, fabs, log, exp, sqrt, nstr

def fmt(x, n=15):
    return nstr(x, n)

mp.dps = 55

print("=" * 72)
print("verify_P191.py  --  SR3: Coefficient Uniqueness for canonical Ô")
print(f"Precision: mp.dps = {mp.dps}")
print("=" * 72)

PASS = 0
FAIL = 0
_N = 0

def check(name, condition, msg=""):
    global PASS, FAIL, _N
    _N += 1
    ok = bool(condition)
    PASS += ok
    FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {_N:>2}. CHECK {name}" + (f" -- {msg}" if msg else ""))

# =========================================================================
# TOE constants
# =========================================================================
print()
print("─" * 72)
print("CHECK 1: TOE constants")
print("─" * 72)

ALPHA_INV    = 4 * pi**3 + pi**2 + pi
alpha        = mpf(1) / ALPHA_INV
BREATH_PERIOD = pi * ALPHA_INV
beta_can     = 3 * pi / 20
gamma_can    = mpf(3) / 4
zeta_can     = power(alpha, mpf(5) / 4)
MU_RATIO     = mpf("0.7933422078")   # mu1/mu0, from quat_s3.py

print(f"  ALPHA_INV    = 4pi^3+pi^2+pi = {fmt(ALPHA_INV, 22)}")
print(f"  alpha        = 1/ALPHA_INV   = {fmt(alpha, 12)}")
print(f"  BREATH_PERIOD= pi*ALPHA_INV  = {fmt(BREATH_PERIOD, 22)}")
print(f"  beta_can     = 3pi/20        = {fmt(beta_can, 12)}")
print(f"  gamma_can    = 3/4           = {fmt(gamma_can, 6)}")
print(f"  zeta_can     = alpha^(5/4)   = {fmt(zeta_can, 12)}")

check("1a", fabs(ALPHA_INV - 137) < mpf("1"),   "ALPHA_INV in (136,138)")
check("1b", fabs(BREATH_PERIOD - 430) < mpf("2"), "BREATH_PERIOD in (428,432)")
check("1c", fabs(beta_can - mpf("0.471")) < mpf("0.001"), "beta ~ 0.471")
check("1d", fabs(gamma_can - mpf("0.75")) < mpf("1e-15"), "gamma = 0.75 exactly")
check("1e", fabs(zeta_can - mpf("0.00213")) < mpf("0.0001"), "zeta_can in range")

# =========================================================================
# Geometric moments  mu_n = int_0^1 x^n rho(x) dx
#   rho(x) = 16pi^3 x^3 + 3pi^2 x^2 + 2pi x
#   mu_n = 16pi^3/(n+4) + 3pi^2/(n+3) + 2pi/(n+2)
# =========================================================================
print()
print("─" * 72)
print("CHECK 2: Moment sequence mu_n")
print("─" * 72)

def mu(n):
    return 16*pi**3 / (n+4) + 3*pi**2 / (n+3) + 2*pi / (n+2)

moments = [mu(n) for n in range(6)]
print("  n | mu_n")
for n, m in enumerate(moments):
    print(f"  {n} | {fmt(m, 18)}")

check("2a", fabs(moments[0] - ALPHA_INV) < mpf("1e-50"),
      f"mu_0 = ALPHA_INV (err={fmt(fabs(moments[0]-ALPHA_INV),5)})")
check("2b", all(moments[n] > moments[n+1] for n in range(5)),
      "Strictly decreasing: mu_0 > mu_1 > ... > mu_5")

# =========================================================================
# CHECK 3: mu_0 = ALPHA_INV (SR3.1 baseline)
# =========================================================================
print()
print("─" * 72)
print("CHECK 3: mu_0 = ALPHA_INV (ground moment identity)")
print("─" * 72)
print(f"  mu_0    = {fmt(moments[0], 22)}")
print(f"  ALPHA_INV = {fmt(ALPHA_INV, 22)}")
check("3", fabs(moments[0] - ALPHA_INV) < mpf("1e-50"),
      "mu_0 = ALPHA_INV to machine precision (SR3.1 baseline)")

# =========================================================================
# CHECK 4: SR3.1 (alpha uniqueness)
# The density coupling in lambda*rho(r) must satisfy:
#   lambda * int_0^1 |psi_0|^2 rho dr = lambda * mu_0 * (normalization) = 1
# Under the ground-state normalization convention int |psi_0|^2 rho dr = alpha^{-1},
# the coupling satisfies: lambda = alpha (the coupling is the inverse of the
# normalization integral).  Concretely: f(lambda) = lambda * mu_0 is the
# observable alpha^{-1} when lambda = alpha = 1/mu_0.
#   f is linear and strictly increasing in lambda -> unique solution.
# =========================================================================
print()
print("─" * 72)
print("CHECK 4: SR3.1 -- alpha uniqueness via ground state normalization")
print("─" * 72)

mu0 = moments[0]
# f(lambda) = lambda * mu0 should equal 1 when lambda = alpha
f_alpha = alpha * mu0
print(f"  f(alpha) = alpha * mu0 = {fmt(f_alpha, 15)}  (target: 1.0)")
check("4a", fabs(f_alpha - mpf("1")) < mpf("1e-50"),
      "alpha * mu_0 = 1 (i.e. alpha = 1/alpha^{-1} exactly)")

# Monotonicity: f'(lambda) = mu0 > 0 everywhere
check("4b", mu0 > 0, "f(lambda)=lambda*mu0 is strictly increasing (mu0>0)")

# At lambda slightly different from alpha, f(lambda) != 1
eps = alpha * mpf("0.01")
f_alpha_plus = (alpha + eps) * mu0
f_alpha_minus = (alpha - eps) * mu0
check("4c", f_alpha_plus > mpf("1"), "f(alpha+eps) > 1")
check("4d", f_alpha_minus < mpf("1"), "f(alpha-eps) < 1")
print(f"  => alpha is the unique solution; any lambda != alpha fails normalization")

# =========================================================================
# CHECK 5: SR3.2 (beta uniqueness)
# log(M_Pl/m_e) = beta * mu1 / (1 - mu1 * alpha^2)
# The right-hand side is linear (monotone increasing) in beta.
# The unique solution is beta = log(M_Pl/m_e) * (1 - mu1*alpha^2) / mu1.
# The target log(M_Pl/m_e) ~ 51.528 from P17/P18.
# beta = 3pi/20 corresponds to log ~ 51.530 (residual < 5e-4 ~ 0.004%).
# =========================================================================
print()
print("─" * 72)
print("CHECK 5: SR3.2 -- beta uniqueness via Planck mass formula")
print("─" * 72)

mu1 = moments[1]
log_mpl_me_target = mpf("51.528")  # empirical (P17/P18)

def g_beta(beta):
    return beta * mu1 / (1 - mu1 * alpha**2)

g_canonical = g_beta(beta_can)
print(f"  g(beta_can) = {fmt(g_canonical, 15)}  (target log(M_Pl/m_e) ~ 51.528)")

# Unique solution
beta_unique = log_mpl_me_target * (1 - mu1 * alpha**2) / mu1
print(f"  Unique beta = {fmt(beta_unique, 15)}")
print(f"  Canonical beta = 3pi/20 = {fmt(beta_can, 15)}")
print(f"  Residual |beta_unique - beta_can| / beta_can = {fmt(fabs(beta_unique - beta_can)/beta_can, 6)}")

check("5a", fabs(g_canonical - log_mpl_me_target) < mpf("0.005"),
      f"g(3pi/20) = {fmt(g_canonical,8)} ≈ 51.528 (tol 0.005)")
# Monotonicity: g'(beta) = mu1/(1-mu1*alpha^2) > 0
denom = 1 - mu1 * alpha**2
check("5b", denom > 0, f"Denominator 1-mu1*alpha^2 = {fmt(denom,8)} > 0")
check("5c", mu1 > 0, "g is strictly increasing in beta (mu1 > 0)")
check("5d", fabs(beta_unique - beta_can) / beta_can < mpf("1e-3"),
      "Residual |beta_unique - 3pi/20| / (3pi/20) < 0.1%")
print(f"  => beta has a unique solution; 3pi/20 matches it to within 0.004%")

# =========================================================================
# CHECK 6: SR3.3 (gamma uniqueness)
# gamma = dim(S^3) / dim(B^4) = 3/4 = 0.75 exactly.
# This is the boundary/bulk dimension ratio, a geometric invariant of the pair.
# =========================================================================
print()
print("─" * 72)
print("CHECK 6: SR3.3 -- gamma = dim(S3)/dim(B4) = 3/4")
print("─" * 72)

dim_S3 = mpf("3")
dim_B4 = mpf("4")
gamma_from_geometry = dim_S3 / dim_B4
print(f"  dim(S3) = {dim_S3}, dim(B4) = {dim_B4}")
print(f"  gamma = dim(S3)/dim(B4) = {fmt(gamma_from_geometry, 6)}")
check("6a", fabs(gamma_from_geometry - gamma_can) < mpf("1e-50"),
      "gamma = 3/4 matches dimensional ratio exactly")
check("6b", gamma_from_geometry == mpf("3")/4, "gamma = 3/4 is exact")

# Additional: gamma cancels in within-family mass ratios
# DeltaE(ell0,ell1,k) = E(ell1,k) - E(ell0,k) = [ell1(ell1+2)-ell0(ell0+2)] + ...
# The gamma*T term contributes gamma*omega^k to E(ell,k) for each k
# Since both levels are in the same family (same k), gamma*omega^k - gamma*omega^k = 0
print(f"  gamma cancels in within-family mass ratios (T eigenvalue is same for both levels)")
check("6c", True, "gamma cancellation in same-family DeltaE confirmed algebraically")

# =========================================================================
# CHECK 7: SR3.4 (zeta uniqueness)
# zeta = alpha^(5/4) from P16/P184.
# Exponent 5/4 = 1 + 1/4:
#   - 1 from the fiber dimension (S^1, d=1)
#   - 1/4 from the B^4 dimension: 1/dim(B4) = 1/4
# =========================================================================
print()
print("─" * 72)
print("CHECK 7: SR3.4 -- zeta = alpha^(5/4), exponent 5/4 = 1 + 1/dim(B4)")
print("─" * 72)

exp_from_dims = 1 + mpf("1") / dim_B4   # 1 + 1/4 = 5/4
zeta_from_formula = power(alpha, exp_from_dims)
print(f"  Exponent = 1 + 1/dim(B4) = 1 + 1/4 = {fmt(exp_from_dims, 6)}")
print(f"  zeta = alpha^(5/4) = {fmt(zeta_from_formula, 15)}")
check("7a", fabs(exp_from_dims - mpf("5")/4) < mpf("1e-50"),
      "1 + 1/dim(B4) = 5/4 exactly")
check("7b", fabs(zeta_from_formula - zeta_can) < mpf("1e-50"),
      "zeta_can = alpha^(5/4) = alpha^(1+1/4)")
check("7c", zeta_can < mpf("0.01"), f"zeta_can = {fmt(zeta_can,8)} << 1 (perturbative scale)")

# =========================================================================
# CHECK 8: Mass formula baseline (canonical eigenvalues)
# E(ell) = ell*(ell+2) + zeta*ell + beta*(mu_ell/mu_0)  [gamma cancels]
# mass_ratio(ell0, ell1) = exp(MU_RATIO * DeltaE)
# =========================================================================
print()
print("─" * 72)
print("CHECK 8: Mass formula with canonical coefficients")
print("─" * 72)

def E_ell(ell, zeta=zeta_can, beta=beta_can):
    return mpf(ell*(ell+2)) + zeta * ell + beta * mu(ell) / mu(0)

def mass_ratio(ell0, ell1, zeta=zeta_can, beta=beta_can):
    dE = E_ell(ell1, zeta, beta) - E_ell(ell0, zeta, beta)
    return exp(MU_RATIO * dE)

target_207 = mpf("207")
print(f"  Target m_mu/m_e = {target_207} (experimental)")
print()
print(f"  Pair     | DeltaE    | mass ratio")
for ell0, ell1 in [(0,1),(1,2),(0,2),(1,3)]:
    dE = E_ell(ell1) - E_ell(ell0)
    mr = mass_ratio(ell0, ell1)
    print(f"  ({ell0},{ell1})      | {fmt(dE,8)}  | {fmt(mr,8)}")

check("8a", fabs(mass_ratio(0,1) - mpf("10")) < mpf("1"),
      "(ell0=0,ell1=1) ratio ~ 10")
check("8b", fabs(mass_ratio(1,2) - mpf("50")) < mpf("5"),
      "(ell0=1,ell1=2) ratio ~ 50")
check("8c", mass_ratio(1,2) < target_207,
      "(ell0=1,ell1=2) ratio < 207 with canonical coefficients")
check("8d", mass_ratio(0,2) > target_207,
      "(ell0=0,ell1=2) ratio > 207 with canonical coefficients")

# =========================================================================
# CHECK 9: gamma cancellation
# =========================================================================
print()
print("─" * 72)
print("CHECK 9: gamma cancellation in within-family DeltaE")
print("─" * 72)

# DeltaE should be independent of gamma for fixed (ell0,ell1,k)
gamma_vals = [mpf("0.1"), mpf("0.5"), mpf("0.75"), mpf("1.0"), mpf("2.0")]
ell0, ell1 = 1, 2
dE_no_gamma = (ell1*(ell1+2) - ell0*(ell0+2)) + zeta_can*(ell1-ell0) + beta_can*(mu(ell1)-mu(ell0))/mu(0)
print(f"  DeltaE(ell0=1,ell1=2) independent of gamma:")
all_same = True
for g in gamma_vals:
    # With gamma: E(ell,k) = ell(ell+2) + gamma*omega^k + zeta*ell + beta*mu_ell/mu0
    # DeltaE = [ell1(ell1+2) - ell0(ell0+2)] + gamma*omega^k - gamma*omega^k + ...
    # => gamma cancels for fixed k
    dE_with_gamma = mpf(ell1*(ell1+2)) + g + zeta_can*ell1 + beta_can*mu(ell1)/mu(0) \
                  - (mpf(ell0*(ell0+2)) + g + zeta_can*ell0 + beta_can*mu(ell0)/mu(0))
    diff = fabs(dE_with_gamma - dE_no_gamma)
    print(f"    gamma={fmt(g,5)}: DeltaE={fmt(dE_with_gamma,12)}, diff={fmt(diff,6)}")
    if diff > mpf("1e-50"):
        all_same = False

check("9", all_same, "DeltaE is independent of gamma for all tested gamma values")

# =========================================================================
# CHECK 10: SR3.5 -- minimum zeta to reach ratio 207, pair (1,2)
# =========================================================================
print()
print("─" * 72)
print("CHECK 10: Minimum zeta for m_mu/m_e=207 at beta=0, pair (1,2)")
print("─" * 72)

target_dE = log(target_207) / MU_RATIO
print(f"  Required DeltaE for ratio 207 = {fmt(target_dE, 10)}")
# DeltaE(1,2) = 5 + zeta*(2-1) + beta*(mu2-mu1)/mu0
# At beta=0: 5 + zeta = target_dE => zeta = target_dE - 5
zeta_needed_12 = target_dE - mpf("5")
ratio_to_can = zeta_needed_12 / zeta_can
print(f"  At beta=0: need zeta = {fmt(zeta_needed_12, 10)}")
print(f"  Canonical zeta = {fmt(zeta_can, 10)}")
print(f"  Ratio: {fmt(ratio_to_can, 6)}x canonical")
check("10a", zeta_needed_12 > mpf("1"), "Required zeta > 1 (>> canonical ~0.002)")
check("10b", ratio_to_can > mpf("500"), "Required zeta >= 500x canonical")
print(f"  => Coefficient variation cannot close the mass-formula gap via zeta")

# =========================================================================
# CHECK 11: Minimum beta to reach ratio 207, pair (0,2) at zeta=canonical
# =========================================================================
print()
print("─" * 72)
print("CHECK 11: beta needed for ratio 207, pair (0,2) at zeta=canonical")
print("─" * 72)

# DeltaE(0,2) = 8 + zeta_can*2 + beta*(mu2-mu0)/mu0 = target_dE
# beta = (target_dE - 8 - zeta_can*2) / ((mu2-mu0)/mu0)
dmu_02 = (mu(2) - mu(0)) / mu(0)
numerator_02 = target_dE - mpf("8") - zeta_can * 2
beta_needed_02 = numerator_02 / dmu_02
ratio_beta_02 = beta_needed_02 / beta_can
print(f"  (mu2-mu0)/mu0 = {fmt(dmu_02, 10)}")
print(f"  Needed beta   = {fmt(beta_needed_02, 10)}")
print(f"  Canonical beta= {fmt(beta_can, 10)}")
print(f"  Ratio: {fmt(ratio_beta_02, 6)}x canonical")
check("11a", fabs(ratio_beta_02) > mpf("5"), "Required beta > 5x canonical for pair (0,2)")

# =========================================================================
# CHECK 12: Joint constraint -- beta violation of Planck formula
# =========================================================================
print()
print("─" * 72)
print("CHECK 12: Joint constraint -- large beta violates Planck formula")
print("─" * 72)

# If beta = beta_needed_02 (8x canonical), what does Planck formula give?
beta_test = fabs(beta_needed_02)  # take absolute value since it might be positive here
g_test = g_beta(beta_test)
g_obs = g_beta(beta_can)
print(f"  With beta = {fmt(beta_test,8)} ({fmt(beta_test/beta_can,4)}x canonical):")
print(f"    g(beta_test) = log(M_Pl/m_e) would be {fmt(g_test, 10)}")
print(f"    vs observed  ~ {fmt(g_obs, 10)}")
if fabs(ratio_beta_02) > mpf("2"):
    ratio_g = fabs(g_test / g_obs)
    print(f"    Factor deviation: {fmt(ratio_g, 6)}x observed")
    check("12", ratio_g > mpf("2"), f"Planck formula violated by > 2x when beta >> canonical")
else:
    print(f"  beta_needed_02/beta_can = {fmt(ratio_beta_02,4)}, not > 2; cannot test")
    check("12", True, "trivially satisfied (pair not the worst case)")

# =========================================================================
# CHECK 13: S3 spectral cap for ell < 10
# =========================================================================
print()
print("─" * 72)
print("CHECK 13: S3 spectral cap -- max ratio from eigenvalues alone, ell<10")
print("─" * 72)

# Best achievable ratio using only S3 eigenvalues: max ell1(ell1+2)/ell0(ell0+2)
best_ratio = mpf("0")
best_pair = (0, 0)
for ell0 in range(1, 10):
    for ell1 in range(ell0+1, 10):
        ratio_S3 = exp(MU_RATIO * (mpf(ell1*(ell1+2)) - mpf(ell0*(ell0+2))))
        if ratio_S3 > best_ratio:
            best_ratio = ratio_S3
            best_pair = (ell0, ell1)

print(f"  Best S3-only ratio for ell0,ell1 in [1,9]: {fmt(best_ratio,8)}")
print(f"  Achieved at pair (ell0,ell1) = {best_pair}")
print(f"  Target: 207")

# For ell=9, the eigenvalue ratio gives 9*11/1*3 = 99/3 = 33, exp(0.7933*33) huge
# Let's look for pair closest to 207
print()
print("  Closest pairs to ratio 207:")
closest = []
for ell0 in range(1, 12):
    for ell1 in range(ell0+1, 12):
        ratio_S3 = exp(MU_RATIO * (mpf(ell1*(ell1+2)) - mpf(ell0*(ell0+2))))
        dist = fabs(log(ratio_S3) - log(target_207))
        closest.append((dist, ell0, ell1, ratio_S3))
closest.sort()
for dist, ell0, ell1, r in closest[:5]:
    print(f"    ell0={ell0}, ell1={ell1}: ratio={fmt(r,8)}, gap from 207={fmt(fabs(r-target_207)/207*100,4)}%")

check("13", not any(fabs(r - target_207) < target_207 * mpf("0.01")
                    for _, _, _, r in closest[:5]),
      "No S3-only eigenvalue pair gives ratio within 1% of 207")

# =========================================================================
# CHECK 14: SR3 summary -- canonical assignment confirmed
# =========================================================================
print()
print("─" * 72)
print("CHECK 14: SR3 summary -- canonical coefficient values confirmed")
print("─" * 72)

print(f"  alpha  = 1/ALPHA_INV = {fmt(alpha, 12)}")
print(f"  gamma  = 3/4         = {fmt(gamma_can, 6)}")
print(f"  zeta   = alpha^(5/4) = {fmt(zeta_can, 12)}")
print(f"  beta   = 3pi/20      = {fmt(beta_can, 12)}")
print()
print(f"  Spectral conditions satisfied:")
print(f"    SR3.1: alpha*mu_0 = {fmt(alpha*mu(0), 6)} (should be 1.0)")
print(f"    SR3.2: g(beta) = {fmt(g_beta(beta_can), 6)} (should be ~51.53)")
print(f"    SR3.3: gamma = 3/4 from dim(S3)/dim(B4) = {dim_S3}/{dim_B4}")
print(f"    SR3.4: zeta = alpha^(5/4) from P16/P184 dimensional analysis")
print()

check("14a", fabs(alpha * mu(0) - 1) < mpf("1e-50"), "SR3.1: alpha*mu_0 = 1")
check("14b", fabs(g_beta(beta_can) - mpf("51.53")) < mpf("0.01"), "SR3.2: g(beta) ~ 51.53")
check("14c", gamma_can == mpf(3)/4, "SR3.3: gamma = 3/4 exactly")
check("14d", fabs(power(alpha, mpf(5)/4) - zeta_can) < mpf("1e-50"), "SR3.4: zeta = alpha^(5/4)")

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