# -*- coding: ascii -*-
# VERIF_pLC_B12.py
# Verification autonome de : sup_{k>=0} B(2**k alpha) >= C+1 pour C = 1..8 (et 9, 10, 11 en option).
# Python pur. Aucune dependance externe (ni numpy ni scipy). Aucun fichier intermediaire.
# Aucun caractere de disjonction binaire dans ce fichier ; les puissances s'ecrivent **.
# Parametres FIXES ci-dessous : lancer sans argument (double-clic ou bouton Executer).
#
# David Niedbala Giraudin, ORCID 0009-0009-1526-1178
# Technical assistance: Claude (Anthropic)

from math import gcd, isqrt
from fractions import Fraction
import time, sys
from array import array
from bisect import bisect_left

# ----------------------------------------------------------------------------- PARAMETRES
MAX_C        = 8      # verifier C = 1..MAX_C.  Couts mesures, machine a 2 coeurs :
                      #   C <= 8 : 14 s         C <= 9  : 1 min 30
                      #   C <= 10 : 6 min 15    C <= 11 : 34 min, pic 31,4 M candidats,
                      # moins de 4 Go.  MAX_C = 11 reconstruit la table du Theoreme 1.
BUDGET_SEC   = 7200   # budget total ; les valeurs de C trop lourdes sont SAUTEES.
MAX_CAND     = 40000000 # garde-fou memoire : un niveau au-dela est saute (annonce avant).
FAIRE_TEMOINS = True  # rejouer deux temoins explicites (arithmetique exacte, rapide).
# -----------------------------------------------------------------------------

# valeurs de controle (calculees par ce programme ; C <= 5 conformes a Vitorino-Vukusic)
ATTENDU_N = {1: 1, 2: 1, 3: 3, 4: 5, 5: 5, 6: 8, 7: 15, 8: 21, 9: 29, 10: 30, 11: 31}
ATTENDU_D = {3: [17], 4: [17], 5: [17, 33, 41], 6: [17, 33, 41, 1241]}

