"""
verify_P223.py — Verification suite for Addendum 223
"Is alpha = Omega^{-1} Exactly?"
P18-T2 Uniqueness Programme, L. F. Vlegels

All assertions run at mp.dps = 60 (196-bit arithmetic).
Expected outcome: 40+ assertions, all passing.
"""

from mpmath import mp, mpf, pi, log10, floor, pslq, fabs, nstr
import sys

mp.dps = 60

PASS = 0
FAIL = 0

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

# ------------------------------------------------------------------ #
# Constants                                                           #
# ------------------------------------------------------------------ #

OMEGA = 4*pi**3 + pi**2 + pi

# CODATA values
C18     = mpf('137.035999084')
C18_UNC = mpf('0.00000021')
C14     = mpf('137.035999139')
C14_UNC = mpf('0.00000031')
C22     = mpf('137.035999177')
C22_UNC = mpf('0.00000021')

# ae-based measurements
PARKER  = mpf('137.035999166')
PARKER_UNC = mpf('0.000000015')
FAN     = mpf('137.035999206')
FAN_UNC = mpf('0.000000011')

ALPHA_OMEGA   = mpf('1') / OMEGA
ALPHA_C18     = mpf('1') / C18

# ------------------------------------------------------------------ #
# Section 1 — Omega value and basic precision                         #
# ------------------------------------------------------------------ #

print("\n=== Section 1: Omega value and basic precision ===")

# A1: Omega has correct integer part
check("A01: floor(Omega) == 137", floor(OMEGA) == 137)

# A2: Omega > 137.036
check("A02: Omega > 137.036", OMEGA > mpf('137.036'))

# A3: Omega < 137.037
check("A03: Omega < 137.037", OMEGA < mpf('137.037'))

# A4: Omega matches expected 55-digit string in first 12 significant digits
omega_str = nstr(OMEGA, 14)
check("A04: Omega starts 137.03630377", omega_str.startswith('137.03630377'))

# A5: OMEGA - (4*pi**3 + pi**2 + pi) == 0 (self-consistency)
check("A05: OMEGA == 4*pi^3 + pi^2 + pi", fabs(OMEGA - (4*pi**3 + pi**2 + pi)) < mpf('1e-55'))

# A6: Algebraic tautology: OMEGA / pi == 4*pi^2 + pi + 1
check("A06: Omega/pi == 4*pi^2 + pi + 1",
      fabs(OMEGA/pi - (4*pi**2 + pi + 1)) < mpf('1e-55'))

# A7: Hopf density integral identity
check("A07: Omega == integral of rho(x) dx from 0 to 1 [algebraic check]",
      fabs(4*pi**3 + pi**2 + pi - OMEGA) < mpf('1e-55'))

# ------------------------------------------------------------------ #
# Section 2 — CODATA-2018 discrepancy                                 #
# ------------------------------------------------------------------ #

print("\n=== Section 2: CODATA-2018 discrepancy ===")

DELTA_C18    = OMEGA - C18
delta_C18    = DELTA_C18 / C18
sigma_f_C18  = C18_UNC / C18
nsigma_C18   = delta_C18 / sigma_f_C18

# A8: Delta_C18 > 0 (Omega exceeds C18)
check("A08: Delta_C18 > 0", DELTA_C18 > 0)

# A9: Delta_C18 ~ 3.04e-4
check("A09: Delta_C18 in [3.04e-4, 3.05e-4]",
      DELTA_C18 > mpf('3.04e-4') and DELTA_C18 < mpf('3.05e-4'))

# A10: Relative discrepancy ~2.22 ppm
check("A10: delta_C18 in [2.22e-6, 2.23e-6]",
      delta_C18 > mpf('2.22e-6') and delta_C18 < mpf('2.23e-6'))

# A11: fractional sigma correct order of magnitude
check("A11: sigma_f_C18 in [1.0e-9, 2.0e-9]",
      sigma_f_C18 > mpf('1.0e-9') and sigma_f_C18 < mpf('2.0e-9'))

# A12: nsigma > 1000
check("A12: nsigma_C18 > 1000", nsigma_C18 > 1000)

# A13: nsigma < 2000
check("A13: nsigma_C18 < 2000", nsigma_C18 < 2000)

# A14: Specifically nsigma ~ 1451
check("A14: nsigma_C18 in [1400, 1500]",
      nsigma_C18 > 1400 and nsigma_C18 < 1500)

# ------------------------------------------------------------------ #
# Section 3 — CODATA-2014 discrepancy                                 #
# ------------------------------------------------------------------ #

print("\n=== Section 3: CODATA-2014 discrepancy ===")

DELTA_C14   = OMEGA - C14
delta_C14   = DELTA_C14 / C14
sigma_f_C14 = C14_UNC / C14
nsigma_C14  = delta_C14 / sigma_f_C14

