"""
verify_P155.py — Addendum P155: Flag-Variety Orbit Average and Gap G1

Computes three integrals:
  I_G2   : G₂ Haar average of B(Ad_g H_s, H_s)²   (analytic, Schur)
  I_SU3  : SU(3) Haar average of B_{g₂}(Ad_g H_s, H_s)²   (analytic + MC)
  TARGET : Ω/9 = π·α⁻¹/9

Copyright: Léon Fernando Vlegels. MIT License.
"""

import numpy as np
import math
import sys

PASS = FAIL = 0
def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond); PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")

# ── Constants ────────────────────────────────────────────────────────────────
ALPHA_INV     = 4*math.pi**3 + math.pi**2 + math.pi   # ≈ 137.036
BREATH_PERIOD = math.pi * ALPHA_INV                    # Ω ≈ 430.512
TARGET        = BREATH_PERIOD / 9                      # Ω/9 ≈ 47.835
DISCRETE      = 48                                     # B_{G₂}(H_s,H_s) from P151
G1            = (432 - BREATH_PERIOD) / BREATH_PERIOD  # ≈ 3.456e-3

print("=" * 65)
print("P155 FLAG-VARIETY ORBIT AVERAGE VERIFICATION")
print("=" * 65)
print(f"ALPHA_INV        = {ALPHA_INV:.10f}")
print(f"BREATH_PERIOD    = {BREATH_PERIOD:.10f}")
print(f"TARGET (Ω/9)     = {TARGET:.10f}")
print(f"DISCRETE         = {DISCRETE}")
print(f"G1 = (432-Ω)/Ω  = {G1:.8e}  ({G1*100:.6f}%)")
print()

# ── Integral I: G₂ Haar average (analytic) ───────────────────────────────────
# By Schur orthogonality for the adjoint representation of G₂ (dim=14, real):
#   ∫_{G₂} B(Ad_g H_s, H_s)² dg = B(H_s,H_s)² / dim(G₂) = 48²/14
print("-" * 65)
print("INTEGRAL I: G₂ HAAR AVERAGE  (Schur, analytic)")
print("-" * 65)
I_G2 = float(DISCRETE**2) / 14.0
print(f"I_G₂ = B(H_s,H_s)² / dim(G₂) = {DISCRETE}²/14 = {I_G2:.6f}")
print(f"TARGET                                          = {TARGET:.6f}")
print(f"ratio I_G₂ / TARGET = {I_G2/TARGET:.6f}  (overshoots by {(I_G2/TARGET-1)*100:.2f}%)")
print()

# ── Analytic setup for SU(3) Haar average ────────────────────────────────────
# B_{su3}(H_s,H_s) = Σ_{α∈Φ(A₂)} α(H_s)² :
#   A₂ roots ±(eᵢ−eⱼ), projections onto H_s=diag(1,−1,0): ±2, ±1, ∓1
A2_root_projs = [2, -2, 1, -1, -1, 1]   # 6 roots of A₂
B_su3_Hs = sum(x**2 for x in A2_root_projs)   # = 12

# Schur for adjoint of SU(3) (dim=8, real orthogonal):
#   ∫_{SU(3)} f(g)² dg = B_{su3}(H_s,H_s)² / dim(su(3))
#   where f(g) = B_{su3}(Ad_g H_s, H_s)
# But we need B_{g₂}, not B_{su3}.  The functional f(g) = B_{su3}(Ad_g H_s, H_s) is
# directly measured by the formula B_{su3}(h, H_s) = 6(h₁ − h₂) (derived in §2 of the paper).
# This is still the su(3) Killing form.  The G₂ Killing form restricted to su(3) is
# B_{g₂}|_{su(3)×su(3)} = c · B_{su3};  the Monte Carlo measures B_{su3}.
I_SU3_analytic = float(B_su3_Hs)**2 / 8.0   # = 144/8 = 18

print("-" * 65)
print("INTEGRAL II: SU(3) HAAR AVERAGE  (Schur, analytic)")
print("-" * 65)
print(f"A₂ root projections onto H_s = diag(1,−1,0): {A2_root_projs}")
print(f"B_{{su3}}(H_s,H_s) = Σα(H_s)² = {B_su3_Hs}")
print(f"Schur formula: I_SU3 = B²/dim(su3) = {B_su3_Hs}²/8 = {I_SU3_analytic:.6f}")
print(f"TARGET                                          = {TARGET:.6f}")
print(f"ratio I_SU3 / TARGET = {I_SU3_analytic/TARGET:.6f}  (undershoots by {(1 - I_SU3_analytic/TARGET)*100:.2f}%)")
print()

# Structural identity: DISCRETE = B_{su3}² / |Φ+(A₂)|
B_sq_over_pos_roots = float(B_su3_Hs)**2 / 3.0   # 144/3 = 48
print(f"Structural identity: B_{{su3}}(H_s,H_s)²/|Φ⁺(A₂)| = {B_su3_Hs}²/3 = {B_sq_over_pos_roots:.1f}")
check(1, f"DISCRETE = {DISCRETE}  (match: B_su3(H_s,H_s)\u00b2/|\u03a6\u207a(A\u2082)| == DISCRETE)",
      B_sq_over_pos_roots == DISCRETE)
print()

