#!/usr/bin/env python3 """lss-evaluate.py - recompute the closed form of the two-loop kite family and verify it against the reference value shipped beside this script. NOTE (2026-09-27): the closed form is restated as Eq. (5.4) of "The loop expansion of galaxy clustering in redshift space" (Schwartz, Ivanov, Mishra-Sharma; Paper I of the series now on the site). It was first obtained in the earlier single manuscript that Papers I and II supersede, where it was called Theorem 1; the labels "Theorem 1" and "the paper" below refer to that statement, which Paper I carries unchanged. The object is the last open master of the two-loop massless kite family, K(1, nu2, nu3, 1, nu5) at d = 3 in the K(nu1..nu5) convention of arXiv:1708.08130 (open since Kotikov 1996 / Grozin 2012). The paper's closure theorem (Theorem 1) writes it, per block j in {1,2,3} of the Gegenbauer/DIS representation (arXiv:2303.09203 eqs. 25-27), as Ihat_j = Sa_j*Q0_j(inf) + 2*Sb_j*Q1_j(inf) - W_j K = kpref * pref * M0 * sum_j t00_j * Ihat_j where the head pieces Sa, Sb, Q0(inf), Q1(inf) are single hypergeometric pFq(1) sums, W_2 = 0 exactly, and W_1, W_3 are the paper's two new transcendentals (single-sum wedge definitions). The raw wedge sum is two-scale and defeats sequence acceleration, so this script evaluates the single-scale sum V_j = sum_k Psihat_j(k) [ (1-2k)(Q0_j(inf)-Q0_j(k-1)) + 2(Q1_j(inf)-Q1_j(k-1)) ] and uses the exact Abel-summation identity W_j = Sa_j*Q0_j(inf) + 2*Sb_j*Q1_j(inf) - V_j (so Ihat_j = V_j). Every series term is a Gaussian rational multiple of the first term, generated by exact term-ratio recurrences; Levin-u acceleration supplies the limits. Modes ----- default Evaluate every piece of Theorem 1 at the certified point B1 (nu = 1, 3/5+21i/100, 9/20-33i/100, 1, 2/5+17i/100), assemble K, and compare against lss-reference-B1.json (beside this script; 307 certified digits, computed by the independent Levin-DIS row organization -- a different summation organization from the one used here). At the default depth the assembly saturates the reference (measured agreement 307.4 digits, i.e. every certified digit); the run PASSES only if it agrees to at least GATE_DIGITS = 250 digits, a bar nothing but a genuine defect can miss at this depth. --check Rediscover Theorem 1 as an integer relation: build the 9-vector { Kn(reference); the six t00-weighted head products; X1 = t00_1*W_1; X3 = t00_3*W_3 } and require PSLQ to return the coefficient vector (1,-1,-2,-1,-2,-1,-2,1,1) -- up to overall sign, which PSLQ does not fix -- on the real leg at TWO working precisions (140 and 190 digits) and on the imaginary leg (140 digits), with two controls that MUST fail: (a) random-target null -- the reference slot replaced by a deterministic pseudo-random number of the same magnitude must yield NO relation at the same height cap; (b) perturbation control -- the reference perturbed by one part in 1e-40 must LOSE the relation. A control that succeeds where it must fail exits 3. The reference slot comes from the vendored certified value, not from this script's own assembly, so the relation ties two independent summation organizations together; --check always runs at the default depth (--depth is refused there, exit 5). --audit HELD. The per-graph audit of the paper's exact two-loop matter table (140 rows, checked against the AFLSZ Planck kernel bank) needs that table and bank to be public; they are not yet released. Until the data release ships, this flag reports the hold and exits 4. --depth N WD Override the default series budget N and working precision WD (digits) in the default mode, e.g. --depth 120 400 for a fast partial run (measured: 106.5 agreed digits in 2-5 s depending on contention). Below the default depth the 250-digit gate fails by construction -- useful for watching the digit ladder, not for certification. Measured ladder on the build server: N=40/WD=150 34.0 d, N=80/250 70.3 d, N=120/400 106.5 d, N=160/550 142.8 d, N=260/900 233.7 d, N=400/1300 (default) 307.4 d. Exit codes ---------- 0 pass (default gate met; --check relation found, stable, controls fail) 2 verification failure (agreement below the gate, or the --check relation not found / not stable across precisions) 3 control anomaly (a must-fail control passed) 4 held (--audit requested before the data release it needs) 5 usage or environment error (bad arguments, missing mpmath, missing or byte-damaged lss-reference-B1.json) Determinism and timing ---------------------- stdout is byte-deterministic: identical on every run of the same mode at the same depth (no timestamps, no wall times; the random-target null uses a fixed seed). Timing goes to stderr. Measured with /usr/bin/time on the build server (2026-09-03, one core of a heavily loaded 96-core machine): default depth wall 56 s lightly contended and 147 s heavily contended, max RSS 34 MB; --check wall 94-143 s; the first stdout line is printed before mpmath loads (well under a second). The work is one single-thread mpmath process with a 34 MB footprint; on an ordinary laptop core expect a few minutes at most. Python 3.8+ with mpmath (built against 1.3.0); nothing else, no network. """ import hashlib import json import os import sys FIRST_LINE = ('lss-evaluate: two-loop kite closure (Theorem 1, Abel-split ' 'W1/W3) at point B1') print(FIRST_LINE, flush=True) try: import mpmath as mp from mpmath import mpf, mpc except ImportError: print('ERROR: this script needs mpmath (pip install mpmath)') sys.exit(5) HERE = os.path.dirname(os.path.abspath(__file__)) REF_FILE = os.path.join(HERE, 'lss-reference-B1.json') REF_SHA256 = '79c5342baa01acbc86640ac197a17d12df9ff1ef2afdda44af2330c3d0e6f61b' # Default budget: series terms N and working digits WD. Measured # 2026-09-03 on the build server (/usr/bin/time): (400, 1300) agrees with the # certified reference to 307.4 digits -- saturating all 307 certified # digits -- in 56-147 s wall depending on machine contention. The gate # sits below the measured agreement by # a margin only a genuine defect could eat (the arithmetic itself is # deterministic; the margin covers mpmath-version differences in the # Levin-u internals, not noise), and far above anything a wrong assembly # could reach. N_DEFAULT, WD_DEFAULT = 400, 1300 GATE_DIGITS = 250 # Expected PSLQ vector for --check (Theorem 1 as an integer relation). RELATION = (1, -1, -2, -1, -2, -1, -2, 1, 1) PSLQ_MAXCOEFF = 10**6 D = mpf(3) # spacetime dimension G = mp.gamma def a0(x, d): return G(d / 2 - x) / G(x) def digits(x, y): if x == y: return 999.0 return float(-mp.log10(abs(x - y) / (abs(x) + abs(y)) * 2)) # ---- factor tables of the three DIS blocks: (Gnk, Rnk, Gk, Rk), lists of # constants c with term Gamma(c+n+k) (Gnk), 1/Gamma(c+n+k) (Rnk), # Gamma(c-k) (Gk), 1/Gamma(c-k) (Rk). def tables(A, d): A1, A2, A3, A4, A5 = A A15, A34, A125 = A1 + A5, A3 + A4, A1 + A2 + A5 A145, A345, A1345 = A1 + A4 + A5, A3 + A4 + A5, A1 + A3 + A4 + A5 A12345 = A1 + A2 + A3 + A4 + A5 h = d / 2 I1 = ([A4, h - A5, d - A125], [d - A15, A34, h], [A15 - h, h - A34], [h - A4, A5, A125 - h]) I2 = ([A145 - h, A1, h - A2], [A15, A1345 - h, h], [h - A15, d - A1345], [d - A145, h - A1, A2]) I3 = ([h - A3, d - A345, 3 * h - A12345], [d - A34, 3 * h - A1345, h], [A34 - h, A1345 - d], [A3, A345 - h, A12345 - d]) return [I1, I2, I3] def t00(tab): Gnk, Rnk, Gk, Rk = tab v = mpc(1) for c in Gnk: v *= G(c) for c in Rnk: v *= mp.rgamma(c) for c in Gk: v *= G(c) for c in Rk: v *= mp.rgamma(c) return v def nratio(tab, n): """Phihat(n+1)/Phihat(n).""" Gnk, Rnk, Gk, Rk = tab r = mpc(1) for c in Gnk: r *= (c + n) for c in Rnk: r /= (c + n) return r def rpsi(tab, k): """Psihat(k+1)/Psihat(k); exact zero when a numerator factor vanishes.""" Gnk, Rnk, Gk, Rk = tab r = mpf(-1) / (k + 1) for c in tab[3]: r *= (c - k - 1) for c in tab[2]: r /= (c - k - 1) return r def levsum(gen, N): """Levin-u limit of the series from gen; (value, nterms, truncated).""" L = mp.mp.levin(method='levin', variant='u') S, s = [], mpc(0) best, beste = None, mp.inf for i in range(N): t, dead = next(gen) if dead: # the series truncates exactly -> s is EXACT return s, i, True s += t if S or s != 0: # start the table at the first nonzero psum S.append(s) if len(S) >= 12 and (i & 3) == 3: v, e = L.update_psum(S) if e < beste and mp.isfinite(v): best, beste = v, e return (best if best is not None else s), N, False def gen_phi(tab, weight): t, m = mpc(1), 0 while True: yield weight(m) * t, (t == 0) t = t * nratio(tab, m) m += 1 def gen_psi(tab, weight): t, k = mpc(1), 0 while True: yield weight(k) * t, (t == 0) t = t * rpsi(tab, k) k += 1 def gen_vee(tab, Q0inf, Q1inf): """Single-scale wedge partner V_j (see header); Abel: W_j = Sa*Q0inf + 2*Sb*Q1inf - V_j exactly, and Ihat_j = V_j.""" psi, phi = mpc(1), mpc(1) q0, q1 = mpc(0), mpc(0) # Q0(k-1), Q1(k-1) k = 0 while True: yield psi * ((1 - 2 * k) * (Q0inf - q0) + 2 * (Q1inf - q1)), (psi == 0) q0 += phi q1 += k * phi phi = phi * nratio(tab, k) psi = psi * rpsi(tab, k) k += 1 def load_reference(): """Vendored certified reference; byte-pinned, exit 5 on any damage.""" try: raw = open(REF_FILE, 'rb').read() except OSError: print('ERROR: lss-reference-B1.json not found beside this script') sys.exit(5) sha = hashlib.sha256(raw).hexdigest() if sha != REF_SHA256: print('ERROR: lss-reference-B1.json sha256 mismatch') print(' expected ' + REF_SHA256) print(' got ' + sha) sys.exit(5) R = json.loads(raw) return R def b1_nus(): """The certified point B1, exact rationals at current precision.""" one = mpf(1) return [mpc(one), mpc(mpf(3) / 5, mpf(21) / 100), mpc(mpf(9) / 20, mpf(-33) / 100), mpc(one), mpc(mpf(2) / 5, mpf(17) / 100)] def compute_pieces(N, wd): """All Theorem-1 pieces at B1: returns dict P and prefactors.""" import time mp.mp.dps = wd nus = b1_nus() d = D A = [d / 2 - nus[3], d / 2 - nus[1], d / 2 - nus[0], d / 2 - nus[2], d / 2 - nus[4]] tabs = tables([mpc(x) for x in A], d) t0s = [t00(tab) for tab in tabs] P = {'t00': t0s} t_start = time.time() for j in (0, 1, 2): s = str(j + 1) P['Q0_' + s], _, _ = levsum(gen_phi(tabs[j], lambda m: 1), N) P['Q1_' + s], _, _ = levsum(gen_phi(tabs[j], lambda m: m), N) P['Sa_' + s], _, _ = levsum(gen_psi(tabs[j], lambda k: 1 - 2 * k), N) P['Sb_' + s], _, _ = levsum(gen_psi(tabs[j], lambda k: 1), N) P['V_' + s], _, tr = levsum( gen_vee(tabs[j], P['Q0_' + s], P['Q1_' + s]), N) # exact Abel identity P['W_' + s] = (P['Sa_' + s] * P['Q0_' + s] + 2 * P['Sb_' + s] * P['Q1_' + s] - P['V_' + s]) print(' block %s summed (series budget %d terms each)' % (s, N), flush=True) print(' block %s wall %.1f s' % (s, time.time() - t_start), file=sys.stderr, flush=True) Kn = sum(t0s[j] * P['V_%d' % (j + 1)] for j in range(3)) M0 = (d / 2 - 1) * G(d - 2) pref = (mp.pi**d * G(d / 2 - 1) / G(d - 2) * a0(A[0], d) * a0(A[3], d) / a0(sum(A) - d, d)) kpref = mpf(4)**(-d) * mp.pi**(-2 * d) for x in nus: kpref *= a0(x, d) kpref *= a0(3 * d / 2 - sum(nus), d) P['Kn'] = Kn P['K'] = kpref * pref * M0 * Kn P['norm'] = kpref * pref * M0 return P def get_ref_value(R): return mpc(mp.mpf(R['value']['re']), mp.mpf(R['value']['im'])) def run_default(N, wd): R = load_reference() print('reference: lss-reference-B1.json sha256 MATCH ' '(%d certified digits, engine: %s)' % (R['value']['certified_digits'], R['value']['engine']), flush=True) print('computing Theorem 1 pieces at B1, series budget N=%d, ' 'working precision %d digits' % (N, wd), flush=True) P = compute_pieces(N, wd) ref = get_ref_value(R) d_W2 = abs(P['W_2']) agree = digits(P['K'], ref) print('', flush=True) print(' W_2 (exactly zero in Theorem 1) |computed| = %s' % mp.nstr(d_W2, 3), flush=True) print(' K assembled = %s' % mp.nstr(P['K'], 20), flush=True) print(' K reference = %s' % mp.nstr(ref, 20), flush=True) print(' agreement: %.1f digits (gate: >= %d)' % (min(agree, 999.0), GATE_DIGITS), flush=True) ok = agree >= GATE_DIGITS and d_W2 < mpf(10)**(-GATE_DIGITS) print('', flush=True) if ok: print('OVERALL: PASS -- the closure reproduces the independent ' '307-digit reference at this depth.', flush=True) return 0 print('OVERALL: FAIL -- agreement below the gate.', flush=True) return 2 def pslq_leg(vec_real, dps, tag): """One PSLQ leg at the given precision; returns the relation or None.""" with mp.workdps(dps): v = [+x for x in vec_real] # round into this precision scale = max(abs(x) for x in v) v = [x / scale for x in v] rel = mp.pslq(v, tol=mpf(10)**(12 - dps), maxcoeff=PSLQ_MAXCOEFF, maxsteps=200000) print(' %s: %s' % (tag, rel), flush=True) return rel def norm_sign(rel): if rel is None: return None t = tuple(rel) return t if t[0] >= 0 else tuple(-c for c in t) def run_check(N, wd): R = load_reference() print('reference: lss-reference-B1.json sha256 MATCH', flush=True) print('building the 9-vector {Kn_ref; t00-weighted heads; X1; X3} at ' 'B1, N=%d, %d working digits' % (N, wd), flush=True) P = compute_pieces(N, wd) ref = get_ref_value(R) Kn_ref = ref / P['norm'] basket = [Kn_ref] for j in (1, 2, 3): tj = P['t00'][j - 1] basket.append(tj * P['Sa_%d' % j] * P['Q0_%d' % j]) basket.append(tj * P['Sb_%d' % j] * P['Q1_%d' % j]) basket.append(P['t00'][0] * P['W_1']) basket.append(P['t00'][2] * P['W_3']) # Two-precision real leg + imaginary leg. The assembled value at the # default depth carries 307.4 agreed digits (measured); the legs sit # far below the constituent accuracy. legs = [('Re leg, 140 digits', [x.real for x in basket], 140), ('Re leg, 190 digits', [x.real for x in basket], 190), ('Im leg, 140 digits', [x.imag for x in basket], 140)] found = [] print('', flush=True) print('positive legs (must each return the Theorem-1 vector ' '%s, up to overall sign):' % (RELATION,), flush=True) for tag, v, dps in legs: found.append(norm_sign(pslq_leg(v, dps, tag))) ok_pos = all(f == RELATION for f in found) stable = len(set(found)) == 1 and found[0] is not None if not ok_pos or not stable: print('', flush=True) print('OVERALL: FAIL -- the integer relation was not recovered ' 'or is not stable across precisions.', flush=True) return 2 print('', flush=True) print('controls (each MUST fail to find the relation):', flush=True) # (a) random-target null: fixed seed -> deterministic output. import random rng = random.Random(20260903) with mp.workdps(wd): rnd = mpf(rng.random()) * abs(Kn_ref.real) null_v = [rnd] + [x.real for x in basket[1:]] rel_null = pslq_leg(null_v, 140, 'random-target null, 140 digits') # (b) perturbation control: reference off by 1e-40. with mp.workdps(wd): pert = Kn_ref.real * (1 + mpf(10)**(-40)) pert_v = [pert] + [x.real for x in basket[1:]] rel_pert = pslq_leg(pert_v, 140, 'perturbed reference (1e-40), 140 digits') bad = [] if norm_sign(rel_null) == RELATION: bad.append('random-target null recovered the relation') if norm_sign(rel_pert) == RELATION: bad.append('perturbed reference kept the relation') print('', flush=True) if bad: for b in bad: print('CONTROL ANOMALY: ' + b, flush=True) print('OVERALL: FAIL -- a must-fail control passed.', flush=True) return 3 print(' random-target null: no relation at height <= %d ' '(failed as required)' % PSLQ_MAXCOEFF, flush=True) print(' perturbation control: relation LOST ' '(failed as required)', flush=True) print('', flush=True) print('OVERALL: PASS -- PSLQ recovers (1,-1,-2,-1,-2,-1,-2,1,1) on ' 'both real legs and the imaginary leg, and both controls fail ' 'as required.', flush=True) return 0 def run_audit(): print('', flush=True) print('--audit is HELD: the per-graph audit of the paper\'s exact ' 'two-loop matter table (140 rows, checked against the AFLSZ ' 'Planck kernel bank) needs that table and kernel bank to be ' 'public, and they are not yet released. The evaluator and its ' '--check ship now; the table audit ships with the data ' 'release.', flush=True) return 4 def main(argv): N, wd = N_DEFAULT, WD_DEFAULT args = list(argv) depth_given = '--depth' in args if depth_given: i = args.index('--depth') try: N, wd = int(args[i + 1]), int(args[i + 2]) except (IndexError, ValueError): print('usage: lss-evaluate.py [--check | --audit] ' '[--depth N WD]') return 5 del args[i:i + 3] if N < 20 or wd < 60: print('ERROR: --depth needs N >= 20 and WD >= 60') return 5 if args == []: return run_default(N, wd) if args == ['--check']: if depth_given: print('ERROR: --check always runs at the default depth; ' 'drop --depth') return 5 return run_check(N, wd) if args == ['--audit']: return run_audit() print('usage: lss-evaluate.py [--check | --audit] [--depth N WD]') return 5 if __name__ == '__main__': sys.exit(main(sys.argv[1:]))