back to table · edit · history · where entries came from · files · download
4145 bytes, as of the version from 2026-10-28 00:41 (current). Recorded here, not run.
"""Mehta integrals for Gaussian ensembles -- numberdb.org/T566.
Run with SageMath (or a compatible passagemath installation):
sage -pip install numberdb
sage -python generate.py # verify every stored row
sage -python generate.py --fill # write entries to this draft
sage -python checks.py # independent integral checks
--fill requires NUMBERDB_API_KEY. It only writes entries to the existing
private draft; it does not offer, review or publish the table. To pass a key
through a wrapper, set NUMBERDB_KEY_FROM_STDIN=1 and pipe the key on stdin;
NUMBERDB_PUBLISH=1 selects --fill for wrappers that do not forward arguments.
The full labelled-coordinate integral has weight exp(-sum(x_i^2)/2).
Gamma factors are reduced exactly to q * sqrt(2)^s * sqrt(pi)^p before any
ball arithmetic. For n>=1 the power p is positive, so no selected row is
rational. The selection is beta=1,2,4 and n=1,...,50; see argument-selection.md.
A 96-bit guard is measured by checks.py, which also verifies the serialized
100-digit intervals against independent exact Gaussian moment integrals.
"""
import os
import sys
from math import factorial
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
MAX_N = 50
DYSON_EXPONENTS = (1, 2, 4)
WORKING_GUARD = 96
def half_gamma(twice_argument):
"""Return (q,p) for Gamma(k/2)=q*sqrt(pi)^p, k a positive integer."""
k = int(twice_argument)
if k < 1:
raise ValueError('Gamma argument must be positive')
if k % 2 == 0:
return QQ(factorial(k // 2 - 1)), 0
m = (k - 1) // 2
return QQ(factorial(2 * m)) / QQ(4 ** m * factorial(m)), 1
def exact_components(beta, n):
"""Return q,s,p for Z=q*sqrt(2)^s*sqrt(pi)^p with s=0 or 1."""
beta, n = int(beta), int(n)
if beta not in DYSON_EXPONENTS or n < 1:
raise ValueError('Expected beta in (1,2,4) and n>=1')
q = QQ(2) ** (n // 2)
p = n
denominator, dp = half_gamma(2 + beta)
for j in range(1, n + 1):
numerator, np = half_gamma(2 + j * beta)
q *= numerator / denominator
p += np - dp
assert q > 0 and p > 0
return q, n % 2, p
def evaluate_components(parts, precision):
q, s, p = parts
field = RealBallField(precision)
result = field(q) * field(2).sqrt() ** s * field.pi().sqrt() ** p
if not result.is_finite() or not result > 0:
raise ArithmeticError('Nonpositive or nonfinite integral enclosure')
return result
class MehtaIntegrals(numberdb.Generator):
table = 'T566'
parameters = ('beta', 'n')
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'checks.py', 'check.py', 'argument-selection.md',
'independent-checks.md', 'numerical-evidence.json', 'skill-evidence.json')
def enumerate(self, max_n=MAX_N):
for beta in DYSON_EXPONENTS:
for n in range(1, int(max_n) + 1):
yield {'beta': str(beta), 'n': str(n)}
def value(self, params, digits):
parts = exact_components(params['beta'], params['n'])
result = evaluate_components(parts, numberdb.bits(digits, losing=WORKING_GUARD))
if int(result.accuracy()) < numberdb.bits(digits):
raise ArithmeticError('Insufficient relative accuracy')
return result
def authenticated_client():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
key = sys.stdin.read().strip()
if '=' in key and key.split('=', 1)[0].isidentifier():
key = key.split('=', 1)[1].strip().strip('\"\'')
return numberdb.Client(api_key=key)
return numberdb.Client()
if __name__ == '__main__':
generator = MehtaIntegrals()
client = authenticated_client()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
outcome = generator.publish(client=client, message='Fill 150 Mehta integrals certified by independent Gaussian moment checks')
print(outcome)
report = generator.verify(sample=None, client=client)
print(report)
sys.exit(0 if report.ok else 1)