# ── Monte Carlo: SU(3) Haar average ──────────────────────────────────────────
def random_su3(rng):
    """Haar-distributed SU(3) via QR of complex Gaussian (Mezzadri 2007)."""
    z = (rng.standard_normal((3, 3)) + 1j * rng.standard_normal((3, 3))) / math.sqrt(2)
    q, r = np.linalg.qr(z)
    # Fix diagonal of R to be real positive — ensures Haar measure
    d = np.diag(r)
    q = q * (d / np.abs(d))[np.newaxis, :]
    # Project to SU(3): remove determinant phase
    det = np.linalg.det(q)
    q[:, 0] /= det
    return q

H_s_mat = np.diag([1.0, -1.0, 0.0])
N_MC = 200_000
rng = np.random.default_rng(42)

print("-" * 65)
print(f"INTEGRAL II: SU(3) HAAR AVERAGE  (Monte Carlo, N={N_MC:,})")
print("-" * 65)

vals = np.empty(N_MC)
for i in range(N_MC):
    g = random_su3(rng)
    h = g @ H_s_mat @ g.conj().T
    h_diag = np.real(np.diag(h))
    # B_{su3}(Ad_g H_s, H_s) = 6(h₁ − h₂)  (see §2)
    B_val = 6.0 * (h_diag[0] - h_diag[1])
    vals[i] = B_val**2

I_SU3_mc  = float(np.mean(vals))
std_err   = float(np.std(vals)) / math.sqrt(N_MC)
print(f"I_SU3 (MC)       = {I_SU3_mc:.6f}  ±  {std_err:.6f}  (1σ)")
print(f"I_SU3 (analytic) = {I_SU3_analytic:.6f}")
print(f"MC vs analytic discrepancy: {abs(I_SU3_mc - I_SU3_analytic):.4f}  ({abs(I_SU3_mc-I_SU3_analytic)/I_SU3_analytic*100:.3f}%)")
print(f"TARGET (Ω/9)     = {TARGET:.6f}")
print(f"ratio I_SU3_MC / TARGET = {I_SU3_mc/TARGET:.6f}")
print(f"gap I_SU3 − TARGET      = {I_SU3_mc - TARGET:.4f}")
print()

# ── Summary table ─────────────────────────────────────────────────────────────
print("=" * 65)
print("SUMMARY TABLE")
print("=" * 65)
header = f"{'Quantity':<38} {'Value':>10}  {'Ratio/TARGET':>13}"
print(header)
print("-" * 65)
rows = [
    ("Discrete B_{G₂}(H_s,H_s)  [P151]",    DISCRETE,          DISCRETE/TARGET),
    ("G₂ Haar avg  48²/14        [analytic]", I_G2,              I_G2/TARGET),
    ("SU(3) Haar avg  12²/8      [analytic]", I_SU3_analytic,    I_SU3_analytic/TARGET),
    ("SU(3) Haar avg             [MC]",       I_SU3_mc,          I_SU3_mc/TARGET),
    ("TARGET  Ω/9",                           TARGET,            1.0),
]
for label, val, ratio in rows:
    print(f"  {label:<36} {val:>10.4f}  {ratio:>13.6f}")
print()

# ── Interpolation ──────────────────────────────────────────────────────────────
print("-" * 65)
print("INTERPOLATION: I(t) = (1−t)·I_SU3 + t·DISCRETE  →  TARGET")
print("-" * 65)
t_interp = (TARGET - I_SU3_analytic) / (DISCRETE - I_SU3_analytic)
print(f"  t* = (Ω/9 − 18) / (48 − 18)")
print(f"     = ({TARGET:.6f} − 18) / 30")
print(f"     = {TARGET - 18:.8f} / 30")
print(f"     = {t_interp:.10f}")
print(f"  1 − t* = {1-t_interp:.10f}")
print(f"  (1−t*) as multiple of G1: {(1-t_interp)/G1:.6f}")
print(f"  → no recognised geometric constant; (1−t*) ≈ {(432-BREATH_PERIOD)/30:.6f} = (432−Ω)/30")
print()

# ── G1 status ──────────────────────────────────────────────────────────────────
print("=" * 65)
print("G1 CLOSING STATUS")
print("=" * 65)
gap_pct = abs(I_SU3_analytic - TARGET) / TARGET * 100
print(f"  G1                            = {G1*100:.6f}%")
print(f"  |I_SU3 − TARGET| / TARGET     = {gap_pct:.4f}%")
print(f"  Ratio (orbit gap) / G1        = {gap_pct/(G1*100):.2f}×")
print()
print("CONCLUSION: The SU(3) orbit average I_SU3 = 18 undershoots")
print(f"  TARGET = {TARGET:.4f} by {gap_pct:.2f}%, which is {gap_pct/(G1*100):.1f}× G1.")
print("  The DISCRETE value 48 (P151) remains the unique approach")
print(f"  to Ω/9 within G1 = {G1*100:.4f}%.  G1 is NOT closed by orbit averaging.")
print()
print("STRUCTURAL NOTE:")
print(f"  I_SU3 × dim(su(3)) / |Φ⁺(A₂)| = 18 × 8/3 = {18*8/3:.1f} = DISCRETE")
print(f"  Equivalently: DISCRETE = B_su3(H_s,H_s)² / |Φ⁺(A₂)| = 144/3 = 48")
print("  The discrete value is recoverable from the Schur second moment")
print("  by dividing by the number of positive roots of A₂ — a new")
print("  representation of the P151 value, but not a derivation of G1.")

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