back to table · edit · history · where entries came from · files · download
7365 bytes, as of the version from 2026-10-28 15:11 (current). Recorded here, not run.
"""Reproducible checks for Gram points -- numberdb.org/T583.
sage -python checks.py # local computation
sage -python checks.py --stored # read the draft; authenticate
The campaign's agents/sage.sh can instead set NUMBERDB_KEY_FROM_STDIN=1;
that selects stored mode and reads its credential only from stdin. Output is
JSON evidence. The independent mpmath implementation is used at 150 and 200
places. The Arb endpoint and phase-identity tests use fresh 600-bit fields.
No website writes are performed. The main generator is needed beside this file.
"""
import sys
import os
import json
import time
import hashlib
import pathlib
import platform
import importlib.metadata
from fractions import Fraction
from decimal import Decimal
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.rings.complex_arb import ComplexBallField
from generate import GramPoints, client_from_environment
import mpmath as mp
R = RealBallField(600)
C = ComplexBallField(600)
def theta(t):
return C(QQ(1)/4,t/2).log_gamma().imag()-t*R.pi().log()/2
def derivative(t):
return (C(QQ(1)/4,t/2).psi().real()-R.pi().log())/2
def rational(text):
f = Fraction(text)
return QQ(f.numerator)/f.denominator
def decimal_interval(text):
value = Decimal(text.rstrip('?'))
assert value.is_finite()
ulp = rational(str(Decimal(1).scaleb(value.as_tuple().exponent)))
center = rational(str(value))
return center-ulp,center+ulp,ulp
def certificate(n,left,right):
fleft = theta(R(left))-n*R.pi()
fright = theta(R(right))-n*R.pi()
return (fleft.is_finite() and fright.is_finite()
and fleft < 0 < fright)
def mp_rational(q):
return mp.mpf(int(q.numerator()))/int(q.denominator())
def check_rows(rows):
# Known-answer controls precede every root or phase check.
assert theta(R(0)).is_zero()
control = theta(R(1))
assert control.is_finite() and -R(2) < control < -R(1)
assert abs(control-rational('-1.767547952812290388302216')) < R('1e-24')
assert C(1).log_gamma().is_zero()
mp.mp.dps = 150
assert mp.siegeltheta(0) == 0
assert abs(mp.siegeltheta(1)-mp.mpf('-1.767547952812290388302216')) < mp.mpf('1e-24')
assert abs(mp.grampoint(0)-mp.mpf('17.84559954041086081682634')) < mp.mpf('1e-23')
assert derivative(R(7)) > 0 and theta(R(7)) < -R.pi()
assert list(rows) == [str(n) for n in range(-1,201)]
failures = []
evidence = []
previous_right = QQ(7)
minimum_bits = None
for key,text in rows.items():
n = int(key)
left,right,ulp = decimal_interval(text)
assert len(Decimal(text.rstrip('?')).as_tuple().digits) == 100
assert left > previous_right
previous_right = right
assert certificate(n,left,right), ('stored decimal bracket',n)
assert not certificate(n,left+10*ulp,right+10*ulp)
assert not certificate(n,left-10*ulp,right-10*ulp)
interval = R(left).union(R(right))
d = derivative(interval)
assert d.is_finite() and d > 0
minimum_bits = min(int(interval.accuracy()),minimum_bits or 10**9)
# Independent implementation at two precisions, checked at full
# stored interval precision, with additional 140-digit agreement.
with mp.workdps(150):
g150 = mp.grampoint(n)
assert mp_rational(left) < g150 < mp_rational(right)
saved150 = mp.nstr(g150,150)
with mp.workdps(200):
g200 = mp.grampoint(n)
assert mp_rational(left) < g200 < mp_rational(right)
assert abs(g200-mp.mpf(saved150)) < mp.mpf('1e-140')
# Gamma itself gives a separate phase computation, without log-gamma.
gamma = C(QQ(1)/4,interval/2).gamma()
unit_gamma = gamma/abs(gamma)
parity = 1 if n%2==0 else -1
expected = parity*C(0,interval*R.pi().log()/2).exp()
difference = unit_gamma-expected
for part in [difference.real(),difference.imag()]:
assert part.is_finite() and part.contains_zero()
assert abs(part) < R('1e-90')
zeta = C(QQ(1)/2,interval).zeta()
assert zeta.real().is_finite() and zeta.imag().is_finite()
assert zeta.imag().contains_zero() and abs(zeta.imag()) < R('1e-90')
assert zeta.real().rad_as_ball() < R('1e-90')
good = zeta.real() > 0
assert good or zeta.real() < 0
if not good:
failures.append(n)
with mp.workdps(150):
zvalue = mp.siegelz(mp.grampoint(n))
assert abs(zvalue-parity*mp.zeta(mp.mpf('0.5')+1j*mp.grampoint(n))) < mp.mpf('1e-140')
sage_z = parity*zeta.real()
assert mp_rational(sage_z.lower().exact_rational()) < zvalue < mp_rational(sage_z.upper().exact_rational())
evidence.append({'n':n,'number':text,'endpoint_signs':[-1,1],
'derivative_positive':True,'gram_sign':1 if good else -1})
return {'entries':len(rows),'endpoint_certificates':len(rows),
'perturbed_brackets_rejected':2*len(rows),
'mpmath_150_and_200_checks':len(rows),
'gamma_phase_checks':len(rows),'zeta_real_and_hardy_z_checks':len(rows),
'minimum_decimal_interval_accuracy_bits':minimum_bits,
'gram_sign_failures':failures,'rows':evidence,
'controls':{'theta_zero':True,'theta_one_nonzero':to_text(control,20),
'loggamma_one':True,'mpmath_g0_reference':True,
'theta_prime_seven':to_text(derivative(R(7)),20),
'theta_seven_less_than_minus_pi':True}}
def main(stored=False):
start=time.monotonic()
gen=GramPoints()
if stored:
client=client_from_environment()
doc=numberdb.table('T583',client=client)
block=doc['Numbers']
if isinstance(block,list):
rows={str(e['params']['n']):e['number'] for e in block}
else:
rows={str(k):v['number'] if isinstance(v,dict) else v for k,v in block.items()}
else:
balls={p['n']:gen.value(p,100) for p in gen.enumerate()}
rows={k:to_text(v,100) for k,v in balls.items()}
report=check_rows(rows)
if not stored:
report['minimum_generator_accuracy_bits']=min(int(v.accuracy()) for v in balls.values())
report.update(table='T583',mode='stored' if stored else 'local',
skill_manifest='agents/runs/20261028T143558Z-build-skills/manifest.json',
skill_sha256='e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c',
versions={'Python':platform.python_version(),
'numberdb':importlib.metadata.version('numberdb'),
'mpmath':mp.__version__,
'SageMath':__import__('sage.version',fromlist=['version']).version},
source_hashes={name:hashlib.sha256(pathlib.Path(__file__).with_name(name).read_bytes()).hexdigest()
for name in ['generate.py','checks.py','argument-selection.md','certification.md','skill-manifest.json']},
elapsed_seconds=round(time.monotonic()-start,3))
print(json.dumps(report,indent=2))
if __name__=='__main__':
main('--stored' in sys.argv or os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1')