# -*- coding: utf-8 -*-
# verify_ball_N1840.py  --  verification rigoureuse, independante, du contre-exemple de degre 29 440
# au probleme de la valeur moyenne de Smale (K = 1).
#
# A LANCER EN PYTHON PUR (pas dans Sage) :  py verify_ball_N1840.py
# Prerequis, une seule fois :  py -m pip install python-flint numpy
# Aucun argument : tous les reglages sont ci-dessous. Pas de scipy.
# Les deux fichiers polynomial_N1840.json et centres_N1840.csv doivent etre dans le MEME dossier.
#
# Ce que fait le programme :
#  1. auto-test de l'arithmetique ;
#  2. lit le polynome EXACT (coefficients rationnels) et controle sa coherence (degre 16N, lambda = N*gamma,
#     phase exactement unimodulaire, coefficient dominant non nul) ;
#  3. lit les 29 439 centres et calcule leurs distances mutuelles (disjonction des boites) ;
#  4. pour chaque point critique de la tranche : Newton en haute precision, test de Krawczyk (existence ET
#     unicite d'une racine de P' dans une boite), puis minoration RIGOUREUSE de
#     |Q| = |P(c)-P(b)| / (|c-b| |P'(b)|) sur toute la boite, en arithmetique de boules a 1280 bits.
#  Le contre-exemple exige min |Q| > 1 sur les 29 439 = deg P' points critiques.
#  Un passage complet demande plusieurs heures : le programme traite une TRANCHE [DEBUT, FIN) et s'arrete
#  avant TEMPS_MAX secondes. Les centres sont classes du plus serre au plus large : la tranche 0-300
#  contient les cas les plus durs.
# Credits : construction de Adamczewski et Colbrook, arXiv:2610.36158 ; parametres re-optimises par
# David Niedbala Giraudin.

DEBUT = 0           # premier point critique de la tranche (0 = le plus serre)
FIN = 29439         # fin de tranche (exclue) ; 29439 = tout
TEMPS_MAX = 10**9    # secondes
PREC = 1280         # bits de precision
JSON = "polynomial_N1840.json"
CENTRES = "centres_N1840.csv"

import sys, os, time
print("VERIF_SMALE_1840 : demarrage", flush=True)
T0 = time.time()
try:
    import flint
except ImportError:
    print("ERREUR : python-flint manque. Tapez :  py -m pip install python-flint", flush=True); sys.exit(1)
import numpy as np
import json
from fractions import Fraction
from flint import acb, arb, acb_poly, ctx

def A(x):            # conversion exacte : Fraction / complex(float) -> acb exact
    if isinstance(x, (complex, np.complexfloating)):
        return acb(arb(Fraction(float(x.real)).numerator)/arb(Fraction(float(x.real)).denominator),
                   arb(Fraction(float(x.imag)).numerator)/arb(Fraction(float(x.imag)).denominator))
    if isinstance(x, (float, np.floating)):
        f = Fraction(float(x)); return acb(arb(f.numerator)/arb(f.denominator))
    if isinstance(x, Fraction): return acb(arb(x.numerator)/arb(x.denominator))
    return acb(x)

