back to table · edit · history · where entries came from · files · download
10921 bytes, as of the version from 2026-10-27 23:22 (current). Recorded here, not run.
"""Independent checks of T564 using eigenvalue integrals, not Selberg products.
Run: sage -python checks.py
Set NUMBERDB_CHECK_STORED=1 to recheck all database values independently.
NUMBERDB_KEY_FROM_STDIN=1 reads the draft's API key from stdin.
Andreief's determinant identity gives beta=2 as n! times a moment determinant.
The confluent Vandermonde and de Bruijn identities give beta=4 and beta=1 as
n! times Pfaffians. All matrix operations are written out over QQ. For beta=1,
the substitution (x,y)=(st,s(1-t)) reduces skew moments to incomplete beta
integrals on [0,1/2]. At integer parameters these are binomial sums. At
half-integer parameters the incomplete-beta recurrences start at pi/2; the
pi coefficient cancels from every skew moment exactly. An odd-size bordering
moment retains sqrt(pi) exactly when r is even. This shares neither the gamma
product nor gamma-reduction code with generate.py.
"""
import os
import sys
import json
import hashlib
import platform
from math import factorial, comb
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from generate import LaguerreSelberg, reduced_integral
from check import exactness, measure, agrees_with, names_its_rings
from numberdb._write import to_text
def determinant(a):
a = [row[:] for row in a]
ans = 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]
ans = -ans
p = a[k][k]
ans *= p
for i in range(k+1, len(a)):
f = a[i][k] / p
for j in range(k+1, len(a)):
a[i][j] -= f*a[k][j]
return ans
def pfaffian(a):
a = [row[:] for row in a]
n = len(a)
assert n % 2 == 0
ans = QQ(1)
for k in range(0, n, 2):
pivot = next((j for j in range(k+1, n) if a[k][j]), None)
if pivot is None:
return QQ(0)
if pivot != k+1:
a[k+1], a[pivot] = a[pivot], a[k+1]
for row in a:
row[k+1], row[pivot] = row[pivot], row[k+1]
ans = -ans
p = a[k][k+1]
ans *= p
for i in range(k+2, n):
for j in range(i+1, n):
a[i][j] += (a[k+1][i]*a[k][j] - a[k][i]*a[k+1][j])/p
a[j][i] = -a[i][j]
return ans
def half_moment(k):
"""Gamma(k+1/2)/sqrt(pi), from the integration-by-parts recurrence."""
ans = QQ(1)
for j in range(k):
ans *= QQ(2*j+1)/2
return ans
def incomplete_half(p_int, q_int):
"""B_(1/2)(p_int+1/2,q_int+1/2) as constant + coefficient*pi."""
const, coeff = QQ(0), QQ(1)/2
p, q = QQ(1)/2, QQ(1)/2
for _ in range(p_int):
const = (p*const - QQ(1)/ZZ(2)**int(p+q))/(p+q)
coeff = p*coeff/(p+q)
p += 1
for _ in range(q_int):
const = (q*const + QQ(1)/ZZ(2)**int(p+q))/(p+q)
coeff = q*coeff/(p+q)
q += 1
return const, coeff
def skew_moment(r, i, j):
if r % 2:
p, q = (r+1)//2+i, (r+1)//2+j
regularized = QQ(sum(comb(p+q-1,k) for k in range(p,p+q)))/ZZ(2)**(p+q-1)
return QQ(factorial(p-1)*factorial(q-1))*(2*regularized-1)
p, q = r//2+i, r//2+j
const, coeff = incomplete_half(p,q)
gamma_sum = ZZ(factorial(p+q))
assert 2*gamma_sum*coeff == half_moment(p)*half_moment(q)
return 2*gamma_sum*const
def independent_integral(key):
beta, n, r = key
if beta == 2:
a = r+1
matrix = [[QQ(factorial(a+i+j-1)) for j in range(n)] for i in range(n)]
return ZZ(factorial(n))*determinant(matrix), 0
if beta == 4:
a = 2*(r+1)
matrix = [[QQ((j-i)*factorial(a+i+j-2)) for j in range(2*n)] for i in range(2*n)]
return ZZ(factorial(n))*pfaffian(matrix), 0
size = n + n%2
matrix = [[QQ(0) for _ in range(size)] for _ in range(size)]
for i in range(n):
for j in range(i+1,n):
matrix[i][j] = skew_moment(r,i,j)
matrix[j][i] = -matrix[i][j]
exponent = int(n % 2 and r % 2 == 0)
if n % 2:
for i in range(n):
m = half_moment(r//2+i) if r%2 == 0 else QQ(factorial((r+1)//2+i-1))
matrix[i][n] = m
matrix[n][i] = -m
return ZZ(factorial(n))*pfaffian(matrix), exponent
def ball_from_pair(pair, bits=640):
q, e = pair
field = RealBallField(bits)
return field(q)*field.pi().sqrt()**e
def claimed_interval(text):
"""Exact endpoints for the database's decimal +/- one last place rule."""
s = str(text).lower()
if '.' not in s and 'e' not in s:
q = QQ(s)
return q,q
mantissa, exponent = (s.split('e') + ['0'])[:2] if 'e' in s else (s,'0')
decimals = len(mantissa.split('.')[1]) if '.' in mantissa else 0
f = Fraction(s)
center = QQ(f.numerator)/f.denominator
unit = QQ(10)**(int(exponent)-decimals)
return center-unit, center+unit
def relative_bits(b):
radius = b.rad().exact_rational()
midpoint = abs(b.mid().exact_rational())
assert radius > 0 and midpoint > radius
k = int((midpoint/radius).numerator().nbits() - (midpoint/radius).denominator().nbits())
while radius*ZZ(2)**k > midpoint:
k -= 1
return k
def assert_contained(text, pair):
lo, hi = claimed_interval(text)
if pair[1] == 0:
assert lo == hi == pair[0], (text,pair)
else:
b = ball_from_pair(pair)
assert b.is_finite() and relative_bits(b) > 580
field = b.parent()
assert b > field(lo) and b < field(hi), text
def controls():
assert determinant([[QQ(2),QQ(1)],[QQ(3),QQ(4)]]) == 5
assert pfaffian([[QQ(0),QQ(2),QQ(3),QQ(5)],[-QQ(2),QQ(0),QQ(7),QQ(11)],[-QQ(3),-QQ(7),QQ(0),QQ(13)],[-QQ(5),-QQ(11),-QQ(13),QQ(0)]]) == 28
# Integrate e^-x-y (x-y)^2 by expanding into three separable moments.
assert independent_integral((2,2,0)) == (2,0)
# Integrate e^-x-y |x-y| using two ordered triangular domains.
assert independent_integral((1,2,1)) == (1,0)
# x=t^2 transforms this scalar integral to twice the Gaussian half integral.
assert independent_integral((1,1,0)) == (1,1)
assert independent_integral((4,1,0)) == (1,0)
assert independent_integral((4,2,0)) == (72,0)
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
controls()
gen = LaguerreSelberg()
keys = [tuple(int(p[k]) for k in gen.parameters) for p in gen.enumerate()]
actual = {key: reduced_integral(*key) for key in keys}
other = {key: independent_integral(key) for key in keys}
assert not agrees_with(actual, other.__getitem__)
assert not exactness(actual)
assert not names_its_rings(os.path.join(os.path.dirname(__file__), 'generate.py'))
scalar = recurrence = decimal_count = 0
raw, retained = {}, []
for key in keys:
beta,n,r = key
q,e = actual[key]
assert e == int(beta==1 and n%2==1 and r%2==0)
if n == 1:
assert actual[key] == other[key]
scalar += 1
# All three classes: increasing r by 2 raises a by beta.
factor = QQ(1)
a = QQ(beta*(r+1))/2
c = QQ(beta)/2
for j in range(n):
for k in range(beta):
factor *= a+j*c+k
q2,e2 = independent_integral((beta,n,r+2))
assert (q*factor,e) == (q2,e2)
recurrence += 1
params = dict(zip(gen.parameters,map(str,key)))
val = gen.value(params,100)
written = to_text(val, digits=100)
raw[','.join(map(str,key))] = written
assert_contained(written,other[key])
if e:
assert val.is_finite()
retained.append(int(relative_bits(val)))
decimal_count += 1
stored_count = 0
stored_recurrence = 0
stored_recurrence_pairs = 0
if os.environ.get('NUMBERDB_CHECK_STORED') == '1':
from numberdb._generate import _stored_entries
document = numberdb.table('T564')
entries = _stored_entries(document)
assert set(entries) == set(raw)
for identity,written in entries.items():
key = tuple(map(int,identity.split(',')))
assert_contained(written, other[key])
stored_count += 1
beta,n,r = key
factor = QQ(1)
for j in range(n):
for k in range(beta):
factor *= QQ(beta*(r+1+j))/2+k
lo,hi = claimed_interval(written)
target = ','.join(map(str,(beta,n,r+2)))
if target in entries:
tlo,thi = claimed_interval(entries[target])
assert lo*factor <= thi and tlo <= hi*factor
stored_recurrence_pairs += 1
else:
q2,e2 = independent_integral((beta,n,r+2))
if e2 == 0:
assert lo*factor == hi*factor == q2
else:
b = ball_from_pair((q2,e2))
assert b > b.parent()(lo*factor) and b < b.parent()(hi*factor)
stored_recurrence += 1
report = gen.verify(sample=None)
assert report.ok, str(report)
print('VERIFY',str(report))
evidence = {
'table':'T564','skill_manifest':'agents/runs/20261027T230956Z-build-skills/manifest.json',
'skill_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c',
'entries':len(keys),'exact_entries':len(keys)-decimal_count,'ball_entries':decimal_count,
'independent_moment_checks':len(keys),'r_plus_two_recurrence_checks':recurrence,
'scalar_checks':scalar,'stored_independent_checks':stored_count,
'stored_recurrence_checks':stored_recurrence,
'stored_recurrence_pairs':stored_recurrence_pairs,
'serialized_interval_containment_checks':len(keys),
'working_bits':numberdb.bits(100,losing=96),
'minimum_retained_relative_bits':min(retained),
'independent_ball_bits':640,'measurement':measure(raw),
'controls':'determinant 5; Pfaffian 28; Z(2,2,0)=2; Z(1,2,1)=1; Z(1,1,0)=sqrt(pi); Z(4,1,0)=1; Z(4,2,0)=72',
'python':platform.python_version(),
'source_sha256':{name:hashlib.sha256(open(os.path.join(os.path.dirname(__file__),name),'rb').read()).hexdigest() for name in ('generate.py','checks.py','check.py')},
'methods':{'beta=1':'exact skew-moment Pfaffian with incomplete beta at 1/2',
'beta=2':'exact Hankel moment determinant',
'beta=4':'exact confluent-Vandermonde Pfaffian'},
}
print('EVIDENCE '+json.dumps(evidence,sort_keys=True))
print('VALUES '+json.dumps(raw,sort_keys=False))
if __name__ == '__main__':
main()