back to table · edit · history · where entries came from · files · download
9344 bytes, as of the version from 2026-10-25 14:25 (current). Recorded here, not run.
"""Independent checks for Steklov eigenvalues -- numberdb.org/T489.
Run: sage -python checks.py
For a private draft set NUMBERDB_API_KEY, then:
sage -python checks.py --stored
The generator must be beside this file. In stored mode every decimal is read
from the API, parsed as its exact last-place interval and checked against an
independent 180-digit computation. No write or review action is performed.
Controls: disk u=x has Rayleigh quotient 1; on the unit square,
u=(x-1/2)(y-1/2) has quotient 2. Positive rectangle modes are recomputed using
mpmath's trigonometric Robin quotients at 140 and 180 digits, independently
from the generator's pole-free Arb determinants. A complete separate branch
enumeration is sorted before comparison. Energy integrals use analytic
antiderivatives of squared factors and derivatives. Annulus eigenvalues are
independently obtained from the weighted boundary Dirichlet-to-Neumann matrix,
then tested on r^m + B r^-m with the inner normal -d/dr.
The square separation parameters are also checked against the beam roots in
T486. Its saved 100-digit entries are compared at their actual last-place
intervals; no assumption about equality of rounded midpoints is made.
"""
import json, os, sys
from decimal import Decimal
import mpmath as mp
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 generate as gen
def decimal_rational(s):
d=Decimal(str(s));t=d.as_tuple()
n=int(''.join(map(str,t.digits)) or '0')*(-1 if t.sign else 1)
return QQ(n)*QQ(10)**t.exponent
def written_interval(s):
s=str(s)
if '.' not in s and 'e' not in s.lower():
q=QQ(s);return q,q
q=decimal_rational(s);ulp=QQ(10)**Decimal(s).as_tuple().exponent
return q-ulp,q+ulp
def inside_written(s,x):
lo,hi=written_interval(s)
q=decimal_rational(mp.nstr(x,175))
assert lo<=q<=hi,(s,mp.nstr(x,110))
def mp_roots(a,b,dps):
with mp.workdps(dps):
out=[(mp.mpf(0),None)]
if a==b:out.append((1/a,None))
for ta,hb in [(a,b),(b,a)]:
for trig in ('sin','cos'):
for hyp in ('cosh','sinh'):
for k in range(6):
if trig=='sin' and hyp=='sinh' and k==0 and hb<=ta:continue
left=k*mp.pi/ta if trig=='sin' else (k+mp.mpf('.5'))*mp.pi/ta
right=left+mp.pi/(2*ta)
def f(x):
t=x*mp.cot(ta*x) if trig=='sin' else -x*mp.tan(ta*x)
h=x*mp.tanh(hb*x) if hyp=='cosh' else x*mp.coth(hb*x)
return t-h
# These quotients decrease, with opposite endpoint signs.
lo=left+mp.mpf('1e-6');hi=right-mp.mpf('1e-6')
assert f(lo)>0 and f(hi)<0
# mpmath Anderson/bisection is independent of the Arb code.
alpha=mp.findroot(f,(lo,hi),solver='anderson',tol=mp.mpf(10)**(-dps+15),maxsteps=1200)
assert left<alpha<right
sigma=alpha*(mp.tanh(hb*alpha) if hyp=='cosh' else mp.coth(hb*alpha))
out.append((+sigma,(+alpha,ta,hb,trig,hyp,k)))
out.sort(key=lambda row:row[0])
return out[:20]
def squares_integrals(alpha,length,kind):
# Integral of F^2 and (F')^2 on [-length,length].
x=alpha*length
if kind=='cos':return length+mp.sin(2*x)/(2*alpha),alpha**2*(length-mp.sin(2*x)/(2*alpha))
if kind=='sin':return length-mp.sin(2*x)/(2*alpha),alpha**2*(length+mp.sin(2*x)/(2*alpha))
if kind=='cosh':return length+mp.sinh(2*x)/(2*alpha),alpha**2*(-length+mp.sinh(2*x)/(2*alpha))
return -length+mp.sinh(2*x)/(2*alpha),alpha**2*(length+mp.sinh(2*x)/(2*alpha))
def quotient(alpha,a,b,trig,hyp):
t,dt=squares_integrals(alpha,a,trig);h,dh=squares_integrals(alpha,b,hyp)
tv=getattr(mp,trig)(alpha*a);hv=getattr(mp,hyp)(alpha*b)
return (dt*h+t*dh)/(2*tv**2*h+2*hv**2*t)
def annulus_independent(dps):
with mp.workdps(dps):
result=[(mp.mpf(0),None),(mp.mpf(3)/(2*mp.log(2)),('radial',))]
for m in range(1,25):
# Boundary values at radii 1 and 2 give the Dirichlet-to-Neumann
# matrix [[m*coth(mL), -m*csch(mL)],
# [-m*csch(mL)/2, m*coth(mL)/2]].
# Multiplying boundary values by sqrt(radius) symmetrizes it.
l=m*mp.log(2)
matrix=mp.matrix([[m*mp.coth(l),-m*mp.csch(l)/mp.sqrt(2)],
[-m*mp.csch(l)/mp.sqrt(2),m*mp.coth(l)/2]])
vals=mp.eigsy(matrix,eigvals_only=True)
for j in range(2):result.extend([(+vals[j],(m,j))]*2)
result.sort(key=lambda row:row[0]);return result[:20]
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
stored='--stored' in sys.argv or os.environ.get('NUMBERDB_CHECK_STORED')=='1'
# The wrapper forwards only selected flags; a tiny runner may call main.
table=numberdb.table(gen.TABLE) if stored else None
mp.mp.dps=180
# Known nonzero controls are performed before the method is trusted.
disk_control=mp.quad(lambda t:mp.cos(t)**2,[0,mp.pi,2*mp.pi])
assert abs(mp.pi/disk_control-1)<mp.mpf('1e-170')
a=mp.mpf('.5');area_energy=2*(2*a)*(2*a**3/3)
boundary_mass=4*a*a*(2*a**3/3)
assert area_energy/boundary_mass==2
print('Controls: disk u=x gives 1; square bilinear mode gives 2.')
independently={}
for domain,a,b in [('unit-square',mp.mpf('.5'),mp.mpf('.5')),('rectangle-1x2',mp.mpf('.5'),mp.mpf(1))]:
low=mp_roots(a,b,140);high=mp_roots(a,b,180)
for (x,_),(y,mode) in zip(low,high):
assert abs(x-y)<mp.mpf('1e-120')
if mode:
alpha,ta,hb,trig,hyp,k=mode
assert abs(quotient(alpha,ta,hb,trig,hyp)-y)<mp.mpf('1e-155')
independently[domain]=high
for domain,a,b in [('unit-square',mp.mpf(1),mp.mpf(1)),('rectangle-1x2',mp.mpf(1),mp.mpf(2))]:
scaled=mp_roots(a,b,160)
for (x,_),(y,_) in zip(independently[domain],scaled):
assert abs(x-2*y)<mp.mpf('1e-140')
print('Scaling independently checked on all 40 rectangle modes.')
alpha=mp.findroot(lambda x:mp.cos(x)*mp.cosh(x)+1,(1,2))
assert abs(alpha*mp.tanh(alpha/2)-independently['unit-square'][1][0])<mp.mpf('1e-160')
print('Displayed mpmath snippet checked.')
independently['unit-disk']=[(mp.mpf(n//2),None) for n in range(1,21)]
low=annulus_independent(140);high=annulus_independent(180)
for (x,_),(y,mode) in zip(low,high):
assert abs(x-y)<mp.mpf('1e-120')
if mode and mode[0]!='radial':
m,j=mode
beta=(m+y)/(m-y)
# Inner boundary: -u'(1)=sigma*u(1).
assert abs((-m+m*beta)-y*(1+beta))<mp.mpf('1e-155')
# Outer boundary: +u'(2)=sigma*u(2).
residual=m*(mp.mpf(2)**(m-1)-beta*mp.mpf(2)**(-m-1))-y*(mp.mpf(2)**m+beta*mp.mpf(2)**(-m))
assert abs(residual)<mp.mpf('1e-145')
radial=mp.mpf(3)/(2*mp.log(2))
# u=log(r)-1/sigma meets both radial boundary conditions.
assert abs(mp.mpf('.5')-radial*(mp.log(2)-1/radial))<mp.mpf('1e-170')
independently['annulus-1-2']=high
minimum_digits=10000;count=0;written={}
for domain in gen.DOMAINS:
rows=gen.spectrum(domain)
for n,(row,(x,_)) in enumerate(zip(rows,independently[domain]),1):
if stored:
r=table['Numbers'][domain][str(n)]
text=r['number'] if isinstance(r,dict) else r
else:text=to_text(row['number'],100)
inside_written(text,x);count+=1;written[domain+','+str(n)]=str(text)
if hasattr(row['number'],'rad') and row['number'].rad():
v=row['number'];assert v.is_finite()
rel=QQ(v.rad().exact_rational())/abs(QQ(v.mid().exact_rational()))
digits=-mp.log10(mp.mpf(str(rel.numerator()))/int(rel.denominator()))
minimum_digits=min(minimum_digits,float(digits))
assert digits>110
print(domain,'first values:',[mp.nstr(v,15) for v,_ in independently[domain][:6]])
# External family comparison: read the actual beam table, not its generator.
beams=numberdb.table('T486')['Numbers']
compared=0
for row in gen.spectrum('unit-square'):
if 'alpha' not in row:continue
alpha=row['alpha']; found=0
for boundary in ['clamped-free','clamped-clamped']:
for r in beams[boundary].values():
s=r['number'] if isinstance(r,dict) else r
lo,hi=written_interval(s)
B=alpha.parent()
if alpha.overlaps(B(lo).union(B(hi))):found+=1
assert found==1
compared+=1
print('Independent decimal interval checks:',count)
print('Square separation parameters matched to T486:',compared)
print('Minimum relative enclosure digits:',round(minimum_digits,3))
print('Stored mode:',stored)
print('All controls, precision comparisons, energy and boundary checks passed.')
print('WRITTEN_JSON',json.dumps(written,sort_keys=True))
if __name__=='__main__':main()