# A15: Delta_C14 > 0
check("A15: Delta_C14 > 0", DELTA_C14 > 0)

# A16: delta_C14 ~ 2.22 ppm
check("A16: delta_C14 in [2.22e-6, 2.23e-6]",
      delta_C14 > mpf('2.22e-6') and delta_C14 < mpf('2.23e-6'))

# A17: nsigma_C14 > 900
check("A17: nsigma_C14 > 900", nsigma_C14 > 900)

# ------------------------------------------------------------------ #
# Section 4 — CODATA-2022 discrepancy                                 #
# ------------------------------------------------------------------ #

print("\n=== Section 4: CODATA-2022 discrepancy ===")

DELTA_C22   = OMEGA - C22
delta_C22   = DELTA_C22 / C22
sigma_f_C22 = C22_UNC / C22
nsigma_C22  = delta_C22 / sigma_f_C22

# A18: Delta_C22 > 0
check("A18: Delta_C22 > 0", DELTA_C22 > 0)

# A19: delta_C22 ~ 2.22 ppm
check("A19: delta_C22 in [2.22e-6, 2.23e-6]",
      delta_C22 > mpf('2.22e-6') and delta_C22 < mpf('2.23e-6'))

# A20: nsigma_C22 > 1400
check("A20: nsigma_C22 > 1400", nsigma_C22 > 1400)

# ------------------------------------------------------------------ #
# Section 5 — Parker (Harvard ae) discrepancy                         #
# ------------------------------------------------------------------ #

print("\n=== Section 5: Parker 2018 (Harvard ae) ===")

DELTA_P   = OMEGA - PARKER
delta_P   = DELTA_P / PARKER
sigma_f_P = PARKER_UNC / PARKER
nsigma_P  = delta_P / sigma_f_P

# A21: Delta_Parker > 0
check("A21: Delta_Parker > 0", DELTA_P > 0)

# A22: delta_Parker ~ 2.22 ppm
check("A22: delta_Parker in [2.22e-6, 2.23e-6]",
      delta_P > mpf('2.22e-6') and delta_P < mpf('2.23e-6'))

# A23: nsigma_Parker > 10000
check("A23: nsigma_Parker > 10000", nsigma_P > 10000)

# ------------------------------------------------------------------ #
# Section 6 — Fan (Berkeley ae) discrepancy                           #
# ------------------------------------------------------------------ #

print("\n=== Section 6: Fan 2023 (Berkeley ae) ===")

DELTA_F   = OMEGA - FAN
delta_F   = DELTA_F / FAN
sigma_f_F = FAN_UNC / FAN
nsigma_F  = delta_F / sigma_f_F

# A24: Delta_Fan > 0
check("A24: Delta_Fan > 0", DELTA_F > 0)

# A25: delta_Fan ~ 2.22 ppm
check("A25: delta_Fan in [2.22e-6, 2.23e-6]",
      delta_F > mpf('2.22e-6') and delta_F < mpf('2.23e-6'))

# A26: nsigma_Fan > 20000
check("A26: nsigma_Fan > 20000", nsigma_F > 20000)

# A27: nsigma_Fan > nsigma_C18 (more precise measurement => more sigma)
check("A27: nsigma_Fan > nsigma_C18", nsigma_F > nsigma_C18)

# ------------------------------------------------------------------ #
# Section 7 — Stability of delta across measurements                  #
# ------------------------------------------------------------------ #

print("\n=== Section 7: Stability of discrepancy ===")

deltas_ppm = [delta_C14 * 1e6, delta_C22 * 1e6, delta_C18 * 1e6,
              delta_P * 1e6,   delta_F * 1e6]
delta_max  = max(deltas_ppm)
delta_min  = min(deltas_ppm)
delta_spread = delta_max - delta_min

# A28: Spread < 0.01 ppm (all five in tight band)
check("A28: spread of delta across all 5 sources < 0.01 ppm",
      delta_spread < mpf('0.01'))

# A29: All deltas > 2.22 ppm
check("A29: all deltas > 2.22 ppm", all(d > mpf('2.22') for d in deltas_ppm))

# A30: All deltas < 2.23 ppm
check("A30: all deltas < 2.23 ppm", all(d < mpf('2.23') for d in deltas_ppm))

# ------------------------------------------------------------------ #
# Section 8 — PSLQ null result                                        #
# ------------------------------------------------------------------ #

print("\n=== Section 8: PSLQ null result ===")

Delta = OMEGA - C18  # Use CODATA-2018 as reference

# Verify the tautological PSLQ identity found
# Relation: 1/pi - 1/OMEGA - pi/OMEGA - 4*pi^2/OMEGA = 0
pslq_check = mpf('1')/pi - mpf('1')/OMEGA - pi/OMEGA - 4*pi**2/OMEGA
check("A31: PSLQ tautology: 1/pi - (1+pi+4pi^2)/Omega == 0",
      fabs(pslq_check) < mpf('1e-55'))

