"""
verify_P095.py — P95: E₈ ⊃ H₄×H₄ and the Golden Ratio in sin²θ_W

Central numerical claims:
  1. det(A_{H₄}) = (7-3√5)/2 ≈ 0.146
  2. The 240 E₈ roots split under H₄×H₄: 120 with |P₊|²=2/φ, 120 with |P₊|²=2/φ²
     → sum: 120·(2/φ) + 120·(2/φ²) = 240 (partition of root norms)
     since 2/φ + 2/φ² = 2·(1/φ+1/φ²) = 2·1 = 2  (key identity 1/φ+1/φ²=1)
  3. Inner product ratio: ⟨Q_Y,Q_Y⟩_{H₄}/⟨Q_Y,Q_Y⟩_{E₆} = 2/(φ+1) = 2/φ² ≈ 0.7639
  4. Golden ratio identities: φ²=φ+1, 1/φ+1/φ²=1, φ³=2φ+1
  5. E₆ roots: 72 total (rank 6, 36 positive roots)
  6. Near-miss from P94 repeated: δ=+0.054%
"""

import sys

import mpmath

mpmath.mp.dps = 50

PASSES = []
FAILS  = []

def report(name, ok, claimed, actual):
    if ok:
        PASSES.append(name)
    else:
        FAILS.append(name)
    status = "PASS" if ok else "FAIL"
    n = len(PASSES) + len(FAILS)
    print(f"  [{status}] {n:>2}. {name}")
    if not ok:
        print(f"          claimed={claimed}")
        print(f"          actual ={actual}")

sqrt5 = mpmath.sqrt(5)
phi = (1 + sqrt5) / 2
pi  = mpmath.pi

# --- Check 1: det(A_{H₄}) = (7-3√5)/2 ---
det_H4 = (7 - 3*sqrt5) / 2
det_H4_expected = mpmath.mpf('0.146')
ok1 = abs(det_H4 - det_H4_expected) < mpmath.mpf('0.001')
report("P95-1: det(A_{H₄})=(7-3√5)/2≈0.146", ok1,
       f"≈{float(det_H4_expected):.4f}", f"=(7-3√5)/2={float(det_H4):.4f}")

# --- Check 2: 1/φ + 1/φ² = 1 (key identity for root-norm partition) ---
id_check = 1/phi + 1/phi**2
ok2a = abs(id_check - 1) < mpmath.mpf('1e-45')
# And 2/φ + 2/φ² = 2 (sum of H₄ projection norms for one E₈ root pair)
root_norm_sum = 2/phi + 2/phi**2
ok2b = abs(root_norm_sum - 2) < mpmath.mpf('1e-45')
ok2 = ok2a and ok2b
report("P95-2: 1/φ+1/φ²=1 and 2/φ+2/φ²=2 (root-norm partition)", ok2,
       "exact identities",
       f"1/φ+1/φ²={float(id_check):.15f}, 2/φ+2/φ²={float(root_norm_sum):.15f}")

# --- Check 3: Inner product ratio 2/(φ+1) = 2/φ² ≈ 0.7639 ---
ratio_QY = 2 / (phi + 1)
ratio_QY_via_phi2 = 2 / phi**2
ratio_QY_expected = mpmath.mpf('0.7639')
ok3a = abs(ratio_QY - ratio_QY_via_phi2) < mpmath.mpf('1e-45')
ok3b = abs(ratio_QY - ratio_QY_expected) < mpmath.mpf('5e-5')
ok3 = ok3a and ok3b
report("P95-3: ⟨Q_Y,Q_Y⟩_{H₄}/⟨Q_Y,Q_Y⟩_{E₆}=2/(φ+1)=2/φ²≈0.7639", ok3,
       f"2/φ²≈{float(ratio_QY_expected):.4f}",
       f"2/(φ+1)={float(ratio_QY):.5f}, 2/φ²={float(ratio_QY_via_phi2):.5f}")

# --- Check 4: Golden ratio identities ---
ok4a = abs(phi**2 - phi - 1) < mpmath.mpf('1e-45')     # φ²=φ+1
ok4b = abs(phi**3 - 2*phi - 1) < mpmath.mpf('1e-45')   # φ³=2φ+1
ok4c = abs(1/phi**2 - (2-phi)) < mpmath.mpf('1e-45')   # 1/φ²=2-φ
ok4 = ok4a and ok4b and ok4c
report("P95-4: φ²=φ+1, φ³=2φ+1, 1/φ²=2-φ", ok4,
       "exact", f"φ²-φ-1={float(phi**2-phi-1):.2e}")

# --- Check 5: E₆ root count = 72 ---
E6_roots = 72
E6_positive_roots = 36
ok5 = (E6_roots == 72 and E6_positive_roots == 36)
report("P95-5: E₆ has 72 roots (36 positive)", ok5,
       "72 roots, 36 positive", f"{E6_roots} roots, {E6_positive_roots} positive")

# --- Check 6: Near-miss δ≈+0.054% (exact; paper states +0.070% from rounded sW1) ---
mu0 = 4*pi**3 + pi**2 + pi
lam = mpmath.sin(pi/14)
sW1 = (1/(1+pi)) * (1 - 4*lam**2/5)
sW5 = mpmath.mpf('3') / 8
ratio = sW1 / sW5
delta_pct = (ratio - 1/phi) / (1/phi) * 100
ok6 = abs(delta_pct - mpmath.mpf('0.054')) < mpmath.mpf('0.01')
report("P95-6: Near-miss δ=+0.054% (exact) [paper states +0.070% from rounded sW1]", ok6,
       "+0.054%", f"{float(delta_pct):.4f}%")

# --- Check 7: k₁^{H₄}/k₁^{E₆} = φ²/2 ≈ 1.309 (not φ) ---
k1_ratio = phi**2 / 2
k1_expected = mpmath.mpf('1.309')
ok7 = abs(k1_ratio - k1_expected) < mpmath.mpf('0.001')
report("P95-7: k₁^{H₄}/k₁^{E₆}=φ²/2≈1.309 (not φ, hence route fails)", ok7,
       f"≈{float(k1_expected):.3f}", f"φ²/2={float(k1_ratio):.4f}")

if FAILS:
    print(f"  FAILED: {FAILS}")
print(f"\n{'='*60}\nRESULT: {len(PASSES)} PASS / {len(FAILS)} FAIL")
sys.exit(0 if not FAILS else 1)