#!/usr/bin/env python3
"""
verify_P172.py — Numerical verification for Addendum 172:
    The CM Uniqueness Theorem

Proves that tau=i (D=-4) is the unique imaginary quadratic CM point where
j(tau_D) is a power of J_short = 12.

All floating-point computations use mpmath at mp.dps=55 (55 decimal digits).
All integer assertions use exact Python integers.

Assertions:
  1.  All nine Stark-Heegner j-values as exact integers (from the classical list)
  2.  Powers of 12 in range: 12^k for k=1..17 bracketing max |j|
  3.  For each of the nine j-values: is it ±12^k for any k in 1..100?
      Assert that only j(D=-4) = 1728 passes; all others fail.
  4.  1728 = 12^3 = J_short^3 (exact)
  5.  -3375 = -15^3        (not a power of 12: 3375 = 3^3 * 5^3 has factor 5)
  6.  8000 = 20^3          (not a power of 12: 8000 = 2^6 * 5^3 has factor 5)
  7.  -32768 = -2^15       (not a power of 12: 2^15 has no factor 3)
  8.  -884736 = -96^3      (not a power of 12: exponent 2k=15 non-integer)
  9.  -884736000 = -960^3  (not a power of 12: divisible by 5^3)
  10. -147197952000 = -5280^3  (not a power of 12: divisible by 5^3 and 11^3)
  11. -262537412640768000 = -640320^3  (not a power of 12: div. by 5^3, 23^3, 29^3)
  12. Numerical: j(i) = 1728 via q-expansion to at least 50 significant figures
  13. The set {j-values} ∩ {±12^k : k >= 1} = {1728}  (uniqueness)
  14. Bonus: e^{pi*sqrt(163)} is within 10^{-10} of -j(D=-163) + 744 = 262537412640768744
             (the Ramanujan near-integer)

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

import sys
from mpmath import (
    mp, mpf, exp, pi, sqrt, log, log10, fabs, nstr, power
)

mp.dps = 55   # 55 decimal digits throughout

# ── TOE constant ──────────────────────────────────────────────────────────────

J_SHORT = 12

# ── The nine Stark-Heegner j-values (exact rational integers) ─────────────────

# Sources: classical CM theory (Weber, Silverman, Zagier).
# Each j-value is exact — no approximation.
STARK_HEEGNER = {
    -3:   0,
    -4:   1728,
    -7:   -3375,
    -8:   8000,
    -11:  -32768,
    -19:  -884736,
    -43:  -884736000,
    -67:  -147197952000,
    -163: -262537412640768000,
}

# ── q-expansion of j(i) ───────────────────────────────────────────────────────

def j_at_i_qexpansion(nterms=500):
    """
    Compute j(i) via the q-expansion at tau = i.

    tau = i  =>  q = exp(2*pi*i*tau) = exp(2*pi*i*i) = exp(-2*pi)  (real, positive)

    j(tau) = E4(tau)^3 / Delta(tau), where:
        E4(tau)    = 1 + 240 * sum_{n>=1} sigma_3(n) * q^n
        Delta(tau) = q * prod_{n>=1} (1 - q^n)^24

    Returns a real mpmath float.
    """
    # q = exp(-2*pi) at tau = i
    q = exp(-2 * pi)

    # sigma_3(n) = sum of cubes of divisors of n
    def sigma3(n):
        return sum(d**3 for d in range(1, n+1) if n % d == 0)

    # E4 via q-expansion
    E4 = mpf(1)
    for n in range(1, nterms + 1):
        E4 += 240 * sigma3(n) * q**n

    # Delta via eta^24: Delta = q * prod_{n>=1}(1 - q^n)^24
    log_prod = mpf(0)
    for n in range(1, nterms + 1):
        log_prod += 24 * log(1 - q**n)
    Delta = q * exp(log_prod)

    return E4**3 / Delta


# ── Helpers ───────────────────────────────────────────────────────────────────

PASS_COUNT = 0
FAIL_COUNT = 0


def check(label, condition, detail=""):
    global PASS_COUNT, FAIL_COUNT
    if condition:
        PASS_COUNT += 1
    else:
        FAIL_COUNT += 1
    n = PASS_COUNT + FAIL_COUNT
    print(f"  [{'PASS' if condition else 'FAIL'}] {n:>2}. {label}")
    if detail:
        print(f"             {detail}")
    return condition


def is_power_of_12(n, max_k=100):
    """
    Check whether n == ±12^k for some k in 1..max_k.
    Returns (True, sign, k) or (False, None, None).
    n must be an integer.
    """
    for sign in (1, -1):
        val = sign * n
        if val <= 0:
            continue
        pw = 12
        for k in range(1, max_k + 1):
            if pw == val:
                return True, sign, k
            if pw > val:
                break
            pw *= 12
    return False, None, None


def prime_factors(n):
    """Return set of prime factors of |n|."""
    n = abs(n)
    factors = set()
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors.add(d)
            n //= d
        d += 1
    if n > 1:
        factors.add(n)
    return factors


# ═══════════════════════════════════════════════════════════════════════════════
print("=" * 72)
print("verify_P172.py — The CM Uniqueness Theorem")
print(f"  J_short = {J_SHORT}")
print(f"  mpmath precision: mp.dps = {mp.dps}")
print("=" * 72)


# ────────────────────────────────────────────────────────────────────────────
# Assertion 1: The nine Stark-Heegner j-values are correct integers
# ────────────────────────────────────────────────────────────────────────────
print("\n── Assertion 1: Stark-Heegner j-values (exact integer check) ───────────")

# Cross-check auxiliary cube expressions for each j-value
CUBE_CHECKS = {
    -3:   (None, None),          # j = 0 (no cube form)
    -4:   (12, 3),               # 1728 = 12^3
    -7:   (-15, 3),              # -3375 = -15^3
    -8:   (20, 3),               # 8000 = 20^3
    -11:  (-32, None),           # -32768 = -2^15 (not a simple cube; use 2^15 check)
    -19:  (-96, 3),              # -884736 = -96^3
    -43:  (-960, 3),             # -884736000 = -960^3
    -67:  (-5280, 3),            # -147197952000 = -5280^3
    -163: (-640320, 3),          # -262537412640768000 = -640320^3
}

for D, j in sorted(STARK_HEEGNER.items()):
    base, exp_val = CUBE_CHECKS[D]
    if base is None:
        check(f"D={D:4d}: j = {j}", j == 0, f"j = 0 (E4 vanishes at rho)")
    elif exp_val == 3:
        expected = base**3
        check(f"D={D:4d}: j = {j} = {base}^3",
              j == expected,
              f"{base}^3 = {expected}")
    else:
        # D = -11: j = -2^15
        check(f"D={D:4d}: j = {j} = -2^15",
              j == -(2**15),
              f"-2^15 = {-(2**15)}")


# ────────────────────────────────────────────────────────────────────────────
# Assertion 2: Powers of 12 bracketing the largest |j|
# ────────────────────────────────────────────────────────────────────────────
print("\n── Assertion 2: Powers-of-12 table up to 12^17 ─────────────────────────")
max_abs_j = max(abs(j) for j in STARK_HEEGNER.values())
print(f"  max |j| = {max_abs_j}  (D=-163)")

pw = 1
k_below = None
k_above = None
for k in range(1, 20):
    pw *= 12
    if pw <= max_abs_j:
        k_below = k
    if pw > max_abs_j and k_above is None:
        k_above = k
        break

check(f"12^{k_below} <= max|j| < 12^{k_above}",
      12**k_below <= max_abs_j < 12**k_above,
      f"12^{k_below} = {12**k_below}, max|j| = {max_abs_j}, 12^{k_above} = {12**k_above}")
check("max k to check is at most 17",
      k_above <= 17,
      f"k_above = {k_above}")


# ────────────────────────────────────────────────────────────────────────────
# Assertion 3: Only D=-4 gives j = ±12^k  (main uniqueness check)
# ────────────────────────────────────────────────────────────────────────────
print("\n── Assertion 3: Uniqueness check — only D=-4 gives j = ±12^k ───────────")
power_of_12_hits = []
for D, j in sorted(STARK_HEEGNER.items()):
    hit, sign, k = is_power_of_12(j)
    if hit:
        power_of_12_hits.append((D, j, sign, k))
        check(f"D={D:4d}: j={j} IS ±12^k  →  k={k}, sign={sign}",
              True, f"{sign}*12^{k} = {sign * 12**k}")
    else:
        check(f"D={D:4d}: j={j:25d} is NOT ±12^k",
              True, "excluded")

check("Exactly one hit: D=-4",
      len(power_of_12_hits) == 1 and power_of_12_hits[0][0] == -4,
      f"hits = {power_of_12_hits}")


# ────────────────────────────────────────────────────────────────────────────
# Assertions 4–11: Individual factorisation checks (per-discriminant)
# ────────────────────────────────────────────────────────────────────────────
print("\n── Assertions 4–11: Individual factorisation checks ─────────────────────")

# 4. 1728 = 12^3 = J_short^3
check("1728 = 12^3 = J_short^3 (exact)",
      STARK_HEEGNER[-4] == J_SHORT**3,
      f"J_short^3 = {J_SHORT**3}")

# 5. -3375 = -15^3; 3375 = 3^3 * 5^3, so factor 5 disqualifies it
j7 = STARK_HEEGNER[-7]
check("-3375 = -15^3  (factored)",
      j7 == -(15**3), f"-15^3 = {-(15**3)}")
check("-3375: prime factors of |j| include 5 → not a power of 12",
      5 in prime_factors(j7),
      f"prime factors of 3375: {sorted(prime_factors(j7))}")

# 6. 8000 = 20^3; 8000 = 2^6 * 5^3, factor 5 disqualifies it
j8 = STARK_HEEGNER[-8]
check("8000 = 20^3  (factored)",
      j8 == 20**3, f"20^3 = {20**3}")
check("8000: prime factors include 5 → not a power of 12",
      5 in prime_factors(j8),
      f"prime factors of 8000: {sorted(prime_factors(j8))}")

# 7. -32768 = -2^15; no factor of 3, so not a power of 12
j11 = STARK_HEEGNER[-11]
check("-32768 = -2^15  (exact)",
      j11 == -(2**15), f"-2^15 = {-(2**15)}")
check("-32768: prime factors of |j| = {2} (no factor 3) → not a power of 12",
      prime_factors(j11) == {2},
      f"prime factors of 32768: {sorted(prime_factors(j11))}")

# 8. -884736 = -96^3; 96 = 2^5 * 3, so 96^3 = 2^15 * 3^3;
#    for 12^k = 2^{2k}*3^k = 2^15*3^3 we'd need 2k=15 (non-integer)
j19 = STARK_HEEGNER[-19]
check("-884736 = -96^3  (factored)",
      j19 == -(96**3), f"-96^3 = {-(96**3)}")
# 96^3 = 2^15 * 3^3; 12^k = 2^{2k}*3^k; needs 2k=15 → k=7.5 (not integer)
check("-884736: for 12^k = 884736 need 2k=15 → non-integer k; not a power of 12",
      884736 == 2**15 * 3**3,
      f"884736 = 2^15 * 3^3 = {2**15 * 3**3}")

# 9. -884736000 = -960^3; 960 = 2^6 * 3 * 5, so 960^3 = 2^18*3^3*5^3
j43 = STARK_HEEGNER[-43]
check("-884736000 = -960^3  (factored)",
      j43 == -(960**3), f"-960^3 = {-(960**3)}")
check("-884736000: prime factors include 5 → not a power of 12",
      5 in prime_factors(j43),
      f"prime factors of 960: {sorted(prime_factors(960))}")

# 10. -147197952000 = -5280^3; 5280 = 2^5 * 3 * 5 * 11, so includes 5, 11
j67 = STARK_HEEGNER[-67]
check("-147197952000 = -5280^3  (factored)",
      j67 == -(5280**3), f"-5280^3 = {-(5280**3)}")
check("-147197952000: prime factors include 5 and 11 → not a power of 12",
      5 in prime_factors(j67) and 11 in prime_factors(j67),
      f"prime factors of 5280: {sorted(prime_factors(5280))}")

# 11. -262537412640768000 = -640320^3; 640320 = 2^6 * 3 * 5 * 23 * 29
j163 = STARK_HEEGNER[-163]
check("-262537412640768000 = -640320^3  (factored)",
      j163 == -(640320**3), f"-640320^3 = {-(640320**3)}")
check("-262537412640768000: prime factors include 5, 23, 29 → not a power of 12",
      all(p in prime_factors(j163) for p in (5, 23, 29)),
      f"prime factors of 640320: {sorted(prime_factors(640320))}")


# ────────────────────────────────────────────────────────────────────────────
# Assertion 12: Numerical j(i) = 1728 to ≥50 significant figures
# ────────────────────────────────────────────────────────────────────────────
print("\n── Assertion 12: Numerical j(i) = 1728 via q-expansion (50 sig figs) ───")
j_num = j_at_i_qexpansion(nterms=500)
abs_err = fabs(j_num - mpf(1728))
sig_figs = int(-log10(abs_err / mpf(1728) + mpf('1e-200')))
print(f"  j(i) computed  = {nstr(j_num, 25)}")
print(f"  absolute error = {nstr(abs_err, 5)}")
print(f"  significant figures ≥ {sig_figs}")
check("j(i) = 1728 to at least 50 significant figures",
      sig_figs >= 50,
      f"sig figs ≈ {sig_figs}")


# ────────────────────────────────────────────────────────────────────────────
# Assertion 13: Set intersection {j-values} ∩ {±12^k : k≥1} = {1728}
# ────────────────────────────────────────────────────────────────────────────
print("\n── Assertion 13: {j-values} ∩ {±12^k : k≥1} = {1728} ──────────────────")
powers_of_12 = set()
pw = 12
for k in range(1, 101):
    powers_of_12.add(pw)
    powers_of_12.add(-pw)
    pw *= 12

j_values = set(STARK_HEEGNER.values())
intersection = j_values & powers_of_12
check("{j-values} ∩ {±12^k : k=1..100} = {1728}",
      intersection == {1728},
      f"intersection = {intersection}")


# ────────────────────────────────────────────────────────────────────────────
# Assertion 14: Ramanujan near-integer — e^{pi*sqrt(163)} ≈ -j(D=-163) + 744
# ────────────────────────────────────────────────────────────────────────────
print("\n── Assertion 14: Ramanujan near-integer e^{π√163} ──────────────────────")

exp_pi_sqrt_163 = exp(pi * sqrt(mpf(163)))
# Exact prediction from CM theory: -j(D=-163) + 744 = 262537412640768000 + 744
ramanujan_integer = mpf(-STARK_HEEGNER[-163]) + 744  # = 262537412640768744
error = fabs(exp_pi_sqrt_163 - ramanujan_integer)

print(f"  e^{{pi*sqrt(163)}}  = {nstr(exp_pi_sqrt_163, 30)}")
print(f"  -j(D=-163) + 744   = {int(ramanujan_integer)}  (exact integer)")
print(f"  error              = {nstr(error, 5)}")

check("e^{pi*sqrt(163)} is within 10^{-10} of -j(D=-163)+744",
      error < mpf('1e-10'),
      f"|error| = {nstr(error, 4)}")

check("|j(D=-163)| = 262537412640768000 (exact)",
      abs(STARK_HEEGNER[-163]) == 262537412640768000,
      f"value = {abs(STARK_HEEGNER[-163])}")

# The near-integer: e^{pi*sqrt(163)} is within 10^{-12} of 262537412640768744
near_integer_value = 262537412640768744
near_int_error = fabs(exp_pi_sqrt_163 - mpf(near_integer_value))
print(f"  |e^{{pi*sqrt(163)}} - {near_integer_value}| = {nstr(near_int_error, 5)}")
check("e^{pi*sqrt(163)} within 10^{-10} of 262537412640768744",
      near_int_error < mpf('1e-10'),
      f"|error| = {nstr(near_int_error, 5)}")


# ── Summary ───────────────────────────────────────────────────────────────────
print()
print("=" * 72)
if FAIL_COUNT == 0:
    print("Theorem confirmed:")
    print("  Among all imaginary quadratic CM points tau_D,")
    print("  the unique point where j(tau_D) = J_short^k (k>=1) is")
    print(f"  D = -4 (tau = i), giving j(i) = 1728 = {J_SHORT}^3 = J_short^3.")

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