back to table · edit · history · where entries came from · files · download
9933 bytes, as of the version from 2026-11-09 23:28. Recorded here, not run.
"""Independent checks for Log-sine-cosine integrals, numberdb.org/T699.
Run: sage -python checks.py [--stored]
Uses mpmath quadrature, independently expanded beta coefficients, and a
harmonic central-binomial series. --stored also verifies every read-back
interval and the companion-table identity. Credentials use the ordinary
numberdb configuration; NUMBERDB_KEY_FROM_STDIN=1 is supported for the runner.
The final EVIDENCE_JSON line is the complete numerical record.
"""
import os
import sys
import json
import platform
from decimal import Decimal
from math import comb, factorial
from pathlib import Path
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.env import SAGE_VERSION
import mpmath as mp
from generate import LogSineCosine, compute
def mpq(q):
return mp.mpf(str(q.numerator()))/mp.mpf(str(q.denominator()))
def interval(text):
d=Decimal(text); a,b=d.as_integer_ratio(); c=QQ(a)/b
ulp=QQ(10)**d.as_tuple().exponent
return c-ulp,c+ulp
def inside(text,value):
lo,hi=interval(text)
if hasattr(value,'lower'):
return bool(value.lower()>lo and value.upper()<hi)
return mpq(lo)<value<mpq(hi)
def quad(m,n,t,dps):
"""Direct trigonometric quadrature; no series or coefficient code shared."""
with mp.workdps(dps):
t=mpq(QQ(t)); theta=mp.pi*t
def f(v):
if v==0 or v==mp.pi:
return mp.mpf(0)
return -mp.log(2*mp.sin(v/2))**(m-1)*mp.log(abs(2*mp.cos(v/2)))**(n-1)
splits=sorted(set([mp.mpf(0),theta]+[x for x in [mp.pi/3,mp.pi/2,2*mp.pi/3] if x<theta]))
return +mp.quad(f,splits,method='tanh-sinh')
def beta_moment(m,n,t,dps=150):
"""Differentiate the incomplete beta integral, then set x=r^2."""
with mp.workdps(dps):
z=mp.sin(mp.pi*mpq(QQ(t))/2)
def integrand(r):
if r==0:return mp.mpf(0)
return -2*mp.log(2*r)**(m-1)*(mp.log(2)+mp.log1p(-r*r)/2)**(n-1)/mp.sqrt(1-r*r)
return +mp.quad(integrand,[0,z/2,z])
def beta_coefficients(degree):
"""Coefficients of exp(log(G/pi)), independent of the generator's series."""
L={}
for d in range(2,degree+1):
c=(-1)**d*mp.zeta(d)/(mp.mpf(d)*2**d)
for a in range(d+1):
b=d-a
L[a,b]=c*((2**d-1 if a==0 or b==0 else 0)-comb(d,a))
G={(0,0):mp.mpf(1)}
for d in range(1,degree+1):
for a in range(d+1):
b=d-a
G[a,b]=sum((mp.mpf(i+j)/d*l*G.get((a-i,b-j),mp.mpf(0))
for (i,j),l in L.items() if i<=a and j<=b),mp.mpf(0))
return G
def harmonic_series(dps=180):
with mp.workdps(dps):
# h=H_(2j), c=binomial(2j,j)/16^j; ratios use exact integer factors.
h=mp.mpf(0); c=mp.mpf(1); total=mp.mpf(0)
for j in range(1,2000):
c*=mp.mpf(2*j*(2*j-1))/(16*j*j)
h+=mp.mpf(1)/(2*j-1)+mp.mpf(1)/(2*j)
term=c*h/(2*j+1)**3; total+=term
if abs(term)<mp.mpf(10)**(-dps-10):
return +total,j
raise AssertionError('series did not converge')
def records(doc):
block=doc['Numbers']
if isinstance(block,list):
return {tuple(str(r['params'][p]) for p in ('m','n','t')):r['number'] for r in block}
groups=doc.get('Display properties',{}).get('group parameters',[['m'],['n'],['t']])
result={}
def walk(node,depth,params):
if depth==len(groups):
result[tuple(params[p] for p in ('m','n','t'))]=node['number'] if isinstance(node,dict) else node
return
for key,value in node.items():
parts=[x.strip() for x in key.split(',')]
assert len(parts)==len(groups[depth])
walk(value,depth+1,dict(params,**dict(zip(groups[depth],parts))))
walk(block,0,{})
return result
def boundary_checks():
# Pure-power relations across every selected angular value and order prefix.
count=0
with mp.workdps(150):
for order in range(2,9):
full=quad(order,1,'1',140)
for t in ['1/6','1/3','1/2','5/6','1']:
direct=quad(1,order,t,140)
complement=1-QQ(t)
reflected=quad(order,1,str(complement),140) if complement else mp.mpf(0)
assert abs(direct-(full-reflected))<mp.mpf('1e-115')*max(1,abs(direct))
base=quad(order,1,t,140)
certified=compute(order,1,t,120)
assert abs(base-mp.mpf(str(certified.mid())))<mp.mpf('1e-115')*max(1,abs(base))
count+=1
return count
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
key=sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].isidentifier():key=key.split('=',1)[1].strip().strip('\"\'')
os.environ['NUMBERDB_API_KEY']=key
stored='--stored' in sys.argv or os.environ.get('T699_CHECK_STORED')=='1'
mp.mp.dps=190
# Known nonzero controls and a deliberately wrong answer, before measurements.
assert abs(mp.quad(lambda x:x*x,[0,1])-mp.mpf(1)/3)<mp.mpf('1e-185')
control=quad(1,1,'1/2',130)
assert abs(control+mp.pi/2)<mp.mpf('1e-125')
assert abs(quad(2,2,'1',130)-mp.pi**3/24)<mp.mpf('1e-120')
assert not inside(to_text(compute(2,2,'1',100),100),mp.pi**3/24+mp.mpf('1e-95'))
print('Known-answer and wrong-answer controls passed',flush=True)
gen=LogSineCosine(); params=list(gen.enumerate())
assert len(params)==48 and len({tuple(p.values()) for p in params})==48
values={}; evidence=[]; minbits=100000
for p in params:
key=tuple(p[k] for k in ('m','n','t'))
m,n,t=int(p['m']),int(p['n']),p['t']
low=gen.value(p,100); high=gen.value(p,140); text=to_text(low,100)
minbits=min(minbits,int(low.accuracy()))
assert inside(text,low) and inside(text,high)
q130=quad(m,n,t,130); q180=quad(m,n,t,180)
assert inside(text,q130) and inside(text,q180),key
relative=abs(q130-q180)/max(abs(q180),mp.mpf('1e-100'))
assert relative<mp.mpf('1e-115'),(key,relative)
values[key]=text
evidence.append({'params':p,'number':text,'accuracy_bits':int(low.accuracy()),
'quad_130':mp.nstr(q130,140),'quad_180':mp.nstr(q180,190),
'quadrature_relative_difference':mp.nstr(relative,8)})
print('checked',','.join(key),'bits',low.accuracy(),flush=True)
if stored:
readback=records(numberdb.table('T699'))
assert set(readback)==set(values)
assert readback==values
assert list(readback)==list(values)
values=readback
print('Read-back entries match all 48 certified decimal intervals',flush=True)
incomplete_count=0
for key,text in values.items():
if key[2]=='1':continue
assert inside(text,beta_moment(int(key[0]),int(key[1]),key[2])),key
incomplete_count+=1
G=beta_coefficients(8); endpoint_count=0; midpoint_count=0
for key,text in values.items():
m,n,t=int(key[0]),int(key[1]),key[2]
if t!='1':continue
expected=-mp.pi*factorial(m-1)*factorial(n-1)*G[m-1,n-1]
assert inside(text,expected),key
endpoint_count+=1
a,b=values[str(m),str(n),'1/2'],values[str(n),str(m),'1/2']
alo,ahi=interval(a); blo,bhi=interval(b); flo,fhi=interval(text)
assert mpq(alo+blo)<expected<mpq(ahi+bhi)
assert max(alo+blo,flo)<min(ahi+bhi,fhi)
midpoint_count+=1
# Published exceptional-angle reflection and harmonic-binomial identity.
full32=-mp.pi*mp.zeta(3)/4
reflection=mp.mpf(values['2','3','1/6'])+mp.mpf(values['3','2','5/6'])
assert abs(reflection-full32)<mp.mpf('3e-99')
H,terms=harmonic_series()
# A separate direct integral supplies Ls4(5*pi/6) before database read-back.
ls4=quad(4,1,'5/6',180)
rhs=2*ls4+mp.mpf(values['3','2','1/3'])-2*mp.mpf(values['3','2','1/6'])+4*mp.mpf(values['3','2','5/6'])-2*mp.pi*mp.zeta(3)
assert abs(H-rhs)<mp.mpf('1e-97')
companion=False
if stored:
doc=numberdb.table('T697'); z=doc['Numbers']['4']['0']['5/6']
if isinstance(z,dict):z=z['number']
lo,hi=interval(z); assert mpq(lo)<ls4<mpq(hi)
terms_rows=[(2,z),(1,values['3','2','1/3']),(-2,values['3','2','1/6']),(4,values['3','2','5/6'])]
lower=upper=mp.mpf(0)
for c,s in terms_rows:
a,b=interval(s); a,b=mpq(a),mpq(b)
lower+=c*(a if c>0 else b);upper+=c*(b if c>0 else a)
offset=2*mp.pi*mp.zeta(3)
assert lower-offset<H<upper-offset
companion=True
report=gen.verify(sample=None);print(report,flush=True);assert report.ok
boundary_count=boundary_checks()
manifest=json.loads(Path(__file__).with_name('skill-manifest.json').read_text())
result={'table':'T699','skill_manifest':'agents/runs/20261109T231810Z-build-skills/manifest.json',
'archived_skill_versions':manifest['versions'],'entries':48,'stored_checked':stored,
'working_bits':525,'series_degree':557,'minimum_accuracy_bits':minbits,
'quadrature_precisions':[130,180],'beta_endpoint_checks':endpoint_count,'incomplete_beta_moment_checks':incomplete_count,
'midpoint_reflection_checks':midpoint_count,'pure_power_boundary_checks':boundary_count,'exceptional_reflection_checks':1,
'harmonic_series_terms':terms,'harmonic_identity_residual':mp.nstr(H-rhs,12),
'companion_T697_interval_identity_checked':companion,
'known_answer_controls_passed':True,'wrong_answer_control_rejected':True,
'versions':{'python':platform.python_version(),'sage':SAGE_VERSION,'mpmath':mp.__version__,'numberdb':numberdb.__version__},
'values':evidence}
print('EVIDENCE_JSON '+json.dumps(result),flush=True)
if __name__=='__main__':main()