class BallFamily:
    """Les parametres sont pris EXACTEMENT comme les flottants (double) stockes : le polynome certifie est celui-la."""
    def __init__(self, F, prec=192, json_path=None):
        ctx.prec = prec
        self.exact = None
        if json_path is not None:
            d = json.load(open(json_path)); self.exact = d
            fq = lambda s_: Fraction(s_)
            ra = lambda s_: arb(Fraction(s_).numerator)/arb(Fraction(s_).denominator)
            def cq(x):
                if 're_numerator' in x: return acb(arb(int(x['re_numerator']))/arb(int(x['denominator'])), arb(int(x['im_numerator']))/arb(int(x['denominator'])))
                return acb(ra(x['re']), ra(x['im']))
            self._cq = cq
            def arr(lst, blocks):
                M = [[acb(0)]*16 for _ in range(blocks)]
                for e in lst: M[e['ell']][e['k']] = acb(ra(e['re']), ra(e['im']))
                return M
            self.exR0 = arr(d['correction_coefficients'], int(d['correction_blocks']))
            self.exX = arr(d['extra_correction_coefficients'], int(d['extra_correction_blocks']))
            self.exE = arr(d['tail_coefficients'], int(d['tail_blocks']))
            self.exeta = [cq(x) for x in d['free_coefficients']]
        self.F = F; self.m = 16; self.N = F.N; self.L = F.L; self.tb = F.tail_blocks
        if self.exact is not None:
            d = self.exact; self.Aa = A(Fraction(d['A'])); self.r = self._cq(d['r']); self.a = self._cq(d['a']); self.bc = self._cq(d['b'])
            self.xi = [self._cq(x) for x in d['correction_centers']]
        else:
            self.Aa = A(F.A); self.r = A(F.r); self.a = A(F.a); self.bc = A(F.bc); self.xi = [A(x) for x in F.xi]
        self.mu = (self.xi[0] + self.xi[1])/2
        self.C = 1 + self.Aa*self.Aa; self.S = self.Aa + 1/self.Aa; rb = self.r.conjugate(); I = acb(0, 1)
        self.d0 = self.C + rb**16*(-self.Aa*self.Aa + I*self.Aa); self.d1 = I*self.Aa - rb**16*(1 + I*self.Aa)
        self.alpha = -self.C*self.d1/(I*self.Aa*self.d0)
        # polynomes A_k(v) = sum_j c_{k,j} v**j, k = 1..15
        self.Ak = []
        for k in range(1, 16):
            nu = 1 - arb(k)/16; b = acb(1); cc = acb(1); co = [acb(1)]
            for j in range(1, self.L + 1):
                b = b*(nu - j + 1)/j; cc = self.alpha*cc + b; co.append(cc)
            p = acb_poly(co); self.Ak.append((p, p.derivative(), p.derivative().derivative()))
        self.wk = [(1 - abs(self.r)**2)*self.C/self.d0*rb**(k-1)*self.bc**k for k in range(1, 16)]
        # corrections en polynome de z (degre modere) : R0, X, E
        self.zq = acb_poly([I*self.Aa/self.S] + [0]*15 + [1/self.S])     # q(z)
        def corr_poly(M):
            out = acb_poly([0]); qpow = acb_poly([1]); zs = acb_poly([0, acb(5)/6])
            rows = M if isinstance(M, list) else [[A(complex(M[l, k])) for k in range(16)] for l in range(M.shape[0])]
            for ell in range(len(rows)):
                w = acb_poly([0]); zk = acb_poly([1])
                for k in range(16):
                    w = w + zk*rows[ell][k]
                    zk = zk*zs
                out = out + w*qpow; qpow = qpow*self.zq
            return out
        Z = acb_poly([0, 1]); u = Z - self.a
        K = (Z - self.xi[0])**2*(Z - self.xi[1])**2
        T = acb_poly([0]); t = Z - self.mu
        eta = self.exeta if self.exact is not None else [A(complex(x)) for x in F.eta]
        for j in range(3, -1, -1): T = T*t + eta[j]
        R0 = self.exR0 if self.exact is not None else F.R0
        Xc = self.exX if self.exact is not None else F.Xc
        Ec = self.exE if self.exact is not None else F.Ec
        self.Pcorr = u**2*(K*corr_poly(R0) + T) + u**2*corr_poly(Xc)      # Delta0 + X
        self.Epoly = corr_poly(Ec)                                          # E(z)
        if self.exact is not None:
            d = self.exact
            self.gamma = A(Fraction(d['high_power_coefficient'])); self.omega = self._cq(d['tail_phase']); self.base = self._cq(d['normalization_base'])
            assert int(d['N']) == F.N and int(d['L']) == F.L and int(d['tail_blocks']) == F.tail_blocks
        else:
            self.gamma = A(complex(F.gamma)); self.omega = A(complex(F.omega)); self.base = A(complex(F.base))
        self.ha = None; self.ha = self.h_jet(self.a)[0]
    def h_jet(self, z):
        I = acb(0, 1); q = (z**16 + I*self.Aa)/self.S; dq = 16*z**15/self.S; ddq = 240*z**14/self.S
        v = I*q; dv = I*dq; ddv = I*ddq
        x = self.alpha*v; L = self.L; xL = x**L
        G = (1 - xL)/(1 - x)
        n = -L*x**(L-1)*(1 - x) + (1 - xL)
        Gp = n/(1 - x)**2
        Gpp = (-L*(L-1)*x**(L-2)*(1 - x)**2 + 2*n)/(1 - x)**3
        dx = self.alpha*dv; ddx = self.alpha*ddv
        Gj = (G, Gp*dx, Gpp*dx*dx + Gp*ddx)
        rb15 = self.r.conjugate()**15*(1 + I*self.Aa)
        J = (-self.r*self.C*(1 + v) + rb15*z**16, -self.r*self.C*dv + rb15*16*z**15, -self.r*self.C*ddv + rb15*240*z**14)
        f0 = (J[0]*Gj[0])/self.d0; f1 = (J[1]*Gj[0] + J[0]*Gj[1])/self.d0; f2 = (J[2]*Gj[0] + 2*J[1]*Gj[1] + J[0]*Gj[2])/self.d0
        for k in range(1, 16):
            p, dp, ddp = self.Ak[k-1]
            a0 = p(v); a1 = dp(v)*dv; a2 = ddp(v)*dv*dv + dp(v)*ddv
            z0 = z**k; z1 = k*z**(k-1); z2 = k*(k-1)*z**(k-2) if k >= 2 else acb(0)
            w = self.wk[k-1]
            f0 += w*z0*a0; f1 += w*(z1*a0 + z0*a1); f2 += w*(z2*a0 + 2*z1*a1 + z0*a2)
        return f0, f1, f2
    def P_jet(self, z):
        h0, h1, h2 = self.h_jet(z); u = z - self.a
        f0 = u*h0/self.ha; f1 = (h0 + u*h1)/self.ha; f2 = (2*h1 + u*h2)/self.ha
        pc = self.Pcorr; dpc = pc.derivative(); ddpc = dpc.derivative()
        c0, c1, c2 = pc(z), dpc(z), ddpc(z)
        I = acb(0, 1); q = (z**16 + I*self.Aa)/self.S; dq = 16*z**15/self.S; ddq = 240*z**14/self.S
        N, tb = self.N, self.tb
        E = self.Epoly; e0, e1, e2 = E(z), E.derivative()(z), E.derivative().derivative()(z)
        # Psi = gamma q**tb + E ; H = omega q**(N-tb) Psi
        s0 = self.gamma*q**tb + e0
        s1 = self.gamma*tb*q**(tb-1)*dq + e1
        s2 = self.gamma*(tb*(tb-1)*q**(tb-2)*dq*dq + tb*q**(tb-1)*ddq) + e2
        M = N - tb
        g0 = q**M; g1 = M*q**(M-1)*dq; g2 = M*(M-1)*q**(M-2)*dq*dq + M*q**(M-1)*ddq
        H0 = self.omega*g0*s0; H1 = self.omega*(g1*s0 + g0*s1); H2 = self.omega*(g2*s0 + 2*g1*s1 + g0*s2)
        return f0 + c0 + H0, f1 + c1 + H1, f2 + c2 + H2


