#!/usr/bin/env python3
"""
verify_P233.py — Verifier for Addendum 233: SR3-ζ Hopf Degree Coupling Proposition

Checks (35 total, mp.dps=60):
  Section 1  — TOE constants and ζ value                         (C01–C06)
  Section 2  — Dimensional formula p = 1 + dim(S¹)/dim(B⁴)      (C07–C14)
  Section 3  — Formula equivalences (P16/P18 forms)              (C15–C20)
  Section 4  — Alternative Hopf fibers: excluded exponents        (C21–C26)
  Section 5  — Killing norm matching (P184)                       (C27–C30)
  Section 6  — Mass-ratio gap consistency (P196 Thm A)           (C31–C33)
  Section 7  — SR3 status cross-checks                           (C34–C35)

All arithmetic uses mpmath at dps=60.
Copyright: Léon Fernando Vlegels. License: MIT. 2026-05-23.
"""

from mpmath import mp, mpf, pi, log, fabs, exp, power
mp.dps = 60

# ── assertion harness ─────────────────────────────────────────────────────────
PASS = 0
FAIL = 0

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

# ── TOE constants ─────────────────────────────────────────────────────────────
OMEGA   = 4*pi**3 + pi**2 + pi        # alpha^{-1}  (P03 identity)
ALPHA   = 1 / OMEGA                    # fine-structure constant
BETA    = 3*pi / 20                    # moment coupling
GAMMA   = mpf('3') / 4                 # layer-cycle coefficient

p_hopf   = mpf('5') / mpf('4')        # HDCP exponent
zeta_can = ALPHA ** p_hopf             # zeta = alpha^{5/4}

dim_S1  = mpf('1')                     # dim Hopf fiber S^1
dim_S2  = mpf('2')                     # dim Hopf base S^2
dim_S3  = mpf('3')                     # dim boundary sphere S^3 = dB^4
dim_B4  = mpf('4')                     # dim bulk B^4

TARGET_RATIO = mpf('206.7682830')      # experimental m_mu/m_e (CODATA 2018)

# ── mass-formula helpers (P196) ───────────────────────────────────────────────
MU0 = 16*pi**3/4 + 3*pi**2/3 + 2*pi/2          # = OMEGA = alpha^{-1}

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

MU1 = mu(1)
MU2 = mu(2)
MU3 = mu(3)
R10 = MU1 / MU0    # ratio used in mass formula

def s3_eig(ell):
    ell = mpf(ell)
    return ell*(ell + 2)

def delta_E(p, ell0, ell1):
    p = mpf(p)
    return (s3_eig(ell1) - s3_eig(ell0)) + ALPHA**p*(ell1 - ell0) + BETA*(mu(ell1) - mu(ell0))/MU0

def m_ratio(p, ell0, ell1):
    return exp(R10 * delta_E(p, ell0, ell1))

# =============================================================================
print("\n=== Section 1: TOE constants and zeta value ===")

check("C01  alpha^-1 = 4pi^3+pi^2+pi approx 137.036",
      fabs(OMEGA - mpf('137.036')) < mpf('0.001'))

check("C02  alpha = 1/OMEGA approx 7.297e-3",
      fabs(ALPHA - mpf('7.2974e-3')) < mpf('1e-6'))

check("C03  p_hopf = 5/4 exactly",
      fabs(p_hopf - mpf('1.25')) < mpf('1e-55'))

check("C04  zeta_can = alpha^{5/4} approx 2.133e-3",
      fabs(zeta_can - mpf('2.133e-3')) < mpf('5e-6'))

check("C05  zeta_can > 0  (positive, as required for coupling)",
      zeta_can > 0)

check("C06  zeta_can < alpha  (zeta < alpha since p > 1 and alpha < 1)",
      zeta_can < ALPHA)

# =============================================================================
print("\n=== Section 2: Dimensional formula p = 1 + dim(S^1)/dim(B^4) ===")

check("C07  dim(S^1) = 1",
      fabs(dim_S1 - 1) < mpf('1e-55'))

check("C08  dim(B^4) = 4",
      fabs(dim_B4 - 4) < mpf('1e-55'))

