back to table · edit · history · where entries came from · files · download
9143 bytes, as of the version from 2026-11-09 23:51 (current). Recorded here, not run.
"""Independent checks of Generalized log-sine integrals -- numberdb.org/T697.
sage -python checks.py # check all computed entries, JSON to stdout
sage -python checks.py --stored # additionally check the authenticated draft
The numerical quadratures and finite polylogarithm/Gamma computations share no
numerical code with the series generator. The generator supplies its enumeration
and values only. Instructions: skill-manifest.json, originally archived at
agents/runs/20261109T220207Z-build-skills/manifest.json.
"""
import os
import sys
import json
from decimal import Decimal
from fractions import Fraction
from math import factorial
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.rational_field import QQ
import mpmath as mp
from generate import LogSine, evaluate, WORKING_GUARD_BITS, TAIL_GUARD_DIGITS
MANIFEST = 'agents/runs/20261109T220207Z-build-skills/manifest.json'
SKILL_SHA256 = '886390324dd1971e381a43cf3b46b822ab0bc67d62f31d8481e29ec79556746b'
def quadrature(n,k,t,dps):
with mp.workdps(dps):
theta=mp.pi*mp.mpf(t)
if n-1-k==0:
return -mp.quad(lambda v:v**k,[0,theta/2,theta])
return -mp.quad(lambda v:v**k*mp.log(2*mp.sin(v/2))**(n-1-k),[0,theta/2,theta])
def gamma_coefficients(R):
c=[R(0),R(0)]+[(-1)**j*(1-R(2)**(1-j))*R(j).zeta()/j for j in range(2,10)]
f=[R(1)]
for m in range(1,10):
f.append(sum((j*c[j]*f[m-j] for j in range(1,m+1)),R(0))/m)
return [-R.pi()*factorial(m)*f[m] for m in range(10)]
def one_log(k,t,C):
z=(C.pi()*C(QQ(t))*C.gen(0)).exp()
theta=C.pi()*C(QQ(t)); i=C.gen(0)
out=sum(((-1)**j*factorial(k)//factorial(k-j)*theta**(k-j)*i**(-j-1)*z.polylog(j+2)
for j in range(k+1)),C(0))
out-=(-1)**k*factorial(k)*i**(-k-1)*C(k+2).zeta()
return out.real()
def interval(text):
# Exact decimal parsing, including its specified one-last-place uncertainty.
d=Decimal(str(text)); center=Fraction(d)
if '.' not in str(text) and 'e' not in str(text).lower():return QQ(center),QQ(center)
radius=Fraction(10)**d.as_tuple().exponent
return QQ(center-radius),QQ(center+radius)
def contains(text,ball):
lo,hi=interval(text)
if not hasattr(ball,'lower'):return lo<=ball<=hi
return lo<=ball.lower() and ball.upper()<=hi
def records(document):
block=document['Numbers']
if isinstance(block,list):return block
out=[]
for n,ks in block.items():
for k,ts in ks.items():
for t,v in ts.items():
row={'params':{'n':str(n),'k':str(k),'t':str(t)}}
row.update(v if isinstance(v,dict) else {'number':v});out.append(row)
return out
def main(stored=False):
mp.mp.dps=180
R=RealBallField(700);C=ComplexBallField(700)
# Controls precede measurements, including a deliberate wrong-answer rejection.
for t in ('1/3','1/2','2/3','5/6','1'):
for n in (1,4,10):
val=quadrature(n,n-1,t,130);theta=mp.pi*mp.mpf(t)
assert abs(val+theta**n/n)<mp.mpf('1e-120')*max(1,abs(val))
assert not contains('1.00000000000000000000',R(2))
rawzero=evaluate(2,0,'1',140,exact_zero=False)
assert rawzero.contains_zero() and R(rawzero.rad())<R(10)**-170
gamma=gamma_coefficients(R)
rows=[];values={};high={};one_count=0;gamma_count=0
worst_quad=mp.mpf(0);worst_precision=mp.mpf(0)
for index,p in enumerate(LogSine().enumerate()):
n,k,t=int(p['n']),int(p['k']),p['t']; key=(n,k,t)
v=evaluate(n,k,t,100); vh=evaluate(n,k,t,140)
values[key]=v;high[key]=vh
text=to_text(v,100)
assert contains(text,vh),(key,'high precision outside decimal')
if hasattr(v,'overlaps'):
assert v.is_finite() and vh.is_finite() and v.overlaps(vh)
mid=mp.mpf(str(vh.mid()))
else:mid=mp.mpf(int(vh))
q1=quadrature(n,k,t,130);q2=quadrature(n,k,t,170)
scale=max(1,abs(mid))
e=abs(mid-q2)/scale;ep=abs(q1-q2)/scale
assert e<mp.mpf('1e-120') and ep<mp.mpf('1e-115'),(key,mp.nstr(e),mp.nstr(ep))
worst_quad=max(worst_quad,e);worst_precision=max(worst_precision,ep)
if k==0 and t=='1':
expected=gamma[n-1]
assert expected.is_finite() and (expected==0 if key==(2,0,'1') else expected.accuracy()>600)
assert expected.contains_zero() if key==(2,0,'1') else vh.overlaps(expected)
gamma_count+=1
if n-1-k==1:
expected=one_log(k,t,C)
assert expected.is_finite()
if key==(2,0,'1'):assert expected.contains_zero() and R(expected.rad())<R(10)**-190
else:assert expected.accuracy()>600 and vh.overlaps(expected),(key,'polylog moment')
one_count+=1
rows.append({'params':p,'number':text,'accuracy_bits':int(v.accuracy()) if hasattr(v,'accuracy') else None,
'radius':str(v.rad()) if hasattr(v,'rad') else '0',
'quadrature_130':mp.nstr(q1,125),'quadrature_170':mp.nstr(q2,165),
'relative_quadrature_difference':mp.nstr(e,8)})
if (index+1)%30==0:print('checked',index+1,'entries',file=sys.stderr,flush=True)
# Known low-order values, independently supplied by Borwein–Straub Example 1.
assert high[(3,0,'1')].overlaps(-R.pi()**3/12)
assert high[(4,0,'1')].overlaps(3*R.pi()*R(3).zeta()/2)
# Specific companion relationships, using their current 100-digit intervals.
companions={tid:numberdb.table(tid) for tid in ('T175','T681','T685')}
def val(d,*keys):
v=d['Numbers']
for key in keys:v=v[str(key)]
return v['number'] if isinstance(v,dict) else v
companion_count=0
for t in ('1/3','1/2','2/3'):
assert contains(val(companions['T175'],'2',t),high[(2,0,t)]);companion_count+=1
for r in (2,3):
assert contains(val(companions['T681'],'1',r),-high[(r+1,0,'1')]/R.pi());companion_count+=1
for r in range(1,10):
expected=(high[(r+1,0,'1/3')]-high[(r+1,0,'1')])/R.pi()
assert contains(val(companions['T685'],'3',r),expected);companion_count+=1
stored_count=0
stored_identity_checks=0
if stored:
document=numberdb.table('T697');rs=records(document)
assert [r['params'] for r in rs]==list(LogSine().enumerate())
back={}
for row in rs:
p=row['params'];key=(int(p['n']),int(p['k']),p['t'])
assert contains(row['number'],high[key]),key
# Read-back values are separately compared to the independent integral.
lo,hi=interval(row['number']);q=mp.mpf(next(r['quadrature_170'] for r in rows if r['params']==p))
if key!=(2,0,'1'):
assert mp.mpf(str(lo.numerator()))/int(lo.denominator())<q<mp.mpf(str(hi.numerator()))/int(hi.denominator())
else:assert str(row['number'])=='0'
center=(R(lo)+R(hi))/2
back[key]=center.add_error(((R(hi)-R(lo))/2).upper())
stored_count+=1
for key,v in back.items():
n,k,t=key
if k==0 and t=='1':
assert v.overlaps(gamma[n-1]);stored_identity_checks+=1
if n-1-k==1:
assert v.overlaps(one_log(k,t,C));stored_identity_checks+=1
if (n-1-k)%2==0:assert v<0
for r in (2,3):
lo,hi=interval(val(companions['T681'],'1',r))
v=-back[(r+1,0,'1')]/R.pi()
assert v.lower()<=hi and lo<=v.upper();stored_identity_checks+=1
for r in range(1,10):
lo,hi=interval(val(companions['T685'],'3',r))
v=(back[(r+1,0,'1/3')]-back[(r+1,0,'1')])/R.pi()
assert v.is_finite() and v.accuracy()>280
assert v.lower()<=hi and lo<=v.upper();stored_identity_checks+=1
evidence={'table':'T697','skill_manifest':MANIFEST,'skill_sha256':{'local':SKILL_SHA256,'live':SKILL_SHA256},
'entries':len(rows),'minimum_accuracy_bits':min(r['accuracy_bits'] for r in rows if r['accuracy_bits'] is not None),
'working_bits':numberdb.bits(100,losing=WORKING_GUARD_BITS),'series_degree':numberdb.bits(100+TAIL_GUARD_DIGITS,losing=0),
'higher_precision_digits':140,'quadrature_working_digits':[130,170],
'maximum_scaled_quad_error':mp.nstr(worst_quad,10),'maximum_scaled_quad_precision_difference':mp.nstr(worst_precision,10),
'controls':15,'negative_control_rejected':True,'exact_zero_raw_radius':str(rawzero.rad()),
'gamma_checks':gamma_count,'finite_polylog_checks':one_count,'companion_checks':companion_count,
'stored_intervals_checked':stored_count,'stored_identity_checks':stored_identity_checks,'versions':{'numberdb':numberdb.__version__,'mpmath':mp.__version__},'rows':rows}
print(json.dumps(evidence,indent=2))
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
main('--stored' in sys.argv or os.environ.get('LOGSINE_STORED')=='1')