class Majorant:
    """Majoration rigoureuse de |P| sur le disque |zeta - z| <= R, par inegalite triangulaire terme a terme."""
    def __init__(self, B):
        self.B = B; F = B.F
        self.absA = [[abs(x) for x in B.Ak[k][0].coeffs()] for k in range(15)]   # |c_{k,j}| (arb)
        def absmat(M):
            if isinstance(M, list): return [[abs(x) for x in row] for row in M]
            return [[abs(A(complex(M[l, k]))) for k in range(16)] for l in range(M.shape[0])]
        if B.exact is not None:
            self.aR0, self.aX, self.aE = absmat(B.exR0), absmat(B.exX), absmat(B.exE); self.aeta = [abs(x) for x in B.exeta]
        else:
            self.aR0, self.aX, self.aE = absmat(F.R0), absmat(F.Xc), absmat(F.Ec); self.aeta = [abs(A(complex(x))) for x in F.eta]
        self.awk = [abs(w) for w in B.wk]
        self.aalpha = abs(B.alpha); self.ar = abs(B.r); self.aA = abs(B.Aa); self.S = abs(B.S); self.C = abs(B.C)
        self.ad0 = abs(B.d0); self.aha = abs(B.ha); self.aa = abs(B.a); self.aom = abs(B.omega); self.agam = abs(B.gamma)
        self.amu = abs(B.mu); self.axi = [abs(x) for x in B.xi]
        self.a1iA = abs(1 + acb(0, 1)*B.Aa)
    def bound(self, z, R):
        B = self.B; rz = abs(z) + R                                   # |zeta| <= rz (arb, borne sup via .upper())
        rz = arb(rz.upper())
        I = acb(0, 1)
        qz = abs((z**16 + I*B.Aa)/B.S)                             # |q(z)| au centre (boule exacte)
        Qb = arb((qz + 16*rz**15*R/self.S).upper())                # |q(zeta)| <= |q(z)| + sup|q'| R
        x = self.aalpha*Qb
        L = B.L
        Gb = (1 + x**L)/(1 - x) if (x.upper() < 1) else arb(L)*(1 + x**L)    # |G| = |sum_{j<L} (alpha v)**j|
        Jb = self.ar*self.C*(1 + Qb) + self.ar**15*self.a1iA*rz**16
        hb = Jb*Gb/self.ad0
        for k in range(15):
            s = arb(0); p = arb(1)
            for cj in self.absA[k]: s += cj*p; p *= Qb
            hb += self.awk[k]*rz**(k+1)*s
        fb = (rz + self.aa)*hb/self.aha
        zs = arb(5)/6*rz
        def blk(M):
            s = arb(0); qp = arb(1)
            for row in M:
                w = arb(0); zp = arb(1)
                for x_ in row: w += x_*zp; zp *= zs
                s += w*qp; qp *= Qb
            return s
        Kb = (rz + self.axi[0])**2*(rz + self.axi[1])**2
        Tb = arb(0); tp = arb(1)
        for e in self.aeta: Tb += e*tp; tp *= (rz + self.amu)
        cb = (rz + self.aa)**2*(Kb*blk(self.aR0) + Tb) + (rz + self.aa)**2*blk(self.aX)
        tb = B.tb; M = B.N - tb
        Hb = self.aom*Qb**M*(self.agam*Qb**tb + blk(self.aE))
        return arb((fb + cb + Hb).upper())



