#!/usr/bin/env python3 r"""Anisotropic simple-cubic lattice Green function (Zucker's final problem) -- compute W_S(alpha,beta,gamma; w) at runtime at any admissible kinematic point, via the certified order-5 Picard-Fuchs / moment-recurrence route, checked live against an independent Bessel-Laplace quadrature. W_S(alpha,beta,gamma;w) = (1/pi^3) int_{[0,pi]^3} d^3k / (w - alpha cos k1 - beta cos k2 - gamma cos k3) = (1/w) sum_{m>=0} p_{2m}(A,B,C) / w^{2m}, A = alpha^2, B = beta^2, C = gamma^2. Self-contained (mpmath + standard library only). NOTHING numerical about the result is stored: every value printed is computed live, and the ONLY embedded data are exact integers/rationals: * the moment formula (exact, from the definition): p_{2m} = binom(2m,m)/4^m * sum_{i+j+k=m} [m!/(i!j!k!)]^2 A^i B^j C^k (only p_0..p_12 are ever computed this way at eval time -- they are the transport seeds -- plus further indices in --check); * the order-5 Picard-Fuchs operator as its equivalent (r,s)=(6,5) moment recurrence sum_{j=0..6} c_j(n;A,B,C) p_{2(n+j)} = 0, with c_j exact integer/rational polynomials, embedded verbatim below (L5_COEFFS). Provenance: built by 81-point exact-rational interpolation, independently confirmed by a GF(p) symbolic nullspace (kernel dim 1 at two primes), and CERTIFIED to annihilate the moment sequence identically in (A,B,C) by a creative-telescoping certificate (ct_certificate_symbolic.json, CERT_RESULT_2026-07-05.md) -- a theorem, not a fit. The leading coefficient c_6(n) = (n+4)(n+5)(n+6)^3 is CONSTANT in (A,B,C) and positive for n >= 0, so the forward transport is defined at every rate point and every n >= 0. Evaluation route (the "bootstrapped" route -- no DE solving, no quadrature): 1. seeds p_0..p_12 exact in Q from the moment formula; 2. exact Fraction transport p_{2m}, m <= N, by the certified recurrence (this is DE/PF-operator transport in the spectral variable: the recurrence IS the PF operator written on series coefficients); 3. W_S = (1/w) sum_m p_{2m}/w^{2m}, summed in mpmath at dps + guard (all terms exact rationals; the series is positive and monotone for real w > alpha+beta+gamma, so no cancellation), N accepted ONLY by the runtime-certified tail bound (see the fail-closed block below; the closed-form N(dps) formula is a starting seed, never the answer). Live oracle (independent route, shares no arithmetic with the above): W_S = int_0^inf e^{-w t} I_0(alpha t) I_0(beta t) I_0(gamma t) dt by tanh-sinh quadrature (mpmath), implemented here from the formula. Agreement at the reference points: 140.3-142.3 digits at five (w; rates) points not used in the construction, two rate families. Retained literals: NONE. There is no stored digit string anywhere in this file; every check below is live-vs-live at the precision you request. (Note: with no stored strings there is no fast-start cache and no separate --recompute path -- the compute path IS the only path.) Fail-closed certification: * Series truncation: the N(dps) formula in nterms() is DEMOTED to a starting seed. Acceptance is by the PROVEN a-priori tail bound p_{2m} = E[(alpha cos k1 + beta cos k2 + gamma cos k3)^{2m}] <= rho^{2m}, rho = sqrt A + sqrt B + sqrt C, so rel. tail <= q^{N+1}/(1-q), q = (rho/w)^2 < 1 (domain-enforced), plus the (N+3)*10^-(dps+19) positive-sum rounding term; W_series accepts only when this bound < 10^-(dps+TAIL_GUARD), escalates N x1.5 by EXACT continuation of the same recurrence, and RAISES at the cap (64x seed) rather than print an uncertified value. * Oracle Laplace truncation: I_0(x) <= cosh x <= e^x gives the proven tail bound e^(-gap*tmax)/gap * w for the [0,tmax] cutoff; enforced < 10^-(dps+TMAX_GUARD) at runtime (tmax x1.5 to 64x seed, then raise). * Oracle quadrature: accepted only when two evaluations at successive working precisions (genuinely different tanh-sinh node sets) agree to < 10^-(dps+ORA_GUARD) relative; the COARSER member of the first agreeing pair is returned (healthy runs: bit-identical to the seed evaluation); RuntimeError if the ladder never agrees. (Double-refinement semantics, reimplemented inline.) * Live agreement RAISES: in point mode the series-vs-oracle agreement is now a raising check, >= max(15, dps - 8) digits required (measured healthy: dps+16.8 worst near-threshold, dps+20..21 generic; >= 10^24 headroom -- a healthy run can never fire it). * Crank test: --check Check D recomputes the series value at dps and dps+40 (different certified depths N) and requires agreement >= dps. Interface: --dps N working precision (default 50); result precision grows with N (series + oracle are both recomputed at N). --point A B C W kinematic point: A,B,C = SQUARED hop rates (positive rationals, e.g. 2 3 5 for rates sqrt2,sqrt3,sqrt5), W = spectral variable (rational, > sqrt A+sqrt B+sqrt C). Default: 2 3 5 8 -- the paper's generic point (alpha,beta,gamma)=(sqrt2,sqrt3,sqrt5) at w=8. --check check mode: (i) exact-transport check -- recurrence- transported moments == direct-formula moments, Fraction equality, at three further rate points + a deep index; (ii) two-precision check -- series vs live oracle at dps AND 2*dps, agreement must be >= max(30, dps-12) and must GROW by >= 0.4*dps under doubling; (iii) kinematic points from a second rate family, checked live against the quadrature; (iv) crank: series value at dps vs dps+40 (independent certified depths) must agree to >= dps d. --no-oracle point mode: skip the live oracle (series value only; the certified series bound still applies). Measured walls (reference box, single core): transport N=400 ~0.07 s; oracle ~1 s at dps 60, ~3 s at dps 120 (x2 since 07-05 pm: every oracle call carries its two-depth verification pass); full --check --dps 60 still well under 2 min. Honest limitations (what this evaluator does NOT cover): * Domain: real w with w/(alpha+beta+gamma) >= 1.02 (enforced; the moment series converges for |w| > alpha+beta+gamma, but cost grows like 1/log(w/threshold) and the live oracle needs the Laplace integral to converge). The THRESHOLD point w* = alpha+beta+gamma -- where the Polya return probability lives -- is NOT computable by this script: the paper evaluates it by ODE-Richardson transport against a Hankel-tail-corrected quadrature, a 37-digit cross-check; the 140-digit agreements are the interior ones reproduced here. * Complex w, and w on/inside the branch cut [-(a+b+g), a+b+g], not covered. * A,B,C must be positive rationals (squared hop rates; the rates themselves may be irrational, as at the certified point). Degenerate symmetric slices (cubic A=B=C, tetragonal B=C) are fine -- the certified operator annihilates the moments for ALL (A,B,C), the minimal order merely drops there (5 -> 4 -> 3, the transcendental-rank ladder). Changelog: 2026-07-05 authored: exact recurrence transport + live Bessel-Laplace tanh-sinh check; L5_COEFFS embedded verbatim (certified 2026-07-05). 2026-07-05 hardening: every truncation (series length N, oracle Laplace cutoff tmax, tanh-sinh depth) demoted to a seed inside a refine-until-certified-bound fail-closed loop; proven tail bounds enforced at runtime; point-mode agreement promoted to a raising check; crank Check D (dps vs dps+40) wired into --check; certified BOUND lines printed. Values at default dps unchanged. """ import argparse import time from fractions import Fraction as F from math import comb import mpmath as mp # =========================================================================== # The certified (r,s)=(6,5) moment recurrence sum_j c_j(n;A,B,C) p_{2(n+j)}=0 # c_j(n) = sum_k L5_COEFFS[(j,k)] n^k; c_j homogeneous of degree 6-j in # (A,B,C); normalized c_{6,5}=1. Embedded verbatim. # =========================================================================== L5_COEFFS = { (0,0): '-315*(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)**3/16', (0,1): '-1371*(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)**3/16', (0,2): '-261*(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)**3/2', (0,3): '-177*(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)**3/2', (0,4): '-27*(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)**3', (0,5): '-3*(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)**3', (1,0): '35*(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)*(103*A**3 - 103*A**2*B - 103*A**2*C - 103*A*B**2 + 642*A*B*C - 103*A*C**2 + 103*B**3 - 103*B**2*C - 103*B*C**2 + 103*C**3)/8', (1,1): '(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)*(27541*A**3 - 27541*A**2*B - 27541*A**2*C - 27541*A*B**2 + 195114*A*B*C - 27541*A*C**2 + 27541*B**3 - 27541*B**2*C - 27541*B*C**2 + 27541*C**3)/24', (1,2): '(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)*(2303*A**3 - 2303*A**2*B - 2303*A**2*C - 2303*A*B**2 + 17970*A*B*C - 2303*A*C**2 + 2303*B**3 - 2303*B**2*C - 2303*B*C**2 + 2303*C**3)/2', (1,3): '4*(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)*(425*A**3 - 425*A**2*B - 425*A**2*C - 425*A*B**2 + 3561*A*B*C - 425*A*C**2 + 425*B**3 - 425*B**2*C - 425*B*C**2 + 425*C**3)/3', (1,4): '8*(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)*(17*A**3 - 17*A**2*B - 17*A**2*C - 17*A*B**2 + 150*A*B*C - 17*A*C**2 + 17*B**3 - 17*B**2*C - 17*B*C**2 + 17*C**3)', (1,5): '2*(A**2 - 2*A*B - 2*A*C + B**2 - 2*B*C + C**2)*(19*A**3 - 19*A**2*B - 19*A**2*C - 19*A*B**2 + 174*A*B*C - 19*A*C**2 + 19*B**3 - 19*B**2*C - 19*B*C**2 + 19*C**3)/3', (2,0): '-35*(1313*A**4 + 3852*A**3*B + 3852*A**3*C - 10330*A**2*B**2 - 1468*A**2*B*C - 10330*A**2*C**2 + 3852*A*B**3 - 1468*A*B**2*C - 1468*A*B*C**2 + 3852*A*C**3 + 1313*B**4 + 3852*B**3*C - 10330*B**2*C**2 + 3852*B*C**3 + 1313*C**4)/48', (2,1): '-(126817*A**4 + 103788*A**3*B + 103788*A**3*C - 461210*A**2*B**2 - 187692*A**2*B*C - 461210*A**2*C**2 + 103788*A*B**3 - 187692*A*B**2*C - 187692*A*B*C**2 + 103788*A*C**3 + 126817*B**4 + 103788*B**3*C - 461210*B**2*C**2 + 103788*B*C**3 + 126817*C**4)/48', (2,2): '-(7589*A**4 - 756*A**3*B - 756*A**3*C - 13666*A**2*B**2 - 13936*A**2*B*C - 13666*A**2*C**2 - 756*A*B**3 - 13936*A*B**2*C - 13936*A*B*C**2 - 756*A*C**3 + 7589*B**4 - 756*B**3*C - 13666*B**2*C**2 - 756*B*C**3 + 7589*C**4)/3', (2,3): '-(3375*A**4 - 1956*A**3*B - 1956*A**3*C - 2838*A**2*B**2 - 7348*A**2*B*C - 2838*A**2*C**2 - 1956*A*B**3 - 7348*A*B**2*C - 7348*A*B*C**2 - 1956*A*C**3 + 3375*B**4 - 1956*B**3*C - 2838*B**2*C**2 - 1956*B*C**3 + 3375*C**4)/3', (2,4): '-239*A**4 + 204*A**3*B + 204*A**3*C + 70*A**2*B**2 + 596*A**2*B*C + 70*A**2*C**2 + 204*A*B**3 + 596*A*B**2*C + 596*A*B*C**2 + 204*A*C**3 - 239*B**4 + 204*B**3*C + 70*B**2*C**2 + 204*B*C**3 - 239*C**4', (2,5): '-(59*A**4 - 60*A**3*B - 60*A**3*C + 2*A**2*B**2 - 164*A**2*B*C + 2*A**2*C**2 - 60*A*B**3 - 164*A*B**2*C - 164*A*B*C**2 - 60*A*C**3 + 59*B**4 - 60*B**3*C + 2*B**2*C**2 - 60*B*C**3 + 59*C**4)/3', (3,0): '-147*(179*A**3 - 515*A**2*B - 515*A**2*C - 515*A*B**2 + 3858*A*B*C - 515*A*C**2 + 179*B**3 - 515*B**2*C - 515*B*C**2 + 179*C**3)/8', (3,1): '-7*(276*A**3 - 1522*A**2*B - 1522*A**2*C - 1522*A*B**2 + 12489*A*B*C - 1522*A*C**2 + 276*B**3 - 1522*B**2*C - 1522*B*C**2 + 276*C**3)', (3,2): '(1309*A**3 + 13419*A**2*B + 13419*A**2*C + 13419*A*B**2 - 134580*A*B*C + 13419*A*C**2 + 1309*B**3 + 13419*B**2*C + 13419*B*C**2 + 1309*C**3)/3', (3,3): '4*(422*A**3 + 610*A**2*B + 610*A**2*C + 610*A*B**2 - 9081*A*B*C + 610*A*C**2 + 422*B**3 + 610*B**2*C + 610*B*C**2 + 422*C**3)/3', (3,4): '48*(3*A**3 + A**2*B + A**2*C + A*B**2 - 36*A*B*C + A*C**2 + 3*B**3 + B**2*C + B*C**2 + 3*C**3)', (3,5): '4*(9*A**3 - A**2*B - A**2*C - A*B**2 - 78*A*B*C - A*C**2 + 9*B**3 - B**2*C - B*C**2 + 9*C**3)/3', (4,0): '9*(1345*A**2 + 714*A*B + 714*A*C + 1345*B**2 + 714*B*C + 1345*C**2)', (4,1): '(46945*A**2 + 25706*A*B + 25706*A*C + 46945*B**2 + 25706*B*C + 46945*C**2)/4', (4,2): '(25853*A**2 + 16342*A*B + 16342*A*C + 25853*B**2 + 16342*B*C + 25853*C**2)/6', (4,3): '3*(477*A**2 + 422*A*B + 422*A*C + 477*B**2 + 422*B*C + 477*C**2)/2', (4,4): '47*A**2 + 82*A*B + 82*A*C + 47*B**2 + 82*B*C + 47*C**2', (4,5): '(A**2 + 14*A*B + 14*A*C + B**2 + 14*B*C + C**2)/3', (5,0): '-12610*(A + B + C)', (5,1): '-73967*(A + B + C)/6', (5,2): '-9569*(A + B + C)/2', (5,3): '-2764*(A + B + C)/3', (5,4): '-88*(A + B + C)', (5,5): '-10*(A + B + C)/3', (6,0): '4320', (6,1): '4104', (6,2): '1548', (6,3): '290', (6,4): '27', (6,5): '1', } R_ORD, S_DEG = 6, 5 # --------------------------------------------------------------------------- # Guard rails. Calibration (measured at 12 points x 2 precisions covering # the default point, the second rate family, the cubic case, large w, and # near-threshold r=1.0220): # * seed tail bound vs tol 10^-(dps+TAIL_GUARD): margin 10^6.4..10^7.5 # generic, 10^1.0 at the enforced domain edge r=1.022 -- the seed N0 # never escalates in-domain; escalation is benign (adds exact terms). # * oracle two-depth agreement: measured 10^-(dps+25.9..27.8) everywhere # vs tol 10^-(dps+ORA_GUARD=8): >= 10^17.9 headroom, never fires healthy. # * live agreement: measured dps+16.8 (near-threshold) .. dps+21.3 vs # raising bar max(15, dps-AGREE_MARGIN=8): >= 10^24.8 headroom. # --------------------------------------------------------------------------- TAIL_GUARD = 10 # series: certified bound must beat 10^-(dps+TAIL_GUARD) SER_CAP_MULT = 64 # series: escalation cap N <= SER_CAP_MULT * N0 TMAX_GUARD = 12 # oracle: proven Laplace-tail bound < 10^-(dps+TMAX_GUARD) ORA_GUARD = 8 # oracle: two-depth agreement < 10^-(dps+ORA_GUARD) ORA_STEP = 15 # oracle: working-dps ladder increment per rung ORA_MAXITER = 6 # oracle: ladder rungs before the fail-closed raise AGREE_MARGIN = 8 # point mode: RAISE if agreement < max(15, dps-AGREE_MARGIN) def digits(a, b): """-log10 |a-b|/|b| : measured agreement in decimal digits. NB: call under a workdps context at least as precise as the operands (mpmath arithmetic rounds to the ambient precision).""" if a == b: return float('inf') return float(-mp.log10(abs((a - b) / b))) def frac_mp(fr): return mp.mpf(fr.numerator) / mp.mpf(fr.denominator) # =========================================================================== # Exact ingredients: moments from the definition, recurrence transport # =========================================================================== def moments_direct(A, B, C, N): """p_{2m}, m=0..N, exact Fractions, straight from the moment formula (the definition-side algebra; used only for seeds + the --check comparisons).""" A, B, C = F(A), F(B), F(C) fac = [1] for i in range(1, N + 1): fac.append(fac[-1] * i) out = [] for m in range(N + 1): s = F(0) fm = fac[m] for i in range(m + 1): Ai = A ** i for j in range(m - i + 1): k = m - i - j tri = fm // (fac[i] * fac[j] * fac[k]) s += F(tri * tri) * Ai * B ** j * C ** k out.append(F(comb(2 * m, m), 4 ** m) * s) return out def l5_at(A, B, C): """Evaluate the exact symbolic coefficients at rational (A,B,C). The strings involve only +,-,*,**,/ on integers and A,B,C; with A,B,C bound to Fractions every operation stays exact in Q.""" ns = {'A': F(A), 'B': F(B), 'C': F(C), '__builtins__': {}} return {jk: F(eval(expr, ns)) for jk, expr in L5_COEFFS.items()} def _l5_cpoly(A, B, C): """Coefficient polynomials c_j(n) as Horner lists, exact in Q.""" c = l5_at(A, B, C) return [[c[(j, k)] for k in range(S_DEG + 1)] for j in range(R_ORD + 1)] def _cj(cpoly, j, n): v = F(0) for k in range(S_DEG, -1, -1): # Horner in n v = v * n + cpoly[j][k] return v def _extend_transport(a, cpoly, M): """EXACT continuation of the certified recurrence: extends the transported prefix a[0..len(a)-1] to a[0..M] with the SAME recurrence in the same exact Fraction arithmetic. Escalation never recomputes differently -- a fresh transport() to M yields identical Fractions. c_6(n) is a constant positive integer polynomial, so every step is defined for all n >= 0 (checked anyway, fail-closed).""" if len(a) < R_ORD + 1: raise RuntimeError("lgf _extend_transport: need >= %d seed moments, " "got %d" % (R_ORD + 1, len(a))) a = list(a) # start one step back: re-derive the last prefix entry from the # recurrence as well (for a fresh 7-seed transport this is the original # n=0 overwrite of a[6], which keeps the recurrence-vs-seed overlap # inside Check A's checked set; for a continuation it recomputes an # identical Fraction at negligible cost) for n in range(len(a) - R_ORD - 1, M - R_ORD + 1): c6 = _cj(cpoly, R_ORD, n) if not c6 > 0: raise RuntimeError("lgf _extend_transport: leading coefficient " "c_6(%d) = %s not positive -- recurrence data " "corrupted" % (n, c6)) v = -sum(_cj(cpoly, j, n) * a[n + j] for j in range(R_ORD)) / c6 if n + R_ORD < len(a): a[n + R_ORD] = v else: a.append(v) return a[:M + 1] def transport(A, B, C, M, seeds=None): """p_{2m}, m=0..M, by exact forward transport of the certified recurrence from the 7 directly-computed seeds.""" a = list(seeds if seeds is not None else moments_direct(A, B, C, R_ORD)) if len(a) != R_ORD + 1: raise RuntimeError("lgf transport: need exactly %d seeds, got %d" % (R_ORD + 1, len(a))) return _extend_transport(a, _l5_cpoly(A, B, C), M) # =========================================================================== # The evaluator: f(point, dps), plus the live independent oracle # =========================================================================== def nterms(A, B, C, w, dps): """STARTING SEED N0 for the truncation order (DEMOTED -- this closed-form guess is a speed knob, never the answer; acceptance lives in W_series's certified-bound loop). Tail of sum p_{2m} xi^m is geometric with ratio (rho/w)^2, rho = sqrt A + sqrt B + sqrt C (nearest singularity of the PF operator = physical threshold). Also enforces the domain (w/rho >= 1.02) and the cost guard -- both fail-closed raises, never silent truncations.""" with mp.workdps(30): rho = mp.sqrt(A) + mp.sqrt(B) + mp.sqrt(C) ratio = mp.mpf(F(w).numerator) / F(w).denominator / rho if ratio < mp.mpf('1.02'): raise ValueError( "w/(alpha+beta+gamma) = %s < 1.02: point too close to (or on/" "below) the threshold w* = sqrt A+sqrt B+sqrt C. The interior " "series route does not cover it; see the limitations block." % mp.nstr(ratio, 8)) n = int((dps + 12) * mp.ln(10) / (2 * mp.ln(ratio))) + 20 if n > 20000: raise ValueError("would need %d series terms; refusing (cost guard). " "Move w away from threshold or lower --dps." % n) return n def _series_cert_bound(A, B, C, w, N, dps): """PROVEN relative error bound for the N-truncated series (a-priori). p_{2m} = E[(alpha cos k1 + beta cos k2 + gamma cos k3)^{2m}] <= rho^{2m} with rho = sqrt A + sqrt B + sqrt C, so the tail obeys sum_{m>N} p_{2m}/w^{2m} <= q^{N+1}/(1-q), q = (rho/w)^2 < 1, and since W = (1/w)(p_0 + ...) >= 1/w the same number bounds the RELATIVE truncation error. Add the positive-monotone-sum rounding term (N+3)*10^-(dps+19) (workdps(dps+20) half-ulps, l1). Evaluated at workdps(40) with an upward slack factor 1+1e-30 that dominates the bound's own evaluation rounding (~1e-38 rel); TAIL_GUARD absorbs it.""" with mp.workdps(40): wq = F(w) q = ((mp.sqrt(A) + mp.sqrt(B) + mp.sqrt(C)) * wq.denominator / wq.numerator) ** 2 q *= 1 + mp.mpf('1e-30') if not q < 1: raise RuntimeError("lgf _series_cert_bound: q=(rho/w)^2 = %s >= 1 " "at (A,B,C;w)=(%s,%s,%s;%s) -- outside the " "series domain" % (mp.nstr(q, 8), A, B, C, w)) return q ** (N + 1) / (1 - q) + (N + 3) * mp.mpf(10) ** (-(dps + 19)) def W_series(A, B, C, w, dps, mom=None): """W_S by the certified route. All series terms are exact rationals; mpmath enters only in the final positive monotone sum. The truncation N sits in a refine-until-bound loop -- accepted only when the PROVEN bound of _series_cert_bound beats 10^-(dps+TAIL_GUARD); escalation x1.5 by EXACT continuation of the same recurrence (_extend_transport); RuntimeError at the cap, never an uncertified value. Returns (val, N, mom, certified_bound).""" N0 = nterms(A, B, C, w, dps) # seed only (speed knob) ncap = SER_CAP_MULT * N0 tol = mp.mpf(10) ** (-(dps + TAIL_GUARD)) N = N0 while True: bound = _series_cert_bound(A, B, C, w, N, dps) if bound < tol: break if N >= ncap: raise RuntimeError( "lgf W_series: certified tail bound %s >= tol %s at " "(A,B,C;w)=(%s,%s,%s;%s), dps=%d, N=%d (seed N0=%d, cap %d) " "-- refusing to print an uncertified value" % (mp.nstr(bound, 3), mp.nstr(tol, 3), A, B, C, w, dps, N, N0, ncap)) N = min(ncap, max(int(1.5 * N), N + 8)) if mom is None: mom = transport(A, B, C, N) elif len(mom) <= N: # EXACT continuation mom = _extend_transport(mom, _l5_cpoly(A, B, C), N) w = F(w) with mp.workdps(dps + 15 + 5): xi = 1 / (frac_mp(w) ** 2) s = mp.mpf(0) x = mp.mpf(1) for m in range(N + 1): s += frac_mp(mom[m]) * x x *= xi val = s / frac_mp(w) return val, N, mom, bound def W_oracle(A, B, C, w, dps, full=False): """Independent live oracle: Bessel-Laplace representation W_S = int_0^inf e^{-w t} I_0(a t) I_0(b t) I_0(c t) dt, a=sqrt A etc., by tanh-sinh quadrature (reimplemented from the formula; shares no arithmetic with the moment/recurrence route). Fail-closed double-refinement quadrature: * [0,tmax] truncation: the seed tmax formula is a starting guess; the PROVEN tail bound e^(-gap*tmax)/gap * w (via I_0(x) <= e^x and W >= 1/w) is enforced < 10^-(dps+TMAX_GUARD) at runtime, tmax x1.5 up to 64x seed, then RuntimeError. * quadrature depth: accepted only when evaluations at two successive working precisions (wd, wd+ORA_STEP -- genuinely different node sets) agree to < 10^-(dps+ORA_GUARD) relative. The COARSER member of the first agreeing pair is returned, so a healthy seed evaluation is bit-identical to the seed route; the ladder raises after ORA_MAXITER rungs. (Measured healthy agreement: 10^-(dps+25.9..27.8); tol leaves >= 10^17.9 headroom.) Returns the value, or (value, cert dict) when full=True.""" Aq, Bq, Cq, wq = F(A), F(B), F(C), F(w) with mp.workdps(dps + 25): a, b, c = mp.sqrt(Aq), mp.sqrt(Bq), mp.sqrt(Cq) ww = frac_mp(wq) gap = ww - (a + b + c) if not gap > 0: raise ValueError("lgf W_oracle: w = %s <= threshold rho = %s -- " "Laplace integral diverges" % (wq, mp.nstr(a + b + c, 8))) tmax0 = (dps + 25) * mp.ln(10) / gap + 10 # seed only ttol = mp.mpf(10) ** (-(dps + TMAX_GUARD)) tmax = tmax0 while True: tail_rel = mp.e ** (-gap * tmax) / gap * ww # PROVEN, I0<=e^x if tail_rel < ttol: break if tmax > 64 * tmax0: raise RuntimeError( "lgf W_oracle: certified Laplace-tail bound %s >= tol %s " "at (A,B,C;w)=(%s,%s,%s;%s), dps=%d, tmax=%s (seed %s, " "cap 64x)" % (mp.nstr(tail_rel, 3), mp.nstr(ttol, 3), Aq, Bq, Cq, wq, dps, mp.nstr(tmax, 8), mp.nstr(tmax0, 8))) tmax *= mp.mpf('1.5') pts = [p for p in [0, mp.mpf('1e-3'), 1, 10, 100] if p < tmax] + [tmax] f = lambda t: (mp.exp(-ww * t) * mp.besseli(0, a * t) * mp.besseli(0, b * t) * mp.besseli(0, c * t)) vprev, eprev = mp.quad(f, pts, method='tanh-sinh', error=True, maxdegree=10) wd_prev = dps + 25 tol = mp.mpf(10) ** (-(dps + ORA_GUARD)) for rung in range(1, ORA_MAXITER + 1): wd = dps + 25 + ORA_STEP * rung with mp.workdps(wd): a2, b2, c2 = mp.sqrt(Aq), mp.sqrt(Bq), mp.sqrt(Cq) ww2 = frac_mp(wq) f2 = lambda t: (mp.exp(-ww2 * t) * mp.besseli(0, a2 * t) * mp.besseli(0, b2 * t) * mp.besseli(0, c2 * t)) pts2 = ([p for p in [0, mp.mpf('1e-3'), 1, 10, 100] if p < tmax] + [tmax]) # same certified cutoff vcur, ecur = mp.quad(f2, pts2, method='tanh-sinh', error=True, maxdegree=10) agree = abs(vprev - vcur) / abs(vcur) if agree < tol: cert = {'agree': agree, 'tol': tol, 'wd': wd_prev, 'wd_next': wd, 'tail_rel': tail_rel, 'tmax': tmax, 'quad_err_est': ecur} return (vprev, cert) if full else vprev vprev, eprev, wd_prev = vcur, ecur, wd raise RuntimeError( "lgf W_oracle: two-depth agreement %s >= tol %s after %d rungs " "(working dps %d..%d) at (A,B,C;w)=(%s,%s,%s;%s), dps=%d -- " "quadrature did not certify; refusing to return a value" % (mp.nstr(agree, 3), mp.nstr(tol, 3), ORA_MAXITER, dps + 25, wd, Aq, Bq, Cq, wq, dps)) def evaluate(A, B, C, w, dps=50, oracle=True): """f(kinematic point, dps): returns dict with the series value, the live oracle value and the measured agreement. A,B,C,w positive rationals, w > sqrt A + sqrt B + sqrt C (see nterms for the enforced margin).""" A, B, C, w = F(A), F(B), F(C), F(w) if min(A, B, C) <= 0 or w <= 0: raise ValueError("need A,B,C,w positive rationals") t0 = time.time() val, N, _, cert_bound = W_series(A, B, C, w, dps) res = {'W': val, 'N': N, 'cert_bound': cert_bound, 't_series_s': time.time() - t0, 'oracle': None} if oracle: t0 = time.time() vo, ocert = W_oracle(A, B, C, w, dps, full=True) with mp.workdps(dps + 25): d = digits(vo, val) bar_live = max(15.0, dps - AGREE_MARGIN) if not d >= bar_live: # RAISING check; catches nan raise RuntimeError( "lgf live-agreement gate: certified series vs certified " "Bessel-Laplace oracle agree to only %.2f d < required %.1f d " "at (A,B,C;w)=(%s,%s,%s;%s), dps=%d (measured healthy: " ">= dps+16.8; embedded data or arithmetic corrupted) -- " "series bound %s, oracle agreement %s" % (d, bar_live, A, B, C, w, dps, mp.nstr(cert_bound, 3), mp.nstr(ocert['agree'], 3))) res['oracle'] = {'W_oracle': vo, 'digits': d, 'dps': dps, 'cert': ocert, 't_oracle_s': time.time() - t0} return res # =========================================================================== # Gate mode # =========================================================================== DEF_POINT = (F(2), F(3), F(5), F(8)) # the paper's generic point, w=8 HELD_OUT = [(F(2), F(3), F(7), F(9)), # second rate family (F(1), F(4), F(9), F(8))] # integer rates (1,2,3) def check(dps, point=DEF_POINT): T0 = time.time() A, B, C, w = point dps2 = 2 * dps bar1 = max(30.0, dps - 12) print("Gate mode: point (A,B,C;w)=(%s,%s,%s;%s), dps=%d, doubling at %d" % (A, B, C, w, dps, dps2)) print(" bar: exact-transport gates must be EXACT (Fraction ==); live " "agreement >= %.0f d at dps,\n and must grow by >= %.0f d " "under dps doubling" % (bar1, 0.4 * dps)) ok = True print("\n== Check A: exact transport (recurrence vs definition, further points) ==") for (Ah, Bh, Ch) in [(A, B, C)] + [(p[0], p[1], p[2]) for p in HELD_OUT]: t0 = time.time() tr = transport(Ah, Bh, Ch, 60) di = moments_direct(Ah, Bh, Ch, 60) eq = sum(1 for m in range(61) if tr[m] == di[m]) good = (eq == 61) ok &= good print(" (A,B,C)=(%s,%s,%s): %d/61 transported moments == direct " "(m=7..60 held out) [%.1f s] %s" % (Ah, Bh, Ch, eq, time.time() - t0, "PASS" if good else "FAIL")) t0 = time.time() deep = 150 trd = transport(A, B, C, deep)[deep] did = moments_direct(A, B, C, deep)[deep] good = (trd == did) ok &= good print(" deep spot m=%d at (%s,%s,%s): exact match = %s [%.1f s] %s" % (deep, A, B, C, trd == did, time.time() - t0, "PASS" if good else "FAIL")) print("\n== Check B: two-precision live-quadrature check at the point ==") agr = {} for dd in (dps, dps2): t0 = time.time() r = evaluate(A, B, C, w, dps=dd) agr[dd] = r['oracle']['digits'] print(" dps=%4d: N=%4d terms, W = %s" % (dd, r['N'], mp.nstr(r['W'], min(dd, 40)))) print(" live oracle agreement = %7.2f d " "[series %.1f s, oracle %.1f s]" % (agr[dd], r['t_series_s'], r['oracle']['t_oracle_s'])) g1 = agr[dps] >= bar1 g2 = agr[dps2] >= agr[dps] + 0.4 * dps ok &= g1 and g2 print(" agreement >= %.0f d at dps=%d: %s; growth %.1f d " "(>= %.0f required): %s" % (bar1, dps, "PASS" if g1 else "FAIL", agr[dps2] - agr[dps], 0.4 * dps, "PASS" if g2 else "FAIL")) print("\n== Check C: second-family kinematic points (live, dps=%d) ==" % dps) for (Ah, Bh, Ch, wh) in HELD_OUT: r = evaluate(Ah, Bh, Ch, wh, dps=dps) d = r['oracle']['digits'] good = d >= bar1 ok &= good print(" (A,B,C;w)=(%s,%s,%s;%s): W = %s" % (Ah, Bh, Ch, wh, mp.nstr(r['W'], min(dps, 40)))) print(" agreement = %7.2f d (bar %.0f d) %s" % (d, bar1, "PASS" if good else "FAIL")) print("\n== Check D: crank, series at dps=%d vs dps+40=%d " "(independent certified depths) ==" % (dps, dps + 40)) for (Ah, Bh, Ch, wh) in [(A, B, C, w)] + HELD_OUT: t0 = time.time() vD, ND, _, bD = W_series(Ah, Bh, Ch, wh, dps) vD40, ND40, _, bD40 = W_series(Ah, Bh, Ch, wh, dps + 40) with mp.workdps(dps + 60): dcr = digits(vD, vD40) good = dcr >= dps ok &= good print(" (A,B,C;w)=(%s,%s,%s;%s): N %d -> %d, D/D+40 agreement " "%.2f d (>= %d required) [%.1f s] %s" % (Ah, Bh, Ch, wh, ND, ND40, dcr, dps, time.time() - t0, "PASS" if good else "FAIL")) print("\nOVERALL: %s [total wall %.1f s]" % ("PASS" if ok else "FAIL", time.time() - T0)) return ok def main(): ap = argparse.ArgumentParser( description="Anisotropic sc-lattice Green function W_S via the " "certified order-5 PF/moment-recurrence route, checked " "live against an independent Bessel-Laplace quadrature " "(see module docstring for domain limits).") ap.add_argument('--dps', type=int, default=50, help='working precision (default 50)') ap.add_argument('--point', nargs=4, metavar=('A', 'B', 'C', 'W'), help='squared rates A B C (positive rationals) and ' 'spectral variable W (default: 2 3 5 8)') ap.add_argument('--check', action='store_true', help='run the exact-transport, two-precision, and ' 'second-family live-quadrature checks') ap.add_argument('--no-oracle', action='store_true', help='point mode: skip the live oracle') args = ap.parse_args() pt = tuple(F(x) for x in args.point) if args.point else DEF_POINT if args.check: raise SystemExit(0 if check(args.dps, pt) else 1) A, B, C, w = pt T0 = time.time() print("Point (A,B,C;w)=(%s,%s,%s;%s), rates (sqrtA,sqrtB,sqrtC), dps=%d" % (A, B, C, w, args.dps)) r = evaluate(A, B, C, w, dps=args.dps, oracle=not args.no_oracle) print(" W_S =", mp.nstr(r['W'], args.dps)) print(" series: %d exact-rational terms via certified transport, %.1f s" % (r['N'], r['t_series_s'])) print(" certified BOUND: series truncation+rounding rel. error <= %s " "(tol 1e-%d)" % (mp.nstr(r['cert_bound'], 3), args.dps + TAIL_GUARD)) if r['oracle']: print(" live Bessel-Laplace oracle (dps=%d): agreement = %.2f d " " [%.1f s]" % (r['oracle']['dps'], r['oracle']['digits'], r['oracle']['t_oracle_s'])) oc = r['oracle']['cert'] print(" certified BOUND: oracle two-depth agreement %s (tol 1e-%d, " "wd %d->%d), Laplace-tail %s; raising gate >= %.0f d: PASS" % (mp.nstr(oc['agree'], 3), args.dps + ORA_GUARD, oc['wd'], oc['wd_next'], mp.nstr(oc['tail_rel'], 3), max(15.0, args.dps - AGREE_MARGIN))) print("Total wall time: %.1f s" % (time.time() - T0)) if __name__ == '__main__': main()