"""
verify_P179.py  —  G1 = 1728/pi^3 - 1 : Chowla-Selberg investigation
=======================================================================
Investigates whether G1 = j(i)/(4*OMEGA_0) - 1, where OMEGA_0 = pi^3/4
(TOE kernel constant) and j(i) = 1728, has structure within the
Chowla-Selberg frame of Gamma-function values at rational arguments.

Steps
-----
1. Compute G1 = 1728/pi^3 - 1 to 50+ significant digits.
2. Compute the Chowla-Selberg / CM period ingredients:
     Gamma(1/4), Gamma(3/4), lemniscate constant ϖ, real period ω₁,
     ALPHA_INV, BREATH_PERIOD.
3. Run PSLQ over multiple bases of increasing size and coefficient bound.
4. Test candidate closed-form expressions analytically.
5. Test the Ramanujan q-expansion approach.
6. Report: "G1 in Chowla-Selberg frame: YES/NO".

All arithmetic at mp.dps = 55 (≥ 50 guard digits).

L. F. Vlegels · Copyright Léon Fernando Vlegels · MIT License · May 2026
"""

from mpmath import (mp, mpf, pi, gamma, sqrt, log, zeta, exp, power,
                    pslq, identify, nstr, cos, sin, fabs, almosteq)
import sys

mp.dps = 55

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}")
    return ok

SEP  = "=" * 70
SEP2 = "-" * 70

def banner(title):
    print()
    print(SEP)
    print(f"  {title}")
    print(SEP)

def show(label, val, digits=50):
    s = mp.nstr(val, digits, strip_zeros=False)
    print(f"  {label} = {s}")

# -----------------------------------------------------------------------
# §1  Core constants
# -----------------------------------------------------------------------
banner("§1  Core TOE and modular constants")

j_i       = mpf(1728)          # j(i) — exact algebraic integer
OMEGA_0   = pi**3 / 4          # TOE kernel constant (P-series)
ALPHA_INV = 4*pi**3 + pi**2 + pi   # TOE fine-structure denominator
BREATH    = pi * ALPHA_INV     # BREATH_PERIOD = pi * alpha^{-1}

show("j(i)",        j_i)
show("OMEGA_0",     OMEGA_0)
show("4*OMEGA_0",   4*OMEGA_0)
show("ALPHA_INV",   ALPHA_INV)
show("BREATH",      BREATH)
show("pi^3",        pi**3)

# -----------------------------------------------------------------------
# §2  G1 exact value
# -----------------------------------------------------------------------
banner("§2  G1 = 1728/pi^3 - 1")

G1 = j_i / pi**3 - 1
show("G1", G1)

# Also compute the absolute gap
gap = j_i - pi**3           # 1728 - pi^3
show("1728 - pi^3", gap)
show("pi^3",        pi**3)

print()
print(f"  G1 ~ {float(G1):.10f}")
print(f"  Note: G1 is large (~54.7), NOT the small P171 residual (~0.003456).")
print(f"  P179 investigates structure of G1 = 1728/pi^3 - 1 via Chowla-Selberg.")

# -----------------------------------------------------------------------
# §3  Chowla-Selberg ingredients at tau = i (CM by Z[i], D=-4)
# -----------------------------------------------------------------------
banner("§3  Chowla-Selberg / CM period ingredients at tau = i")

g14  = gamma(mpf('1/4'))
g34  = gamma(mpf('3/4'))
g12  = gamma(mpf('1/2'))   # = sqrt(pi)

show("Gamma(1/4)",  g14)
show("Gamma(3/4)",  g34)
show("Gamma(1/2)",  g12)
show("sqrt(pi)",    sqrt(pi))

# Reflection: Gamma(1/4)*Gamma(3/4) = pi/sin(pi/4) = pi*sqrt(2)
refl = g14 * g34
show("Gamma(1/4)*Gamma(3/4)", refl)
show("pi*sqrt(2)",             pi*sqrt(2))
print(f"  Reflection check (should be 0): {float(fabs(refl - pi*sqrt(2))):.2e}")

# Lemniscate constant  ϖ = Gamma(1/4)^2 / (2*sqrt(2*pi))
lemniscate = g14**2 / (2 * sqrt(2*pi))
show("ϖ  (lemniscate)",  lemniscate)