check("C09  dim(S^1)/dim(B^4) = 1/4",
      fabs(dim_S1/dim_B4 - mpf('1')/4) < mpf('1e-55'))

p_formula = 1 + dim_S1/dim_B4
check("C10  p_formula = 1 + dim(S^1)/dim(B^4) = 5/4",
      fabs(p_formula - p_hopf) < mpf('1e-55'))

check("C11  p_formula expressed as 5/4 (exact rational)",
      fabs(p_formula - mpf('5')/4) < mpf('1e-55'))

check("C12  p_formula > 1  (radial contribution is the base)",
      p_formula > 1)

check("C13  p_formula - 1 = 1/4  (fiber contribution alone)",
      fabs(p_formula - 1 - mpf('1')/4) < mpf('1e-55'))

check("C14  1/dim(B^4) = 0.25 exactly",
      fabs(1/dim_B4 - mpf('0.25')) < mpf('1e-55'))

# =============================================================================
print("\n=== Section 3: Formula equivalences (P16/P18 forms) ===")

# P16 form: 1 + dim(S^3)/(dim(B^4)*3)
p_P16 = 1 + dim_S3 / (dim_B4 * 3)
check("C15  P16 form: 1 + dim(S^3)/(dim(B^4)*3) = 1 + 3/12 = 5/4",
      fabs(p_P16 - p_hopf) < mpf('1e-55'))

# P18 form: 1 + 3/(4*3)
p_P18 = 1 + mpf('3') / (4*3)
check("C16  P18 form: 1 + 3/(4*3) = 5/4",
      fabs(p_P18 - p_hopf) < mpf('1e-55'))

check("C17  P16 form = P18 form = HDCP form (all three identical)",
      fabs(p_P16 - p_P18) < mpf('1e-55') and fabs(p_P16 - p_formula) < mpf('1e-55'))

# Hopf decomposition: dim(S^3) = dim(S^1) + dim(S^2)
check("C18  Hopf decomposition: dim(S^3) = dim(S^1) + dim(S^2) = 1 + 2 = 3",
      fabs(dim_S1 + dim_S2 - dim_S3) < mpf('1e-55'))

# Factor reduction: dim(S^3)/(dim(B^4)*dim(S^3)) = dim(S^1)/dim(B^4)
p_fiber_P16 = dim_S3 / (dim_B4 * dim_S3)  # = 1/dim(B4)
check("C19  dim(S^3)/(dim(B^4)*dim(S^3)) = 1/dim(B^4) = dim(S^1)/dim(B^4)",
      fabs(p_fiber_P16 - dim_S1/dim_B4) < mpf('1e-55'))

# Fiber fraction: S^1 is 1/3 of S^3 weight * 3/4 bulk weight = 1/4
fiber_of_S3   = dim_S1 / dim_S3        # = 1/3
boundary_of_B4 = dim_S3 / dim_B4      # = 3/4
fiber_of_B4   = fiber_of_S3 * boundary_of_B4
check("C20  Fiber weight: (dim(S^1)/dim(S^3))*(dim(S^3)/dim(B^4)) = 1/3 * 3/4 = 1/4",
      fabs(fiber_of_B4 - dim_S1/dim_B4) < mpf('1e-55'))

# =============================================================================
print("\n=== Section 4: Alternative Hopf fibers -- excluded exponents ===")

# S^3 fiber (hypothetical): p = 1 + dim(S^3)/dim(B^4) = 1 + 3/4 = 7/4
p_S3_fiber = 1 + dim_S3/dim_B4
check("C21  S^3 fiber gives p = 7/4 != 5/4  (excluded by P18 assembly choice)",
      fabs(p_S3_fiber - mpf('7')/4) < mpf('1e-55') and fabs(p_S3_fiber - p_hopf) > mpf('0.4'))

# S^2 fiber (hypothetical): p = 1 + dim(S^2)/dim(B^4) = 1 + 2/4 = 3/2
p_S2_fiber = 1 + dim_S2/dim_B4
check("C22  S^2 fiber gives p = 3/2 != 5/4  (not the Hopf fiber)",
      fabs(p_S2_fiber - mpf('3')/2) < mpf('1e-55') and fabs(p_S2_fiber - p_hopf) > mpf('0.2'))

