"""Independent checks for T693, runnable as: sage -python checks.py. With NUMBERDB_STORED=1, read T693 and test the written intervals as well. The separate formal-series computation uses reciprocal/exponential convolutions, not the generator's Bernoulli routine or logarithmic-derivative recurrence. The Kummer differential equation and direct Arb U computations are additional checks external to the defining finite sum. See derivation.md for limitations. """ import json import os import sys from math import factorial from functools import reduce import numberdb.sage as numberdb from sage.rings.rational_field import QQ from sage.rings.real_arb import RealBallField from sage.rings.real_mpfi import RealIntervalField from sage.rings.complex_arb import ComplexBallField from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.polynomial.laurent_polynomial_ring import LaurentPolynomialRing from sage.env import SAGE_VERSION from numberdb._write import to_text from check import exactness, measure, agrees_with, names_its_rings import generate N = 60 S = PolynomialRing(QQ, 'z') z = S.gen() def independent_polynomials(limit): # Invert (exp(u)-1)/u using exact convolution, avoiding bernoulli(). inv = [QQ.one()] for k in range(1, limit+2): inv.append(-sum(inv[k-r]/QQ(factorial(r+1)) for r in range(1,k+1))) lam = [QQ.zero()]+inv[2:limit+2] # Sum (-z lambda(u))^m/m! by truncated Cauchy products. power = [QQ.one()]+[QQ.zero()]*limit p = [S.one()]+[S.zero()]*limit for m in range(1,limit+1): power = [sum(power[k-r]*lam[r] for r in range(1,k+1)) for k in range(limit+1)] for k in range(m,limit+1): p[k] += power[k]*(-z)**m/QQ(factorial(m)) c=[] R=PolynomialRing(QQ,'x'); x=R.gen() for ell in range(limit+1): v=S.zero() for k in range(ell+1): j=ell-k # H(nu,j)=(-1)^j (1/2-nu)_j (1/2+nu)_j/j!. poch=QQ.one() for r in range(j): poch *= (QQ(1)/2-(k+1)+r)*(QQ(1)/2+(k+1)+r) v += QQ((-1)**j)*poch/QQ(factorial(j)*4**j)*z**k*p[k] assert all(v[r]==0 for r in range(1,v.degree()+1,2)) c.append(R([v[2*m]*(-4)**m for m in range(ell+1)])) return c,inv def written_bounds(text): from decimal import Decimal if '/' in text: return (QQ(text),QQ(text)) d=Decimal(text); n,dn=d.as_integer_ratio(); center=QQ(n)/QQ(dn) if '.' not in text and 'e' not in text.lower(): return center,center e=d.as_tuple().exponent ulp=QQ(10)**e return center-ulp,center+ulp def encloses(text,v): lo,hi=written_bounds(text) if isinstance(v,type(QQ.one())):return lo<=v<=hi return bool(v.is_finite() and lo < QQ(v.lower()) and QQ(v.upper()) < hi) def kummer_residuals(polys): L=LaurentPolynomialRing(QQ,'t');t=L.gen() a=[sum(L(c)*QQ(-1)**m/QQ(4)**m*t**(4*m-ell) for m,c in enumerate(p.list())) for ell,p in enumerate(polys)] return [4*a[l+1].derivative()-a[l].derivative(2)+(t*t+QQ(3)/4/t**2)*a[l] for l in range(len(a)-1)] def run(): if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1': numberdb.configure(api_key=sys.stdin.read().strip()) out={'table':'T693','sage':SAGE_VERSION,'numberdb':numberdb.__version__, 'skill_manifest':'agents/runs/20261109T154709Z-build-skills/manifest.json', 'skill_sha256':'886390324dd1971e381a43cf3b46b822ab0bc67d62f31d8481e29ec79556746b'} ps=generate.coefficient_polynomials(N) ips,inv=independent_polynomials(N) assert inv[:4]==[QQ(1),-QQ(1)/2,QQ(1)/12,QQ(0)] assert encloses('1.25',QQ(5)/4) and not encloses('1.25',QQ(3)/2) assert not exactness(ps) and not exactness(ips) assert not names_its_rings(generate.__file__) assert not agrees_with(dict(enumerate(ps)),lambda l:ips[l]) bad=list(ps[:2]);bad[1]=bad[1]+1 assert any(kummer_residuals(bad)) assert all(r==0 for r in kummer_residuals(ps)) expected=[ [1], [QQ(3)/16,QQ(1)/3],[-QQ(15)/512,QQ(5)/16,QQ(1)/18], [QQ(105)/8192,QQ(35)/512,QQ(23)/160,QQ(1)/162], [-QQ(4725)/524288,-QQ(105)/8192,QQ(973)/5120,QQ(137)/4320,QQ(1)/1944], [QQ(72765)/8388608,QQ(3465)/524288,QQ(2849)/16384,QQ(17365)/193536,QQ(229)/51840,QQ(1)/29160] ] assert all(ps[i]==generate.R(p) for i,p in enumerate(expected)) const=QQ.one() for ell,p in enumerate(ps): if ell:const*=QQ(4-(2*ell-1)**2)/QQ(16*ell) assert p.degree()==ell and p[ell]==QQ(1)/QQ(factorial(ell)*3**ell) and p[0]==const R=RealBallField(700); pi2=R.pi()**2 main={i:generate.evaluate(p) for i,p in enumerate(ps)} other={i:ips[i](pi2) for i in range(N+1)} texts={i:to_text(v,digits=100) for i,v in main.items()} for i in range(N+1): assert other[i].is_finite() and (i==0 or other[i].accuracy()>650) assert encloses(texts[i],main[i]) and (i==0 or encloses(texts[i],other[i])) neg=[i for i in range(N+1) if main[i]<0] assert neg[0]==50 and all(main[i]>0 for i in range(50)) cancellation=[] for p,v in zip(ps,main.values()): if p.degree()>0: abs_sum=sum(abs(c)*pi2**k for k,c in enumerate(p.list())) cancellation.append(float((abs_sum/abs(R(v))).log()/R(10).log())) out.update(controls='reciprocal series, exact interval inclusion, corrupted coefficient and nonzero U identity passed', entries=N+1,exact_polynomial_agreements=N+1,printed_formula_agreements=6,kummer_ode_residuals=60, minimum_accuracy_bits=min(int(v.accuracy()) for i,v in main.items() if i), working_bits=numberdb.bits(100,losing=generate.WORKING_GUARD), minimum_independent_accuracy_bits=min(int(other[i].accuracy()) for i in range(1,N+1)), max_cancellation_decimal_digits=max(cancellation),negative_orders=neg, raw_ball_measurement=measure(main),written_measurement=measure(texts), values={str(k):v for k,v in texts.items()},polynomials={str(i):str(p) for i,p in enumerate(ps)}) # Corpus ingredients read through numberdb, never altered. bern=numberdb.table('T13')['Numbers'] hank=numberdb.table('T569')['Numbers'] def entries(block): if isinstance(block,list): return [(next(iter(r['params'].values())),r['number']) for r in block] return [(k,v['number'] if isinstance(v,dict) else v) for k,v in block.items()] bc=0 for k,v in entries(bern): k=int(k) if 0<=k<=N+1: assert QQ(str(v))==inv[k]*QQ(factorial(k));bc+=1 hc=0 Hring=PolynomialRing(QQ,'v');v=Hring.gen() for j,expr in entries(hank): j=int(j) if j>N:continue poly=Hring(str(expr)) for k in range(N-j+1): h=QQ.one() for r in range(1,j+1):h*=QQ(4*(k+1)**2-(2*r-1)**2)/QQ(4*r) assert 2**j*poly(k+1)==h;hc+=1 out['bernoulli_ingredient_checks']=bc;out['hankel_ingredient_checks']=hc # Direct complex Tricomi values validate amplitude, phase, and truncation. C=ComplexBallField(1200);rb=RealBallField(1200);zz=2*C.pi()*C.gen(0) control=C(2).hypergeometric_U(1,2) assert control.real().is_finite() and control.imag().is_finite() and (QQ(1)/2 in control.real()) and (0 in control.imag()) assert control.real().rad() < QQ(2)**-1000 and control.imag().rad() < QQ(2)**-1000 out['tricomi_control']={'real':str(control.real()),'imaginary':str(control.imag())} used=main if os.environ.get('NUMBERDB_STORED')=='1': pairs=entries(numberdb.table('T693')['Numbers']) assert [int(k) for k,v in pairs]==list(range(N+1)) used={} for k,text in pairs: k=int(k); assert encloses(str(text),main[k]) if k: assert encloses(str(text),other[k]) lo,hi=written_bounds(str(text)) used[k]=rb(RealIntervalField(1200)(lo,hi)) assert used[k].is_finite() and (k==0 or used[k].accuracy()>320) assert [i for i in range(N+1) if used[i]<0]==neg assert str(pairs[0][1])=='1' out['stored_interval_checks']=len(pairs) out['stored_sign_checks']=len(pairs) tests=[] for n in [25,100,400,1600]: val=(-4*C(factorial(n-1))*zz.hypergeometric_U(n,0)).real() assert val.is_finite() and val.accuracy()>500 r=rb(n);pi=rb.pi();phase=2*(pi*r).sqrt() amplitude=rb(2)**(QQ(9)/4)*pi**(QQ(3)/4)*r**(-QQ(3)/4)*(-phase).exp() target=val/amplitude for L in [0,1,3,5,10,20,40,60]: approx=sum(rb(used[l])*(2*pi*r)**(-QQ(l)/2)*(phase+pi*l/4+3*pi/8).sin() for l in range(L+1)) error=abs(target-approx) tests.append({'n':n,'L':L,'normalized_error':str(error),'error_mid':float(error.mid()),'U_accuracy_bits':int(val.accuracy())}) out['tricomi_tests']=tests for n in [25,100,400,1600]: errors={row['L']:row['error_mid'] for row in tests if row['n']==n} assert errors[5]