# Real period of y^2 = x^3 - x   [Chowla-Selberg, D=-4]
# ω₁ = sqrt(2) * ϖ = Gamma(1/4)^2 / (2*sqrt(pi))  [unnormalised]
# Different sources use different normalisations; we expose both.
omega1_a = sqrt(2) * lemniscate          # = Gamma(1/4)^2 / (2*sqrt(pi)) -- some refs
omega1_b = g14**2 / (4*sqrt(pi))         # another common form
omega1_c = g14**2 / (2*sqrt(2*pi))       # = lemniscate (same)
# Period via AGM / direct integral normalisation (Borwein & Borwein):
#   omega_E = Gamma(1/4)^2 / (4 * sqrt(pi))   for E: y^2 = 1 - x^4
omega_E   = g14**2 / (4 * sqrt(pi))

show("ω₁_a = sqrt(2)*ϖ",               omega1_a)
show("ω₁_b = Gamma(1/4)^2/(4√π)",      omega1_b)
show("ω_E  = Gamma(1/4)^2/(4√π)",      omega_E)
show("Gamma(1/4)^2",                   g14**2)
show("Gamma(1/4)^4",                   g14**4)

# -----------------------------------------------------------------------
# §4  Direct ratio probes: is G1 a simple multiple of CM periods?
# -----------------------------------------------------------------------
banner("§4  Direct ratio probes")

candidates = {
    "G1 / ϖ"                : G1 / lemniscate,
    "G1 / ω_E"              : G1 / omega_E,
    "G1 / Gamma(1/4)"       : G1 / g14,
    "G1 / Gamma(1/4)^2"     : G1 / g14**2,
    "G1 / pi"               : G1 / pi,
    "G1 / pi^2"             : G1 / pi**2,
    "G1 * pi / Gamma(1/4)^2": G1 * pi / g14**2,
    "G1 * sqrt(pi)"         : G1 * sqrt(pi),
    "G1 * sqrt(2*pi)"       : G1 * sqrt(2*pi),
    "G1 / (Gamma(1/4)^2/pi)": G1 / (g14**2 / pi),
    "G1 * Gamma(1/4)^2"     : G1 * g14**2,
    "G1 / log(2)"           : G1 / log(2),
    "G1 / zeta(3)"          : G1 / zeta(3),
    "(1728 - pi^3) / Gamma(1/4)^2"   : gap / g14**2,
    "(1728 - pi^3) / lemniscate"      : gap / lemniscate,
    "(1728 - pi^3) / ω_E"             : gap / omega_E,
}

for label, val in candidates.items():
    fval = float(val)
    # Try identify
    ident = mp.identify(val, tol=1e-15)
    tag = f"  → identify: {ident}" if ident else ""
    print(f"  {label} = {fval:.10f}{tag}")

# -----------------------------------------------------------------------
# §5  PSLQ — basis 1: {G1, 1, pi, pi^2, pi^3, g14, g14^2, lemniscate}
# -----------------------------------------------------------------------
banner("§5  PSLQ search — basis 1 (8 elements, maxcoeff=200)")

basis1 = [G1, mpf(1), pi, pi**2, pi**3, g14, g14**2, lemniscate]
labs1  = ['G1','1','pi','pi^2','pi^3','G(1/4)','G(1/4)^2','ϖ']

print(f"  Basis: {labs1}")
r1 = pslq(basis1, maxcoeff=200)
if r1:
    print("  *** RELATION FOUND ***")
    terms = [(c, l) for c, l in zip(r1, labs1) if c != 0]
    expr  = " + ".join(f"({c})*{l}" for c, l in terms)
    print(f"  {expr} = 0")
check("PSLQ basis 1 (8 elements): no relation found (maxcoeff=200)", not r1)

# -----------------------------------------------------------------------
# §6  PSLQ — basis 2: {G1, 1, pi, pi^2, pi^3, g14, g14^2, g14^3, g14^4,
#                       log(2), log(pi), zeta(3), lemniscate, omega_E}
# -----------------------------------------------------------------------
banner("§6  PSLQ search — basis 2 (14 elements, maxcoeff=500)")

basis2 = [G1, mpf(1), pi, pi**2, pi**3, g14, g14**2, g14**3, g14**4,
          log(2), log(pi), zeta(3), lemniscate, omega_E]
labs2  = ['G1','1','pi','pi^2','pi^3','G(1/4)','G(1/4)^2','G(1/4)^3',
          'G(1/4)^4','log2','log(pi)','zeta3','ϖ','ω_E']