# p = 1 (radial only, no fiber): excluded
p_radial_only = mpf('1')
check("C23  p = 1 (radial only) != 5/4  (missing fiber contribution 1/4)",
      fabs(p_radial_only - p_hopf) > mpf('0.2'))

# p = gamma = 3/4: the layer coefficient, not an exponent for zeta
check("C24  p != gamma = 3/4  (gamma is the layer coefficient, not the Hopf exponent)",
      fabs(p_hopf - GAMMA) > mpf('0.4'))

# p = 2 (hypothetical B^4 x S^1 product without projection): excluded
p_product = mpf('2')
check("C25  p = 2 != 5/4  (naive product dim=2 is excluded)",
      fabs(p_product - p_hopf) > mpf('0.7'))

# p = 1/2: another candidate excluded
p_half = mpf('1') / 2
check("C26  p = 1/2 != 5/4  (below 1, reverses radial scaling direction)",
      fabs(p_half - p_hopf) > mpf('0.7'))

# =============================================================================
print("\n=== Section 5: Killing norm matching (P184) ===")

KILLING_NORM_SQ = mpf('24')    # |H_{U(1)}|^2 = 24 = 2|Phi_{G_2}| (P182)

check("C27  |H_{U(1)}|^2 = 24 = 2|Phi_{G_2}|  (P182 / P184)",
      fabs(KILLING_NORM_SQ - 24) < mpf('1e-55'))

# |zeta*H_{U(1)}|^2 = zeta^2 * 24 = alpha^{5/2} * 24
killing_op       = zeta_can**2 * KILLING_NORM_SQ
killing_expected = ALPHA**(mpf('5')/2) * 24
check("C28  |zeta*H_{U(1)}|^2 = alpha^{5/2} * 24  (Killing norm with zeta = alpha^{5/4})",
      fabs(killing_op - killing_expected) < mpf('1e-50'))

# Check that zeta^2 = alpha^{5/2}
check("C29  zeta^2 = alpha^{5/2}  (squaring the 5/4 exponent gives 5/2)",
      fabs(zeta_can**2 - ALPHA**(mpf('5')/2)) < mpf('1e-55'))

# Normalization c=1: zeta_can / alpha^{5/4} = 1
c_factor = zeta_can / ALPHA**p_hopf
check("C30  zeta_can / alpha^{5/4} = 1  (normalization c = 1; no extra factor)",
      fabs(c_factor - 1) < mpf('1e-55'))

# =============================================================================
print("\n=== Section 6: Mass-ratio gap consistency (P196 Thm A) ===")

r01 = m_ratio(p_hopf, 0, 1)
r12 = m_ratio(p_hopf, 1, 2)
r23 = m_ratio(p_hopf, 2, 3)

print(f"       pair (0,1) at p=5/4: {float(r01):.4f}")
print(f"       pair (1,2) at p=5/4: {float(r12):.4f}")
print(f"       pair (2,3) at p=5/4: {float(r23):.4f}")

check("C31  pair (0,1) at p=5/4: m_ratio < 206.77  (below muon/electron gap)",
      r01 < TARGET_RATIO)

check("C32  pair (1,2) at p=5/4: m_ratio < 206.77  (below muon/electron gap)",
      r12 < TARGET_RATIO)

check("C33  pair (2,3) at p=5/4: m_ratio > 206.77  (above muon/electron gap)",
      r23 > TARGET_RATIO)

# =============================================================================
print("\n=== Section 7: SR3 status cross-checks ===")

# SR3-alpha: MU0 = OMEGA = alpha^{-1}  (P03 identity)
check("C34  SR3-alpha: MU0 = OMEGA = 4pi^3+pi^2+pi  (P03 identity holds)",
      fabs(MU0 - OMEGA) < mpf('1e-50'))

# SR3-zeta: p = 1 + dim(S^1)/dim(B^4) = 5/4 uniquely  (HDCP)
check("C35  SR3-zeta: p = 1 + dim(S^1)/dim(B^4) = 5/4  (HDCP closes SR3-zeta)",
      fabs(p_formula - mpf('5')/4) < mpf('1e-55'))

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