# ---------------- 1. auto-test ----------------
ctx.prec = 200
_t = (arb(2).sqrt()*arb(2).sqrt() - 2).abs_upper()
assert float(_t) < 1e-50, "auto-test sqrt"
assert pow(7, 5) == 16807 and int(Fraction(3, 4)*4) == 3, "auto-test entiers"
_z = acb(1, 1)*acb(1, -1)
assert abs(float(_z.real.mid()) - 2.0) < 1e-30 and float(abs(_z.imag).upper()) < 1e-30, "auto-test complexes"
print("1. auto-test arithmetique : OK", flush=True)

# ---------------- 2. polynome exact ----------------
HERE = os.path.dirname(os.path.abspath(__file__))
d = json.load(open(os.path.join(HERE, JSON)))
N = int(d['N']); L = int(d['L']); tb = int(d['tail_blocks'])
assert int(d['degree']) == 16*N, "degre"
assert L == N - 2, "L"
assert Fraction(d['lambda']) == N*Fraction(d['high_power_coefficient']), "lambda = N gamma"
ph = d['tail_phase']; assert Fraction(ph['re'])*Fraction(ph['re']) + Fraction(ph['im'])*Fraction(ph['im']) == 1, "phase unimodulaire"
assert Fraction(d['high_power_coefficient']) != 0, "coefficient dominant non nul"
print(f"2. polynome exact lu : N = {N}, degre = {16*N}, deg P' = {16*N-1} ; controles de coherence OK", flush=True)

class _F: pass
F = _F(); F.N = N; F.L = L; F.tail_blocks = tb

# ---------------- 3. centres et disjonction ----------------
pts = []
for line in open(os.path.join(HERE, CENTRES)):
    if line.startswith('#') or line.startswith('re') or not line.strip(): continue
    a_, b_ = line.strip().split(','); pts.append(complex(float(a_), float(b_)))
c = np.array(pts)
assert len(c) == 16*N - 1, "nombre de centres"
H = 1e-3; grid = {}
for i, z in enumerate(c): grid.setdefault((int(np.floor(z.real/H)), int(np.floor(z.imag/H))), []).append(i)
nn = np.full(len(c), np.inf)
for (gx, gy), lst in grid.items():
    voisins = []
    for dx in (-1, 0, 1):
        for dy in (-1, 0, 1): voisins += grid.get((gx+dx, gy+dy), [])
    for i in lst:
        for j in voisins:
            if j != i:
                dd = abs(c[i] - c[j])
                if dd < nn[i]: nn[i] = dd
nn = np.minimum(nn, H)     # tout voisin hors des 9 cellules est a distance >= H
print(f"3. {len(c)} centres lus ; distance minimale entre deux centres = {nn.min():.3e}", flush=True)