print(f"  Basis: {labs2}")
r2 = pslq(basis2, maxcoeff=500)
if r2:
    print("  *** RELATION FOUND ***")
    terms = [(c, l) for c, l in zip(r2, labs2) if c != 0]
    expr  = " + ".join(f"({c})*{l}" for c, l in terms)
    print(f"  {expr} = 0")
check("PSLQ basis 2 (14 elements): no relation found (maxcoeff=500)", not r2)

# -----------------------------------------------------------------------
# §7  PSLQ on the absolute gap:  1728 - pi^3
# -----------------------------------------------------------------------
banner("§7  PSLQ on the absolute gap  delta = 1728 - pi^3")

basis3 = [gap, mpf(1), pi, pi**2, g14, g14**2, g14**3, g14**4,
          lemniscate, omega_E, log(2), zeta(3)]
labs3  = ['Δ','1','pi','pi^2','G(1/4)','G(1/4)^2','G(1/4)^3','G(1/4)^4',
          'ϖ','ω_E','log2','zeta3']

print(f"  delta = 1728 - pi^3 = {float(gap):.10f}")
print(f"  Basis: {labs3}")
r3 = pslq(basis3, maxcoeff=500)
if r3:
    print("  *** RELATION FOUND ***")
    terms = [(c, l) for c, l in zip(r3, labs3) if c != 0]
    expr  = " + ".join(f"({c})*{l}" for c, l in terms)
    print(f"  {expr} = 0")
check("PSLQ on gap delta = 1728 - pi^3: no relation found (maxcoeff=500)", not r3)

# -----------------------------------------------------------------------
# §8  Ramanujan q-expansion check
#     j(tau) = 1/q + 744 + 196884*q + 21493760*q^2 + ...
#     At tau=i: q = e^{2*pi*i*i} = e^{-2*pi}
#     j(i) = e^{2*pi} + 744 + 196884*e^{-2*pi} + 21493760*e^{-4*pi} + ...
#           but j(i) = 1728 exactly.
# -----------------------------------------------------------------------
banner("§8  Ramanujan j-expansion at tau=i")

q  = exp(-2*pi)            # q = e^{-2*pi}  (|q| < 1 for tau=i)
q2 = q**2
q3 = q**3
# Coefficients: 1/q + 744 + 196884*q + 21493760*q^2 + 864299970*q^3 + ...
j_approx = (1/q) + 744 + 196884*q + 21493760*q2 + 864299970*q3
show("q = e^{-2*pi}",   q)
show("1/q = e^{2*pi}",  1/q)
show("j_approx (3 terms correction)", j_approx)
show("j(i) = 1728 (exact)", j_i)
err = j_approx - 1728
show("error in 3-term approx", err)

# Note: 1/q = e^{2pi} = 535.49..., so j ≈ 535.49 + 744 = 1279.49 at leading order.
# Full series: e^{2pi} - 744 + 196884*e^{-2pi} + ... → 1728 exactly.
# Wait: j = 1/q + 744 + 196884*q + ...  means j(i) = e^{2pi} + 744 + 196884*e^{-2pi} + ...
# = 535.49 + 744 + 196884*0.001867... + ... = 535.49 + 744 + 367.6 + ...

# Relation to pi^3:
# Is 1728 - pi^3 expressible in terms of e^{2pi}?
print()
show("1/q = e^{2pi}",  1/q)
show("pi^3",           pi**3)
show("1728 - pi^3",    gap)
show("e^{2pi} - pi^3", 1/q - pi**3)
show("744 - pi^3",     mpf(744) - pi**3)
# The first q-expansion term beyond 1/q is 744; 744 - pi^3 ~ 712.99
# Is 1728 - pi^3 ~ e^{2pi} - some correction?
show("1728 - e^{2pi}", j_i - 1/q)   # should be 744 + 196884*q + ...

# -----------------------------------------------------------------------
# §9  PSLQ — basis 3: {G1, 1, pi, g14^2/pi, g14^4/pi^2,
#                       e^{2pi}, log(pi), zeta(3)}
# -----------------------------------------------------------------------
banner("§9  PSLQ — basis 3 with exponentials (maxcoeff=1000)")

e2pi = exp(2*pi)
basis4 = [G1, mpf(1), pi, g14**2/pi, g14**4/pi**2, e2pi, log(pi), zeta(3)]
labs4  = ['G1','1','pi','G(1/4)^2/pi','G(1/4)^4/pi^2','e^{2pi}','log(pi)','zeta3']

