# -*- coding: utf-8 -*-
# Verificateur exact : p_c(site, Z2) <= 0.6463
# Python pur (aucun module externe), arithmetique entiere exacte. Lancer tel quel, sans argument.
# Auteur : David Niedbala Giraudin (Fino), ORCID 0009-0009-1526-1178, avec assistance IA (Claude, Anthropic).
#
# Chaine verifiee (methode de substitution par regions a sommets fantomes, regions 4x4) :
#   etape 1 : la loi de region de H a 0.5755 domine stochastiquement celle de LG a 0.50001
#   etape 2 : la loi de region de S a 0.6463 domine stochastiquement celle de H a 0.5755
# LG = reseau carre avec les deux diagonales dans une face sur deux (graphe adjoint de Z2, seuil site 1/2),
# H  = reseau carre avec les diagonales des faces (impair, impair) seulement, S = reseau carre.
from fractions import Fraction as Fr
from collections import deque, defaultdict
import time, sys, gc

# --- Reglages (a modifier si la machine manque de memoire) ---
FAIRE_AUTOTEST = True
FAIRE_ETAPE_1 = True    # pic memoire estime ~110 Mo (mesure sur la machine de developpement)
FAIRE_ETAPE_2 = True    # pic memoire estime ~60 Mo
# Si l'execution s'arrete a l'etape 1, mettre FAIRE_ETAPE_1 = False, relancer pour l'etape 2,
# puis l'inverse : le certificat complet est la conjonction des deux lignes VRAI.

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

T0 = time.time()

def cls_face(i, j, ox, oy):
    a, b = (i + ox) % 2, (j + oy) % 2
    if a == 1 and b == 1: return 'white'
    if a == 0 and b == 0: return 'black'
    return None

def box_region(W, H, diag, off):
    ox, oy = off
    V = [(x, y) for x in range(W) for y in range(H)]; Vs = set(V); E = []
    for (x, y) in V:
        for d in ((1, 0), (0, 1)):
            w = (x + d[0], y + d[1])
            if w in Vs: E.append(((x, y), w))
    for i in range(W - 1):
        for j in range(H - 1):
            if cls_face(i, j, ox, oy) in diag:
                E.append(((i, j), (i + 1, j + 1))); E.append(((i + 1, j), (i, j + 1)))
    return V, E

def ghosts(V, diag, off):
    ox, oy = off; Vs = set(V); G = []
    for (x, y) in V:
        for d in ((1, 0), (-1, 0), (0, 1), (0, -1)):
            w = (x + d[0], y + d[1])
            if w not in Vs: G.append(((x, y), ('S', (x, y), w)))
    for (x, y) in V:
        for dx in (1, -1):
            for dy in (1, -1):
                w = (x + dx, y + dy)
                if w in Vs: continue
                if cls_face(min(x, w[0]), min(y, w[1]), ox, oy) in diag:
                    G.append(((x, y), ('D', (x, y), w)))
    return G

def law_poly(V, E, G):
    n = len(V); idx = {v: i for i, v in enumerate(V)}
    gid = {}; gv = []
    for v, g in G:
        if g not in gid: gid[g] = len(gid)
        gv.append((idx[v], gid[g]))
    ng = len(gid); out = {}
    Ei = [(idx[a], idx[b]) for a, b in E]
    for mask in range(1 << n):
        par = list(range(n))
        def f(a):
            while par[a] != a:
                par[a] = par[par[a]]; a = par[a]
            return a
        for a, b in Ei:
            if (mask >> a) & 1 and (mask >> b) & 1:
                ra, rb = f(a), f(b)
                if ra != rb: par[ra] = rb
        lab = [None] * ng
        for vi, g in gv:
            lab[g] = ('o', f(vi)) if (mask >> vi) & 1 else ('c', g)
        m = {}; key = []
        for l in lab:
            if l not in m: m[l] = len(m)
            key.append(m[l])
        key = tuple(key)
        if key not in out: out[key] = [0] * (n + 1)
        out[key][bin(mask).count('1')] += 1
    return out, n, gid

