#!/usr/bin/env python3
"""
PROVENANCE: PROOF

PST Computation 136 — Gauge anomaly cancellation and the N_c = 3 condition
==========================================================================
Witness computation for the anomaly-cancellation consistency check
(sec:gauge-sm). With the one-generation Standard-Model hypercharge assignment
(convention Q = T_3 + Y/2)
    Y(Q_L)=+1/3, Y(u_R)=+4/3, Y(d_R)=-2/3, Y(L_L)=-1, Y(e_R)=-2,
this evaluates all six triangle-anomaly coefficients as exact functions of the
colour rank N_c and shows:
  * at N_c = 3 all six vanish exactly (the theory is gauge-consistent);
  * the [U(1)_Y]^3 coefficient equals -2 N_c + 6 and the [SU(2)]^2 U(1)_Y
    coefficient equals (N_c - 3)/6, so BOTH vanish iff N_c = 3.
Hence anomaly freedom, given the SM hypercharge pattern, forces N_c = 3, the
dimension of the M_3(C) colour factor of A_F.

This is exact rational arithmetic, hence PROVENANCE PROOF. It is a CONSISTENCY
CHECK given the SM hypercharges (which PST inherits, see sec:gauge-sm), not a
first-principles derivation of the hypercharges or of anomaly cancellation.

Left-handed Weyl content of one generation (right-handed fields entered as
left-handed conjugates: colour 3 -> 3bar, Y -> -Y):
    field  colour  weak  Y        su3(+1 for 3, -1 for 3bar, 0 singlet)
    Q_L    N_c     2     +1/3     +1
    u_R^c  N_c     1     -4/3     -1
    d_R^c  N_c     1     +2/3     -1
    L_L    1       2     -1        0
    e_R^c  1       1     +2        0
"""
from fractions import Fraction as F

T = F(1, 2)   # Dynkin index of the fundamental of SU(N)


def content(Nc):
    # (n_colour, n_weak, Y, su3_sign)
    return [
        ("Q_L", Nc, 2, F(1, 3), +1),
        ("uRc", Nc, 1, F(-4, 3), -1),
        ("dRc", Nc, 1, F(2, 3), -1),
        ("L_L", 1, 2, F(-1), 0),
        ("eRc", 1, 1, F(2), 0),
    ]


def anomalies(Nc):
    fl = content(Nc)
    grav_U1 = sum(nc * nw * Y for _, nc, nw, Y, _ in fl)                    # [grav]^2 U(1)
    U1_cubed = sum(nc * nw * Y ** 3 for _, nc, nw, Y, _ in fl)             # [U(1)]^3
    SU2sq_U1 = sum(nc * T * Y for _, nc, nw, Y, _ in fl if nw == 2)        # [SU(2)]^2 U(1)
    SU3sq_U1 = sum(nw * T * Y for _, nc, nw, Y, s in fl if s != 0)         # [SU(3)]^2 U(1)
    SU3_cubed = sum(nw * s for _, nc, nw, Y, s in fl if s != 0)            # [SU(3)]^3
    SU2_cubed = F(0)                                                       # SU(2) has no cubic anomaly
    return {
        "[grav]^2 U(1)": grav_U1,
        "[U(1)]^3": U1_cubed,
        "[SU(2)]^2 U(1)": SU2sq_U1,
        "[SU(3)]^2 U(1)": SU3sq_U1,
        "[SU(3)]^3": SU3_cubed,
        "[SU(2)]^3": SU2_cubed,
    }


print("Anomaly coefficients as functions of N_c:")
print(f"  {'N_c':>3} | " + " | ".join(f"{k:>14}" for k in anomalies(3)))
for Nc in (1, 2, 3, 4, 5):
    a = anomalies(Nc)
    print(f"  {Nc:>3} | " + " | ".join(f"{str(v):>14}" for v in a.values()))

# exact N_c-dependence of the two constraining conditions
u1_3 = [anomalies(n)["[U(1)]^3"] for n in (0, 1)]
su2_1 = [anomalies(n)["[SU(2)]^2 U(1)"] for n in (0, 1)]
slope_u1 = u1_3[1] - u1_3[0]
slope_su2 = su2_1[1] - su2_1[0]
print(f"\n[U(1)]^3(N_c)      = {slope_u1} N_c + {u1_3[0]}   (= -2 N_c + 6 ; zero iff N_c=3)")
print(f"[SU(2)]^2 U(1)(N_c) = {slope_su2} N_c + {su2_1[0]}   (= (N_c-3)/6 ; zero iff N_c=3)")

at3 = anomalies(3)
all_zero_at_3 = all(v == 0 for v in at3.values())
forces_3 = (anomalies(3)["[U(1)]^3"] == 0 and anomalies(2)["[U(1)]^3"] != 0
            and anomalies(4)["[U(1)]^3"] != 0)
print("\nAll six anomalies vanish at N_c=3:", all_zero_at_3)
print("Anomaly freedom forces N_c=3 (cubic & SU(2)^2 U(1) nonzero for N_c!=3):", forces_3)
print("\nRESULT:", "PASS" if (all_zero_at_3 and forces_3) else "FAIL")