""" VERIF_KORENBLUM_0446.py - verifier for the bound c2 >= 223/500 = 0.446 (Korenblum's constant of the Bergman space A2). Author: David Niedbala Giraudin (ORCID 0009-0009-1526-1178). Technical assistance: Claude (Anthropic). Ancillary file of the note "The Moebius pseudodistance of the annulus and an improved lower bound for Korenblum's constant" (Hexagon deposit). Run with plain Python 3 (numpy only; no scipy, no Sage). No argument, no external file. Running time: a few seconds. What the script checks, in rigorous arithmetic: Wikstrom's dual moment criterion (arXiv:2607.17748, Lemma 3.1) for mu = sum_j w_j delta_{t_j} + (Lebesgue measure on [R, 1]), R = 99/100, with beta(t) = 1 - 1/G(t), G = (1 + sqrt(1 + S(t)**2)) / 2, where S(t) = sinh(D/2), D = hyperbolic diameter (Moebius pseudodistance) of the circle |z| = sqrt(t) in the annulus c < |z| < 1, given by the product with positive factors (Theorem 2 of the note) S(t) = P(c) * 2 sqrt(t)/(1-t) * prod_n (1+c**(2n-1) t)(1+c**(2n-1)/t) / ((1-c**(2n) t)(1-c**(2n)/t)), P(c) = prod_n (1+c**(2n))**2 / (1-c**(2n-1))**2. One needs, for every k >= 1: A_k = sum_j w_j t_j**k - R**(k+1)/(k+1) > 0 B_k = R**(k+1)/(k+1) + int_R**1 t**k / G(t) dt - sum_j w_j beta(t_j) t_j**k > 0. Arithmetic: * S, G, beta: exact fixed-point integers (2**-PREC) with directed rounding, plus an explicit bound for the remainder of the infinite product. * moment sums: IEEE double floats (numpy), all terms positive, explicit relative error bound, widened by a safety factor. * integral: dyadic grid, midpoint rule (a lower bound for the convex t**k, Jensen's inequality), 1/G bounded below by its value at the right endpoint (S increasing for t > c: Theorem 3(iii) of the note). * calibration: the seven upper bounds for beta published by Wikstrom (section 5, c = 0.4263, Wang's constant) are reproduced first. Output messages are in French. The last line reads "RESULTAT : critere certifie en c = 223/500 = 0.446 ; donc c2 >= 0.446". """ import math, time, sys from fractions import Fraction as Fr import numpy as np T0 = time.time() def say(*a): print(*a, flush=True) def fail(msg): say("ECHEC :", msg) sys.exit(1) # ---------------- auto-test arithmetique ---------------- if not (2 ** 10 == 1024 and 7 // 2 == 3 and (-7) // 2 == -4 and math.isqrt(10 ** 20) == 10 ** 10 and float(Fr(1, 4)) == 0.25 and np.float64(3.0) * np.float64(0.5) == 1.5): fail("auto-test arithmetique (environnement non standard, par exemple preparseur Sage)") say("auto-test arithmetique : OK") # ---------------- donnees exactes ---------------- C = Fr(223, 500) R = Fr(99, 100) S2 = C * C ATOMS = [ (Fr(49729, 250000), Fr(42261186256165887, 100000000000000000)), (Fr(60895666661691, 100000000000000), Fr(3213347492292841, 25000000000000000)), (Fr(61180781079179, 100000000000000), Fr(5117313952474759, 25000000000000000)), (Fr(4187499456063, 5000000000000), Fr(30953435822431, 2500000000000000)), (Fr(41934234956647, 50000000000000), Fr(3274174589692609, 25000000000000000)), (Fr(18639525158413, 20000000000000), Fr(4984739013621153, 100000000000000000)), (Fr(93247222724023, 100000000000000), Fr(180033676111197, 25000000000000000)), (Fr(12117035664079, 12500000000000), Fr(107784329378477, 10000000000000000)), (Fr(4847931161323, 5000000000000), Fr(1227418071154803, 100000000000000000)), (Fr(98448736353123, 100000000000000), Fr(408308945472583, 50000000000000000)), (Fr(98482420908651, 100000000000000), Fr(14560557107887, 20000000000000000)), (Fr(49469854898371, 50000000000000), Fr(9044125245481, 6250000000000000)), (Fr(98955114839119, 100000000000000), Fr(213851621217, 6250000000000000)), (Fr(99, 100), Fr(16824391869377, 50000000000000000)), ] TJ = [a for a, b in ATOMS]; WJ = [b for a, b in ATOMS] if not (TJ[0] == S2 and TJ[-1] == R and sum(WJ) == R and all(w > 0 for w in WJ) and all(S2 <= t <= R for t in TJ) and all(x < y for x, y in zip(TJ, TJ[1:]))): fail("donnees des atomes") say(f"c = {C} = {float(C)}, R = {R}, {len(TJ)} atomes, masse totale 1 : OK") PREC = 200 ONE = 1 << PREC NTERMS = 30 # termes du produit ; reste borne explicitement def fx_dn(q): # Fraction -> entier virgule fixe, arrondi vers le bas return (q.numerator * ONE) // q.denominator def fx_up(q): return -((-q.numerator * ONE) // q.denominator) def mul_dn(a, b): return (a * b) >> PREC def mul_up(a, b): return -((-a * b) >> PREC) def div_dn(a, b): return (a << PREC) // b def div_up(a, b): return -((-a << PREC) // b) def sqrt_dn(a): return math.isqrt(a << PREC) def sqrt_up(a): s = math.isqrt(a << PREC) return s if s * s == (a << PREC) else s + 1 # puissances de c, encadrees CP_DN = [ONE]; CP_UP = [ONE] for j in range(1, 2 * NTERMS + 2): CP_DN.append(mul_dn(CP_DN[-1], fx_dn(C))) CP_UP.append(mul_up(CP_UP[-1], fx_up(C))) # constante P(c) majoree P_UP = ONE for n in range(1, NTERMS + 1): num = (ONE + CP_UP[2 * n]); num = mul_up(num, num) den = (ONE - CP_UP[2 * n - 1]); den = mul_dn(den, den) P_UP = div_up(mul_up(P_UP, num), den) # reste du produit infini (n > NTERMS), valable pour t >= c**2 : # log(facteur n) <= 2c**(2n) + 2c**(2n-1)/(1-c) + c**(2n-1)(1 + 1/c**2) + 2c**(2n-2)/(1-c**2) # somme geometrique en n de raison c**2 ; on majore exp(E) par 1 + 2E (E < 1). cq = C r2 = cq * cq first = NTERMS + 1 geo = 1 / (1 - r2) Erest = (2 * cq ** (2 * first) + 2 * cq ** (2 * first - 1) / (1 - cq) + cq ** (2 * first - 1) * (1 + 1 / r2) + 2 * cq ** (2 * first - 2) / (1 - r2)) * geo if not Erest < Fr(1, 10 ** 15): fail("reste du produit") TAIL_UP = fx_up(1 + 2 * Erest) def S_upper(t): """majorant de S(t) (t Fraction, c**2 < t < 1), en virgule fixe""" tu = fx_up(t); td = fx_dn(t) itu = fx_up(1 / t) v = div_up(mul_up(2 * ONE, sqrt_up(tu)), fx_dn(1 - t)) v = mul_up(v, P_UP) for n in range(1, NTERMS + 1): a = CP_UP[2 * n - 1]; ad = CP_DN[2 * n]; au = CP_UP[2 * n] num = mul_up(ONE + mul_up(a, tu), ONE + mul_up(a, itu)) den = mul_dn(ONE - mul_up(au, tu), ONE - mul_up(au, itu)) v = div_up(mul_up(v, num), den) return mul_up(v, TAIL_UP) def invG_lower(Su): """minorant de 1/G = 2/(1 + sqrt(1 + S**2)) a partir d'un majorant de S""" root = sqrt_up(ONE + mul_up(Su, Su)) return div_dn(2 * ONE, ONE + root) def to_float_dn(x): # entier virgule fixe -> flottant minorant f = float(Fr(x, ONE)) return math.nextafter(f, 0.0) def to_float_up(x): f = float(Fr(x, ONE)) return math.nextafter(f, math.inf) def frac_dn(q): return math.nextafter(float(q), 0.0) def frac_up(q): return math.nextafter(float(q), math.inf) # ---------------- 0. calibration sur les valeurs publiees de Wikstrom (c = 0.4263) ---------------- # avec la constante de Wang G_W = (1 + sqrt(1 + 4 S**2))/2 ; controle informatif en flottants def S_float(c, t, N=60): v = 2 * math.sqrt(t) / (1 - t) for n in range(1, N + 1): v *= (1 + c ** (2 * n)) ** 2 / (1 - c ** (2 * n - 1)) ** 2 v *= (1 + c ** (2 * n - 1) * t) * (1 + c ** (2 * n - 1) / t) / ((1 - c ** (2 * n) * t) * (1 - c ** (2 * n) / t)) return v PUB_T = [0.5786106116843371, 0.8163605367760075, 0.9200761835179418, 0.9633238309788337, 0.9813028389745606, 0.9881242650409469, 0.99] PUB_B = [0.989685929536115, 0.995533783756722, 0.998127951092143, 0.999157174953505, 0.999573932287887, 0.999730257446456, 0.999773067105237] dev = 0.0 for t, b in zip(PUB_T, PUB_B): s = S_float(0.4263, t) bw = 1 - 2 / (1 + math.sqrt(1 + 4 * s * s)) dev = max(dev, abs(bw - b)) say(f"calibration (Wikstrom, c = 0.4263, constante de Wang) : ecart maximal {dev:.2e} sur 7 valeurs publiees") if dev > 1e-9: fail("calibration") # ---------------- 1. beta aux noeuds ---------------- BETA_UP = [1.0] # convention beta = 1 en t = c**2 for t in TJ[1:]: g = invG_lower(S_upper(t)) BETA_UP.append(to_float_up(ONE - g)) say("beta majore aux noeuds :", ", ".join(f"{b:.10f}" for b in BETA_UP)) # ---------------- 2. grille dyadique de la queue [R, U] ---------------- SH = 51 # points multiples de 2**-51 : milieux exacts en double DEN = 1 << SH U = Fr(99999, 100000) X = [-((-R.numerator * DEN) // R.denominator)] # premier point >= R XU = (U.numerator * DEN) // U.denominator EPSN, EPSD = 1, 5000 # pas geometrique en 1 - t while X[-1] < XU: d = DEN - X[-1] step = max(1, (d * EPSN) // EPSD) X.append(min(X[-1] + step, XU)) NP = len(X) - 1 say(f"queue : {NP} morceaux dyadiques sur [{X[0] / DEN:.12f}, {X[-1] / DEN:.12f}]") G_LO = np.empty(NP); H = np.empty(NP); MID = np.empty(NP) for i in range(NP): G_LO[i] = to_float_dn(invG_lower(S_upper(Fr(X[i + 1], DEN)))) # 1/G decroissant : bord droit H[i] = (X[i + 1] - X[i]) / DEN # exact MID[i] = (X[i] + X[i + 1]) / (2 * DEN) # exact if i % 5000 == 0: say(f" {i} / {NP} morceaux encadres [{time.time() - T0:.0f} s]") if not (np.all(np.diff(G_LO) <= 0) and np.all(G_LO > 0)): fail("1/G doit etre positif et decroissant sur la grille") say(f" 1/G minore sur les {NP} morceaux [{time.time() - T0:.0f} s]") u = 2.0 ** -53 def lo(x, m): # minorant de la valeur exacte d'un calcul de m operations positives return x * (1 - 2.2 * m * u) def hi(x, m): return x * (1 + 2.2 * m * u) T_DN = np.array([frac_dn(t) for t in TJ]); T_UP = np.array([frac_up(t) for t in TJ]) W_DN = np.array([frac_dn(w) for w in WJ]); W_UP = np.array([frac_up(w) for w in WJ]) BU = np.array(BETA_UP) R_DN = frac_dn(R); R_UP = frac_up(R) NA = len(TJ) # ---------------- 3. premiere famille A_k ---------------- wl = WJ[-1] K1 = int(R / wl) + 1 if not wl * (K1 + 1) >= R: fail("argument du dernier atome") say(f"A_k : controle direct pour 1 <= k <= {K1} ; au-dela, w_dernier (k+1) >= R suffit") pw = np.ones(NA); pr = R_UP; minA = (math.inf, 0) for k in range(1, K1 + 1): pw = pw * T_DN; pr = pr * R_UP left = lo(float(np.sum(W_DN * pw)), k + NA + 2) right = hi(pr / (k + 1), k + 3) if not left > right: fail(f"A_{k}") rel = (left - right) / right if rel < minA[0]: minA = (rel, k) say(f" A_k > 0 pour tout k ; marge relative minimale {minA[0]:.4e} en k = {minA[1]} [{time.time() - T0:.0f} s]") # ---------------- 4. seconde famille B_k ---------------- A_ = Fr(998, 1000); B_ = Fr(999, 1000) iA = next(i for i in range(len(X)) if X[i] >= A_ * DEN) iB = max(i for i in range(len(X)) if X[i] <= B_ * DEN) a_f = X[iA] / DEN; b_f = X[iB] / DEN m_lo = to_float_dn(invG_lower(S_upper(Fr(X[iB], DEN)))) say(f" queue : sur [{a_f:.9f}, {b_f:.9f}], 1/G >= {m_lo:.7e}") pg = G_LO * H * MID # puissance k = 1 pw = np.ones(NA); pd = R_DN; minB = (math.inf, 0); k = 1; K2 = None while True: if k > 1: pg = pg * MID pw = pw * T_UP; pd = pd * R_DN integ = lo(float(np.sum(pg)), k + NP + 3) first_ = lo(pd / (k + 1), k + 3) atoms_ = hi(float(np.sum(W_UP * BU * pw)), k + NA + 3) Bk = first_ + integ - atoms_ if not Bk > 0: fail(f"B_{k}") rel = Bk / first_ if rel < minB[0]: minB = (rel, k) if k % 200 == 0: # minorant de m (b-a) a**k par multiplication repetee ak = 1.0 for _ in range(k): ak *= a_f tailL = lo(m_lo * (b_f - a_f) * ak, k + 3) tailR = hi(float(np.sum(W_UP * pw)), k + NA + 3) if tailL > tailR: K2 = k break say(f" B_k > 0 jusqu'a k = {k} [{time.time() - T0:.0f} s]") k += 1 if not (a_f > float(R)): fail("argument de queue") say(f" B_k > 0 pour 1 <= k <= {K2} ; marge relative minimale {minB[0]:.4e} en k = {minB[1]}") say(f" pour k >= {K2} : m (b-a) a**k > somme w_j t_j**k a k = {K2} ; a chaque pas le membre de gauche") say(f" est multiplie par a = {a_f:.6f} et chaque terme de droite par au plus R = 0.99 < a") say(f"RESULTAT : critere certifie en c = {C} = {float(C)} ; donc c2 >= 0.446 [{time.time() - T0:.0f} s]")