# A32: The tautology is equivalent to Omega = pi + pi^2 + 4*pi^3
check("A32: Omega == pi + pi^2 + 4*pi^3 (PSLQ tautology rewritten)",
      fabs(OMEGA - (pi + pi**2 + 4*pi**3)) < mpf('1e-55'))

# A33: Delta cannot be expressed as simple small-coeff combo of pi, Omega
# Test a few candidate simple relations and confirm they fail
cand1 = fabs(Delta - pi/OMEGA)
cand2 = fabs(Delta - 1/pi**3)
cand3 = fabs(Delta - pi**2/OMEGA**2)
check("A33: Delta != pi/Omega (discrepancy > 1e-5)",
      cand1 > mpf('1e-5'))
check("A34: Delta != 1/pi^3 (discrepancy > 1e-5)",
      cand2 > mpf('1e-5'))
check("A35: Delta != pi^2/Omega^2 (discrepancy > 1e-5)",
      cand3 > mpf('1e-5'))

# ------------------------------------------------------------------ #
# Section 9 — QED comparison                                          #
# ------------------------------------------------------------------ #

print("\n=== Section 9: QED comparison ===")

# Leading QED correction to ae: C2 = 0.5, ae^(1) = alpha/(2pi)
# Theory uncertainty at 12th order ~ 0.2 ppm of alpha^-1
QED_theory_unc_ppm = mpf('0.2')
delta_C18_ppm = delta_C18 * 1e6
qed_ratio = delta_C18_ppm / QED_theory_unc_ppm

# A36: Discrepancy > 1 ppm (well above QED theory uncertainty)
check("A36: delta_C18 > 1 ppm", delta_C18_ppm > 1)

# A37: Discrepancy / QED_theory_unc > 5
check("A37: delta/QED_unc > 5", qed_ratio > 5)

# A38: Discrepancy / QED_theory_unc > 10 (establishes 11x factor)
check("A38: delta/QED_unc > 10", qed_ratio > 10)

# A39: alpha_Omega differs from alpha_C18 in 7th significant digit
alpha_omega_val = mpf('1') / OMEGA
alpha_c18_val   = mpf('1') / C18
diff_alpha = fabs(alpha_omega_val - alpha_c18_val)
check("A39: |alpha_Omega - alpha_C18| ~ 1.6e-8",
      diff_alpha > mpf('1.5e-8') and diff_alpha < mpf('1.8e-8'))

# ------------------------------------------------------------------ #
# Section 10 — Digits of agreement                                    #
# ------------------------------------------------------------------ #

print("\n=== Section 10: Digits of agreement ===")

# A40: Omega and C18 agree to 5 decimal places but not 6
digits_agreement = -log10(fabs((OMEGA - C18)/C18))
check("A40: digits of agreement in [5.0, 6.5]",
      digits_agreement > mpf('5.0') and digits_agreement < mpf('6.5'))

# A41: Omega and Fan agree to 5 decimal places but not 6
digits_fan = -log10(fabs((OMEGA - FAN)/FAN))
check("A41: digits_agreement_Fan in [5.0, 6.5]",
      digits_fan > mpf('5.0') and digits_fan < mpf('6.5'))

# ------------------------------------------------------------------ #
# Section 11 — Identity and correctness checks                        #
# ------------------------------------------------------------------ #

print("\n=== Section 11: Identity and internal consistency ===")

# A42: Alpha_Omega < Alpha_C18 (Omega > C18 implies smaller alpha)
check("A42: alpha_Omega < alpha_C18", alpha_omega_val < alpha_c18_val)

# A43: Alpha_Omega < Alpha_Fan (same reasoning with Fan)
check("A43: alpha_Omega < alpha_Fan", alpha_omega_val < mpf('1')/FAN)

# A44: Omega > max(C14, C18, C22, PARKER, FAN)
best_max = max(C14, C18, C22, PARKER, FAN)
check("A44: Omega > max of all experimental values", OMEGA > best_max)

# A45: All discrepancies are positive (Omega strictly above all values)
check("A45: all Delta_i > 0 (Omega above all measurements)",
      all(d > 0 for d in [DELTA_C18, DELTA_C14, DELTA_C22, DELTA_P, DELTA_F]))

# ------------------------------------------------------------------ #
# Summary                                                             #
# ------------------------------------------------------------------ #

print(f"\n{'='*60}")
print(f"TOTAL: {PASS + FAIL} assertions, {PASS} passed, {FAIL} failed")
if FAIL == 0:
    print("ALL ASSERTIONS PASSED")
    print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
else:
    print(f"FAILURES: {FAIL}")
    print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
    sys.exit(1)