back to table · edit · history · where entries came from · files · download
3866 bytes, as of the version from 2026-10-27 22:51 (current). Recorded here, not run.
"""Selberg integrals for classical Jacobi ensembles -- numberdb.org/T563.
Run with SageMath or compatible passagemath and the numberdb package:
sage -pip install numberdb # install the client once
sage -python generate.py # verify every stored row
sage -python generate.py --fill # write entries to this draft
Set NUMBERDB_API_KEY for access to the draft. --fill never offers or publishes
it. In the campaign, pipe the key with NUMBERDB_KEY_FROM_STDIN=1 and set
NUMBERDB_PUBLISH=1 to fill through agents/sage.sh.
The full-cube Selberg product is reduced to q*pi^k using exact integer and
half-integer gamma identities. k=0 returns an exact rational. Other rows
return real balls at numberdb.bits(digits, losing=96); the build checks
measure retained accuracy and verify the actual serialized intervals.
Selection and conventions are explained in argument-selection.md.
"""
import os
import sys
from functools import lru_cache
from math import factorial
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
WORKING_GUARD = 96
PARAMETERS = ('beta', 'n', 'r', 's')
def selection():
"""Union of small-matrix, near-square, and fixed-aspect-ratio families."""
triples = set()
for n in range(1, 7):
triples.update((n, r, s) for s in range(5) for r in range(s + 1))
for n in range(1, 21):
triples.update((n, r, s) for s in range(3) for r in range(s + 1))
for n in range(1, 13):
triples.update([(n, 0, n), (n, n, n), (n, n, 2*n)])
for beta in (1, 2, 4):
for n, r, s in sorted(triples):
yield dict(zip(PARAMETERS, (beta, n, r, s)))
@lru_cache(None)
def gamma_half(twice_argument):
"""Gamma(t/2) as a rational times sqrt(pi)^e, for positive integer t."""
t = int(twice_argument)
if t < 1:
raise ValueError('gamma argument must be positive')
if t % 2 == 0:
return QQ(factorial(t//2 - 1)), 0
m = (t - 1)//2
return QQ(factorial(2*m)) / (ZZ(4)**m * ZZ(factorial(m))), 1
def exact_form(beta, n, r, s):
"""Return q,k such that Z=q*pi^k, allowing either order of r,s."""
beta, n, r, s = map(int, (beta, n, r, s))
if beta not in (1, 2, 4) or n < 1 or min(r, s) < 0:
raise ValueError('outside the classical rectangular Jacobi domain')
q, e = QQ(1), 0
for j in range(n):
numerators = (beta*(r+1+j), beta*(s+1+j), 2+beta*(j+1))
denominators = (beta*(r+s+n+j+1), 2+beta)
for arg in numerators:
coeff, power = gamma_half(arg)
q *= coeff
e += power
for arg in denominators:
coeff, power = gamma_half(arg)
q /= coeff
e -= power
assert e % 2 == 0
return q, e//2
class JacobiSelberg(numberdb.Generator):
table = 'T563'
parameters = PARAMETERS
type = 'R'
digits = 100
rigour = 'proven'
def enumerate(self):
yield from selection()
def value(self, params, digits):
q, k = exact_form(*(params[name] for name in PARAMETERS))
if k == 0:
return q
field = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
return field(q) * field.pi()**k
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
generator = JacobiSelberg()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Fill checked full-cube Jacobi normalizations',
assisted_by=os.environ.get('NUMBERDB_ASSISTED_BY', 'codex-cli')))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)