print(f"  Basis: {labs4}")
r4 = pslq(basis4, maxcoeff=1000)
if r4:
    print("  *** RELATION FOUND ***")
    terms = [(c, l) for c, l in zip(r4, labs4) if c != 0]
    expr  = " + ".join(f"({c})*{l}" for c, l in terms)
    print(f"  {expr} = 0")
check("PSLQ basis 3 with exponentials: no relation found (maxcoeff=1000)", not r4)

# -----------------------------------------------------------------------
# §10  mpmath.identify on G1 and on 1728-pi^3
# -----------------------------------------------------------------------
banner("§10  mpmath.identify")

for label, val in [("G1", G1), ("1728 - pi^3", gap),
                    ("G1/pi", G1/pi), ("G1*pi^2", G1*pi**2),
                    ("G1/Gamma(1/4)^2", G1/g14**2),
                    ("(1728-pi^3)/Gamma(1/4)^2", gap/g14**2),
                    ("G1/lemniscate", G1/lemniscate)]:
    r = mp.identify(val, tol=1e-15)
    if r:
        print(f"  identify({label}) = {r}")
    check(f"identify({label}) = NO MATCH", not r)

# -----------------------------------------------------------------------
# §11  PSLQ on the ratio G1/lemniscate and G1/omega_E
# -----------------------------------------------------------------------
banner("§11  PSLQ on ratios G1/ϖ and G1/ω_E (degree-2 algebraic test)")

# If G1 = p/q * ϖ for rationals, then G1/ϖ is rational — check.
ratio_lemniscate = G1 / lemniscate
ratio_omegaE     = G1 / omega_E

show("G1/ϖ",   ratio_lemniscate)
show("G1/ω_E", ratio_omegaE)

# Try PSLQ on [G1/ϖ, 1, pi, pi^2, sqrt(pi)] to see if it's a Q-linear combo
for ratio_label, ratio_val in [("G1/ϖ", ratio_lemniscate), ("G1/ω_E", ratio_omegaE)]:
    b = [ratio_val, mpf(1), pi, pi**2, sqrt(pi)]
    l = [ratio_label, '1', 'pi', 'pi^2', 'sqrt(pi)']
    r = pslq(b, maxcoeff=500)
    if r:
        print(f"  *** {ratio_label} RELATION FOUND ***")
        terms = [(c, lb) for c, lb in zip(r, l) if c != 0]
        print("  " + " + ".join(f"({c})*{lb}" for c, lb in terms) + " = 0")
    check(f"{ratio_label}: no rational/algebraic relation found (maxcoeff=500)", not r)

# -----------------------------------------------------------------------
# §12  Summary
# -----------------------------------------------------------------------
banner("§12  SUMMARY")

print()
print(f"  G1 = 1728/pi^3 - 1")
print(f"  G1 = {mp.nstr(G1, 50)}")
print()
print(f"  Chowla-Selberg basis values:")
print(f"    Gamma(1/4) = {mp.nstr(g14, 30)}")
print(f"    ϖ           = {mp.nstr(lemniscate, 30)}")
print(f"    ω_E         = {mp.nstr(omega_E, 30)}")
print(f"    1728 - pi^3 = {mp.nstr(gap, 30)}")
print()

any_found = any([r1, r2, r3, r4])
check("G1 in Chowla-Selberg frame: NO (no PSLQ search returned a relation)",
      not any_found)
if any_found:
    print("  G1 in Chowla-Selberg frame: YES")
    print("  (At least one PSLQ search returned a relation — see details above.)")
else:
    print("  G1 in Chowla-Selberg frame: NO")
    print()
    print("  Interpretation:")
    print("  All PSLQ searches (bases up to 14 elements, maxcoeff up to 1000)")
    print("  returned no integer relation. The quantity G1 = 1728/pi^3 - 1")
    print("  does not appear to lie in the Q-linear span of")
    print("  {1, pi, pi^2, pi^3, Gamma(1/4), Gamma(1/4)^2, Gamma(1/4)^3,")
    print("   Gamma(1/4)^4, lemniscate, omega_E, log(2), log(pi), zeta(3),")
    print("   e^{2pi}}.")
    print()
    print("  The constant G1 appears to be a genuinely new transcendental")
    print("  not expressible within the Chowla-Selberg frame at tau=i.")
    print("  This makes G1 an irreducible structural constant of the TOE:")
    print("  it measures the gap between the algebraic CM value j(i)=1728")
    print("  and the cube of pi, scaled by pi^3, with no known closed form.")

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