# ------------------------------------------------------- matrices de Raney, niveau L / R
def redrun(M):
    """reduction de Raney par PLAGES entieres ; renvoie (matrice reduite, plages emises)."""
    a, b, c, d = M
    runs = []
    while True:
        if a >= c and b >= d:
            if c == 0:
                if d == 0:
                    return (a, b, c, d), runs
                k = b // d
            elif d == 0:
                k = a // c
            else:
                k = min(a // c, b // d)
            if k == 0:
                return (a, b, c, d), runs
            a -= k * c
            b -= k * d
            runs.append((0, k))
        elif c >= a and d >= b:
            if a == 0:
                if b == 0:
                    return (a, b, c, d), runs
                k = d // b
            elif b == 0:
                k = c // a
            else:
                k = min(c // a, d // b)
            if k == 0:
                return (a, b, c, d), runs
            c -= k * a
            d -= k * b
            runs.append((1, k))
        else:
            return (a, b, c, d), runs

def mulrun(M, x, a):
    """M multipliee a droite par L**a (x=0) ou R**a (x=1)."""
    p, q, r, s = M
    if x == 0:
        return (p, p * a + q, r, r * a + s)
    return (p + q * a, q, r + s * a, s)

def matmul(N, M):
    a, b, c, d = N
    p, q, r, s = M
    return (a * p + b * r, a * q + b * s, c * p + d * r, c * q + d * s)

def closure(M0):
    M0 = redrun(M0)[0]
    seen = {M0}
    st = [M0]
    while st:
        M = st.pop()
        for x in (0, 1):
            N = redrun(mulrun(M, x, 1))[0]
            if N not in seen:
                seen.add(N)
                st.append(N)
    return sorted(seen)

IDENT = (1, 0, 0, 1)
S1 = closure((2, 0, 0, 1))
S2 = closure((4, 0, 0, 1))

def candidats(M1, M2):
    """M_e = reduit(N . M_{e-1}) avec N dans S1, filtre contre M_{e-2} par S2."""
    out = set(redrun(matmul(N, M1))[0] for N in S1)
    if M2 is not None:
        out &= set(redrun(matmul(N, M2))[0] for N in S2)
    return sorted(out)

# ------------------------------------------------------------------- auto-test arithmetique
def autotest():
    ok = True
    if len(S1) != 4 or len(S2) != 12:
        print("AUTO-TEST : tailles des transducteurs 2**1 et 2**2 = %d, %d au lieu de 4, 12"
              % (len(S1), len(S2)))
        ok = False
    for j in range(1, 7):
        n = len(closure((2 ** j, 0, 0, 1)))
        att = [4, 12, 35, 99, 269, 707][j - 1]
        if n != att:
            print("AUTO-TEST : cloture de 2**%d = %d au lieu de %d" % (j, n, att))
            ok = False
    # le transducteur doit rendre la vraie fraction continue
    for (p, q, mot, deb) in [(2, 1, [1, 1, 1, 1, 1, 1], [3]), (1, 2, [3, 1, 4, 1, 5], [1])]:
        M = redrun((p, 0, 0, q))[0]
        em = []
        for i, a in enumerate(mot):
            M, runs = redrun(mulrun(M, i % 2, a))
            em += runs
        plages = []
        for (t, k) in em:
            if plages and plages[-1][0] == t:
                plages[-1][1] += k
            else:
                plages.append([t, k])
        if not plages or plages[0][1] != deb[0]:
            print("AUTO-TEST : premiere plage de sortie %s au lieu de %s" % (plages[:1], deb))
            ok = False
    if not ok:
        print("AUTO-TEST ECHOUE : arret.")
        sys.exit(1)
    print("auto-test arithmetique : OK", flush=True)

# --------------------------------------------------------------------------- noyau K_n
def elaguer(succ, n, W):
    outd = array('i', bytes(4 * n))
    ind  = array('i', bytes(4 * n))
    for s in range(n):
        b = s * W
        for a in range(1, W):
            t = succ[b + a]
            if t >= 0:
                outd[s] += 1
                ind[t] += 1
    start = array('i', bytes(4 * (n + 1)))
    tot = 0
    for s in range(n):
        start[s] = tot
        tot += ind[s]
    start[n] = tot
    pos = array('i', start[:n])
    items = array('i', bytes(4 * tot)) if tot else array('i')
    for s in range(n):
        b = s * W
        for a in range(1, W):
            t = succ[b + a]
            if t >= 0:
                items[pos[t]] = s
                pos[t] += 1
    del pos
    alive = bytearray(b'\x01') * n
    pile = [s for s in range(n) if outd[s] == 0 or ind[s] == 0]
    while pile:
        s = pile.pop()
        if not alive[s]:
            continue
        alive[s] = 0
        b = s * W
        for a in range(1, W):
            t = succ[b + a]
            if t >= 0 and alive[t]:
                ind[t] -= 1
                if ind[t] == 0:
                    pile.append(t)
        for k in range(start[s], start[s + 1]):
            u = items[k]
            if alive[u]:
                outd[u] -= 1
                if outd[u] == 0:
                    pile.append(u)
    return alive

def un_niveau(C, K_par, K_m1, K_m2, KS):
    W = C + 1
    nK = len(K_par)
    allow = []
    need = set()
    for s in range(nK):
        cs = candidats(K_m1[s], K_m2[s])
        allow.append(cs)
        need.update(cs)
    need = sorted(need)
    dm = dict((m, i) for i, m in enumerate(need))
    KK = 2 * C
    E = array('i', b'\xff\xff\xff\xff' * (len(need) * 2 * W * KK))
    for M0 in need:
        d0 = dm[M0]
        for x in (0, 1):
            for a in range(1, W):
                Mr, runs = redrun(mulrun(M0, x, a))
                if Mr not in dm:
                    continue
                dr = dm[Mr]
                for pl0 in (0, 1):
                    for pr0 in range(W):
                        pl, pr = pl0, pr0
                        mort = False
                        for (t, k) in runs:
                            if pr == 0 or t != pl:
                                pl, pr = t, k
                            else:
                                pr += k
                            if pr > C:
                                mort = True
                                break
                        if not mort:
                            E[(((d0 * 2 + pl0) * W + pr0) * KK) + x * C + a - 1] = (dr * 2 + pl) * W + pr
    SEE = 2 * W * len(need) + W
    cpar = array('i')
    ce = array('i')
    keys = array('q')
    for s in range(nK):
        for M in allow[s]:
            d0 = dm[M]
            for pl in (0, 1):
                for pr in range(W):
                    if pr == 0 and pl == 1:
                        continue
                    ee = (d0 * 2 + pl) * W + pr
                    cpar.append(s)
                    ce.append(ee)
                    keys.append(s * SEE + ee)
    n = len(ce)
    if n > MAX_CAND:
        return None, n, 0, 0, None, None, None, None
    succ = array('i', b'\xff\xff\xff\xff' * (n * W))
    for i in range(n):
        s = cpar[i]
        ee = ce[i]
        x = K_par[s]
        b0 = ee * KK
        for a in range(1, W):
            e2 = E[b0 + x * C + a - 1]
            if e2 < 0:
                continue
            t = KS[s * W + a]
            if t < 0:
                continue
            kk = t * SEE + e2
            j = bisect_left(keys, kk)
            if j < n and keys[j] == kk:
                succ[i * W + a] = j
    del E
    alive = elaguer(succ, n, W)
    keep = [s for s in range(n) if alive[s]]
    remap = dict((s, i) for i, s in enumerate(keep))
    nk = len(keep)
    ns = array('i', b'\xff\xff\xff\xff' * (nk * W))
    for s in keep:
        for a in range(1, W):
            t = succ[s * W + a]
            if t >= 0 and alive[t]:
                ns[remap[s] * W + a] = remap[t]
    exces = sum(1 for v in ns if v >= 0) - nk
    npar = [K_par[cpar[s]] for s in keep]
    nm1 = [need[ce[s] // (2 * W)] for s in keep]
    nm2 = [K_m1[cpar[s]] for s in keep]
    return "ok", n, nk, exces, npar, nm1, nm2, ns

def discriminant(mot):
    p, pp, q, qp = 1, 0, 0, 1
    for a in mot:
        p, pp, q, qp = p * a + pp, p, q * a + qp, q
    T = p + qp
    det = (-1) ** len(mot)
    g = gcd(gcd(abs(q), abs(qp - p)), abs(pp))
    if g == 0:
        g = 1
    return (T * T - 4 * det) // (g * g)

def chaine(C, tmax):
    W = C + 1
    K_par = [0, 1]
    K_m1 = [IDENT, IDENT]
    K_m2 = [None, None]
    KS = [-1] * (2 * W)
    for p in (0, 1):
        for a in range(1, W):
            KS[p * W + a] = 1 - p
    for lev in range(1, tmax + 1):
        t0 = time.time()
        st, n, nk, exces, npar, nm1, nm2, ns = un_niveau(C, K_par, K_m1, K_m2, KS)
        if st is None:
            print("   C=%d T=%2d : %d candidats, au-dela du garde-fou memoire -> SAUTE"
                  % (C, lev, n), flush=True)
            return None, None
        print("   C=%d T=%2d  candidats=%9d  noyau=%8d  exces=%7d  (%.1f s)"
              % (C, lev, n, nk, exces, time.time() - t0), flush=True)
        K_par, K_m1, K_m2, KS = npar, nm1, nm2, ns
        if nk == 0:
            return lev, []
        if exces == 0:
            vu = set()
            ds = set()
            for s in range(nk):
                if s in vu:
                    continue
                chem = []
                pos = {}
                cur = s
                while cur not in pos and cur not in vu:
                    pos[cur] = len(chem)
                    lt = [a for a in range(1, W) if KS[cur * W + a] >= 0][0]
                    chem.append((cur, lt))
                    cur = KS[cur * W + lt]
                if cur in pos:
                    ds.add(discriminant([l for (_, l) in chem[pos[cur]:]]))
                for (v, _) in chem:
                    vu.add(v)
            return lev, sorted(ds)
    return None, None

# ---------------- temoins explicites des chaines dont le Lemme 10 de Vitorino-Vukusic
# ---------------- affirme l'existence (ces nombres ne figurent pas dans leur article)
def temoins():
    PREC = 700
    def rac(D):
        return Fraction(isqrt(D * 10 ** (2 * PREC)), 10 ** PREC)
    def fc(x, n):
        o = []
        for _ in range(n):
            a = x.numerator // x.denominator
            o.append(a)
            x = x - a
            if x == 0:
                break
            x = 1 / x
        return o
    for (nom, P0, D, Q0, Battendu, Kmax) in [("D=33", 26513, 33, 32768, 5, 12),
                                             ("D=1241", 224485, 1241, 557056, 6, 10)]:
        beta = (Fraction(P0) + rac(D)) / Q0
        pire = 0
        for k in range(0, Kmax + 1):
            c = fc(beta * (2 ** k), 400)
            pire = max(pire, max(c[60:330]))
        etat = "OK" if pire == Battendu else "ECART"
        print("   temoin %-8s : sup_{k=0..%d} B(2**k beta) = %d  (attendu %d)  %s"
              % (nom, Kmax, pire, Battendu, etat), flush=True)

# ------------------------------------------------------------------------------- main
def main():
    t0 = time.time()
    print("VERIF_pLC_B12 -- verification de sup_{k>=0} B(2**k alpha) >= C+1", flush=True)
    print("", flush=True)
    autotest()
    print("", flush=True)
    if FAIRE_TEMOINS:
        print("temoins explicites (chaines du Lemme 10 de Vitorino-Vukusic) :", flush=True)
        temoins()
        print("", flush=True)
    print("chaines de noyaux (T = nombre de transducteurs ; fenetre k = 0..T) :", flush=True)
    resume = []
    for C in range(1, MAX_C + 1):
        if time.time() - t0 > BUDGET_SEC:
            print("   budget epuise : C = %d et au-dela sont SAUTES" % C, flush=True)
            break
        print("", flush=True)
        N, ds = chaine(C, 60)
        if N is None:
            resume.append((C, None, None))
            continue
        att = ATTENDU_N.get(C)
        marque = "OK" if att is None or N == att else "ECART (attendu %d)" % att
        print("   ==> N(%d) = %d transducteurs   %s" % (C, N, marque), flush=True)
        if ds:
            impairs = [d for d in ds if d % 2 == 1]
            mauvais = [d for d in impairs if d % 8 != 1]
            print("       %d discriminants ; %d impairs, dont %d non congrus a 1 mod 8"
                  % (len(ds), len(impairs), len(mauvais)), flush=True)
            if C in ATTENDU_D:
                m = "OK" if [d for d in ds if d % 2 == 1] == ATTENDU_D[C] else "ECART"
                print("       liste impaire = %s   %s" % (impairs, m), flush=True)
        resume.append((C, N, ds))
    print("", flush=True)
    print("RESUME", flush=True)
    best = 0
    for (C, N, ds) in resume:
        if N is None:
            print("   C = %2d : saute" % C, flush=True)
        else:
            print("   C = %2d : N(C) = %2d  =>  sup_k B(2**k alpha) >= %d" % (C, N, C + 1), flush=True)
            best = max(best, C + 1)
    print("", flush=True)
    print("CONCLUSION VERIFIEE ICI : sup_{k>=0} B(2**k alpha) >= %d" % best, flush=True)
    print("temps total %.1f s" % (time.time() - t0), flush=True)

if __name__ == "__main__":
    main()