#!/usr/bin/env python3 r""" ================================================================================ PROUVEUR AUTONOME -- borne inferieure PROUVEE pour le MINIMUM ESSENTIEL de la hauteur de Zhang-Zagier (constante C_82). ================================================================================ Ne depend que de numpy et du fichier CERTIFICATE_C82_ZhangZagier_FINAL.json. Les parametres sont fixes ci-dessous : lancer le fichier tel quel. ess( h(alpha) + h(1-alpha) ) >= TARGET ATTENTION, PORTEE DE L'ENONCE. C'est un enonce sur le MINIMUM ESSENTIEL, non un enonce ponctuel. Au degre n la borne est diminuee de lambda_0 [ log(n)/n + D.pi/n**2 ] (Proposition 2.2 de la note) ; l'enonce essentiel suit par Northcott. Doche 2001 et Flammang 2018 prouvent des bornes PONCTUELLES, valides hors d'un ensemble exceptionnel fini : leurs enonces sont plus forts en nature, bien que numeriquement plus petits. -------------------------------------------------------------------------------- CE QUE CE SCRIPT ETABLIT, ET CE QU'IL SUPPOSE -------------------------------------------------------------------------------- SUPPOSE (demontre ailleurs, pas ici) : (H1) l'inegalite de certificat C_82 >= inf_z A(z) avec A(z) = g(z) - sum_k c_k log|Q_k(z)| - 2 l0 U_eta(z) + l0 I(eta), g(z) = log+|z| + log+|1-z|, valide des que c_k >= 0, l0 >= 0 et sum_k c_k deg Q_k + 2 l0 <= 2 ; (H2) le principe du minimum : A est surharmonique sur chacune des quatre regions decoupees par |z| = 1 et |z-1| = 1, donc son infimum est atteint sur la reunion des deux cercles unite. Les hypotheses numeriques de (H1) sont RE-VERIFIEES ci-dessous a partir du seul fichier de certificat. ETABLIT ICI, et c'est tout le travail du script : inf sur les deux cercles de A >= TARGET par des bornes d'intervalle rigoureuses et une bissection adaptative. -------------------------------------------------------------------------------- LES TROIS TERMES, ET LE SENS DANS LEQUEL CHACUN DOIT ETRE BORNE -------------------------------------------------------------------------------- Sur un intervalle I = [t-H, t+H] du parametre theta, z(theta) = shift + e(i theta), on a |z(theta) - z(t)| = 2|sin((theta-t)/2)| <= |theta - t| <= H. A = g - sum c log|Q| - 2 l0 U + l0 I. Pour MINORER A il faut donc MINORER g, MAJORER sum c log|Q|, MAJORER U, MINORER I(eta). (a) g MINOREE. log+ est croissante, donc g >= log+( max(0, |z(t)| - H) ) + log+( max(0, |1-z(t)| - H) ). (b) log|Q| MAJOREE. AUCUNE RACINE N'EST CALCULEE. Q est holomorphe, le disque de rayon H est convexe, donc par Taylor avec reste integral |Q(zeta)| <= |Q(z_t)| + |Q'(z_t)| d + ( max_disque |Q''| ) d**2 / 2, ou d majore |zeta - z_t|. La derivee est EVALUEE (Horner sur les coefficients entiers) ; seule la derivee SECONDE est bornee par la somme des valeurs absolues des coefficients, et elle entre avec d**2, donc son defaut est sans effet. Mesure : 93 pour cent de l'erreur totale provient des deux polynomes de degre 2 en z, pour lesquels la borne est EXACTE (derivee constante, reste nul). La constante d'erreur est donc irreductible : seul H la fait baisser, par bissection. DEUX ROUTES, choisies par le degre : - route z : d = H, directement. Pour les bas degres. - route w : w = z - z**2. Les polynomes de Flammang sont des polynomes ENTIERS en w (verifie exactement ci-dessous par reexpansion), de degre moitie. |dw/dtheta| = |1 - 2z|, EVALUE, donc d = (|1-2z_t| + 2H) H. Indispensable aux hauts degres : en z leur borne du second ordre explose (degre double, |z| atteignant 2 sur le cercle |z-1| = 1). (c) U_eta MAJOREE. u_j(rho) = log rho si rho >= r_j, sinon log r_j + (rho**2/r_j**2 - 1)/2 ; croissante et continue, donc U_eta <= sum_j m_j u_j( |z_t - c_j| + H ). (d) I(eta) MINOREE : calculee en long double par blocs de lignes, puis diminuee de 1e-12. -------------------------------------------------------------------------------- BUDGET D'ERREUR D'ARRONDI, EXPLICITE -------------------------------------------------------------------------------- Arithmetique IEEE double, eps = 2**-53 = 1.11e-16 ; on suppose seulement, comme il est standard, que +,-,*,/ sont correctement arrondis et que sqrt, exp et log ont au plus 1 ulp d'erreur. (i) z et les distances |z - w| : erreur relative <= 4 eps, donc absolue <= 1.8e-15 car toutes les distances utiles sont <= 4. La demi-largeur est gonflee de HINF = 1e-13, soit 55 fois cette borne : la monotonie garantit alors que la quantite majoree l'est pour tout theta de I. (ii) evaluations de Horner. L'erreur d'un Horner de degre d en arithmetique COMPLEXE (multiplication complexe : erreur relative <= sqrt(2) gamma(2)) est majoree par gamma(4d) sum_k |a_k| |x|**k avec gamma(n) = n eps/(1 - n eps). Cette borne est calculee et AJOUTEE au majorant de |Q| et de |Q'| : c'est le seul sens sur : si l'evaluation derape pres d'une racine, le majorant reste un majorant. Mesure au minimiseur : 6e-10 en valeur relative sur le pire polynome, donc sans effet. (iii) la somme ponderee de (c), 2157 termes, produit matriciel (BLAS) : erreur <= 2157 eps sum_j m_j |u_j| <= 2157 eps * 7 < 2e-12, multipliee par 2 l0 < 0.18. (iv) I(eta) en long double (eps = 1.08e-19) puis DIMINUEE de 1e-12, ce qui majore de plusieurs ordres l'erreur de la somme sur 2157**2 termes. Total majore par 1.1e-12. On retranche EPS = 1e-9 partout, soit 900 fois cette borne, et EPS ne consomme que 0.2 pour cent de la marge disponible. L'enonce est donc complet, sans reserve d'arrondi. -------------------------------------------------------------------------------- AUTO-TEST -------------------------------------------------------------------------------- Le script REFUSE DE DEMARRER si l'un des controles d'entete echoue. Le plus important : sur un echantillon aleatoire d'intervalles, la minoration doit etre inferieure ou egale a A evaluee directement en des points interieurs. Ce test attrape toute erreur de sens sur l'un des quatre termes. COURSE DE REFERENCE (ce fichier, tel quel, sur le certificat FINAL) : cercle |z| = 1 : 1 143 150 intervalles valides, 14 niveaux cercle |z-1| = 1 : 1 079 510 intervalles valides, 14 niveaux total 2 222 660 intervalles (143 secondes sur la machine de reference). auto-test : 0 violation, pire ecart minoration moins A = -3.98e-08. I(eta) recalcule a +1.34e-18 de la valeur du certificat. David Niedbala Giraudin, ORCID 0009-0009-1526-1178. Technical assistance: Claude (Anthropic). ================================================================================ """ import json import math import os import sys import time import numpy as np # ------------------------------------------------------------------ parametres CERT_FILE = "CERTIFICATE_C82_ZhangZagier_FINAL.json" TARGET = float(0.2497815) # cible a prouver N0 = int(2000) # intervalles initiaux par cercle EPS = float(1e-9) # allocation d'arrondi, retranchee partout HINF = float(1e-13) # gonflement de la demi-largeur MAXLEV = int(60) # niveaux de bissection au plus MAXINT = int(8000000) # intervalles simultanes au plus BLK = int(800) # intervalles par bloc (memoire) DCH = int(2000) # disques par bloc (memoire) WDEG_MIN = int(2) # route w des que le degre en w atteint ce seuil SEED = int(20260812) NTEST = int(60) # points de l auto-test par cercle et par demi-largeur DBL_EPS = float(2.0) ** int(-53) if len(sys.argv) > 1: CERT_FILE = sys.argv[1] if len(sys.argv) > 2: TARGET = float(sys.argv[2]) if len(sys.argv) > 3: N0 = int(sys.argv[3]) T0 = time.time() def say(*a): print(*a, flush=True) # --------------------------------------------------- polynomes a coeff entiers def poly_divmod_int(P, D): """Division euclidienne de P par D, listes DESCENDANTES d'entiers. D[0] doit diviser exactement a chaque etape. Rend (Q, R) ou (None, None).""" P = [int(x) for x in P] dq = len(P) - len(D) if dq < 0: return [0], P Q = [0] * (dq + 1) for i in range(dq + 1): c = P[i] if c % D[0] != 0: return None, None c = c // D[0] Q[i] = c for j in range(len(D)): P[i + j] -= c * D[j] return Q, P[dq + 1:] def strip_lead(P): P = list(P) while len(P) > 1 and P[0] == 0: P = P[1:] return P W_DESC = [-1, 1, 0] # w = z - z**2, en liste descendante def to_w(P_desc): """P(z) = sum_j b_j w**j avec w = z - z**2. Rend b ASCENDANT, ou None. Division synthetique entiere repetee : P = w Q + b_0, puis recursion sur Q.""" b = [] P = strip_lead(P_desc) for _ in range(64): if len(P) <= 1: b.append(int(P[0]) if P else 0) return b Q, R = poly_divmod_int(P, W_DESC) if Q is None: return None while len(R) < 2: R = [0] + R if R[0] != 0: return None # reste non constant : pas un polynome en w b.append(int(R[1])) P = strip_lead(Q) return None def poly_mul_int(a, b): r = [0] * (len(a) + len(b) - 1) for i, x in enumerate(a): for j, y in enumerate(b): r[i + j] += x * y return r def expand_from_w(b_asc): """Reexpansion exacte de sum_j b_j w**j en z, liste descendante.""" acc = [0] powr = [1] for bj in b_asc: term = [bj * c for c in powr] n = max(len(acc), len(term)) acc = [0] * (n - len(acc)) + acc term = [0] * (n - len(term)) + term acc = [x + y for x, y in zip(acc, term)] powr = poly_mul_int(powr, W_DESC) return strip_lead(acc) def horner(coef_desc, x): v = np.zeros(np.shape(x), dtype=complex) for a in coef_desc: v = v * x + float(a) return v def deriv_desc(coef_desc): d = len(coef_desc) - 1 return [coef_desc[i] * (d - i) for i in range(d)] if d >= 1 else [0] def gamma_n(n): n = float(n) return float(1.02) * n * DBL_EPS / max(float(1e-12), float(1.0) - n * DBL_EPS) def abs_coef_sum(coef_desc, rad, power_shift=0): """sum_k |a_k| rad**(k+power_shift), majoree (facteur 1.01 pour l'arrondi).""" d = len(coef_desc) - 1 s = np.zeros(np.shape(rad)) for i, a in enumerate(coef_desc): k = d - i s = s + abs(float(a)) * rad ** float(max(0, k + power_shift)) return float(1.01) * s def second_deriv_bound(coef_desc, rad): """majorant de |Q''| sur le disque de rayon rad : sum k(k-1)|a_k| rad**(k-2).""" d = len(coef_desc) - 1 s = np.zeros(np.shape(rad)) for i, a in enumerate(coef_desc): k = d - i if k >= 2: s = s + float(k * (k - 1)) * abs(float(a)) * rad ** float(k - 2) return float(1.01) * s # ------------------------------------------------------------------- lecture if not os.path.exists(CERT_FILE): say("ERREUR : fichier de certificat introuvable : " + str(CERT_FILE)) sys.exit(2) J = json.load(open(CERT_FILE)) L0 = float(J["lambda0"]) ACT = [(float(p["multiplier"]), [int(x) for x in p["coeffs_desc"]]) for p in J["polynomials"]] CEN = (np.array(J["eta"]["centres_re"], dtype=float) + 1j * np.array(J["eta"]["centres_im"], dtype=float)) RAD = np.array(J["eta"]["radii"], dtype=float) MAS = np.array(J["eta"]["masses"], dtype=float) NDISC = len(CEN) say("=" * 78) say("PROUVEUR C_82 -- minimum essentiel de la hauteur de Zhang-Zagier") say("=" * 78) say("certificat : " + str(CERT_FILE)) say(" " + str(NDISC) + " disques, " + str(len(ACT)) + " polynomes, lambda_0 = " + ("%.12f" % L0)) say(" cible TARGET = " + repr(TARGET) + " (certificat : " + str(J.get("proved_lower_bound")) + ")") # ------------------------------------------- controles des hypotheses de (H1) say("") say("--- hypotheses du certificat, relues depuis le seul fichier ---") ok = True cmin = min(c for c, _ in ACT) budget = sum(c * (len(q) - 1) for c, q in ACT) + 2.0 * L0 msum = float(MAS.sum()) say(" tous les multiplicateurs >= 0 : " + str(cmin >= 0.0) + " (min = " + ("%.3e" % cmin) + ")") say(" lambda_0 >= 0 : " + str(L0 >= 0.0)) say(" budget sum c deg Q + 2 l0 : " + ("%.12f" % budget) + " <= 2 : " + str(budget <= 2.0)) say(" somme des masses : " + ("%.15f" % msum)) from fractions import Fraction msum_exact = sum(Fraction(float(x)) for x in J["eta"]["masses"]) say(" somme EXACTE des masses - 1 : " + ("%.3e" % float(msum_exact - 1))) say(" masses > 0 et rayons > 0 : " + str(bool((MAS > 0).all() and (RAD > 0).all()))) DENS = float((MAS / (math.pi * RAD * RAD)).max()) say(" densite max D = max m/(pi r^2): " + ("%.6f" % DENS) + " D.pi < 599 : " + str(DENS * math.pi < 599.0)) if not ((MAS > 0).all() and (RAD > 0).all() and DENS * math.pi < 599.0 and abs(float(msum_exact - 1)) < 1e-15): ok = False if cmin < 0.0 or L0 < 0.0 or budget > 2.0 or abs(msum - 1.0) > 1e-12: ok = False marge = np.inf for s in range(0, NDISC, 256): D = np.abs(CEN[s:s + 256, None] - CEN[None, :]) S = RAD[s:s + 256, None] + RAD[None, :] for i in range(D.shape[0]): D[i, s + i] = np.inf marge = min(marge, float((D - S).min())) say(" disques disjoints : marge = " + ("%+.3e" % marge) + " " + str(marge > 0.0)) if not (marge > 0.0): ok = False # ------------------------------------------------ decomposition en w, exacte say("") say("--- decomposition en w = z - z**2, exacte, verifiee par reexpansion ---") ROUTES = [] for c, q in ACT: b = to_w(q) route = "z" if b is not None and (len(b) - 1) >= WDEG_MIN: back = expand_from_w(b) if back != strip_lead(q): say(" ECHEC de reexpansion sur " + str(q)) ok = False route = "w" elif b is not None: back = expand_from_w(b) if back != strip_lead(q): say(" ECHEC de reexpansion sur " + str(q)) ok = False if route == "w": bd = list(reversed(b)) # descendant en w ROUTES.append(("w", c, bd, deriv_desc(bd))) else: qd = strip_lead(q) ROUTES.append(("z", c, qd, deriv_desc(qd))) nw = sum(1 for r in ROUTES if r[0] == "w") say(" " + str(nw) + " polynomes par la route w, " + str(len(ROUTES) - nw) + " par la route z") for kind, c, cd, _ in ROUTES: say(" route " + kind + " degre " + str(len(cd) - 1).rjust(2) + " multiplicateur " + ("%.4e" % c)) # ------------------------------------------------------ I(eta), borne basse say("") Cl = CEN.astype(np.clongdouble) Rl = RAD.astype(np.longdouble) Ml = MAS.astype(np.longdouble) acc = np.longdouble(0.0) for s in range(0, NDISC, 256): D = np.abs(Cl[s:s + 256, None] - Cl[None, :]).astype(np.longdouble) K = np.where(D > 0, np.log(np.where(D > 0, D, np.longdouble(1.0))), np.longdouble(0.0)) for i in range(D.shape[0]): K[i, s + i] = np.log(Rl[s + i]) - np.longdouble(0.25) acc = acc + Ml[s:s + 256] @ K @ Ml IE = float(acc) IE_LO = IE - float(1e-12) say("--- I(eta) ---") say(" I(eta) recalcule = " + ("%.15e" % IE)) say(" certificat = " + ("%.15e" % float(J["I_eta"])) + " ecart " + ("%+.3e" % (IE - float(J["I_eta"])))) say(" borne inferieure = " + ("%.15e" % IE_LO)) if not ok: say("") say("AUTO-TEST D'ENTETE ECHOUE -- arret.") sys.exit(3) # ----------------------------------------------------------------- le noyau def u_of(rho, r): return np.where(rho >= r, np.log(np.maximum(rho, float(1e-300))), np.log(r) + (rho * rho / (r * r) - float(1.0)) / float(2.0)) def logQ_upper(zt, Hh): """Majorant de sum_k c_k log|Q_k(z)| pour |z - zt| <= Hh (route z) ou |w - w(zt)| <= (|1-2 zt| + 2 Hh) Hh (route w).""" out = np.zeros(np.shape(zt), dtype=float) azt = np.abs(zt) wt = zt * (float(1.0) - zt) awt = np.abs(wt) jac = np.abs(float(1.0) - float(2.0) * zt) + float(2.0) * Hh dlt = jac * Hh for kind, c, cd, dd in ROUTES: if kind == "z": x, ax, step = zt, azt, Hh else: x, ax, step = wt, awt, dlt qv = np.abs(horner(cd, x)) qd = np.abs(horner(dd, x)) eq = gamma_n(4 * (len(cd) - 1)) * abs_coef_sum(cd, ax) ed = gamma_n(4 * max(1, len(dd) - 1)) * abs_coef_sum(dd, ax) b2 = second_deriv_bound(cd, ax + step) maj = (qv + eq) + step * (qd + ed) + step * step * b2 / float(2.0) out = out + c * np.log(np.maximum(maj, float(1e-300))) return out def U_upper(zt, Hh): U = np.zeros(np.shape(zt), dtype=float) for s in range(0, NDISC, DCH): cc = CEN[s:s + DCH] rr = RAD[s:s + DCH] mo = MAS[s:s + DCH] dd = np.abs(zt[:, None] - cc[None, :]) + Hh[:, None] U = U + u_of(dd, rr[None, :]) @ mo return U def lower_bound(t, H, shift): out = np.empty(np.shape(t), dtype=float) for s in range(0, len(t), BLK): tt = t[s:s + BLK] Hh = H[s:s + BLK] + HINF zt = float(shift) + np.exp(1j * tt) g = (np.log(np.maximum(float(1.0), np.maximum(float(0.0), np.abs(zt) - Hh))) + np.log(np.maximum(float(1.0), np.maximum(float(0.0), np.abs(float(1.0) - zt) - Hh)))) out[s:s + BLK] = (g - logQ_upper(zt, Hh) - float(2.0) * L0 * U_upper(zt, Hh) + L0 * IE_LO - EPS) return out def A_exact(t, shift): """A evaluee directement, pour l'auto-test et pour l'information.""" out = np.empty(np.shape(t), dtype=float) for s in range(0, len(t), BLK): tt = t[s:s + BLK] z = float(shift) + np.exp(1j * tt) a = (np.log(np.maximum(float(1.0), np.abs(z))) + np.log(np.maximum(float(1.0), np.abs(float(1.0) - z)))) for kind, c, cd, _ in ROUTES: x = z if kind == "z" else z * (float(1.0) - z) a = a - c * np.log(np.maximum(np.abs(horner(cd, x)), float(1e-300))) U = np.zeros(np.shape(tt), dtype=float) for u0 in range(0, NDISC, DCH): cc = CEN[u0:u0 + DCH] rr = RAD[u0:u0 + DCH] mo = MAS[u0:u0 + DCH] U = U + u_of(np.abs(z[:, None] - cc[None, :]), rr[None, :]) @ mo out[s:s + BLK] = a - float(2.0) * L0 * U + L0 * IE return out # ------------------------------------------------------------- AUTO-TEST cle say("") say("--- auto-test : la minoration doit etre <= A en tout point de I ---") rng = np.random.default_rng(SEED) bad = 0 worst = -np.inf for shift in (float(0.0), float(1.0)): tc = rng.uniform(float(0.0), float(2.0) * math.pi, NTEST) for Hval in (float(1e-2), float(1e-3), float(1e-5), float(1e-7)): Hs = np.full(tc.shape, Hval) lb = lower_bound(tc, Hs, shift) for frac in (float(0.0), float(0.5), float(-0.5), float(1.0), float(-1.0)): av = A_exact(tc + frac * Hval, shift) d = lb - av worst = max(worst, float(d.max())) bad += int((d > float(1e-12)).sum()) say(" violations (minoration > A au point) : " + str(bad)) say(" pire ecart minoration - A : " + ("%+.3e" % worst) + " (doit etre <= 0)") if bad > 0: say("") say("AUTO-TEST ECHOUE -- la minoration n'en est pas une. Arret.") sys.exit(4) say(" auto-test PASSE") # ------------------------------------------------------------- la bissection def prove_circle(shift): e = np.linspace(float(0.0), float(2.0) * math.pi, N0 + 1) lo, hi = e[:-1].copy(), e[1:].copy() done = 0 for lev in range(MAXLEV): t = (lo + hi) / float(2.0) H = (hi - lo) / float(2.0) b = lower_bound(t, H, shift) good = b >= TARGET done += int(good.sum()) bad_m = ~good if not bad_m.any(): say(" cercle " + str(shift) + " : CERTIFIE. " + str(done) + " intervalles valides, " + str(lev + 1) + " niveaux (" + ("%.0f" % (time.time() - T0)) + "s)") return True lo2, hi2 = lo[bad_m], hi[bad_m] wmax = float((hi2 - lo2).max()) if wmax < float(1e-13): say(" cercle " + str(shift) + " : ECHEC. pire minoration = " + ("%.10f" % float(b[bad_m].max())) + " a largeur " + ("%.2e" % wmax)) return False md = (lo2 + hi2) / float(2.0) lo = np.concatenate([lo2, md]) hi = np.concatenate([md, hi2]) say(" niv " + str(lev).rjust(3) + " : " + str(len(lo)).rjust(9) + " en attente, largeur max " + ("%.2e" % float((hi - lo).max())) + ", pire " + ("%.10f" % float(b[bad_m].max())) + " (" + ("%.0f" % (time.time() - T0)) + "s)") if len(lo) > MAXINT: say(" ABANDON : plus de " + str(MAXINT) + " intervalles.") return False return False say("") say("--- preuve, cercle |z| = 1 puis |z-1| = 1 ---") all_ok = True for shift in (float(0.0), float(1.0)): all_ok = all_ok and prove_circle(shift) say("") if all_ok: bar = " " + "#" * 72 say(bar) say(" ### PROUVE : inf sur les deux cercles de A >= " + repr(TARGET)) say(" ###") say(" ### ess( h(alpha) + h(1-alpha) ) >= " + repr(TARGET)) say(" ### soit zeta(alpha) >= " + ("%.7f" % math.exp(TARGET))) say(" ###") say(" ### MINIMUM ESSENTIEL, non ponctuel. Au degre n la borne est") say(" ### diminuee de lambda_0 [ log(n)/n + D.pi/n**2 ] (Prop. 2.2).") say(" ### Doche 2001 et Flammang 2018 prouvent des bornes PONCTUELLES,") say(" ### plus fortes en nature bien que numeriquement plus petites.") say(" ###") say(" ### Record publie : 0.2487457988") say(" ### (V. Flammang, Int. J. Number Theory 14 (2018) 2663-2671,") say(" ### Theoreme 1, zeta >= 1.282416).") say(bar) else: say(" NON PROUVE a la cible " + repr(TARGET) + ".") say(" temps total " + ("%.0f" % (time.time() - T0)) + "s")