back to table · edit · history · where entries came from · files · download
8921 bytes, as of the version from 2026-10-28 00:24. Recorded here, not run.
"""Independent checks of Mehta integrals -- numberdb.org/T566.
sage -python checks.py # local enumeration, no API or key
sage -python checks.py --stored # every entry read from the private draft
For --stored, set NUMBERDB_API_KEY or use NUMBERDB_KEY_FROM_STDIN=1.
NUMBERDB_CHECK_STORED=1 also selects stored checks. Output is JSON evidence.
The integral evaluator contains no Gamma factors or Mehta products. See
independent-checks.md for the determinant and Pfaffian derivations.
"""
import json
import os
import platform
from pathlib import Path
from functools import lru_cache
from math import factorial
from fractions import Fraction
from decimal import Decimal
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.env import SAGE_VERSION
from numberdb._write import to_text
from check import exactness, measure, names_its_rings, agrees_with
import generate
def determinant(rows):
a = [[QQ(x) for x in row] for row in rows]
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
d = a[k][k]
result *= d
for i in range(k + 1, len(a)):
ratio = a[i][k] / d
for j in range(k + 1, len(a)):
a[i][j] -= ratio * a[k][j]
return result
def pfaffian(rows):
a = [[QQ(x) for x in row] for row in rows]
size = len(a)
assert size % 2 == 0
result = QQ(1)
for k in range(0, size, 2):
pivot = next((j for j in range(k + 1, size) 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]
result = -result
d = a[k][k + 1]
result *= d
for i in range(k + 2, size):
for j in range(i + 1, size):
a[i][j] += (a[k + 1][i] * a[k][j] - a[k][i] * a[k + 1][j]) / d
a[j][i] = -a[i][j]
return result
@lru_cache(None)
def moment(k):
"""Integral x^k exp(-x^2/2) dx divided by sqrt(2*pi)."""
if k < 0 or k % 2:
return QQ(0)
result = QQ(1)
for j in range(1, k, 2):
result *= j
return result
@lru_cache(None)
def skew_moment(i, j):
"""Double sign(y-x) moment divided by sqrt(pi)."""
if (i + j) % 2 == 0:
return QQ(0)
if i % 2:
return -skew_moment(j, i)
k = (j - 1) // 2
result = QQ(0)
# Integral_x^infty y^(2k+1)exp(-y^2/2)dy
# = exp(-x^2/2) sum_l 2^(k-l) k!/l! x^(2l).
for ell in range(k + 1):
coefficient = QQ(2 ** (k - ell) * factorial(k)) / QQ(factorial(ell))
power = i + 2 * ell
result += 2 * coefficient * moment(power) / QQ(2) ** (power // 2)
return result
@lru_cache(None)
def integral_components(beta, n):
"""Independent exact integral as q*sqrt(2)^s*sqrt(pi)^p."""
if beta == 2:
d = determinant([[moment(i+j) for j in range(n)] for i in range(n)])
return QQ(factorial(n)) * d * QQ(2) ** (n//2), n % 2, n
if beta == 4:
p = pfaffian([[(j-i)*moment(i+j-1) for j in range(2*n)] for i in range(2*n)])
return QQ(factorial(n)) * p * QQ(2) ** (n//2), n % 2, n
assert beta == 1
size = n + n % 2
a = [[QQ(0) for _ in range(size)] for _ in range(size)]
for i in range(n):
for j in range(n):
a[i][j] = skew_moment(i,j)
if n % 2:
a[i][n], a[n][i] = moment(i), -moment(i)
return QQ(factorial(n)) * pfaffian(a), n % 2, (n+1)//2
def enclosure(parts, precision=640):
q,s,p = parts
b=RealBallField(precision)
# Evaluate powers of pi separately, including the odd sqrt(pi),
# independently of the generator's repeated sqrt(pi) power.
result=b(q)*b.pi()**(p//2)
if p%2: result *= b.pi().sqrt()
if s: result *= b(2).sqrt()
assert result.is_finite() and result > 0
assert result.accuracy() >= 600
return result
def decimal_bounds(text):
d=Decimal(text)
f=Fraction(d)
mid=QQ(f.numerator)/QQ(f.denominator)
exponent=d.as_tuple().exponent
ulp=QQ(10)**exponent
return mid-ulp,mid+ulp
def contained(ball, text):
lo,hi=decimal_bounds(text)
assert ball.is_finite() and ball.accuracy() >= 380
assert ball.lower().exact_rational() >= lo, text
assert ball.upper().exact_rational() <= hi, text
def controls():
assert determinant([[1,2,3],[0,1,4],[5,6,0]]) == 1
assert determinant([[0,2],[3,4]]) == -6
assert determinant([[1,1],[1,1]]) == 0
assert pfaffian([[0,1,2,3],[-1,0,4,5],[-2,-4,0,6],[-3,-5,-6,0]]) == 8
assert pfaffian([[0,0,1,0],[0,0,0,1],[-1,0,0,0],[0,-1,0,0]]) == -1
assert moment(0)==1 and moment(2)==1 and moment(4)==3 and moment(3)==0
assert skew_moment(0,1)==2 and skew_moment(0,3)==5
for beta in (1,2,4):
assert integral_components(beta,1)==(QQ(1),1,1)
assert integral_components(1,2)==(QQ(4),0,1)
assert integral_components(2,2)==(QQ(4),0,2)
assert integral_components(4,2)==(QQ(24),0,2)
assert integral_components(1,3)==(QQ(6),1,2)
def run(stored=False):
controls()
generator=generate.MehtaIntegrals()
selected=list(generator.enumerate())
parts={(int(p['beta']),int(p['n'])):generate.exact_components(p['beta'],p['n']) for p in selected}
assert not agrees_with(parts,lambda key:integral_components(*key))
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
values={key:generator.value({'beta':str(key[0]),'n':str(key[1])},100) for key in parts}
assert not exactness(values)
texts={key:to_text(value,digits=100) for key,value in values.items()}
for key,text in texts.items():
contained(values[key],text)
contained(enclosure(integral_components(*key)),text)
# Check the dimension recurrence numerically from independently integrated values.
b=RealBallField(640)
recurrences=0
for beta in (1,2,4):
for n in range(2,generate.MAX_N+1):
ratio=(b(2)*b.pi()).sqrt()*b(1+QQ(beta*n)/2).gamma()/b(1+QQ(beta)/2).gamma()
lhs=enclosure(integral_components(beta,n))
rhs=enclosure(integral_components(beta,n-1))*ratio
assert rhs.is_finite() and rhs.accuracy() > 580 and lhs.overlaps(rhs)
recurrences+=1
stored_count=0
if stored:
client=generate.authenticated_client()
doc=numberdb.table(generator.table,client=client)
records=doc['Numbers']
assert isinstance(records,list) and len(records)==len(parts)
expected=[(p['beta'],p['n']) for p in selected]
observed=[(str(r['params']['beta']),str(r['params']['n'])) for r in records]
assert expected==observed
by_key={}
for row in records:
key=(int(row['params']['beta']),int(row['params']['n']))
text=row['number']; assert text==texts[key]
contained(enclosure(integral_components(*key)),text)
by_key[key]=RealBallField(640)(decimal_bounds(text)[0]).union(RealBallField(640)(decimal_bounds(text)[1]))
stored_count+=1
for beta in (1,2,4):
for n in range(2,generate.MAX_N+1):
ratio=(b(2)*b.pi()).sqrt()*b(1+QQ(beta*n)/2).gamma()/b(1+QQ(beta)/2).gamma()
assert by_key[beta,n].overlaps(by_key[beta,n-1]*ratio)
source=Path(__file__).with_name('skill-evidence.json')
evidence={'table':generator.table,'entries':len(parts),'selection':{'beta':[1,2,4],'n_min':1,'n_max':generate.MAX_N},
'controls_passed':True,'independent_exact_component_matches':len(parts),
'generated_enclosures_contained':len(parts),'independent_enclosures_contained':len(parts),
'dimension_recurrences_checked':recurrences,'stored_entries_checked':stored_count,
'working_bits':numberdb.bits(100,losing=generate.WORKING_GUARD),
'minimum_relative_accuracy_bits':min(int(v.accuracy()) for v in values.values()),
'independent_bits':640,'minimum_independent_relative_accuracy_bits':min(int(enclosure(integral_components(*key)).accuracy()) for key in parts),
'serialized_measure':measure(texts),'computed_measure':measure(values),
'exact_rational_rows':0,'all_sqrt_pi_exponents_positive':True,
'versions':{'sage':SAGE_VERSION,'python':platform.python_version(),'numberdb':getattr(numberdb,'__version__','unknown')},
'skill_evidence':json.loads(source.read_text())}
print(json.dumps(evidence,indent=2,default=str))
if __name__=='__main__':
import sys
run('--stored' in sys.argv or os.environ.get('NUMBERDB_CHECK_STORED')=='1')