#!/usr/bin/env python3
"""Exact-integer verification of the sigma0 mass-ratio binding facts.

Replicates GoldenInt.IsPos / PosPair sign test bit-for-bit so every `decide`
in DeltaSpine/MassRatioBinding.lean is guaranteed true before Lean runs.

CODATA 2022: m_mu / m_e = 206.768 2827(46)
  -> certified outer interval [206.7682780, 206.7682874]
     R_lo = 2067682780 / 10^7,  R_hi = 2067682874 / 10^7
"""
from fractions import Fraction
import math

R_LO_NUM, R_HI_NUM, SCALE = 2067682780, 2067682874, 10**7


def phi_pow(n: int):
    """phiPow n = (a, b) meaning a + b*phi. Recurrence: (a,b) -> (b, a+b)."""
    a, b = 1, 0
    for _ in range(n):
        a, b = b, a + b
    return a, b


def pos_pair(s: int, t: int) -> bool:
    """PosPair s t: sign of s + t*sqrt(5) > 0, exactly as in GoldenInt.lean."""
    if s >= 0 and t >= 0 and (s > 0 or t > 0):
        return True
    if s < 0 and t > 0 and s * s < 5 * t * t:
        return True
    if s > 0 and t < 0 and 5 * t * t < s * s:
        return True
    return False


def is_pos(a: int, b: int) -> bool:
    """IsPos (a + b*phi) = PosPair (2a+b) b."""
    return pos_pair(2 * a + b, b)


def rat_lt(p: int, q: int, x) -> bool:
    """RatLt p q x  <=>  p/q < x. ratWitness = (q*a - p, q*b)."""
    a, b = x
    return is_pos(q * a - p, q * b)


def rat_gt(p: int, q: int, x) -> bool:
    """RatGt p q x  <=>  x < p/q. IsPos(ratWitness (-p) q (-x))."""
    a, b = x
    return is_pos(q * (-a) + p, q * (-b))


checks = []

# 1/2: nearest-rung window at first power: phi^11 < R_lo, R_hi < phi^12
checks.append(("window_lower  phi^11 < R_lo",
               rat_gt(R_LO_NUM, SCALE, phi_pow(11))))
checks.append(("window_upper  R_hi < phi^12",
               rat_lt(R_HI_NUM, SCALE, phi_pow(12))))

# 3/4: nearest rung k=11 unique: phi^21 < R_lo^2 and R_hi^2 < phi^23
checks.append(("nearest_lower phi^21 < R_lo^2",
               rat_gt(R_LO_NUM**2, SCALE**2, phi_pow(21))))
checks.append(("nearest_upper R_hi^2 < phi^23",
               rat_lt(R_HI_NUM**2, SCALE**2, phi_pow(23))))

# 5/6: deviation bracket 5/63 < eps < 7/88
#   eps > 5/63 : R_lo^63 > phi^(11*63+5) = phi^698
#   eps < 7/88 : R_hi^88 < phi^(11*88+7) = phi^975
checks.append(("eps_lower  phi^698 < R_lo^63",
               rat_gt(R_LO_NUM**63, SCALE**63, phi_pow(698))))
checks.append(("eps_upper  R_hi^88 < phi^975",
               rat_lt(R_HI_NUM**88, SCALE**88, phi_pow(975))))

for name, ok in checks:
    print(f"{'PASS' if ok else 'FAIL'}  {name}")

# --- context numerics -------------------------------------------------------
phi = (1 + math.sqrt(5)) / 2
lphi = math.log(phi)
eps_lo = math.log(R_LO_NUM / SCALE) / lphi - 11
eps_hi = math.log(R_HI_NUM / SCALE) / lphi - 11
print(f"\neps in [{eps_lo:.9f}, {eps_hi:.9f}]")
print(f"5/63  = {5/63:.9f}   7/88  = {7/88:.9f}")
print(f"1/(4pi)          = {1/(4*math.pi):.9f}")
alpha = 1 / 137.035999177
print(f"1/(4pi) - alpha^2 = {1/(4*math.pi) - alpha**2:.9f}")

# --- can we exclude 1/(4pi)-alpha^2 with a rational? ------------------------
# need m/n with 1/(4pi)-alpha^2 < m/n < eps_lo, then prove eps > m/n.
lo_band = Fraction(1) / (4 * Fraction(math.pi)) - Fraction(alpha) ** 2
hi_band = Fraction(eps_lo)


def simplest_between(lo: Fraction, hi: Fraction) -> Fraction:
    """Simplest fraction strictly inside (lo, hi) via Stern-Brocot."""
    a, b = Fraction(0), Fraction(1, 0) if False else None
    # continued-fraction walk
    lo_n, lo_d = lo.numerator, lo.denominator
    hi_n, hi_d = hi.numerator, hi.denominator

    def walk(ln, ld, hn, hd):
        # find simplest p/q with ln/ld < p/q < hn/hd, all positive
        fl = ln // ld
        if fl + 1 < Fraction(hn, hd):
            return Fraction(fl + 1)
        if ld == 0 or hd == 0:
            return Fraction(fl + 1)
        # recurse on reciprocal of fractional parts
        sub = walk(hd, hn - fl * hd, ld, ln - fl * ld)
        return fl + 1 / sub

    return walk(lo_n, lo_d, hi_n, hi_d)


try:
    f = simplest_between(lo_band, hi_band)
    n = f.denominator
    m = f.numerator
    exp = 11 * n + m
    ok = rat_gt(R_LO_NUM**n, SCALE**n, phi_pow(exp))
    print(f"\nsimplest m/n in (1/4pi-a^2, eps_lo): {m}/{n} = {float(f):.9f}")
    print(f"  requires phi^{exp} < R_lo^{n}: {'PASS' if ok else 'FAIL'}")
except Exception as e:  # noqa: BLE001
    print(f"stern-brocot failed: {e}")

# --- tightest brackets at modest denominators -------------------------------
best_lo, best_hi = Fraction(0), Fraction(1)
for n in range(1, 400):
    m_lo = math.floor(eps_lo * n)  # candidate lower bound m/n < eps_lo
    if Fraction(m_lo, n) > best_lo and m_lo / n < eps_lo:
        # verify exactly against R_lo
        if rat_gt(R_LO_NUM**n, SCALE**n, phi_pow(11 * n + m_lo)):
            best_lo = Fraction(m_lo, n)
            print(f"lower {m_lo}/{n} = {m_lo/n:.9f}  (exp {11*n+m_lo})")
    m_hi = math.ceil(eps_hi * n)  # candidate upper bound m/n > eps_hi
    if Fraction(m_hi, n) < best_hi and m_hi / n > eps_hi:
        if rat_lt(R_HI_NUM**n, SCALE**n, phi_pow(11 * n + m_hi)):
            best_hi = Fraction(m_hi, n)
            print(f"upper {m_hi}/{n} = {m_hi/n:.9f}  (exp {11*n+m_hi})")