back to table · edit · history · where entries came from · files · download
8479 bytes, as of the version from 2026-10-28 04:41. Recorded here, not run.
"""Independent checks for numberdb.org/T571; run: sage -python checks.py.
Use --stored to repeat the checks on values read from NumberDB (a private draft
requires NUMBERDB_API_KEY). No mutation is performed. This checker evaluates
Laguerre integrals by exact moments: a determinant for beta=2, or a de Bruijn
Pfaffian for beta=1. It does not use the generator's gamma product.
For p,q >= 0, I(p,q)=integral_{0<x<y} x^p y^q exp(-x-y) dx dy
= q! sum_{k=0}^q (p+k)!/(k! 2^(p+k+1)). The skew moment is I(p,q)-I(q,p).
Odd ranks are bordered by the single moments p!. See de Bruijn (1955),
On some multiple integrals involving determinants, (1.1), (4.3)-(4.5).
The trace constraint divides the Laguerre integral by (A-1)!, where
A=r*(1+beta*(N-r))+beta*r*(r-1)/2. The orbit volume is independently
assembled from successive real or complex projective-space volumes.
Controls run first: exact determinant/Pfaffian examples and eigenvalue
integrals; all four N=2 volumes come from a disk/ball of radius 1/sqrt(2).
The checks also cover pure-state geometry and the boundary/interior ratio.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from math import factorial
from functools import lru_cache
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import json
import os
import sys
import time
from check import exactness, measure, names_its_rings
from numberdb._write import to_text
from generate import HilbertSchmidtVolumes
def determinant(matrix):
a=[list(map(QQ,row)) for row in matrix]
result=QQ(1)
for k in range(len(a)):
pivot=next((i for i in range(k,len(a)) if a[i][k]),None)
if pivot is None: return QQ(0)
if pivot!=k: a[k],a[pivot]=a[pivot],a[k]; result=-result
t=a[k][k];result*=t
for i in range(k+1,len(a)):
ratio=a[i][k]/t
for j in range(k+1,len(a)):a[i][j]-=ratio*a[k][j]
return result
def pfaffian(matrix):
a=[list(map(QQ,row)) for row in matrix]
result=QQ(1)
while a:
assert len(a)%2==0
pivot=next((j for j in range(1,len(a)) if a[0][j]),None)
if pivot is None:return QQ(0)
if pivot!=1:
a[1],a[pivot]=a[pivot],a[1]
for row in a:row[1],row[pivot]=row[pivot],row[1]
result=-result
t=a[0][1];result*=t
a=[[a[i][j]-(a[0][i]*a[1][j]-a[0][j]*a[1][i])/t
for j in range(2,len(a))] for i in range(2,len(a))]
return result
@lru_cache(None)
def ordered_moment(p,q):
return QQ(factorial(q))*sum((QQ(factorial(p+k))/QQ(factorial(k)*2**(p+k+1))
for k in range(q+1)),QQ(0))
@lru_cache(None)
def laguerre_integral(beta,d,r):
if beta==2:
a=[[factorial(2*d+i+j) for j in range(r)] for i in range(r)]
return QQ(factorial(r))*determinant(a)
size=r+(r%2)
a=[[QQ(0) for j in range(size)] for i in range(size)]
for i in range(r):
for j in range(i+1,r):
a[i][j]=ordered_moment(d+i,d+j)-ordered_moment(d+j,d+i)
a[j][i]=-a[i][j]
if r%2:a[i][r]=QQ(factorial(d+i));a[r][i]=-a[i][r]
return QQ(factorial(r))*pfaffian(a)
def independent_volume(field,N,r,bits=560):
R=RealBallField(bits);beta=1 if field=='real' else 2;d=N-r
A=r*(1+beta*d)+beta*r*(r-1)//2
integral=laguerre_integral(beta,d,r)/QQ(factorial(A-1))
orbit=R(1)
for k in range(d+1,N+1):
if beta==2:orbit*=(2*R.pi())**(k-1)/R(factorial(k-1))
else:orbit*=R(2)**(QQ(k-1)/2)*R.pi()**(QQ(k)/2)/R(QQ(k)/2).gamma()
return R(r).sqrt()*R(integral)*orbit/R(factorial(r))
def endpoints(text):
# Decimal text asserts a radius of one unit in its last place.
x=Decimal(str(text).rstrip('?'));f=Fraction(x)
c=QQ(f.numerator)/QQ(f.denominator)
exponent=x.as_tuple().exponent
ulp=QQ(10)**exponent
return c-ulp,c+ulp
def narrow(value,bits=330):
assert value.is_finite() and value>0
assert value.accuracy()>=bits,(value,value.accuracy())
def contains_ball(text,ball):
lo,hi=endpoints(text)
assert ball.lower()>=lo and ball.upper()<=hi,(text,str(ball))
def run(stored=False):
start=time.monotonic()
# Nontrivial, positive, negative, singular and swapped-pivot controls.
assert determinant([[1,2],[3,4]])==-2
assert determinant([[0,1],[2,3]])==-2
assert determinant([[1,2],[2,4]])==0
a=[[0,2,3,5],[-2,0,7,11],[-3,-7,0,13],[-5,-11,-13,0]]
assert pfaffian(a)==28 and determinant(a)==28**2
assert pfaffian([[0,0,1,0],[0,0,0,1],[-1,0,0,0],[0,-1,0,0]])==-1
assert laguerre_integral(1,0,1)==1
assert laguerre_integral(1,0,2)==1
assert laguerre_integral(2,0,2)==2
print('PASS controls: determinant, Pfaffian, and small eigenvalue integrals')
g=HilbertSchmidtVolumes();params=list(g.enumerate())
assert len(params)==270 and len({tuple(p.values()) for p in params})==270
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
values={tuple(p.values()):g.value(p,g.digits) for p in params}
assert not exactness(values)
if stored:
doc=numberdb.table(g.table)
rows=doc['Numbers']
if isinstance(rows,dict):
texts={(field,int(N),int(r)):value for field,sizes in rows.items()
for N,ranks in sizes.items() for r,value in ranks.items()}
else:
texts={(str(x['params']['field']),int(x['params']['N']),int(x['params']['rank'])):x['number'] for x in rows}
assert len(texts)==270
assert list(texts)==list(values)
else:texts={key:to_text(value,digits=100) for key,value in values.items()}
accuracy=[]
for key,value in values.items():
narrow(value);accuracy.append(int(value.accuracy()))
oracle=independent_volume(*key)
narrow(oracle,480)
assert value.overlaps(oracle),key
contains_ball(texts[key],value)
contains_ball(texts[key],oracle)
print('PASS all 270 exact-moment computations and serialized-interval containments')
if stored:
namespace={}
exec(doc['Programs']['program-sage']['code'],namespace)
for key in values:
snippet=namespace['hs_volume'](*key)
narrow(snippet)
assert snippet.overlaps(independent_volume(*key)),key
print('PASS stored Programs snippet on all 270 arguments')
R=RealBallField(560)
controls={('real',2,1):R(2).sqrt()*R.pi(),('real',2,2):R.pi()/2,
('complex',2,1):2*R.pi(),('complex',2,2):R.pi()*R(2).sqrt()/3}
for key,v in controls.items():contains_ball(texts[key],v)
for field in ('real','complex'):
for N in range(2,17):
if field=='complex':pure=(2*R.pi())**(N-1)/R(factorial(N-1))
else:pure=R(2)**(QQ(N-1)/2)*R.pi()**(QQ(N)/2)/R(QQ(N)/2).gamma()
contains_ball(texts[(field,N,1)],pure)
D=N*N-1 if field=='complex' else N*(N+1)//2-1
ratio=R(D)*R(N*(N-1)).sqrt()
# Interval arithmetic applied to database decimal enclosures.
plo,phi=endpoints(texts[(field,N,N-1)])
vlo,vhi=endpoints(texts[(field,N,N)])
assert plo/vhi<=ratio.lower() and ratio.upper()<=phi/vlo
if field=='complex':
interior=R(N).sqrt()*(2*R.pi())**(N*(N-1)//2)
interior*=R(QQ(__import__('functools').reduce(lambda x,k:x*factorial(k),range(N),1))/QQ(factorial(N*N-1)))
contains_ball(texts[(field,N,N)],interior)
print('PASS four Bloch disk/ball controls, 30 pure-state volumes, 30 boundary ratios, 15 complex full-rank formulas')
report={'table':g.table,'stored_readback':stored,'entries':len(values),'min_accuracy_bits':min(accuracy),
'working_bits':numberdb.bits(100,losing=64),'independent_bits':560,
'serialized_max_characters':max(map(len,texts.values())),
'measure':measure(values),'seconds':round(time.monotonic()-start,3),
'skill_manifest':'agents/runs/20261028T042605Z-build-skills/manifest.json',
'skill_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c'}
print('EVIDENCE '+json.dumps(report,sort_keys=True))
for key in controls:print('CONTROL',key,texts[key])
return report
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
run('--stored' in sys.argv)