"""Broadhurst–Mellit vacuum determinants -- numberdb.org/T648 Run with SageMath: sage -pip install numberdb mpmath # once sage -python generate.py # recompute and verify all three rows sage -python generate.py --compute # numerical evidence, no website write sage -python generate.py --fill # fill draft using NUMBERDB_API_KEY The only selected factor counts are 5,6,8. Unit-scale, order-zero Bessel moments have rows a=0,2. Columns are t and t^3, except at n=8 where the second column is t^3-2*t^5. Read argument-selection.md for the coverage rationale. No modular identity is used to generate these values. Tanh-sinh integration uses Arb Bessel midpoints, not certified integration. Runs at 150 and 200 decimal digits establish agreement; 100 digits are stored. Bessel precision doubles, with a hard limit, until its relative ball radius is below 10^(-working-20). This deals with cancellation in Arb's K0 evaluation. The integration intervals are [0,1,4,16,64,T], T=3*(working+60). At zero every integrand tends to zero. For t>=T, I0(t)<=exp(t) and K0(t)<=sqrt(pi/(2*t))*exp(-t), the latter from cosh(u)>=1+u^2/2. Writing b=n-a, lambda=b-a>0, a moment's omitted tail is at most (pi/2)^(b/2)*T^(r-b/2)*exp(-lambda*T)/(lambda-r/T), provided lambda>r/T. The estimate follows from (t/T)^r<=exp(r*(t-T)/T). Polynomial-column tail bounds follow by the triangle inequality. This controls truncation, but does not turn the quadrature error estimates into certificates. The first two nonzero controls run before the table moments. Additional full-period controls use Zhou's proved two-by-two determinant formulas. All sources controlling selection are in this file; checks.py is independent. """ import numberdb.sage as numberdb from sage.rings.complex_arb import ComplexBallField from sage.env import SAGE_VERSION from functools import lru_cache import mpmath as mp import json import sys SELECTION = (5, 6, 8) WORKING = (150, 200) @lru_cache(maxsize=4) def compute_all(working): with mp.workdps(working + 35): cutoff = 3 * (working + 60) bits = numberdb.bits(working + 40) cache = {} max_bits = [bits] def bessel(t): key = t._mpf_ if key in cache: return cache[key] precision = bits while precision <= 8 * bits: z = ComplexBallField(precision)(mp.nstr(t, working + 32)) ib, kb = z.bessel_I(0).real(), z.bessel_K(0).real() if ib.is_finite() and kb.is_finite(): i, k = mp.mpf(str(ib.mid())), mp.mpf(str(kb.mid())) if (i > 0 and k > 0 and mp.mpf(str(ib.rad())) < i * mp.power(10, -working-20) and mp.mpf(str(kb.rad())) < k * mp.power(10, -working-20)): max_bits[0] = max(max_bits[0], precision) cache[key] = (i, k) return i, k precision *= 2 raise ArithmeticError('Bessel precision exhausted') moments, errors, tails = {}, {}, {} def moment(n, a, r): key = f'{n},{a},{r}' if key in moments: return mp.mpf(moments[key]) def f(t): if not t: return mp.mpf(0) i, k = bessel(t) return i**a * k**(n-a) * t**r with mp.workdps(working): val, err = mp.quad(f, [0,1,4,16,64,cutoff], error=True) assert mp.isfinite(val) and val > 0 nstr = lambda x: mp.nstr(x, working) moments[key], errors[key] = nstr(val), nstr(err) b=n-a; lam=n-2*a bound=(mp.pi/2)**(mp.mpf(b)/2)*mp.mpf(cutoff)**(mp.mpf(r)-mp.mpf(b)/2)*mp.exp(-lam*cutoff)/(mp.mpf(lam)-mp.mpf(r)/cutoff) tails[key] = nstr(bound) return val controls={} for n, expected in ((1,mp.mpf(1)),(2,mp.mpf(1)/2)): discrepancy=abs(moment(n,0,1)/expected-1) assert discrepancy < mp.power(10,-working+10) controls[str(n)] = mp.nstr(discrepancy,8) values, cancellation = {}, {} for n in SELECTION: x,y,u,v = moment(n,0,1),moment(n,0,3),moment(n,2,1),moment(n,2,3) if n==8: y -= 2*moment(n,0,5) v -= 2*moment(n,2,5) value = x*v-y*u assert mp.isfinite(value) and value > 0 values[str(n)] = mp.nstr(value,working) cancellation[str(n)] = mp.nstr(mp.log10((abs(x*v)+abs(y*u))/abs(value)),12) period_controls={} for n, expected in ((5,2*mp.pi**3/mp.sqrt(3**3*5**5)),(6,mp.pi**4/576)): value=moment(n,1,1)*moment(n,2,3)-moment(n,1,3)*moment(n,2,1) discrepancy=abs(value/expected-1) assert discrepancy < mp.power(10,-working+10) period_controls[str(n)] = mp.nstr(discrepancy,8) return {'working_digits':working,'cutoff':cutoff,'bessel_max_bits':max_bits[0], 'bessel_nodes':len(cache),'values':values,'moments':moments, 'quad_error_estimates':errors,'tail_upper_estimates':tails, 'controls_relative_error':controls,'period_controls_relative_error':period_controls, 'determinant_cancellation_digits':cancellation} class VacuumDeterminants(numberdb.Generator): table = 'T648' parameters = ('n',) type = 'R' digits = 100 rigour = 'heuristic (agreement-checked)' files = ('generate.py','checks.py','argument-selection.md','skill-manifest.json','numerical-evidence.json', 'independent-evidence.json','program.py') def enumerate(self): for n in SELECTION: yield {'n':str(n)} def value(self,params,digits): n=str(params['n']) if n not in ('5','6','8') or digits>100: raise ValueError('This generator supports n=5,6,8 and at most 100 digits') return numberdb.agreeing(lambda w:compute_all(w)['values'][n], at=WORKING) if __name__=='__main__': if '--compute' in sys.argv: runs=[compute_all(w) for w in WORKING] with mp.workdps(220): differences={n:mp.nstr(abs(mp.mpf(runs[0]['values'][n])/mp.mpf(runs[1]['values'][n])-1),10) for n in ('5','6','8')} assert all(mp.mpf(v)