# -*- coding: utf-8 -*-
# Verificateur exact : p_c(site, Z2) >= 0.556921  (et >= 0.55692 avec marge plus large)
# Version du 3 novembre 2026 : controle aussi les profondeurs n = 61 et n >= 62 (queue partant de [0,1]).
# Python pur, arithmetique rationnelle exacte (fractions), aucun module externe.
# Aucun argument a taper : lancer tel quel. Duree : quelques secondes.
# Auteur : David Niedbala Giraudin (Fino), ORCID 0009-0009-1526-1178, avec assistance IA (Claude, Anthropic).
from fractions import Fraction as Fr
import itertools, math, sys

def say(*a):
    print(*a, flush=True)

D = 10 ** 40
def dn(x): return Fr(math.floor(x * D), D)
def upr(x): return Fr(math.ceil(x * D), D)

def lfun(X, Y): return (X + 1) / (X + Y + 1)
def ufun(X, Y): return (Y + 1) / (X + Y + 1)

def law(rho, probs, n):
    out = {}
    for x in itertools.product((0, 1), repeat=n):
        tot = Fr(0)
        for s, r in rho.items():
            w = r
            for bit, pr in zip(x, probs[s]):
                w = w * (pr if bit else 1 - pr)
            tot += w
        out[x] = tot
    return out

def upsets(n):
    pts = list(itertools.product((0, 1), repeat=n)); res = []
    for mask in range(1, 2 ** len(pts) - 1):
        S = {pts[i] for i in range(len(pts)) if (mask >> i) & 1}
        if all(tuple(max(a, b) for a, b in zip(x, y)) in S for x in S for y in pts):
            res.append(S)
    return res

U3 = upsets(3)
U4 = upsets(4)

def check(p, W, Wp):
    q = 1 - p
    sP = {'01': p, '10': 1 - q * q, '11': 1 - q * q}
    sM = {'01': 1 - q * q, '10': p, '11': 1 - q * q}
    tgt3 = law({'x': Fr(1)}, {'x': [1 - p] * 3}, 3)
    worst = None
    for psi in (W, Wp):
        Z = q + q * psi + p
        rho = {'01': q / Z, '10': q * psi / Z, '11': p / Z}
        L = law(rho, {s: [sP[s], sP[s], sM[s]] for s in rho}, 3)
        for A in U3:
            m = sum(L[x] for x in A) - sum(tgt3[x] for x in A)
            worst = m if worst is None or m < worst else worst
    Z = 1 - q * q
    rho = {'01': q * p / Z, '10': p * q / Z, '11': p * p / Z}
    L = law(rho, {s: [sP[s], sP[s], sM[s], sM[s]] for s in rho}, 4)
    tgt4 = law({'x': Fr(1)}, {'x': [1 - p] * 4}, 4)
    rootw = min(sum(L[x] for x in A) - sum(tgt4[x] for x in A) for A in U4)
    origin = (1 - q * q) - (1 - p)
    return worst, rootw, origin

def run(p, ordered=True, N=60):
    q = 1 - p; c = q / p
    R0 = (1 / (1 + q), 1 / q); R = [R0]
    for k in range(1, N + 2):
        rmin, rmax = R[k - 1]
        a = dn(lfun(c * rmin * rmin, c * rmax)); b = upr(lfun(c * rmax * rmax, c * rmin))
        R.append((min(R0[0], a), max(R0[1], b)))
    if not ordered:
        R = [R[-1]] * (N + 2)
    def sib(j, t):
        if not ordered or j >= 2: return R[max(j, 1)]
        return R[1] if t == 'P' else R[0]
    def step(j, W, Wp):
        sM = sib(j, 'P'); RP = R[j] if j >= 2 else R[1]
        Pmin = ufun(c * RP[1], c * sM[0] * W); Pmax = ufun(c * RP[0], c * sM[1] * Wp)
        sM2 = sib(j, 'M')
        Mmin = ufun(c * Wp, c * sM2[0] * sM2[0]); Mmax = ufun(c * W, c * sM2[1] * sM2[1])
        return dn(min(Pmin, Mmin)), upr(max(Pmax, Mmax))
    wW, wWp = Fr(1), Fr(0)
    for n in range(1, N + 2):
        Rn = R[n]
        W = dn(ufun(c * Rn[1], c * Rn[0] * Rn[0])); Wp = upr(ufun(c * Rn[0], c * Rn[1] * Rn[1]))
        for j in range(n - 1, 0, -1):
            sM = sib(j, 'P'); RP = R[j] if j >= 2 else R[1]
            Pmin = ufun(c * RP[1], c * sM[0] * W); Pmax = ufun(c * RP[0], c * sM[1] * Wp)
            sM2 = sib(j, 'M')
            Mmin = ufun(c * Wp, c * sM2[0] * sM2[0]); Mmax = ufun(c * W, c * sM2[1] * sM2[1])
            W, Wp = dn(min(Pmin, Mmin)), upr(max(Pmax, Mmax))
        wW = min(wW, W); wWp = max(wWp, Wp)
    if ordered:
        # profondeurs n >= N+2 : au sommet de profondeur n-N-1 on prend seulement psi dans [0,1]
        W, Wp = Fr(0), Fr(1)
        for j in range(N, 0, -1):
            W, Wp = step(j, W, Wp)
        say("  queue (profondeurs >= %d) : [%.15f, %.15f] ; profondeurs <= %d : [%.15f, %.15f]" % (N + 2, float(W), float(Wp), N + 1, float(wW), float(wWp)))
        wW = min(wW, W); wWp = max(wWp, Wp)
    return check(p, wW, wWp)

say("Auto-test : version information totale (ordre fixe) a p=0.4433, marge attendue +1.7e-4 ...")
p = Fr(4433, 10000)
w, r, o = run(p, ordered=False)
say("  marge information totale a p=0.4433 :", float(w), "(attendu positif, ~1.7e-4)")
if not (w > 0):
    say("ECHEC de l'auto-test : arret."); sys.exit(1)

for p in (Fr(11077, 25000), Fr(443079, 1000000)):
    w, r, o = run(p, ordered=True)
    ok = (w > 0) and (r > 0) and (o >= 0)
    say("p =", p, "-> borne 1-p =", float(1 - p))
    say("  marge minimale (18 evenements croissants, deux bornes) :", float(w))
    say("  marge etape racine :", float(r), " ; P(xi(O)=1)-(1-p) :", float(o))
    say("  RESULTAT :", "CERTIFIE" if ok else "ECHEC")
say("Fin.")