back to table · edit · history · where entries came from · files · download
15518 bytes, as of the version from 2026-09-26 18:54 (current). Recorded here, not run.
"""Reproduce T441 Tracy-Widom distribution function values at 100 digits.
Run with Python and the public dependencies `numberdb` and `python-flint`:
$ python generate.py --sample
$ python generate.py --compute-json entries.json --diagnostics-json diagnostics.json
$ python generate.py --check-fredholm --fredholm-json fredholm-controls.json
$ python generate.py
The final command verifies against the website table when network access is
available. Offline checks use `--compute-json`, `--sample`, and
`--check-fredholm`. Publishing is disabled unless `--publish` is explicitly
supplied and `NUMBERDB_API_KEY` is present.
The arguments are unchanged: beta=1,2,4 and every reduced rational s in
[-6,3] with denominator at most four.
The computation follows the Painleve-II formulas already stated in T441. Let
q be the Hastings-McLeod solution, u(s)=int_s^infty q(x) dx, and
M(s)=int_s^infty (x-s) q(x)^2 dx. Then
F1(s) = exp(-M(s)/2) exp(-u(s)/2),
F2(s) = exp(-M(s)),
F4(s) = exp(-M(sqrt(2)s)/2) cosh(u(sqrt(2)s)/2).
The beta=4 formula therefore keeps the table's Bornemann sqrt(2) scaling.
Two Taylor integrations vary working precision, Taylor order, maximum step,
and right boundary: (230 digits, 180 terms, 1/8, 40) and
(270 digits, 230 terms, 1/10, 46). At the right boundary q and q' are Airy
values, and the integrals use Airy tail data rather than being set to zero.
Every target checks J=q'^2-s*q^2-q^4. Every CDF value must agree between the
two runs to relative 1e-110 before 100 digits are written.
Arb is used for multiprecision arithmetic, but step results are replaced by
midpoints. Neither the Taylor remainder nor the nonlinear right-boundary error
is rigorously enclosed. The rigour is heuristic (agreement-checked), not
proven.
"""
import argparse
import functools
import json
import os
from decimal import Decimal, localcontext
from fractions import Fraction
from pathlib import Path
import flint
from flint import acb, arb, arb_mat, arb_poly, arb_series, ctx
import numberdb
DIGITS = 100
CONFIGURATIONS = ((230, 180, 8, 40), (270, 230, 10, 46))
AGREEMENT_TOLERANCE = arb("1e-110")
IDENTITY_TOLERANCE = arb("1e-112")
FREDHOLM_POINTS = (Fraction(-6), Fraction(0), Fraction(3))
def arguments():
return sorted({Fraction(numerator, denominator)
for denominator in range(1, 5)
for numerator in range(-6 * denominator, 3 * denominator + 1)})
def _identity(beta, argument):
return str(beta) + "," + str(Fraction(argument))
def _format_decimal(text, digits=DIGITS):
with localcontext() as decimal_context:
decimal_context.prec = max(280, digits + 80)
return format(Decimal(text), f".{digits - 1}e")
def initial_state(boundary):
value, derivative, _, _ = boundary.airy()
integral = derivative**2 - boundary * value**2
moment = (2 * boundary**2 * value**2 - 2 * boundary * derivative**2 - value * derivative) / 3
tail = acb.integral(lambda point, analytic: point.airy_ai(), boundary, 100,
abs_tol=arb("1e-190"), rel_tol=arb("1e-170"))
if not tail.is_finite() or not tail.imag.contains(0) or not tail.real > 0:
raise ArithmeticError("Invalid Airy tail integral")
if not tail.real.rad() < arb("1e-185"):
raise ArithmeticError("Airy tail integral is insufficiently accurate")
return [value.mid(), derivative.mid(), tail.real.mid(), integral.mid(), moment.mid()]
def taylor_step(position, state, increment, order):
value, derivative, tail, integral, moment = state
series = arb_series([value, derivative], prec=order)
variable = arb_series([position, 1], prec=order)
for precision in range(4, order + 1, 2):
ctx.cap = precision
series = arb_series(series.coeffs(), prec=precision)
series = arb_series([value, derivative], prec=precision) + (
variable * series + 2 * series**3).integral().integral()
ctx.cap = order
polynomial = arb_poly(series.coeffs())
square_integral = (series * series).integral()
result = [
polynomial(increment),
polynomial.derivative()(increment),
tail - polynomial.integral()(increment),
integral - arb_poly(square_integral.coeffs())(increment),
moment - integral * increment + arb_poly(square_integral.integral().coeffs())(increment),
]
if not all(value.is_finite() for value in result):
raise ArithmeticError("Non-finite Taylor result")
return [value.mid() for value in result]
def _cdf_from_state(beta, state):
value, derivative, tail, integral, moment = state
factor = (-moment / 2).exp()
if beta == "1":
result = factor * (-tail / 2).exp()
elif beta == "2":
result = factor**2
elif beta == "4":
result = factor * (tail / 2).cosh()
else:
raise ValueError(f"unknown beta {beta!r}")
if not result.is_finite() or not result > 0 or not result < 1:
raise ArithmeticError(f"Invalid CDF value for beta={beta}")
return result
def _target_set(only=None):
targets = {}
if only is None:
needed = {(_identity("1", argument), "plain", argument, "1") for argument in arguments()}
needed.update({(_identity("2", argument), "plain", argument, "2") for argument in arguments()})
needed.update({(_identity("4", argument), "sqrt2", argument, "4") for argument in arguments()})
else:
needed = set()
for identity in only:
beta, argument_text = identity.split(",", 1)
argument = Fraction(argument_text)
scaling = "sqrt2" if beta == "4" else "plain"
needed.add((identity, scaling, argument, beta))
for identity, scaling, argument, beta in needed:
targets.setdefault((scaling, argument), []).append((identity, beta))
return targets
def compute_configuration(configuration, only=None):
working_digits, order, step_denominator, boundary = configuration
old_precision, old_cap = ctx.prec, ctx.cap
try:
ctx.dps, ctx.cap = working_digits, order
grouped_targets = _target_set(only)
march_targets = []
for scaling, argument in grouped_targets:
point = arb(argument.numerator) / argument.denominator
if scaling == "sqrt2":
point *= arb(2).sqrt()
march_targets.append((point, scaling, argument))
march_targets.sort(key=lambda target: target[0], reverse=True)
position = arb(boundary)
state = initial_state(position)
maximum_step = arb(1) / step_denominator
samples = {}
worst_identity_error = arb(0)
for target, scaling, argument in march_targets:
while position - target > maximum_step:
state = taylor_step(position, state, -maximum_step, order)
position -= maximum_step
increment = (target - position).mid()
if not increment.is_zero():
state = taylor_step(position, state, increment, order)
position = target.mid()
value, derivative, tail, integral, moment = state
identity_error = abs((derivative**2 - position * value**2 - value**4) - integral)
worst_identity_error = max(worst_identity_error, identity_error)
if not identity_error < IDENTITY_TOLERANCE:
raise ArithmeticError(
f"Hamiltonian control failed at {scaling} {argument}: {identity_error}")
for identity, beta in grouped_targets[(scaling, argument)]:
cdf = _cdf_from_state(beta, state)
samples[identity] = cdf.mid().str(working_digits - 5, radius=False)
return samples, {"worst_hamiltonian_error": worst_identity_error.str(12, radius=False)}
finally:
ctx.prec, ctx.cap = old_precision, old_cap
def checked_values(only=None, return_diagnostics=False):
first, first_diag = compute_configuration(CONFIGURATIONS[0], only=only)
second, second_diag = compute_configuration(CONFIGURATIONS[1], only=only)
old_precision = ctx.prec
try:
ctx.dps = 320
if first.keys() != second.keys():
raise ArithmeticError("The two runs cover different arguments")
worst_relative = arb(0)
worst_identity = None
for identity, text in second.items():
value = arb(text)
relative = abs((arb(first[identity]) - value) / value)
if relative > worst_relative:
worst_relative = relative
worst_identity = identity
if not relative < AGREEMENT_TOLERANCE:
raise ArithmeticError(
f"Insufficient numerical agreement at {identity}: {relative}")
values = {identity: _format_decimal(text) for identity, text in second.items()}
diagnostics = {
"configurations": [list(config) for config in CONFIGURATIONS],
"agreement_tolerance": "1e-110",
"identity_tolerance": "1e-112",
"worst_relative_difference": worst_relative.str(12, radius=False),
"worst_relative_difference_entry": worst_identity,
"first_run": first_diag,
"second_run": second_diag,
}
if return_diagnostics:
return values, diagnostics
return values
finally:
ctx.prec = old_precision
@functools.lru_cache(maxsize=1)
def all_checked_values():
return checked_values()
def fredholm_determinants(argument, nodes, length):
quadrature = [arb.legendre_p_root(nodes, index, weight=True) for index in range(nodes)]
abscissas = [length * (point + 1) / 2 for point, weight in quadrature]
roots = [(length * weight / 2).sqrt() for point, weight in quadrature]
kernel = arb_mat(nodes, nodes)
for row in range(nodes):
for column in range(row, nodes):
value, _, _, _ = (abscissas[row] + abscissas[column] + argument).airy()
weighted = roots[row] * value * roots[column]
kernel[row, column] = kernel[column, row] = weighted
determinants = {}
for label, sign in (("minus", -1), ("plus", 1)):
matrix = sign * kernel
for index in range(nodes):
matrix[index, index] += 1
determinant = matrix.det()
if not determinant.is_finite():
raise ArithmeticError(f"Non-finite Fredholm determinant at {argument}")
determinants[label] = determinant
return determinants
def fredholm_cdf_values(argument, nodes=260, length=52):
minus_plus = fredholm_determinants(argument, nodes, length)
minus, plus = minus_plus["minus"], minus_plus["plus"]
return {
"1": minus,
"2": minus * plus,
"4": (minus + plus) / 2,
}
def check_fredholm(nodes=260, length=52):
values = all_checked_values()
old_precision = ctx.prec
try:
ctx.dps = 220
controls = []
for argument in FREDHOLM_POINTS:
point = arb(argument.numerator) / argument.denominator
plain = fredholm_cdf_values(point, nodes=nodes, length=length)
for beta in ("1", "2"):
identity = _identity(beta, argument)
target = arb(values[identity])
relative = abs((plain[beta] - target) / target)
if not relative < arb("1e-95"):
raise ArithmeticError(
f"Fredholm control disagrees at {identity}: {relative}")
controls.append({
"entry": identity,
"nodes": nodes,
"length": length,
"relative_difference": relative.str(12, radius=False),
})
scaled_point = point * arb(2).sqrt()
scaled = fredholm_cdf_values(scaled_point, nodes=nodes, length=length)
identity = _identity("4", argument)
target = arb(values[identity])
relative = abs((scaled["4"] - target) / target)
if not relative < arb("1e-95"):
raise ArithmeticError(f"Fredholm control disagrees at {identity}: {relative}")
controls.append({
"entry": identity,
"nodes": nodes,
"length": length,
"relative_difference": relative.str(12, radius=False),
})
return controls
finally:
ctx.prec = old_precision
class TracyWidomDistributionFunctions(numberdb.Generator):
table = "T441"
parameters = ("beta", "s")
type = "R"
digits = DIGITS
rigour = "heuristic (agreement-checked)"
def enumerate(self):
for beta in ("1", "2", "4"):
for argument in arguments():
yield {"beta": beta, "s": str(argument)}
def value(self, params, digits):
if digits > DIGITS:
raise ValueError("This calibration supports at most 100 significant digits")
identity = _identity(str(params["beta"]), params["s"])
return {"number": all_checked_values()[identity], "digits": DIGITS}
def environment(self):
return {**super().environment(), "python-flint": flint.__version__}
def sample():
identities = ["1,-6", "2,0", "4,-6", "4,3"]
values, diagnostics = checked_values(only=identities, return_diagnostics=True)
return {"values": values, "diagnostics": diagnostics}
def main():
parser = argparse.ArgumentParser()
parser.add_argument("--sample", action="store_true")
parser.add_argument("--compute-json", type=Path)
parser.add_argument("--diagnostics-json", type=Path)
parser.add_argument("--check-fredholm", action="store_true")
parser.add_argument("--fredholm-json", type=Path)
parser.add_argument("--fredholm-nodes", type=int, default=260)
parser.add_argument("--fredholm-length", type=int, default=52)
parser.add_argument("--publish", action="store_true")
options = parser.parse_args()
generator = TracyWidomDistributionFunctions()
if options.sample:
result = sample()
print(json.dumps(result, indent=2))
return
if options.check_fredholm:
controls = check_fredholm(nodes=options.fredholm_nodes, length=options.fredholm_length)
if options.fredholm_json:
options.fredholm_json.write_text(json.dumps(controls, indent=2) + "\n")
print(json.dumps({"fredholm_controls": controls}, indent=2))
return
if options.compute_json:
values, diagnostics = checked_values(return_diagnostics=True)
records = []
for params in generator.enumerate():
identity = _identity(str(params["beta"]), params["s"])
records.append({"params": params, "number": values[identity], "digits": str(DIGITS)})
options.compute_json.write_text(json.dumps(records, indent=2) + "\n")
if options.diagnostics_json:
options.diagnostics_json.write_text(json.dumps(diagnostics, indent=2) + "\n")
print(f"Computed {len(records)} entries at {DIGITS} significant digits")
print(json.dumps(diagnostics, indent=2))
return
if options.publish:
print(generator.publish(
message="Refine Tracy-Widom distribution functions to 100 significant digits",
assisted_by=os.environ.get("NUMBERDB_ASSISTED_BY", "")))
return
report = generator.verify(sample=None)
print(report)
raise SystemExit(0 if report.ok else 1)
if __name__ == "__main__":
main()