"""verify_P165.py — Verification for Addendum 165: The Geometric Identity of B_{F4} = 36.

Checks:
  1. G2 root counts: 6 short (squared length 2), 6 long (squared length 6), total 12
  2. 36 = 6 x 6 = short x long root pair count  (Theorem 1)
  3. 36 = 12 x 3 = |Phi_G2| x |Z/3Z|            (Theorem 2)
  4. Angular check: 36 x 10.0 == 360.0
  5. Angular sub-sector: 360/12 = 30 deg per root direction; 30/3 = 10 deg per monodromy
  6. j(i) = B_G2 x B_F4 = 48 x 36 = 1728        (numerical, B_G2 = 48 from P164)
  7. j(i) = (4 x J_short) x (3 x J_short) with J_short = 12
"""

import sys

from mpmath import mp, mpf, pi, nstr, fabs

mp.dps = 55  # 55 digits for 50 sig figs of safety

PASS = FAIL = 0
N = 0

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

# ── TOE constants ─────────────────────────────────────────────────────────────
ALPHA_INV = 4*pi**3 + pi**2 + pi        # ≈ 137.036
J_short   = mpf(12)
j_i       = J_short**3                  # = 1728
B_G2      = mpf(48)                     # proved P164: |Phi_F4|
B_F4      = mpf(36)                     # this paper

# ── G2 root system in the hyperplane e1+e2+e3=0 inside R^3 ───────────────────
# Short roots (6): ±(e_i - e_j) for i<j — squared length 2
# Long roots  (6): ±(2e_i - e_j - e_k), {i,j,k}={1,2,3} — squared length 6

def sq_len(v):
    return sum(x*x for x in v)

# Build short roots: ±(e_i - e_j) for {i,j} in {{1,2},{1,3},{2,3}}
short_roots = []
pairs = [(0,1), (0,2), (1,2)]   # index pairs for e_i - e_j
for (i, j) in pairs:
    for sign in (+1, -1):
        v = [0, 0, 0]
        v[i] = sign
        v[j] = -sign
        short_roots.append(tuple(v))

# Build long roots: ±(2e_i - e_j - e_k) for each i, {j,k} = {0,1,2}\{i}
long_roots = []
for i in range(3):
    others = [k for k in range(3) if k != i]
    for sign in (+1, -1):
        v = [0, 0, 0]
        v[i]         =  2 * sign
        v[others[0]] = -1 * sign
        v[others[1]] = -1 * sign
        long_roots.append(tuple(v))


# ── 1. Root counts and squared lengths ───────────────────────────────────────
print("=== 1. G2 root counts and squared lengths ===")
print(f"  |Phi_G2_short| = {len(short_roots)}  (expected 6)")
print(f"  |Phi_G2_long|  = {len(long_roots)}   (expected 6)")
print(f"  |Phi_G2|       = {len(short_roots) + len(long_roots)}  (expected 12 = J_short)")

check("|Phi_G2_short| = 6", len(short_roots) == 6)
check("|Phi_G2_long| = 6", len(long_roots) == 6)
check("|Phi_G2| = 12 = J_short", len(short_roots) + len(long_roots) == 12)

short_sq_lens = {sq_len(v) for v in short_roots}
long_sq_lens  = {sq_len(v) for v in long_roots}
print(f"  Short root squared lengths: {short_sq_lens}  (expected {{2}})")
print(f"  Long  root squared lengths: {long_sq_lens}   (expected {{6}})")

check("short root squared lengths = {2}", short_sq_lens == {2})
check("long root squared lengths = {6}", long_sq_lens == {6})

root_length_ratio_sq = 6 / 2
print(f"  Squared length ratio = {root_length_ratio_sq}  (expected 3 = sqrt(3)^2, G2 ratio)")
check("squared length ratio = 3 (G2 ratio)", root_length_ratio_sq == 3)

# Verify all roots lie in the plane e1+e2+e3=0
all_in_plane = all(sum(v) == 0 for v in short_roots + long_roots)
print(f"  All roots in e1+e2+e3=0: {all_in_plane}")
check("all roots in e1+e2+e3=0", all_in_plane)


# ── 2. B_F4 = 36 = 6 x 6 (Theorem 1: root-pair identification) ────────────
print("\n=== 2. B_F4 = 36 = |Phi_short| x |Phi_long| (Theorem 1) ===")
n_short = len(short_roots)
n_long  = len(long_roots)
pair_count = n_short * n_long
print(f"  |Phi_short| x |Phi_long| = {n_short} x {n_long} = {pair_count}  (expected 36)")
check("|Phi_short| x |Phi_long| = 36 (Theorem 1)", pair_count == 36)
check("B_F4 = 36", int(B_F4) == 36)
print(f"  B_F4 = {int(B_F4)} == pair_count = {pair_count}  MATCH")


# ── 3. B_F4 = 36 = 12 x 3 (Theorem 2: monodromy factorisation) ────────────
print("\n=== 3. B_F4 = 36 = |Phi_G2| x |Z/3Z| (Theorem 2) ===")
n_G2_total   = len(short_roots) + len(long_roots)  # 12
n_monodromy  = 3                                    # |Z/3Z|
monodromy_factored = n_G2_total * n_monodromy
print(f"  |Phi_G2| x |Z/3Z| = {n_G2_total} x {n_monodromy} = {monodromy_factored}  (expected 36)")
check("|Phi_G2| x |Z/3Z| = 36 (Theorem 2)", monodromy_factored == 36)
check("B_F4 = monodromy factorisation", int(B_F4) == monodromy_factored)
print(f"  Both factorisations agree: 6x6 = 12x3 = 36")


