back to table · edit · history · where entries came from · files · download
9905 bytes, as of the version from 2026-10-30 08:21 (current). Recorded here, not run.
"""Independent checks for numberdb.org/T621.
sage -python checks.py # pre-fill, no key or table needed
sage -python checks.py --stored # also check the live stored values
Requires generate.py beside this file, mpmath, SageMath, and numberdb.
The corner integral (Aurell--Salomonson, equations 18 and 20) is evaluated
at 150 and 200 decimal digits independently of the finite Dedekind sums.
The small-x series cancels the singular terms algebraically. Its truncation
and the large-x tail have explicit bounds (see derivation.md). Quadrature
on the remaining finite interval is agreement-checked, not certified.
The production certification comes from the finite Arb formula.
Square and triangle controls come from independently separated spectra.
The integral and product routines are adapted from T620 checks.py (downloaded
2026-10-30), adding Neumann axis modes and the polygon determinant ratio.
"""
import sys
import os
import json
import math
import time
from decimal import Decimal
from fractions import Fraction
from functools import lru_cache
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
import mpmath as mp
import generate
from sage.env import SAGE_VERSION
def decimal_interval(text):
d = Decimal(str(text))
center = QQ(Fraction(d).numerator) / Fraction(d).denominator
ulp = QQ(10)**d.as_tuple().exponent
return center-ulp, center+ulp, ulp
def contains_ball(text, ball):
lower, upper, _ = decimal_interval(text)
return bool(ball > lower and ball < upper)
def bounded(ball, minbits=400):
assert ball.is_finite(), 'nonfinite enclosure'
assert ball.accuracy() >= minbits, ('wide enclosure', ball.accuracy())
def spectral_controls(bits):
"""Independent convergent products from the separated square/triangle spectra."""
R = RealBallField(bits)
# log eta(i) from its product. Bound all omitted logarithms by
# rho^(M+1)/[(1-rho)(1-rho^(M+1))].
rho = (-2*R.pi()).exp()
M = 160
sq = -R.pi()/12
for j in range(1, M+1):
sq += (1-rho**j).log()
tail = rho**(M+1)/((1-rho)*(1-rho**(M+1)))
sq += R(0, tail.upper())
sq -= QQ(3)/4*R(2).log()
# Corollary 3.3, Mardby--Rowlett, with side length sqrt(3).
# sigma_1(j)/j <= j. The absolute tail after M is bounded by
# rho^(M+1)*((M+1)-M*rho)/(1-rho)^2.
rho = (-R.pi()*R(3).sqrt()).exp()
tri_sum = R(0)
for j in range(1, M+1):
sigma = sum(d for d in range(1, j+1) if j % d == 0)
tri_sum += QQ(((-1)**j)*sigma)/j*rho**j
tail = rho**(M+1)*((M+1)-M*rho)/(1-rho)**2
tri_sum += R(0, tail.upper())
tri = -(QQ(2)/3*(3*R(3).sqrt()/2).log()
+ R.pi()*R(3).sqrt()/36 + QQ(2)/3*tri_sum)
# Neumann adds the two positive axis sequences to each Dirichlet spectrum.
# Differentiating 2*k^(-s)*zeta(2s) at zero gives log(4*pi^2/k).
# Square side sqrt(2): k=pi^2/2. Triangle side sqrt(3): k=16*pi^2/27.
return {3:tri + (R(27)/4).log(), 4:sq + R(8).log()}
@lru_cache(maxsize=None)
def mp_coefficients(dps):
with mp.workdps(dps+40):
# B_(2k)/(2k)!; direct zeta avoids a second recurrence implementation.
return tuple(2*(-1)**(k+1)*mp.zeta(2*k)/(2*mp.pi)**(2*k)
for k in range(2, 112))
def integral_corner(a, dps):
with mp.workdps(dps+40):
a=mp.mpf(a[0])/a[1]
delta=mp.mpf(1)/8
coeffs=[b*(a**(1-2*k)-a) for k,b in enumerate(mp_coefficients(dps),2)]
def integrand(x):
if x==0:
return mp.mpf(0)
if x <= delta:
bracket=mp.polyval(list(reversed(coeffs)),x*x)*x*x
else:
bracket=(mp.coth(x/(2*a))-a*mp.coth(x/2))/(2*x)-(1/a-a)/12
return bracket/mp.expm1(x)
T=(dps+20)*mp.log(10)
# Known nonzero control runs before trusting quadrature.
control=mp.quad(lambda x: mp.exp(-x), [0,1,4,16,T])
assert abs(control-(1-mp.exp(-T))) < mp.mpf(10)**(-dps)
J=mp.quad(integrand,[0,delta,1,4,16,64,T])
F=(1/a-a)/12*(mp.euler-mp.log(2))-(1/a+3+a)/12*mp.log(a)+J
# Explicit bound for the excluded tail (absolute value).
tail=mp.exp(-T)/(1-mp.exp(-T))*((1+a)*(1+2/mp.expm1(T))/(2*T)+(1/a-a)/12)
return +F,+tail
def mp_polygon(n, dps):
with mp.workdps(dps+40):
a=(n-2,n)
F,tail=integral_corner(a,dps)
h=mp.mpf(1)/n
logR=mp.loggamma(1+h)+mp.loggamma(1-2*h)-mp.loggamma(1-h)
c=mp.mpf(n-1)/(6*(n-2))
L=2*c*logR+mp.log(n)/(3*(n-2))-n*F
area=mp.mpf(n)/2*mp.sin(2*mp.pi/n)
L+=mp.log(area)-mp.mpf(n)/2*mp.log(1-2*h)
return +L,+mp.exp(L),+tail
def main():
start=time.monotonic()
gen=generate.PolygonNeumannDeterminants()
R=RealBallField(750)
controls=spectral_controls(750)
for n,b in controls.items():
bounded(b,680)
l=generate.polygon_logdet(n,750)
assert l.overlaps(b),(n,l,b)
assert not (l+R('0.000001')).overlaps(b)
for n, closed in [(3, R(QQ(27)/4) * R(QQ(1)/3).gamma() / (R.pi()**(QQ(2)/3)*R(3)**(QQ(11)/12))),
(4, R(2)**(QQ(5)/4)*R(QQ(1)/4).gamma()/R.pi()**(QQ(3)/4))]:
bounded(closed,680)
assert generate.polygon_logdet(n,750).exp().overlaps(closed)
print('Separated-spectrum controls and deliberate perturbations passed.',flush=True)
# Integral control: the angle pi has F(1)=0, but also test nonzero F(1/2).
f1,tail=integral_corner((1,1),150)
assert abs(f1)<mp.mpf('1e-140')
fhalf,tail=integral_corner((1,2),150)
with mp.workdps(190):
known=mp.log(2)/12+mp.log(2*mp.pi)/8
assert abs(fhalf-known)<mp.mpf('1e-140')
print('Corner integral controls passed.',flush=True)
stored=None
if '--stored' in sys.argv or os.environ.get('NUMBERDB_CHECK_STORED')=='1':
client=generate.authenticated_client()
companion=numberdb.table('T620',client=client)['Numbers']
stored=numberdb.table('T621',client=client)
records=stored['Numbers']
if isinstance(records,list):
stored={ (str(r['params']['n']),r['params']['convention']):r['number'] for r in records }
else:
stored={(str(n),c):(v['number'] if isinstance(v,dict) else v)
for n,group in records.items() for c,v in group.items()}
assert len(stored)==56
evidence=[]
for n in range(3,31):
lo=generate.polygon_logdet(n,numberdb.bits(100,losing=generate.WORKING_GUARD))
hi=generate.polygon_logdet(n,750)
bounded(lo);bounded(hi,680);assert lo.overlaps(hi)
vals=[]
for convention,ball,high in [('det',lo.exp(),hi.exp()),('logdet',lo,hi)]:
text=to_text(ball,100)
assert contains_ball(text,high),(n,convention,text,high)
_,_,ulp=decimal_interval(text)
assert not contains_ball(text,high+10*ulp)
if stored is not None:
text=stored[str(n),convention]
assert contains_ball(text,high),(n,convention,text,high)
vals.append((convention,text,int(ball.accuracy())))
L150,D150,t150=mp_polygon(n,150)
L200,D200,t200=mp_polygon(n,200)
with mp.workdps(240):
err=max(abs(L150-L200),abs(D150-D200))
assert err < mp.mpf('1e-140'),(n,mp.nstr(err,8))
for convention,text,_ in vals:
expected=D200 if convention=='det' else L200
lower,upper,ulp=decimal_interval(text)
low=mp.mpf(str(lower.numerator()))/mp.mpf(str(lower.denominator()))
high=mp.mpf(str(upper.numerator()))/mp.mpf(str(upper.denominator()))
assert low<expected<high,(n,convention)
row={'n':n,'values':[{'convention':c,'number':t,'accuracy_bits':b} for c,t,b in vals],
'independent_logdet_150':mp.nstr(L150,155),'independent_logdet_200':mp.nstr(L200,205),
'agreement_error':mp.nstr(err,8),'integral_tail_bound_200':mp.nstr(n*t200,8)}
evidence.append(row)
print('ROW '+json.dumps(row),flush=True)
if stored is not None:
for n in range(3,31):
dlow,dhigh,_=decimal_interval(stored[str(n),'det'])
llow,lhigh,_=decimal_interval(stored[str(n),'logdet'])
db=R((dlow+dhigh)/2,(dhigh-dlow)/2)
lb=R((llow+lhigh)/2,(lhigh-llow)/2)
assert db.log().overlaps(lb) and lb.exp().overlaps(db)
if n in controls:
assert lb.overlaps(controls[n])
c=companion[str(n)]
cdlo,cdhi,_=decimal_interval(c['det'])
cllo,clhi,_=decimal_interval(c['logdet'])
cdb=R((cdlo+cdhi)/2,(cdhi-cdlo)/2)
clb=R((cllo+clhi)/2,(clhi-cllo)/2)
alpha=QQ(n-2)/n
area=QQ(n)/2*(2*R.pi()/n).sin()
logr=area.log()-QQ(n)/2*R(alpha).log()
assert (db/cdb).overlaps(logr.exp())
assert (lb-clb).overlaps(logr)
# Algebraic corner heat coefficient and omitted zero mode.
heat=QQ(n)/24*(1/alpha-alpha)
assert heat==QQ(n-1)/(6*(n-2))
assert (heat-1)-heat == -1
print('SUMMARY '+json.dumps({'table':'T621','entries':56,'stored_checked':stored is not None,
'production_bits':numberdb.bits(100,losing=generate.WORKING_GUARD),'comparison_bits':750,
'minimum_accuracy_bits':min(v['accuracy_bits'] for r in evidence for v in r['values']),
'sage':SAGE_VERSION,'numberdb':numberdb.__version__,'mpmath':mp.__version__,'seconds':time.monotonic()-start,
'skill_manifest':'../20261030T080603Z-build-skills/manifest.json',
'skill_sha256':'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199'}),flush=True)
if __name__=='__main__':
main()