# Exact check of the parameter inequalities used in Corollaries 2 and 3 (pure Python, Fractions).
# Auteur : David Niedbala Giraudin (Fino), ORCID 0009-0009-1526-1178, avec assistance IA (Claude, Anthropic).
from fractions import Fraction as F
import sys

def check(w, d):
    """w = varpi, d = delta. Returns (xi1, xi2, xi3) or raises if some inequality fails."""
    lo2 = max(F(1, 4) + 12*w + 4*d, F(2, 5))       # Theorem 1 (48w+16d+4s<1, s = 1/2 - xi2) and xi2 > 2/5
    hi2 = min(F(4, 9) - F(28, 9)*w - F(2, 9)*d, F(55, 136))  # Type III (s3 = 1/2 - xi3) and xi2 < 55/136
    assert lo2 < hi2, "no room for xi2"
    xi2 = (lo2 + hi2)/2; xi3 = xi2
    lo1 = max(F(1, 4) + 7*w + 2*d, (1 - xi2)/2)    # Type I (OpenAI Lemma 2.20) and 2xi1 + xi2 > 1
    hi1 = min((2 - 3*xi2)/2, 4 - 9*xi2)            # 2xi1 + 3xi2 < 2 and xi1 + 9xi2 < 4
    assert lo1 < hi1, "no room for xi1"
    xi1 = (lo1 + hi1)/2
    cond = {
        "72w+24d<1": 72*w + 24*d < 1,
        "68w+14d<1": 68*w + 14*d < 1,
        "Thm1: 48w+16d+4s<1 (s=1/2-xi2)": 48*w + 16*d + 4*(F(1, 2) - xi2) < 1,
        "Thm1: 48w+18d+2s<1 (s=1/2-xi2)": 48*w + 18*d + 2*(F(1, 2) - xi2) < 1,
        "Type III: 1/2-xi3 > 1/18+28w/9+2d/9": F(1, 2) - xi3 > F(1, 18) + F(28, 9)*w + F(2, 9)*d,
        "Type I: xi1 > 1/4+7w+2d": xi1 > F(1, 4) + 7*w + 2*d,
        "2xi1+3xi2<2": 2*xi1 + 3*xi2 < 2,
        "xi2<=xi3": xi2 <= xi3,
        "xi1+9xi2<4": xi1 + 9*xi2 < 4,
        "2xi1+xi2>1": 2*xi1 + xi2 > 1,
        "17xi2<7": 17*xi2 < 7,
        "xi2>4/10 (minorant not trivial)": xi2 > F(4, 10),
        "6(1-2xi2)>1": 6*(1 - 2*xi2) > 1,
        "1/2+2w < 1-xi2": F(1, 2) + 2*w < 1 - xi2,
        "level 1/2+2w < 0.53": F(1, 2) + 2*w < F(53, 100),
        "xi2 < 55/136": xi2 < F(55, 136),
    }
    bad = [k for k, v in cond.items() if not v]
    return xi1, xi2, xi3, bad

def check_corners(w, d, k=10**6):
    """Points close to the ends of the xi2 and xi1 intervals of the note (with the bounds 2/5 and 55/136)."""
    lo2 = max(F(1, 4) + 12*w + 4*d, F(2, 5)); hi2 = min(F(4, 9) - F(28, 9)*w - F(2, 9)*d, F(55, 136))
    assert lo2 < hi2
    bad = []
    for xi2 in (lo2 + (hi2 - lo2)/k, hi2 - (hi2 - lo2)/k):
        lo1 = max(F(1, 4) + 7*w + 2*d, (1 - xi2)/2); hi1 = min(1 - F(3, 2)*xi2, 4 - 9*xi2)
        assert lo1 < hi1
        for xi1 in (lo1 + (hi1 - lo1)/k, hi1 - (hi1 - lo1)/k):
            c = [72*w + 24*d < 1, 68*w + 14*d < 1, 48*w + 16*d + 4*(F(1, 2) - xi2) < 1,
                 48*w + 18*d + 2*(F(1, 2) - xi2) < 1, F(1, 2) - xi2 > F(1, 18) + F(28, 9)*w + F(2, 9)*d,
                 xi1 > F(1, 4) + 7*w + 2*d, 2*xi1 + 3*xi2 < 2, xi1 + 9*xi2 < 4, 2*xi1 + xi2 > 1, 17*xi2 < 7,
                 xi2 > F(2, 5), 6*(1 - 2*xi2) > 1, F(1, 2) + 2*w < 1 - xi2, F(1, 2) + 2*w < F(53, 100), xi2 < F(55, 136)]
            if not all(c): bad.append((float(xi1), float(xi2), [i for i, v in enumerate(c) if not v]))
    return bad

if __name__ == "__main__":
    ok = True
    for w, d in [(F(7, 544) - F(1, 10**6), F(1, 10**8)), (F(7, 544) - F(1, 10**9), F(1, 10**11)),
                 (F(1, 80), F(1, 10**6)), (F(1, 100), F(1, 10**6)), (F(1, 1000), F(1, 10**6))]:
        bad = check_corners(w, d)
        print("corners: varpi =", float(w), " delta =", float(d), " ->", "OK" if not bad else bad)
        ok &= not bad
    for w, d in [(F(7, 544) - F(1, 10**6), F(1, 10**8)), (F(7, 544) - F(1, 10**9), F(1, 10**11)), (F(253, 20000), F(1, 10**8))]:
        xi1, xi2, xi3, bad = check(w, d)
        print("varpi =", float(w), " delta =", float(d), " xi1 =", float(xi1), " xi2 = xi3 =", float(xi2), " ->", "OK" if not bad else bad)
        ok &= not bad
    print("Meeting point: 1/4+12w = 4/9-28w/9 at w =", F(7, 36)/(12 + F(28, 9)), "; xi2 =", F(1, 4) + 12*F(7, 544))
    print("Interval for xi1 at the limit: (", float(F(1, 4) + 7*F(7, 544)), ",", float(4 - 9*F(55, 136)), ")")
    base = F(544, 143)
    print("1/(1/4+7/544) = 544/143 =", float(base))
    for c1 in (F(67, 10**7), F(2, 10**5)):
        print("  with c1 <=", float(c1), ": 544/(143(1-c1)) =", float(base/(1 - c1)))
    assert base/(1 - F(67, 10**7)) < F(380423, 100000) and base/(1 - F(2, 10**5)) < F(380428, 100000)
    print("RESULTAT :", "OK" if ok else "ECHEC")