def d4_maps(W):
    c2 = W - 1   # coordinates doubled around centre
    def mk(k):
        def g(pt):
            x, y = 2 * pt[0] - c2, 2 * pt[1] - c2
            for _ in range(k % 4): x, y = -y, x
            if k >= 4: x = -x
            return ((x + c2) // 2, (y + c2) // 2)
        return g
    return [mk(k) for k in range(8)]

def ghost_perms(gid, W):
    perms = []
    for m in d4_maps(W):
        perm = [None] * len(gid)
        for g, i in gid.items():
            t, v, w = g; perm[i] = gid[(t, m(v), m(w))]
        perms.append(perm)
    return perms

def canon(key):
    mp = {}; out = []
    for k in key:
        if k not in mp: mp[k] = len(mp)
        out.append(mp[k])
    return tuple(out)

def act(key, perm):
    new = [None] * len(key)
    for i, k in enumerate(key): new[perm[i]] = k
    return canon(new)

def eqmask(key):
    m = 0; t = 0; n = len(key)
    for i in range(n):
        for j in range(i + 1, n):
            if key[i] == key[j]: m |= (1 << t)
            t += 1
    return m

def masses(poly, n, p):
    a, b = p.numerator, p.denominator
    pw = [a ** k * (b - a) ** (n - k) for k in range(n + 1)]
    return {k: sum(c[i] * pw[i] for i in range(n + 1)) for k, c in poly.items()}, b ** n

def maxflow(N, edges, s, t):
    g = [[] for _ in range(N)]; to = []; cap = []; rev = []
    for u, v, c in edges:
        e1 = len(to); e2 = e1 + 1
        g[u].append(e1); to.append(v); cap.append(c); rev.append(e2)
        g[v].append(e2); to.append(u); cap.append(0); rev.append(e1)
    flow = 0
    while True:
        lvl = [-1] * N; lvl[s] = 0; dq = deque([s])
        while dq:
            u = dq.popleft()
            for e in g[u]:
                if cap[e] > 0 and lvl[to[e]] < 0:
                    lvl[to[e]] = lvl[u] + 1; dq.append(to[e])
        if lvl[t] < 0: return flow
        it = [0] * N
        while True:
            path = []; node = s; pushed = 0
            while True:
                if node == t:
                    bott = min(cap[e] for e in path)
                    for e in path:
                        cap[e] -= bott; cap[rev[e]] += bott
                    pushed = bott; break
                adv = False
                while it[node] < len(g[node]):
                    e = g[node][it[node]]
                    if cap[e] > 0 and lvl[to[e]] == lvl[node] + 1:
                        path.append(e); node = to[e]; adv = True; break
                    it[node] += 1
                if not adv:
                    if not path: break
                    lvl[node] = -1; e = path.pop(); node = to[rev[e]]; it[node] += 1
            if pushed == 0: break
            flow += pushed

def dominates(W, off, diag_hi, diag_lo, diag_ghost, p_hi, p_lo, label):
    t = time.time()
    V, Ehi = box_region(W, W, diag_hi, off); _, Elo = box_region(W, W, diag_lo, off)
    G = ghosts(V, diag_ghost, off)
    phi, n, gid = law_poly(V, Ehi, G); plo, _, gid2 = law_poly(V, Elo, G)
    perms = ghost_perms(gid, W)
    say("   lois calculees :", len(phi), "et", len(plo), "partitions (", round(time.time() - t), "s )")
    def classes(poly):
        rep = {}; mem = defaultdict(set)
        for k in poly:
            imgs = {act(k, pm) for pm in perms}; r = min(imgs); rep[k] = r; mem[r] |= imgs
        return rep, mem
    rH, mH = classes(phi); rL, mL = classes(plo)
    for k, c in phi.items():
        for pm in perms:
            if phi.get(act(k, pm)) != c: say("ECHEC : loi non invariante"); sys.exit(1)
    for k, c in plo.items():
        for pm in perms:
            if plo.get(act(k, pm)) != c: say("ECHEC : loi non invariante"); sys.exit(1)
    mh, dh = masses(phi, n, p_hi); ml, dl = masses(plo, n, p_lo)
    cH = defaultdict(int); cL = defaultdict(int)
    for k, v in mh.items(): cH[rH[k]] += v * dl
    for k, v in ml.items(): cL[rL[k]] += v * dh
    total = dh * dl
    assert sum(cH.values()) == total and sum(cL.values()) == total
    Ls = [c for c in cL if cL[c] > 0]; Hs = [c for c in cH if cH[c] > 0]
    Hmask = {b: [eqmask(x) for x in mH[b]] for b in Hs}
    iL = {c: 2 + i for i, c in enumerate(Ls)}; iH = {c: 2 + len(Ls) + j for j, c in enumerate(Hs)}
    edges = [(0, iL[c], cL[c]) for c in Ls] + [(iH[c], 1, cH[c]) for c in Hs]
    INF = 4 * total; ne = 0
    for a in Ls:
        ma = eqmask(a); na = ~ma
        for b in Hs:
            if any((ma & ~mb) == 0 for mb in Hmask[b]):
                edges.append((iL[a], iH[b], INF)); ne += 1
    say("   classes :", len(Ls), "/", len(Hs), " aretes d'ordre :", ne, "(", round(time.time() - t), "s )")
    f = maxflow(2 + len(Ls) + len(Hs), edges, 0, 1)
    ok = (f == total)
    say("  ", label, ":", "VRAI" if ok else "FAUX", " (", round(time.time() - t), "s )")
    return ok

if FAIRE_AUTOTEST:
    say("Auto-test (regions 2x2, valeur publiee de Wierman 1995 : 0.679492) ...")
    a = dominates(2, (1, 1), [], ['white'], [], Fr(6795, 10000), Fr(586668, 1000000), "S(0.6795) domine H(0.586668)")
    b = dominates(2, (1, 1), [], ['white'], [], Fr(6794, 10000), Fr(586668, 1000000), "S(0.6794) domine H(0.586668)")
    if not (a and not b):
        say("ECHEC de l'auto-test : arret."); sys.exit(1)
    say("Auto-test reussi.")
s1 = s2 = None
if FAIRE_ETAPE_1:
    say("Etape 1 (regions 4x4, ~110 Mo, ~1 a 5 min selon la machine) ...")
    s1 = dominates(4, (0, 0), ['white'], ['white', 'black'], ['white'], Fr(5755, 10000), Fr(50001, 100000), "H(0.5755) domine LG(0.50001)")
    gc.collect()
if FAIRE_ETAPE_2:
    say("Etape 2 (regions 4x4, ~60 Mo, ~1 a 3 min selon la machine) ...")
    s2 = dominates(4, (1, 1), [], ['white'], [], Fr(6463, 10000), Fr(5755, 10000), "S(0.6463) domine H(0.5755)")
if s1 and s2:
    say("RESULTAT : p_c(site, Z2) <= 0.6463 CERTIFIE (etapes 1 et 2 VRAI)")
elif s1 is False or s2 is False:
    say("RESULTAT : ECHEC")
else:
    say("RESULTAT PARTIEL : etape 1 =", s1, " etape 2 =", s2, " (lancer l'autre etape pour completer)")
say("Duree totale :", round(time.time() - T0), "s")