# ── 4. Angular check: 36 x 10.0 = 360.0 ──────────────────────────────────
print("\n=== 4. Angular tiling: 36 x 10.0 deg = 360.0 deg ===")
deg_per_subsector  = mpf('10')
total_angle        = 36 * deg_per_subsector
print(f"  36 x 10 deg = {float(total_angle)} deg  (expected 360.0)")
check("36 x 10 deg = 360 deg", total_angle == 360)


# ── 5. Angular sub-sector structure ────────────────────────────────────────
print("\n=== 5. Angular sub-sector structure ===")
deg_per_root_dir   = mpf('360') / 12           # 30 deg per root direction
deg_per_subsector2 = deg_per_root_dir / 3       # 10 deg per monodromy subdivision
print(f"  360 / |Phi_G2| = 360 / 12 = {float(deg_per_root_dir)} deg  (expected 30)")
print(f"  30 / |Z/3Z|    = 30  / 3  = {float(deg_per_subsector2)} deg  (expected 10)")
check("360 / |Phi_G2| = 30 deg per root direction", deg_per_root_dir == 30)
check("30 / |Z/3Z| = 10 deg per monodromy sub-sector", deg_per_subsector2 == 10)

# Cross-check: 10 deg = 60 deg / 6 (intra-species period)
deg_per_short_species = mpf('360') / 6          # 60 deg between consecutive short roots
deg_per_sub_from_species = deg_per_short_species / 6
print(f"  360 / 6 = {float(deg_per_short_species)} deg (per short-species period)")
print(f"  60  / 6 = {float(deg_per_sub_from_species)} deg = 10 deg  MATCH")
check("60 / 6 = 10 deg (cross-check via short-species period)", deg_per_sub_from_species == 10)


# ── 6. j(i) = B_G2 x B_F4 = 48 x 36 = 1728 ────────────────────────────────
print("\n=== 6. j(i) = B_G2 x B_F4 = 48 x 36 = 1728 ===")
j_from_factors = B_G2 * B_F4
print(f"  B_G2 x B_F4 = {int(B_G2)} x {int(B_F4)} = {int(j_from_factors)}  (expected 1728)")
check("B_G2 x B_F4 = 48 x 36 = 1728", int(j_from_factors) == 1728)
check("j(i) = J_short^3 = 1728", int(j_from_factors) == int(j_i))
print(f"  j(i) = J_short^3 = 12^3 = {int(j_i)}  MATCH")


# ── 7. j(i) = (4 x J_short) x (3 x J_short) ──────────────────────────────
print("\n=== 7. j(i) = (rank(D4) x J_short) x (|Z/3Z| x J_short) ===")
rank_D4     = mpf(4)
order_Z3    = mpf(3)
B_G2_check  = rank_D4  * J_short    # 4 x 12 = 48
B_F4_check  = order_Z3 * J_short    # 3 x 12 = 36
j_check     = B_G2_check * B_F4_check

print(f"  rank(D4) x J_short  = {int(rank_D4)} x {int(J_short)} = {int(B_G2_check)}"
      f"  (= B_G2 = {int(B_G2)}, match: {int(B_G2_check) == int(B_G2)})")
print(f"  |Z/3Z|   x J_short  = {int(order_Z3)} x {int(J_short)} = {int(B_F4_check)}"
      f"  (= B_F4 = {int(B_F4)}, match: {int(B_F4_check) == int(B_F4)})")
print(f"  j(i) = {int(B_G2_check)} x {int(B_F4_check)} = {int(j_check)}  (expected 1728)")
check("rank(D4) x J_short = 48 = B_G2", int(B_G2_check) == 48)
check("|Z/3Z| x J_short = 36 = B_F4", int(B_F4_check) == 36)
check("j(i) = 48 x 36 = 1728", int(j_check) == 1728)

# Self-consistency: (4 x 12) x (3 x 12) = 12 x 12^2 = 12^3
Js_cubed = J_short**3
print(f"  (4x12) x (3x12) = 12 x 144 = {int(j_check)} = 12^3 = {int(Js_cubed)}"
      f"  MATCH: {int(j_check) == int(Js_cubed)}")
check("(4x12) x (3x12) = 12^3", int(j_check) == int(Js_cubed))


# ── 8. Summary table ──────────────────────────────────────────────────────
print("\n=== Summary ===")
print(f"  G2 short roots:            {len(short_roots)} (squared length 2)")
print(f"  G2 long roots:             {len(long_roots)}  (squared length 6)")
print(f"  G2 total roots = J_short:  {len(short_roots)+len(long_roots)}")
print(f"  B_F4 = 6 x 6 (root pairs): {n_short * n_long}")
print(f"  B_F4 = 12 x 3 (monodromy): {n_G2_total * n_monodromy}")
print(f"  Angular: 36 x 10 deg       = {36*10} deg = full circle")
print(f"  Angular: 30 deg / 3        = {30//3} deg per monodromy sub-sector")
print(f"  B_G2 = rank(D4) x J_short  = {int(rank_D4)} x {int(J_short)} = {int(B_G2_check)}")
print(f"  B_F4 = |Z/3Z|  x J_short  = {int(order_Z3)} x {int(J_short)} = {int(B_F4_check)}")
print(f"  j(i) = B_G2 x B_F4        = {int(B_G2)} x {int(B_F4)} = {int(j_i)}")
print(f"  j(i) = J_short^3           = {int(J_short)}^3 = {int(Js_cubed)}")
print(f"  ALPHA_INV                  = {nstr(ALPHA_INV, 12)}")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)