back to table · edit · history · where entries came from · files · download
8003 bytes, as of the version from 2026-10-26 15:39 (current). Recorded here, not run.
"""Independent origin controls without redoing expensive Taylor integrations.
pip install numberdb python-flint
python probability_controls.py
Nyström determinants and their derivatives check all three laws at zero;
the Airy-kernel resolvent independently checks q(0). Convergence between
two quadrature configurations is heuristic, not an enclosure of quadrature
and infinite-interval errors. No values are changed or approved here.
"""
import argparse
import json
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import numberdb
from flint import arb, arb_mat, ctx
def written_bounds(text):
center = Fraction(Decimal(str(text)))
unit = Fraction(10) ** Decimal(str(text)).as_tuple().exponent
return center - unit, center + unit
def endpoint(value):
mantissa, exponent = value.man_exp()
return Fraction(int(mantissa)) * Fraction(2) ** int(exponent)
def controls(nodes, length, digits):
ctx.dps = digits
quadrature = [arb.legendre_p_root(nodes, index, weight=True) for index in range(nodes)]
points = [(length * (point + 1) / 2).mid() for point, weight in quadrature]
weights = [(length * weight / 2).sqrt().mid() for point, weight in quadrature]
kernel, derivative = arb_mat(nodes, nodes), arb_mat(nodes, nodes)
for row in range(nodes):
for column in range(row, nodes):
value, slope, _, _ = (points[row] + points[column]).airy()
multiplier = weights[row] * weights[column]
kernel[row, column] = kernel[column, row] = (multiplier * value).mid()
derivative[row, column] = derivative[column, row] = (multiplier * slope).mid()
determinants, slopes = [], []
for sign in (-1, 1):
matrix = sign * kernel
for index in range(nodes):
matrix[index, index] += 1
determinant = matrix.det()
determinants.append(determinant)
slopes.append(determinant * matrix.solve(sign * derivative).trace())
minus, plus = determinants
minus_slope, plus_slope = slopes
cdf = {'1': minus, '2': minus * plus, '4': (minus + plus) / 2}
density = {'1': minus_slope, '2': minus_slope * plus + minus * plus_slope,
'4': (minus_slope + plus_slope) / arb(2).sqrt()}
airy = [point.airy()[:2] for point in points]
matrix, rhs = arb_mat(nodes, nodes), arb_mat(nodes, 1)
for row in range(nodes):
rhs[row, 0] = (weights[row] * airy[row][0]).mid()
for column in range(row, nodes):
value = (airy[row][1] ** 2 - points[row] * airy[row][0] ** 2 if row == column else
(airy[row][0] * airy[column][1] - airy[row][1] * airy[column][0]) /
(points[row] - points[column]))
weighted = (weights[row] * weights[column] * value).mid()
matrix[row, column] = matrix[column, row] = -weighted
matrix[row, row] += 1
solution = matrix.solve(rhs)
ai_zero, aip_zero, _, _ = arb(0).airy()
hastings = ai_zero
for column in range(nodes):
value = (ai_zero * airy[column][1] - aip_zero * airy[column][0]) / (-points[column])
hastings += value * weights[column] * solution[column, 0]
return {'T441': cdf, 'T451': density, 'T443': {'q': hastings}}
def number(entry):
return entry['number'] if isinstance(entry, dict) else entry
if __name__ == '__main__':
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--archive-dir', type=Path)
parser.add_argument('--output', type=Path)
options = parser.parse_args()
baseline = {}
for tid in ('T441', 'T443', 'T451'):
baseline[tid] = (json.loads((options.archive_dir / tid / 'baseline.json').read_text())
if options.archive_dir else {'document': numberdb.table(tid)})
records = []
first = None
worst_difference = arb(0)
cache = (json.loads(options.output.read_text()).get('configurations', [])
if options.output and options.output.exists() else [])
for nodes, length, digits in ((170, 36, 170), (210, 44, 210)):
cached = next((record for record in cache if record['nodes'] == nodes and
record['interval_length'] == length and record['working_digits'] == digits), None)
if cached and first is None:
records.append(cached)
first = cached['values']
continue
computed = controls(nodes, length, digits)
values = {}
for tid, laws in computed.items():
values[tid] = {}
for law, value in laws.items():
if not value.is_finite() or not value > 0 or value.rel_accuracy_bits() < 350:
raise ArithmeticError('nonfinite, nonpositive or wide Fredholm control')
entry = (baseline[tid]['document']['Numbers']['0'] if tid == 'T443'
else baseline[tid]['document']['Numbers'][law]['0'])
lower, upper = written_bounds(number(entry))
if not lower < endpoint(value.mid()) < upper:
raise ArithmeticError('independent origin control differs from stored 100 digits: ' + tid)
values[tid][law] = value.mid().str(130, radius=False)
if first is not None:
for tid, laws in values.items():
for law, text in laws.items():
difference = abs((arb(text) - arb(first[tid][law])) / arb(text))
worst_difference = max(worst_difference, difference)
if not difference < arb('1e-102'):
raise ArithmeticError('quadrature settings differ by ' + str(difference) + ' at ' + tid)
first = values
records.append({'nodes': nodes, 'interval_length': length, 'working_digits': digits,
'values': values, 'inside_stored_100_digit_intervals': True})
print('Fredholm controls passed', nodes, length, digits, flush=True)
report = {
'scope': 'seven origin values only, not a full numerical revalidation of 379 entries',
'rigour': 'heuristic (agreement-checked)', 'configurations': records,
'omitted_error_bounds': ['quadrature truncation', 'infinite-interval truncation'],
'independent_of': 'stored Painleve-II Taylor integration',
'convergence_relative_tolerance': '1e-102',
'worst_relative_difference': worst_difference.str(15),
'complete': False}
if options.output:
options.output.write_text(json.dumps(report, indent=2) + '\n')
for tid in ('T441', 'T443', 'T451'):
numbers = baseline[tid]['document']['Numbers']
ctx.prec = 800
if tid == 'T441':
for law in numbers.values():
values = [arb(str(number(law[key]))) for key in sorted(law, key=Fraction)]
if not all(0 < value < 1 for value in values):
raise ArithmeticError('CDF outside probability range')
if not all(left < right for left, right in zip(values, values[1:])):
raise ArithmeticError('CDF not strictly increasing')
elif tid == 'T443':
values = [arb(str(number(numbers[key]))) for key in sorted(numbers, key=Fraction)]
if not all(value > 0 for value in values) or not all(
left > right for left, right in zip(values, values[1:])):
raise ArithmeticError('Hastings-McLeod sign or monotonicity failure')
elif not all(arb(str(number(entry))) > 0 for law in numbers.values() for entry in law.values()):
raise ArithmeticError('density not positive')
print('All retained probability entries pass structural checks', flush=True)
report['complete'] = True
report['all_entries_structurally_checked'] = True
if options.output:
options.output.write_text(json.dumps(report, indent=2) + '\n')
else:
print(json.dumps(report, indent=2))