# ---------------- 4. certification de la tranche ----------------
B = BallFamily(F, PREC, os.path.join(HERE, JSON)); Mj = Majorant(B)
b = B.base; Pb0, Pb1, _ = B.P_jet(b)
print(f"   initialisation en boules : {time.time()-T0:.0f} s", flush=True)
worst = None; fails = []; rmax = 0.0; dmax = 0.0; faits = 0
fin = min(FIN, len(c))
for i in range(DEBUT, fin):
    if time.time() - T0 > TEMPS_MAX:
        print(f"   temps ecoule : arret propre apres {faits} points (reprendre avec DEBUT = {i})", flush=True); break
    z = A(complex(c[i]))
    for _ in range(3):
        p0, p1, p2 = B.P_jet(z); z = z - p1/p2
        z = acb(z.real.mid(), z.imag.mid())
    dmax = max(dmax, abs(complex(float(z.real.mid()), float(z.imag.mid())) - c[i]))
    p0, p1, p2 = B.P_jet(z)
    R = arb(0.3*nn[i])
    Bd = Mj.bound(z, R)
    stp = abs(p1/p2)
    rho0 = arb(2)**int(np.ceil(np.log2(max(1e-40, 1e6*float(stp.upper())))))   # demi-cote dyadique : rayon exact
    Dz = acb(arb(z.real.mid(), rho0), arb(z.imag.mid(), rho0))      # boite effective
    rho = arb(max(Dz.real.rad(), Dz.imag.rad()))      # demi-cote majorant (exact), pour r2
    rin = arb(min(Dz.real.rad(), Dz.imag.rad()))      # demi-cote minorant (exact), pour l'inclusion
    rmax = max(rmax, float(rho.upper()))
    r2 = rho*arb(2).sqrt()
    if not bool(R > r2): print("ERREUR : disque de Cauchy trop petit au point", i, flush=True); sys.exit(1)
    M3 = 6*Bd/(R - r2)**3
    P2D = acb(arb(p2.real.mid(), arb((p2.real.rad() + M3*r2).upper())), arb(p2.imag.mid(), arb((p2.imag.rad() + M3*r2).upper())))
    Kr = z - p1/p2 + (1 - P2D/p2)*(Dz - z)
    inside = bool((Kr.real - z.real).abs_upper() < rin) and bool((Kr.imag - z.imag).abs_upper() < rin)
    dz = Dz - z
    rem = M3*r2*r2*r2/6
    PD = p0 + p1*dz + p2*dz*dz/2
    PD = acb(arb(PD.real.mid(), arb((PD.real.rad() + rem).upper())), arb(PD.imag.mid(), arb((PD.imag.rad() + rem).upper())))
    Qb = (PD - Pb0)/((Dz - b)*Pb1)
    low = Qb.abs_lower()
    if not inside or not bool(low > 1): fails.append(i)      # NaN ou |Q| <= 1 : compte comme echec
    if worst is None or low < worst[0]: worst = (low, i)
    faits += 1
    if faits == 1:
        print(f"   premier point certifie apres {time.time()-T0:.0f} s", flush=True)
    if faits % 10 == 0:
        print(f"   {faits} points ; min |Q| certifie = {float(worst[0]):.15f} ; echecs Krawczyk {len(fails)} ; {time.time()-T0:.0f} s", flush=True)

print("", flush=True)
print("=================== RESULTAT ===================", flush=True)
print(f"tranche traitee : [{DEBUT}, {DEBUT+faits})  soit {faits} points critiques", flush=True)
if faits:
    print("min |Q| certifie sur la tranche >= " + worst[0].lower().str(25, radius=False), flush=True)
    print(f"echecs (Krawczyk ou |Q| non > 1) : {len(fails)}", flush=True)
    print(f"demi-cote max des boites = {rmax:.1e} ; ecart max centre fourni / centre raffine = {dmax:.1e}", flush=True)
    # regle uniforme, valable d'une tranche a l'autre : chaque boite est dans le disque de rayon
    # 1e-9 + sqrt(2)*1e-20 autour de SON centre fourni ; les centres fournis sont a distance >= nn.min().
    sep = (rmax <= 1e-20) and (dmax <= 1e-9) and (nn.min() > 4e-9)
    print(f"boites disjointes (demi-cote <= 1e-20, derive <= 1e-9, centres >= {nn.min():.1e} l'un de l'autre) : {'OUI' if sep else 'NON'}", flush=True)
    ok = (not fails) and bool(worst[0].lower() > 1) and sep
    print("VERDICT TRANCHE : " + ("PASS (min |Q| > 1 certifie)" if ok else "ECHEC"), flush=True)
print(f"duree totale : {time.time()-T0:.0f} s", flush=True)