#!/usr/bin/env python3 """CY3 threshold banana (1,1,1,1,16): evaluate the closed-form connection coefficients c_{5/4}, c_{7/4}. Closed forms established in this work (2026-06-28): c_{5/4} = -1/(sqrt(2*pi)*Gamma(1/4)^2) c_{7/4} = -5*Gamma(1/4)^2/(384*sqrt(2)*pi^(5/2)) c_{5/4}*c_{7/4} = 5/(768*pi^3) (exact) Reference literals below are COPIED from the certified numerical records of this work (Arb ball-certified runs and independent-route gate comparisons); this script only re-evaluates the closed forms with mpmath and prints digit agreement, recomputed live as -log10 of the relative difference. Evaluation interface -------------------- f(dps=320) -> dict {'c54','c74','prod','c32'} of mpmath mpf at dps CLI: python3 cy3-banana-evaluate.py [--dps N] Domain: the shipped closed result is a pair of connection CONSTANTS at the fixed threshold coalescence point sqrt(M5) = 4 = 1+1+1+1 of the masses-squared family (1,1,1,1,16), plus the K3 (1,1,1,9) control constant c_{3/2}. There is no free kinematic variable in this artifact — dps is the only input, and the closed forms evaluate to ARBITRARY precision (more Gamma/pi digits; no stored float caps them). Only the copied oracle STRINGS have finite length; their caps are printed explicitly. The full kinematic period functions of the family are campaign transport objects and are NOT shipped here. Requires only mpmath (pip install mpmath). Runtime: ~1 s at the default dps=320 gate demo including the 2x-dps doubling check; wall time printed. Derivation tiers (2026-09-11; the default path above is unchanged) ------------------------------------------------------------------- python3 cy3-banana-evaluate.py --derive [--bits B] the DIRECT LINEAR decomposition of the transported (1,1,1,1,16) holomorphic period onto the full local Frobenius basis at s = 0, built from the units of direct_linear_extract.py served beside this file (the operator of record read by pin and its exact annihilation of the multinomial-squared series asserted, the exact seed jet with the rigorous tail bound, the ball-arithmetic Taylor transport at B and at B/2 + 100 bits along the path of record, the exact-Q Frobenius blocks, one 6x6 solve). Prints ALL SIX coordinates, each with the two-precision agreement: alpha_1, alpha_2 (the analytic coefficients in the canonical basis -- Phi_1 = t + O(t^3), Phi_2 = t^2 (1 + O(t)), t = e^{-i pi} s; the resonant rho = 1 block converted by x_{s^2} = a_1(0) x_{1,0} + a_1'(0) x_{1,1} + x_{2,0}, x_{s^2 log s} = a_1(0) x_{1,1} + x_{2,1}), the two logarithmic coefficients ell_1, ell_2 (zero), and c_{5/4}, c_{7/4} against the closed forms above (recomputed live) and the Arb literals; then the K3 (1,1,1,9) control in the same code path (A_0 vs sqrt(3) log(24)/(12 pi), A_1 vs -sqrt(3)/(24 pi), A_2 = 0, c_{3/2} vs -sqrt(3)/(36 pi)). B >= 1400, the bit floor of direct_linear_extract.py (default 1400; fewer bits are refused, exit 2). python3 cy3-banana-evaluate.py --hankel [--dps N] [--bits B] the analytic side: threshold_hankel_tail.py --check --family CY3 (the (1,1,1,1,16) tuple: the exact Hankel-tail objects, the closed STRUCTURE of c_{5/4}, c_{7/4} compared exactly with the record, the record strings at their certified lengths, the Frobenius-ratio identity c_{alpha+1}/c_alpha = phi_{alpha,1}), then the Bessel moments I_F = int_0^oo y F dy and M_F = int_0^oo y ln(y) F dy, F = prod_i J_0(m_i y), at two (Y, K, dps) settings ((36, 64, N-12), (46, 90, N); N >= 50, default 60): for (1,1,1,1,16) the vanishing first moment int_0^oo y J_0(y)^4 J_0(4y) dy = 0 and the log-moment identity alpha_1^(s) = int_0^oo y ln(y) J_0(y)^4 J_0(4y) dy against the direct linear solve's alpha_1^(s) (the --derive code path at B bits); for the (1,1,1,9) control I_F = sqrt(3)/(12 pi) and M_F = -sqrt(3)(ln 12 + gamma_E)/(12 pi), with a planted target (ln 13 for ln 12) that must fail by name. Digits of agreement printed at each setting; the bar is 30. The two (Y, K) quadrature settings, (36, 64) and (46, 90), are fixed, so the agreement figures are set by that tail truncation: raising --dps above the default buys no further digits. python3 cy3-banana-evaluate.py --planted direct_linear_extract.py's served planted-operator control: the theta^2 z^1 coefficient of the operator of record shifted by +1, on (1,1,1,1,16) and on (1,1,1,9); the exact annihilation must FAIL by name; exit 3 when it does. The tiers read five files served beside this one, each pinned in VENDORED_PINS by sha256 (a missing file or another sha256 is refused by name, exit 4, before any import): direct_linear_extract.py, threshold_hankel_tail.py, calpha_rings.py, fixtures/direct_linear_fixtures.json, fixtures/hankel_tail_fixtures.json -- the same bytes the threshold-banana bundle serves (one analytic derivation, two bundles). The two modules pin their own inputs again (exit 3 on a mismatch). --derive, --hankel and --planted need python-flint (pip install python-flint) and numpy beside mpmath; the default needs mpmath only. Tier exit codes: 0 pass; 1 a named FAIL; 2 usage; 3 the planted control refused as planted; 4 a vendored file absent or off its pin; 5 python-flint absent. Wall times are printed. """ import argparse import hashlib import importlib.util import math import os import sys import time from fractions import Fraction from mpmath import mp, mpf, sqrt, pi, gamma, log, fabs DEFAULT_DPS = 320 # ---- reference values (copied literal STRINGS, NOT computed here) ---------- # Kept as strings and parsed at the working precision of each use, so no # module-level mpf() call can silently cap them (mp.mpf import-dps footgun). # Route-B Arb ball midpoints (Julia/Arb, z-chart, 1280-bit run of this work; # they agree with the closed forms at the recorded 295.47/295.13 d) REF_C54 = "-0.0303492466882296557301753163706229879005161050829311035232435345318470894858801077932057795752827576466108783708670915808100848637145028098185258535596997066664622708287607349490105140666249302728790072874269914098711256890305915377922537622701370117823999170612808741599469304641894544350236363096393766640875879653" REF_C74 = "-0.006918488931750284187386003874378333926188476030768363894798534336415339266092596593361322224742229001215813306151159038027483581915347906368270661140931980574497977154961383970242892527210559993091610240259155726942669335071958135043526924289292275729419475616152845485861612430740148963024916492260815871702614548" # Route A independent numeric extraction (Python/mpmath transport, s-chart, # dps=80, s_dec=0.25), computed in this work ROUTEA_C54 = "-0.0303492466882296557301753163706229879005161050829311035232435345318470894769" ROUTEA_C74 = "-0.00691848893175028418738600387437833392618847603076836389479853433641533928446" # K3 positive control c_{3/2} = -sqrt(3)/(36*pi), Route A value # (computed in this work, same code path as the CY3 extraction) ROUTEA_C32 = "-0.0153146915394942235975368471753602622610385133435177959110935005636397306519" def agree_digits(a, b): """Matching decimal digits between a and b (relative, live -log10).""" if a == b: return mp.dps return int(-log(fabs((a - b) / b), 10)) def sig_digits(s): """Significant digits carried by a stored decimal literal string.""" return len(s.replace("-", "").replace(".", "").lstrip("0")) def f(dps=DEFAULT_DPS): """Evaluate the closed-form constants of this artifact at precision dps. Returns {'c54', 'c74', 'prod', 'c32'} as mpmath mpf. Domain note: these are constants at the fixed threshold point of the (1,1,1,1,16) family — there is no kinematic argument to vary (see module docstring); any dps is supported, limited only by mpmath. """ with mp.workdps(dps): g14 = gamma(mpf(1) / 4) c54 = -1 / (sqrt(2 * pi) * g14 ** 2) c74 = -5 * g14 ** 2 / (384 * sqrt(2) * pi ** mpf("2.5")) prod = mpf(5) / (768 * pi ** 3) c32 = -sqrt(3) / (36 * pi) return {"c54": +c54, "c74": +c74, "prod": +prod, "c32": +c32} def main(dps=DEFAULT_DPS): t_start = time.time() with mp.workdps(dps): g14 = gamma(mpf(1) / 4) g34 = gamma(mpf(3) / 4) c54 = -1 / (sqrt(2 * pi) * g14 ** 2) c54b = -(g34 ** 2) / (2 * sqrt(2) * pi ** mpf("2.5")) # equivalent form c74 = -5 * g14 ** 2 / (384 * sqrt(2) * pi ** mpf("2.5")) prod = mpf(5) / (768 * pi ** 3) c32 = -sqrt(3) / (36 * pi) ref_c54 = mpf(REF_C54) ref_c74 = mpf(REF_C74) ra_c54 = mpf(ROUTEA_C54) ra_c74 = mpf(ROUTEA_C74) ra_c32 = mpf(ROUTEA_C32) show = lambda x: mp.nstr(x, 60) print("CY3 threshold banana (1,1,1,1,16) — closed-form evaluation") print("=" * 64) print("working precision dps =", dps) print("c_{5/4} closed form :", show(c54)) print(" alt (Gamma(3/4)) :", show(c54b), " | forms agree:", agree_digits(c54, c54b), "d") print(" vs 198d Arb ref :", agree_digits(c54, ref_c54), "d") print(" vs Route A (dps80):", agree_digits(c54, ra_c54), "d", " (recorded cross-route gate: 72.53 d)") print() print("c_{7/4} closed form :", show(c74)) print(" vs 197d Arb ref :", agree_digits(c74, ref_c74), "d") print(" vs Route A (dps80):", agree_digits(c74, ra_c74), "d", " (recorded cross-route gate: 71.02 d)") print() print("product c54*c74 :", show(c54 * c74)) print(" 5/(768 pi^3) :", show(prod), " | identity:", agree_digits(c54 * c74, prod), "d") print() print("K3 control c_{3/2} :", show(c32)) print(" vs Route A (K3) :", agree_digits(c32, ra_c32), "d", " (gate record: 74.8 d)") # ---- dps-doubling check (gate point: c_{5/4}), all recomputed live ---- print() print("dps-doubling check (gate point c_{5/4}; agreements recomputed live)") for d in (dps, 2 * dps): with mp.workdps(d): g14d = gamma(mpf(1) / 4) c54d = -1 / (sqrt(2 * pi) * g14d ** 2) c74d = -5 * g14d ** 2 / (384 * sqrt(2) * pi ** mpf("2.5")) a_id = agree_digits(c54d * c74d, mpf(5) / (768 * pi ** 3)) a_ref = agree_digits(c54d, mpf(REF_C54)) print(" dps=%4d : identity c54*c74 vs 5/(768 pi^3) : %4d d (live, uncapped)" % (d, a_id)) print(" vs stored Arb oracle string : %4d d" % a_ref) print(" stored-oracle-string cap: the copied Arb c_{5/4} literal carries") print(" %d significant digits with recorded accuracy 295.47 d vs the" % sig_digits(REF_C54)) print(" closed form, so the vs-string row saturates near 295 d; the") print(" identity row keeps growing with dps (closed forms are uncapped).") print() print("wall time: %.2f s" % (time.time() - t_start)) # ---- derivation tiers (2026-09-11): --derive / --hankel / --planted --------------------------------------------- # The tiers import five files served BESIDE this one -- the direct linear extraction module, the Hankel-tail module, # its ring table and the two fixtures files -- each pinned here by sha256. A missing file or a file whose sha256 is # not the pin is refused by name (exit 4) before anything is imported. The default path above never touches them. VENDORED_PINS = { "direct_linear_extract.py": "601735c6fa4e1b322eced2bede662eb839f000a02c7b11689e0678cf10418d07", "fixtures/direct_linear_fixtures.json": "c3bb23092d1367f9eb94f0b250c3a93752593b2ea4f391d11e9f04aa28908e22", "threshold_hankel_tail.py": "350d3808cd85efe7f47c17c3bb1f71679e20be54e8c1702ec3fdf83489954bec", "fixtures/hankel_tail_fixtures.json": "3aca127c00cb0dd9e4dad3a441c5c801f3ce5a95ab816a4cbff6e3f2be783997", "calpha_rings.py": "f31b186255d2aa1f41166382ed00881d927eafde7c82ded8c9a30d14a745ce46", } DERIVE_BITS_DEFAULT = 1400 # the bit floor of direct_linear_extract.py (its --direct-linear refuses fewer) HANKEL_DPS_DEFAULT = 60 EXIT_TIER_FAIL, EXIT_TIER_USAGE, EXIT_PLANTED, EXIT_VENDOR, EXIT_DEPENDENCY = 1, 2, 3, 4, 5 MOMENT_BAR_D = 30 # the named bar of the --hankel moment identities (two settings must each clear it) def _here(): return os.path.dirname(os.path.abspath(__file__)) def _sha256_file(path): h = hashlib.sha256() with open(path, "rb") as fh: for chunk in iter(lambda: fh.read(1 << 20), b""): h.update(chunk) return h.hexdigest() def check_vendored(quiet=False): """Every file of VENDORED_PINS present beside this script at its pinned sha256, else exit 4 by name.""" here = _here() for rel, pin in VENDORED_PINS.items(): p = os.path.join(here, rel) if not os.path.isfile(p): print("REFUSED (exit %d): vendored file %s is absent beside %s" % (EXIT_VENDOR, rel, os.path.basename(__file__))) sys.exit(EXIT_VENDOR) got = _sha256_file(p) if got != pin: print("REFUSED (exit %d): vendored file %s sha256 %s... is not the pin %s..." % (EXIT_VENDOR, rel, got[:16], pin[:16])) sys.exit(EXIT_VENDOR) if not quiet: print("vendored files: %d present beside this script at their sha256 pins (%s)" % (len(VENDORED_PINS), ", ".join("%s %s" % (r, h[:16]) for r, h in VENDORED_PINS.items()))) def load_vendored_module(rel): """Import a pinned module from beside this script by explicit path (never through sys.path).""" path = os.path.join(_here(), rel) name = os.path.splitext(os.path.basename(rel))[0] try: spec = importlib.util.spec_from_file_location(name, path) mod = importlib.util.module_from_spec(spec) sys.modules[name] = mod spec.loader.exec_module(mod) except ImportError as exc: print("MISSING DEPENDENCY (exit %d): importing %s needs %s (the direct linear solve runs in python-flint ball " "arithmetic: pip install python-flint)" % (EXIT_DEPENDENCY, rel, exc)) sys.exit(EXIT_DEPENDENCY) return mod def _dig(a, b, cap): """matching decimal digits of a vs b (relative; absolute below 1 when b == 0), capped""" a, b = mp.mpf(a), mp.mpf(b) if b == 0: return float(cap) if a == 0 else min(float(cap), float(-mp.log10(fabs(a)))) if a == b: return float(cap) return min(float(cap), float(-mp.log10(fabs(a - b) / fabs(b)))) def derive_family(D, fixtures, name, bits, frac=0.4, ncheck=200, quiet=False): """The direct linear decomposition of the transported holomorphic period onto the FULL local Frobenius basis at s = 0 (every unipotent block coordinate and every fractional branch), built from direct_linear_extract.py's public units: the exact multinomial-squared series annihilated EXACTLY by the operator of record (fixtures, by pin), the exact seed jet with the rigorous tail bound, the acb Taylor transport along the path of record at BITS and at BITS/2 + 100, the exact-Q Frobenius blocks and branches, one r x r acb solve. Returns a dict; prints progress.""" from flint import acb, arb, acb_mat, fmpq, ctx fam = fixtures["families"][name] msq = tuple(fam["msq"]) Pj = {int(j): [int(v) for v in co] for j, co in fam["operator"]["theta_form_Pj"].items()} por = fam["path_of_record"] sb, lift, sdec = Fraction(str(por["sb"])), Fraction(str(por["lift"])), Fraction(str(por["sdec"][0])) order = max(Pj) degz = max(len(co) - 1 for co in Pj.values()) m = [math.isqrt(M) for M in msq] thr = Fraction(sum(m) ** 2) precs = [bits, bits // 2 + 100] rep = {"family": name, "msq": list(msq), "order": order, "bits": precs, "sb": str(sb), "lift": str(lift), "sdec": str(sdec)} Nser = int(bits * math.log(2) / math.log(float(sb / thr))) + 60 c_def = D.multinomial_squared_series(msq, ncheck + degz + 5) a_alt = [(-1) ** n * c_def[n] for n in range(len(c_def))] res = D.annihilation_residuals(Pj, a_alt, ncheck) nz = sum(1 for x in res if x != 0) rep["annihilation"] = {"ncheck": ncheck, "nonzero_residuals": nz} if not quiet: print(" %s msq=%s: operator of record order %d, z-degree %d (%s %s...); exact annihilation of the " "multinomial-squared series over %d coefficients: %s" % (name, list(msq), order, degz, os.path.basename(fam["operator"]["source"]["path"]), fam["operator"]["source"]["sha256"][:16], ncheck, "residual EXACTLY 0" if nz == 0 else "%d NONZERO residuals" % nz)) if nz: rep["fail"] = "the operator of record does not annihilate the period series" return rep cser = D.extend_by_recurrence(Pj, c_def, max(Nser, ncheck) + 5) R = D.indicial_threshold(Pj, order) exps = D.exponents_from_indicial(R[0]) rep["threshold_exponents"] = {str(k): v for k, v in sorted(exps.items())} p_int = D.s_chart_dform(Pj, order) import numpy as _np zroots = _np.roots(list(reversed([float(cf) for cf in Pj[order]]))) sing_s = [0j] + [complex(-1.0 / zr) for zr in zroots if abs(zr) > 1e-300] if not quiet: print(" threshold exponents (indicial polynomial at s = 0): %s" % rep["threshold_exponents"]) results = {} for prec in precs: ctx.prec = prec T = D.Transport(p_int, prec, sing=sing_s) jet0, tails = D.seed_jet(cser[:Nser + 1], thr, sb, order, prec) q = lambda fr: arb(fmpq(fr.numerator, fr.denominator)) wps = [acb(q(sb), arb(0)), acb(q(sb), q(lift)), acb(q(sdec), q(lift)), acb(q(sdec), arb(0))] t1 = time.time() jet, nsteps = T.path(wps, jet0, frac=frac) t_tr = time.time() - t1 Nthr = int(prec * math.log(2) / -math.log(float(sdec) / 4.0)) + 40 cols, labels, blocks = [], [], {} for rho, mult in sorted(exps.items()): if rho.denominator == 1: aa = D.frobenius_block(R, rho, mult, Nthr, extra=2) J = D.jets_block(aa, rho, mult, sdec, order, prec) blocks[str(rho)] = {"mult": mult, "a_n_eps_head": {str(n): [str(x) for x in aa[n][:mult + 1]] for n in range(min(3, len(aa)))}} for j in range(mult): cols.append(J[j]) labels.append("unip rho=%s j=%d" % (rho, j)) else: aa = D.frobenius_fractional(R, rho, Nthr) cols.append(D.jets_fractional(aa, rho, sdec, order, prec)) labels.append("frac alpha=%s" % rho) A = acb_mat(order, order) for col in range(order): for row in range(order): A[row, col] = cols[col][row] bvec = acb_mat(order, 1) for row in range(order): bvec[row, 0] = jet[row] x = A.solve(bvec) ndig = int(prec * 0.30103) - 8 out = {} for i, lab in enumerate(labels): xi = x[i, 0] up = lambda b: b.abs_upper().str(8, more=True, radius=False) # a rigorous upper bound for |.| of the solve's ball ent = {"x_re": xi.real.str(ndig, radius=False), "x_im": xi.imag.str(ndig, radius=False), "abs_upper": up(abs(xi)), "im_abs_upper": up(xi.imag)} if lab.startswith("frac"): al = Fraction(lab.split("=")[1]) ph = (acb(0, 1) * acb.pi() * acb(fmpq(al.numerator, al.denominator))).exp() cv = xi * ph ent.update({"alpha": str(al), "c_re": cv.real.str(ndig, radius=False), "c_im": cv.imag.str(ndig, radius=False), "c_im_abs_upper": up(cv.imag)}) out[lab] = ent results[prec] = {"labels": labels, "coords": out, "nsteps": nsteps, "transport_s": round(t_tr, 2), "Nser": Nser, "Nthr": Nthr, "seed_tails": ["%.2e" % t for t in tails], "truncation_envelope_max": "%.2e" % getattr(T, "envelope_max", 0.0), "blocks": blocks} if not quiet: print(" %d bits: transport %d steps in %.1f s along the path of record [s_b=%s, +i %s, s_dec=%s] " "(N_series %d, N_threshold %d; seed tail bounds %s; largest per-step truncation envelope %s, midpoint " "arithmetic: an estimate, the two-precision agreement below is the error measure); basis %s" % (prec, nsteps, t_tr, sb, lift, sdec, Nser, Nthr, results[prec]["seed_tails"], results[prec]["truncation_envelope_max"], labels)) rep["results"] = results return rep def _mpc_of(ent, key_re, key_im): return mp.mpc(mp.mpf(ent[key_re]), mp.mpf(ent[key_im])) def canonical_cy3(res): """The (1,1,1,1,16) coordinates in the CANONICAL local basis at s = 0 (each solution zero on the other leading monomials s, s log s, s^2, s^2 log s; s-chart, reached from Im s > 0). The resonant rho = 1 block of multiplicity two [eps^j] s^{1+eps} sum a_n(eps) s^n carries a_1(eps) = a_1(0) + a_1'(0) eps + ...; the rho = 2 block starts at s^2. So: alpha_1 = x_{1,0} (on s), ell_1 = x_{1,1} (on s log s), alpha_2 = a_1(0) x_{1,0} + a_1'(0) x_{1,1} + x_{2,0} (on s^2), ell_2 = a_1(0) x_{1,1} + x_{2,1} (on s^2 log s); c_{5/4}, c_{7/4} = x_alpha e^{i pi alpha} (t = e^{-i pi} s, the lower branch of record).""" co = res["coords"] a1 = [Fraction(v) for v in res["blocks"]["1"]["a_n_eps_head"]["1"]] a10 = mp.mpf(a1[0].numerator) / a1[0].denominator a11 = mp.mpf(a1[1].numerator) / a1[1].denominator x10 = _mpc_of(co["unip rho=1 j=0"], "x_re", "x_im") x11 = _mpc_of(co["unip rho=1 j=1"], "x_re", "x_im") x20 = _mpc_of(co["unip rho=2 j=0"], "x_re", "x_im") x21 = _mpc_of(co["unip rho=2 j=1"], "x_re", "x_im") up = lambda lab, k="abs_upper": mp.mpf(co[lab][k]) # rigorous upper bounds carried from the ball solve return {"alpha1_s": x10, "ell1": x11, "alpha2_s": a10 * x10 + a11 * x11 + x20, "ell2": a10 * x11 + x21, "c54": _mpc_of(co["frac alpha=5/4"], "c_re", "c_im"), "c74": _mpc_of(co["frac alpha=7/4"], "c_re", "c_im"), "a1_0": a1[0], "a1_1": a1[1], "ell1_bound": up("unip rho=1 j=1"), "ell2_bound": fabs(a10) * up("unip rho=1 j=1") + up("unip rho=2 j=1"), "alpha1_im": up("unip rho=1 j=0", "im_abs_upper"), "alpha2_im": fabs(a10) * up("unip rho=1 j=0", "im_abs_upper") + fabs(a11) * up("unip rho=1 j=1", "im_abs_upper") + up("unip rho=2 j=0", "im_abs_upper"), "c54_im": up("frac alpha=5/4", "c_im_abs_upper"), "c74_im": up("frac alpha=7/4", "c_im_abs_upper")} def tbasis_k3(res): """The (1,1,1,9) unipotent triple in the t = e^{-i pi} s basis: B_j = the coordinate on [eps^j] s^{1+eps} sum a_n(eps) s^n, A_2 = -B_2, A_1 = -B_1 + i pi A_2, A_0 = -B_0 + i pi A_1 + (pi^2/2) A_2; c_{3/2} = x_{3/2} e^{3 i pi/2}.""" co = res["coords"] B = [_mpc_of(co["unip rho=1 j=%d" % j], "x_re", "x_im") for j in range(3)] A2 = -B[2] A1 = -B[1] + mp.mpc(0, 1) * pi * A2 A0 = -B[0] + mp.mpc(0, 1) * pi * A1 + (pi ** 2 / 2) * A2 return {"A0": A0, "A1": A1, "A2": A2, "c32": _mpc_of(co["frac alpha=3/2"], "c_re", "c_im"), "A2_bound": mp.mpf(co["unip rho=1 j=2"]["abs_upper"])} def fixture_literals_check(fx_direct, fx_hankel): """The fixtures carry copies of this file's literals labelled 'cy3-banana-evaluate.REF_C54' / '.REF_C74' / '.ROUTEA_C54' / '.ROUTEA_C74' (copied from an earlier state of this file, named there by that state's sha256); the copies must equal this file's literals byte for byte -- checked by VALUE here, whatever this file's sha256.""" want = {"cy3-banana-evaluate.REF_C54": REF_C54, "cy3-banana-evaluate.REF_C74": REF_C74, "cy3-banana-evaluate.ROUTEA_C54": ROUTEA_C54, "cy3-banana-evaluate.ROUTEA_C74": ROUTEA_C74} seen, bad = 0, [] for tag, fx in (("direct_linear_fixtures", fx_direct), ("hankel_tail_fixtures", fx_hankel)): for alpha, lst in fx["families"]["CY3"].get("record_strings", {}).items(): for rs in lst: if rs["label"] in want: seen += 1 if rs["value"] != want[rs["label"]]: bad.append("%s %s [%s]" % (tag, alpha, rs["label"])) print("fixtures' copies of REF_C54 / REF_C74 / ROUTEA_C54 / ROUTEA_C74 vs this file's literals: %d copies, %s" % (seen, "all EQUAL" if not bad and seen else "DIFFER at " + ", ".join(bad) if bad else "none found")) return seen > 0 and not bad def print_derive(rep3, repk, bits): """All six (1,1,1,1,16) coordinates + the (1,1,1,9) control, each with the two-precision agreement; returns fails.""" hi, lo = rep3["bits"] fails = [] dshow = 60 with mp.workdps(int(hi * 0.30103)): H, L = canonical_cy3(rep3["results"][hi]), canonical_cy3(rep3["results"][lo]) cap = lo * 0.30103 vals = f(int(hi * 0.30103)) print() print("(1,1,1,1,16): the six coordinates of the transported holomorphic period on the local Frobenius basis at s = 0") print(" canonical basis: alpha_1 on s, ell_1 on s log s, alpha_2 on s^2, ell_2 on s^2 log s (s-chart, from Im s > 0);") print(" t-basis (t = e^{-i pi} s, Phi_1 = t + O(t^3), Phi_2 = t^2 (1 + O(t))): alpha_1^(t) = -alpha_1^(s), alpha_2^(t) = alpha_2^(s);") print(" rho = 1 block of record: a_1(0) = %s, a_1'(0) = %s (exact)" % (H["a1_0"], H["a1_1"])) for key, lab in (("alpha1_s", "alpha_1^(s) [on s] "), ("alpha2_s", "alpha_2^(s) [on s^2] ")): d2 = _dig(mp.re(L[key]), mp.re(H[key]), cap) imb = H[key.replace("_s", "_im")] print(" %s = %s" % (lab, mp.nstr(mp.re(H[key]), dshow))) print(" |Im| <= %s; two-precision (%d vs %d bits) %.1f d" % (mp.nstr(imb, 3), hi, lo, d2)) if d2 < 100: fails.append("%s two-precision %.1f d < 100" % (key, d2)) print(" alpha_1^(t) = %s ; alpha_2^(t) = %s" % (mp.nstr(-mp.re(H["alpha1_s"]), 30), mp.nstr(mp.re(H["alpha2_s"]), 30))) for key, lab in (("ell1", "ell_1 [on s log s] "), ("ell2", "ell_2 [on s^2 log s]")): bh, bl = H[key + "_bound"], L[key + "_bound"] zh, zl = _dig(bh, 0, hi * 0.30103), _dig(bl, 0, cap) print(" %s : |value| <= %s at %d bits (zero to %.1f d), <= %s at %d bits (zero to %.1f d) -- expected 0 " "(bounds from the ball solve; the transport is midpoint arithmetic, so the %d-bit line is the measure)" % (lab, mp.nstr(bh, 3), hi, zh, mp.nstr(bl, 3), lo, zl, lo)) if min(zh, zl) < 100: fails.append("%s not zero to 100 d (%.1f)" % (key, min(zh, zl))) for key, lab, cf, ref, refname in (("c54", "c_{5/4}", vals["c54"], REF_C54, "REF_C54 (%d-digit Arb literal)" % sig_digits(REF_C54)), ("c74", "c_{7/4}", vals["c74"], REF_C74, "REF_C74 (%d-digit Arb literal)" % sig_digits(REF_C74))): v = mp.re(H[key]) d2 = _dig(mp.re(L[key]), v, cap) dcf = _dig(v, cf, hi * 0.30103) dref = _dig(v, mp.mpf(ref), sig_digits(ref)) print(" %s = %s" % (lab, mp.nstr(v, dshow))) print(" |Im| <= %s; two-precision %.1f d; vs the closed form (this file's f(), live) %.1f d; vs %s %.1f d" % (mp.nstr(H[key + "_im"], 3), d2, dcf, refname, dref)) if d2 < 100 or dcf < 100: fails.append("%s: two-precision %.1f d / closed form %.1f d below 100" % (lab, d2, dcf)) hik, lok = repk["bits"] with mp.workdps(int(hik * 0.30103)): KH, KL = tbasis_k3(repk["results"][hik]), tbasis_k3(repk["results"][lok]) capk = lok * 0.30103 s3 = sqrt(3) A0cf, A1cf, c32cf = s3 * log(24) / (12 * pi), -s3 / (24 * pi), -s3 / (36 * pi) print() print("(1,1,1,9) control, same code path: the unipotent triple in the t-basis and c_{3/2}") for key, lab, cf, cfs in (("A0", "A_0", A0cf, "sqrt(3) log(24)/(12 pi)"), ("A1", "A_1", A1cf, "-sqrt(3)/(24 pi)"), ("c32", "c_{3/2}", c32cf, "-sqrt(3)/(36 pi)")): v = mp.re(KH[key]) d2 = _dig(mp.re(KL[key]), v, capk) dcf = _dig(v, cf, hik * 0.30103) extra = "; vs ROUTEA_C32 %.1f d" % _dig(v, mp.mpf(ROUTEA_C32), sig_digits(ROUTEA_C32)) if key == "c32" else "" print(" %s = %s ; two-precision %.1f d; vs %s %.1f d%s" % (lab, mp.nstr(v, 50), d2, cfs, dcf, extra)) if d2 < 100 or dcf < 100: fails.append("K3 %s: two-precision %.1f d / closed form %.1f d below 100" % (lab, d2, dcf)) za, zb = _dig(KH["A2_bound"], 0, hik * 0.30103), _dig(KL["A2_bound"], 0, capk) print(" A_2 : |value| <= %s at %d bits (zero to %.1f d), <= %s at %d bits (zero to %.1f d) -- the log^2 coefficient, zero by parity" % (mp.nstr(KH["A2_bound"], 3), hik, za, mp.nstr(KL["A2_bound"], 3), lok, zb)) if min(za, zb) < 100: fails.append("K3 A_2 not zero to 100 d (%.1f)" % min(za, zb)) return fails # ---- the Bessel moments behind the unipotent data (the --hankel identities) ------------------------------------ # I_F = int_0^inf y F(y) dy and M_F = int_0^inf y ln(y) F(y) dy, F(y) = prod_i J_0(m_i y): the head int_0^Y by # tanh-sinh on steps of 1/2; the tail int_Y^inf from the K-term Hankel expansion of every J_0 factor, the 2^N sign # patterns sigma (k_sigma = sum sigma_i m_i) integrated in closed form (k = 0) or on the rotated ray y = Y +- i u # (k != 0). The k = 0 group carries conjugate-pair cancellations that are EXACT zeros (parity): a combined tail # coefficient is skipped as zero when it sits at roundoff level RELATIVE TO THE PRE-CANCELLATION MAGNITUDE AT THAT # ORDER j (per order -- not relative to the maximum over all orders, which the factorially growing high orders # dominate and which would drop the whole non-oscillatory tail). Two (Y, K, dps) settings give the digits by agreement. def _pattern_polys(TH, m, K): from itertools import product base = {} for mi in set(m): for sg in (+1, -1): ser = [] for j in range(K + 1): ph = [(1, 0), (0, 1), (-1, 0), (0, -1)][j % 4] if sg < 0: ph = (ph[0], -ph[1]) c = TH.hankel_a(j) / Fraction(mi) ** j ser.append((ph[0] * c, ph[1] * c)) base[(mi, sg)] = ser out = [] for sig in product((+1, -1), repeat=len(m)): ser = [(Fraction(1), Fraction(0))] + [(Fraction(0), Fraction(0))] * K for mi, sg in zip(m, sig): fac = base[(mi, sg)] new = [(Fraction(0), Fraction(0))] * (K + 1) for i, (xa, xb) in enumerate(ser): if xa == 0 and xb == 0: continue for j in range(K + 1 - i): ya, yb = fac[j] na, nb = new[i + j] new[i + j] = (na + xa * ya - xb * yb, nb + xa * yb + xb * ya) ser = new out.append((sum(s * mi for s, mi in zip(sig, m)), -sum(sig), ser)) return out def _tail_integrals(TH, m, Y, K, with_log): N = len(m) prodm = 1 for mi in m: prodm *= mi pref = mp.mpf(2) ** (-N) * (2 / pi) ** (mp.mpf(N) / 2) / sqrt(prodm) Y = mp.mpf(Y) total = mp.mpc(0) groups = {} for k, P, ser in _pattern_polys(TH, m, K): groups.setdefault(k, []).append((P, ser)) for k, lst in groups.items(): comb = [mp.mpc(0)] * (K + 1) comb_abs = [mp.mpf(0)] * (K + 1) # per-order magnitude BEFORE the conjugate-pair cancellation for P, ser in lst: ph = mp.expj(pi * P / 4) for j, (a, b) in enumerate(ser): if a or b: term = ph * (mp.mpf(a.numerator) / a.denominator + mp.mpc(0, 1) * (mp.mpf(b.numerator) / b.denominator)) comb[j] += term comb_abs[j] += fabs(term) if k == 0: for j, cj in enumerate(comb): if fabs(cj) <= comb_abs[j] * mp.mpf(10) ** (-(mp.dps - 10)): continue # an exact zero of the combined tail coefficient at THIS order p = mp.mpf(N) / 2 + j - 1 assert p > 1, "k=0 tail term y^-%s with a nonzero coefficient: the moment diverges" % p if with_log: val = Y ** (1 - p) * (log(Y) / (p - 1) + 1 / (p - 1) ** 2) else: val = Y ** (1 - p) / (p - 1) total += cj * val else: sgn = 1 if k > 0 else -1 kk = abs(k) def g(u, comb=comb, sgn=sgn, k=k, kk=kk): y = Y + sgn * mp.mpc(0, 1) * u poly = mp.mpc(0) yinv = 1 / y pw = mp.mpc(1) for cj in comb: poly += cj * pw pw *= yinv val = y ** (1 - mp.mpf(N) / 2) * poly * mp.expj(k * Y) * mp.exp(-kk * u) * (sgn * mp.mpc(0, 1)) if with_log: val *= log(y) return val Lr = mp.mpf(40) / kk + 2 total += mp.quad(g, mp.linspace(0, Lr * 3, 13)) + mp.quad(g, [Lr * 3, mp.inf]) return pref * total def _head_integral(m, Y, with_log, h=0.5, maxdeg=8): def fint(y): v = y if with_log: v = y * log(y) if y > 0 else mp.mpf(0) for mi in m: v *= mp.besselj(0, mi * y) return v hh = mp.mpf(h) pts = [mp.mpf(i) * hh for i in range(int(Y / h) + 1)] if pts[-1] < Y: pts.append(mp.mpf(Y)) tot = mp.mpf(0) for a_, b_ in zip(pts[:-1], pts[1:]): tot += mp.quad(fint, [a_, b_], maxdegree=maxdeg) return tot def bessel_moment(TH, m, Y, K, dps, with_log): """(value, |Im residual of the tail|, wall s) of I_F (with_log False) or M_F (True) at the setting (Y, K, dps).""" t0 = time.time() with mp.workdps(dps): hd = _head_integral(m, Y, with_log) tl = _tail_integrals(TH, m, Y, K, with_log) return hd + mp.re(tl), fabs(mp.im(tl)), time.time() - t0 def _load_tier_inputs(): check_vendored() D = load_vendored_module("direct_linear_extract.py") TH = load_vendored_module("threshold_hankel_tail.py") fx_d, sha_d = D.load_fixtures(D.FIXTURES_DEFAULT) # the module's own pin (exit 3 on a mismatch) fx_h, sha_h = TH.load_fixtures(TH.FIXTURES_DEFAULT, False) # likewise print("direct_linear_extract.py STAMP %s (fixtures %s...); threshold_hankel_tail.py STAMP %s (fixtures %s...); python-flint %s; mpmath %s" % (D.STAMP, sha_d[:16], TH.STAMP, sha_h[:16], __import__("flint").__version__, __import__("mpmath").__version__)) ok_lit = fixture_literals_check(fx_d, fx_h) return D, TH, fx_d, fx_h, ok_lit def run_derive_tier(bits): t0 = time.time() print("CY3 threshold banana (1,1,1,1,16) -- --derive: the direct linear decomposition at s = 0 (%d and %d bits)" % (bits, bits // 2 + 100)) print("=" * 100) D, TH, fx_d, fx_h, ok_lit = _load_tier_inputs() print("\n-- (1,1,1,1,16)") rep3 = derive_family(D, fx_d, "CY3", bits) print("\n-- (1,1,1,9) control") repk = derive_family(D, fx_d, "K3", bits) fails = [] for r in (rep3, repk): if r.get("fail"): fails.append("%s: %s" % (r["family"], r["fail"])) if not fails: fails = print_derive(rep3, repk, bits) if not ok_lit: fails.append("the fixtures' copies of this file's literals (REF_C54 / REF_C74 / ROUTEA_C54 / ROUTEA_C74) differ from them") print() print("DERIVE VERDICT: %s (bars: every two-precision and closed-form agreement >= 100 d, the log coefficients zero to >= 100 d)" % ("PASS" if not fails else "FAIL")) for fl in fails: print(" FAIL:", fl) print("wall time: %.2f s" % (time.time() - t0)) return 0 if not fails else EXIT_TIER_FAIL def run_hankel_tier(dps, bits): t0 = time.time() if dps < 50: print("usage: --hankel --dps N needs N >= 50 (threshold_hankel_tail.py prints the analytic value beside the record at >= 50 digits)") return EXIT_TIER_USAGE settings = [(36, 64, dps - 12), (46, 90, dps)] print("CY3 threshold banana (1,1,1,1,16) -- --hankel: the Hankel-tail derivation and the Bessel moments (settings (Y, K, dps) = %s)" % settings) print("=" * 100) D, TH, fx_d, fx_h, ok_lit = _load_tier_inputs() print("\n-- threshold_hankel_tail.py --check --family CY3 --dps %d (the (1,1,1,1,16) tuple: the exact tail objects, the closed " "STRUCTURE vs the record, the record strings at their certified lengths, the Frobenius-ratio identity)" % dps) rc_th = TH.main(["--check", "--family", "CY3", "--dps", str(dps)]) print("\n-- alpha_1^(s) by the direct linear solve (the --derive code path, %d bits) for the log-moment identity" % bits) rep3 = derive_family(D, fx_d, "CY3", bits, quiet=True) fails = [] if not rep3.get("fail") else ["CY3: " + rep3["fail"]] alpha1 = a1_2p = None if not fails: hi, lo = rep3["bits"] with mp.workdps(int(hi * 0.30103)): H, L = canonical_cy3(rep3["results"][hi]), canonical_cy3(rep3["results"][lo]) alpha1 = mp.re(H["alpha1_s"]) a1_2p = _dig(mp.re(L["alpha1_s"]), alpha1, lo * 0.30103) print(" alpha_1^(s) = %s (two-precision %d vs %d bits: %.1f d)" % (mp.nstr(alpha1, 50), hi, lo, a1_2p)) print("\n-- the Bessel moments I_F = int_0^oo y F(y) dy and M_F = int_0^oo y ln(y) F(y) dy, F = prod_i J_0(m_i y)") fams = (("K3", (1, 1, 1, 3)), ("CY3", (1, 1, 1, 1, 4))) vals = {} for tag, m in fams: for wl in (False, True): key = "%s_F[%s m=%s]" % ("M" if wl else "I", tag, m) vals[(tag, wl)] = [] for (Y, K, dp) in settings: v, imres, w = bessel_moment(TH, m, Y, K, dp, wl) vals[(tag, wl)].append((v, dp)) with mp.workdps(dp): print(" %s (Y=%d, K=%d, dps=%d): %s (tail Im residual %s) [%.0f s]" % (key, Y, K, dp, mp.nstr(v, 42), mp.nstr(imres, 3), w)) sys.stdout.flush() topd = max(s[2] for s in settings) with mp.workdps(topd + 10): s3 = sqrt(3) tI, tM, tA0 = s3 / (12 * pi), -s3 * (log(12) + mp.euler) / (12 * pi), s3 * log(24) / (12 * pi) print("\n-- identities (digits of agreement at each setting; bar %d d)" % MOMENT_BAR_D) def line(name, pairs, target, absolute=False, bar=MOMENT_BAR_D, planted=False): ds = [] for (v, dp) in pairs: ds.append(_dig(v, target, dp) if not absolute else _dig(fabs(v - target), 0, dp)) verdict = all(d >= bar for d in ds) tag = ("FAIL (as planted)" if not verdict else "UNEXPECTED PASS") if planted else ("PASS" if verdict else "FAIL") print(" %s: %s -> %s" % (name, " / ".join("%.1f d" % d for d in ds), tag)) return verdict ok = [] ok.append(line("(1,1,1,9) I_F vs sqrt(3)/(12 pi) [= -2 A_1 = -3 c_{3/2}]", vals[("K3", False)], tI)) ok.append(line("(1,1,1,9) M_F vs -sqrt(3)(ln 12 + gamma_E)/(12 pi) [<=> A_0 = sqrt(3) log(24)/(12 pi) via A_0 = -M_F - (gamma_E - ln 2) I_F]", vals[("K3", True)], tM)) a0 = [(-vm - (mp.euler - log(2)) * vi, min(dm, di)) for (vm, dm), (vi, di) in zip(vals[("K3", True)], vals[("K3", False)])] ok.append(line("(1,1,1,9) A_0 from the moments vs sqrt(3) log(24)/(12 pi)", a0, tA0)) two = _dig(vals[("K3", True)][0][0], vals[("K3", True)][1][0], settings[0][2]) print(" (1,1,1,9) M_F two-setting agreement: %.1f d" % two) planted_fails = not line("(1,1,1,9) PLANTED: M_F vs -sqrt(3)(ln 13 + gamma_E)/(12 pi) [must FAIL by name]", vals[("K3", True)], -s3 * (log(13) + mp.euler) / (12 * pi), planted=True) ok.append(planted_fails) ok.append(line("(1,1,1,1,16) I_F vs 0 [the s log s coefficient ell_1 = I_F/2 vanishes] (absolute)", vals[("CY3", False)], 0, absolute=True)) if alpha1 is not None: ok.append(line("(1,1,1,1,16) M_F vs alpha_1^(s) of the direct linear solve [the log-moment identity alpha_1^(s) = int_0^oo y ln y J_0(y)^4 J_0(4y) dy]", vals[("CY3", True)], alpha1)) two3 = _dig(vals[("CY3", True)][0][0], vals[("CY3", True)][1][0], settings[0][2]) print(" (1,1,1,1,16) M_F two-setting agreement: %.1f d; M_F = %s" % (two3, mp.nstr(vals[("CY3", True)][1][0], 40))) if rc_th != 0: fails.append("threshold_hankel_tail.py --check --family CY3 returned %d" % rc_th) if not all(ok): fails.append("a moment identity below the bar (see the lines above)") if not ok_lit: fails.append("the fixtures' copies of this file's literals (REF_C54 / REF_C74 / ROUTEA_C54 / ROUTEA_C74) differ from them") print() print("HANKEL VERDICT: %s" % ("PASS" if not fails else "FAIL")) for fl in fails: print(" FAIL:", fl) print("wall time: %.2f s" % (time.time() - t0)) return 0 if not fails else EXIT_TIER_FAIL def run_planted_tier(): """direct_linear_extract.py's served planted-operator control (theta^2, z^1, +1): on the tampered operator the exact annihilation of the multinomial-squared series must FAIL by name -- applied to this bundle's (1,1,1,1,16) operator of record and to the (1,1,1,9) control's; exit 3 when both refuse as planted.""" t0 = time.time() print("CY3 threshold banana (1,1,1,1,16) -- --planted: the planted-operator control (the theta^2 z^1 coefficient of the operator of record + 1)") print("=" * 100) D, TH, fx_d, fx_h, ok_lit = _load_tier_inputs() as_planted = True for name in ("CY3", "K3"): rep, r = D.run_family(name, 600, fx_d, plant_operator=(2, 1, 1), quiet=True) ann = rep.get("annihilation", {}) refused = (r == D.EXIT_PIN and ann.get("verdict") == "FAIL") print(" %s %s planted (2,1,1): exact annihilation over %s coefficients -> %s nonzero residuals: %s" % (name, rep.get("msq"), ann.get("ncheck"), ann.get("nonzero_residuals"), "REFUSED by name (annihilation FAIL, the module's exit %d) as planted" % r if refused else "NOT as planted (rc %d)" % r)) as_planted = as_planted and refused print() print("PLANTED CONTROL: %s" % ("FAILED BY NAME as planted -> exit %d" % EXIT_PLANTED if as_planted else "UNEXPECTED -> exit %d" % EXIT_TIER_FAIL)) print("wall time: %.2f s" % (time.time() - t0)) return EXIT_PLANTED if as_planted else EXIT_TIER_FAIL if __name__ == "__main__": ap = argparse.ArgumentParser( description="CY3 threshold banana (1,1,1,1,16) closed-form connection " "coefficients. The result is a pair of constants at the " "fixed threshold point; --dps is the only free input of the " "default evaluation (no kinematic variable in this artifact). " "Derivation tiers: --derive (the direct linear decomposition at " "s = 0, all six coordinates), --hankel (the Hankel-tail derivation " "and the Bessel-moment identities), --planted (the planted-operator " "control); see the module docstring.") ap.add_argument("--dps", type=int, default=None, help="working decimal precision (default %d; with --hankel the top setting's dps, default %d, >= 50; " "the two (Y, K) quadrature settings are fixed, so a dps above the default buys --hankel no further digits)" % (DEFAULT_DPS, HANKEL_DPS_DEFAULT)) ap.add_argument("--derive", action="store_true", help="the direct linear decomposition of the transported period at s = 0 (all six coordinates + the K3 control)") ap.add_argument("--bits", type=int, default=DERIVE_BITS_DEFAULT, help="ball-arithmetic working precision of --derive / --hankel's direct linear solve (>= %d; default %d)" % (DERIVE_BITS_DEFAULT, DERIVE_BITS_DEFAULT)) ap.add_argument("--hankel", action="store_true", help="the Hankel-tail derivation (threshold_hankel_tail.py --check --family CY3) and the Bessel-moment identities at two fixed (Y, K) settings, (36, 64) and (46, 90)") ap.add_argument("--planted", action="store_true", help="the planted-operator control: the exact annihilation must FAIL by name (exit 3)") args = ap.parse_args() if (args.derive or args.hankel) and args.bits < DERIVE_BITS_DEFAULT: print("usage: --bits must be >= %d (the bit floor of direct_linear_extract.py, whose --direct-linear refuses fewer bits; " "the two-precision pair is B and B/2 + 100): exit %d" % (DERIVE_BITS_DEFAULT, EXIT_TIER_USAGE)) sys.exit(EXIT_TIER_USAGE) if args.planted: sys.exit(run_planted_tier()) if args.derive or args.hankel: rc = 0 if args.derive: rc = run_derive_tier(args.bits) if args.hankel: rc_h = run_hankel_tier(args.dps if args.dps is not None else HANKEL_DPS_DEFAULT, args.bits) rc = rc or rc_h sys.exit(rc) main(args.dps if args.dps is not None else DEFAULT_DPS)