back to table · edit · history · where entries came from · files · download
13533 bytes, as of the version from 2026-10-27 01:06 (current). Recorded here, not run.
"""Independent complete-selection checks, with no website writes.
Run ``python balanced_checks.py`` beside generate.py to check table.json,
or ``python balanced_checks.py --numbers numbers.json`` for a working copy.
Scorer enclosures use certified Hi quadrature on [0,12] and Bi-Hi for Gi;
the omitted Hi tail is bounded by exp(12*x-576)/(pi*(144-x)). Nonzero
polynomial and exponential controls are tested first. All stored decimal
intervals must contain both the generator and independent enclosures.
All values also undergo independent mpmath computations at two precisions.
"""
import argparse
import hashlib
import importlib.metadata
import json
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import mpmath as mp
from flint import acb, arb, ctx
from balanced_values import (
BalancedValues, KINDS, argument_for, as_interval, ball, certified_value,
critical_cells, evaluate, mp_fraction, precision, scorer, sign, widen,
working_precision,
)
def bounds(written):
text = str(written).strip()
if '+/-' in text:
center, radius = (Fraction(part.strip()) for part in text.split('+/-'))
return center - radius, center + radius
if text.startswith('[') and text.endswith(']'):
return tuple(Fraction(part.strip()) for part in text[1:-1].split(','))
if '.' not in text and 'e' not in text.lower():
rational = Fraction(text)
return rational, rational
decimal = Decimal(text)
center = Fraction(decimal)
unit = Fraction(10) ** decimal.as_tuple().exponent
return center - unit, center + unit
def digits_in(written):
text = str(written).strip()
if '+/-' in text or text.startswith('['):
lower, upper = bounds(text)
with mp.workdps(40):
if lower == upper:
return 100
scale = max(abs(lower), abs(upper))
return max(100, int(mp.ceil(mp.log10(mp_fraction(scale / (upper - lower))))) + 2)
if '.' not in text and 'e' not in text.lower():
return 100
return max(100, len(Decimal(text).as_tuple().digits))
def assert_contained(written, interval, label):
lower, upper = bounds(written)
if isinstance(interval, int):
interval_lower = interval_upper = Fraction(interval)
else:
interval_lower, interval_upper = interval.lower, interval.upper
if not lower <= interval_lower <= interval_upper <= upper:
raise ArithmeticError('enclosure outside written interval: ' + label)
def mp_argument(identity):
prefixes = {'ai-prime-zero-': (mp.airyaizero, 1),
'bi-prime-zero-': (mp.airybizero, 1),
'ai-zero-': (mp.airyaizero, 0), 'bi-zero-': (mp.airybizero, 0)}
for prefix, (function, derivative) in prefixes.items():
if identity.startswith(prefix):
return function(int(identity[len(prefix):]), derivative=derivative)
if identity.startswith('gi-critical-'):
index = int(identity[len('gi-critical-'):])
lower, upper = critical_cells()[index]
return mp.findroot(lambda position: mp.diff(mp.scorergi, position),
(mp_fraction(lower), mp_fraction(upper)))
if identity.endswith('/cbrt(3)'):
return int(identity.split('/')[0]) / mp.root(3, 3)
if identity.endswith('*sqrt(2/pi)'):
return int(identity.split('*')[0]) * mp.sqrt(2 / mp.pi)
if identity.startswith('sqrt('):
return mp.sqrt(mp_fraction(Fraction(identity[5:-1])))
if identity.endswith('*pi/2'):
return int(identity[:-5]) * mp.pi / 2
if identity.endswith('*pi'):
return int(identity[:-3]) * mp.pi
if identity.endswith('/sqrt(2)'):
return int(identity[:-8]) / mp.sqrt(2)
return mp_fraction(Fraction(identity))
def mp_value(kind, argument):
if kind in ('ai', 'ai_prime', 'bi', 'bi_prime'):
function = mp.airyai if kind.startswith('ai') else mp.airybi
return function(argument, derivative=int(kind.endswith('prime')))
names = {'gi': 'scorergi', 'hi': 'scorerhi',
'fresnel_s': 'fresnels', 'fresnel_c': 'fresnelc'}
return getattr(mp, names.get(kind, kind))(argument)
def quadrature_hi(argument, digits):
if not argument < 144:
raise ArithmeticError('Hi tail bound requires x < 144')
tolerance = arb(10) ** (-digits - 10)
integral = acb(0)
for left in range(0, 12, 3):
integral += acb.integral(
lambda position, analytic: (argument * position - position ** 3 / 3).exp(),
left, left + 3, rel_tol=tolerance, abs_tol=tolerance)
integral /= arb.pi()
tail = (12 * argument - arb(576)).exp() / (arb.pi() * (144 - argument))
if not integral.imag.contains(0):
raise ArithmeticError('real Hi integral acquired an imaginary part')
result = widen(integral.real, tail)
if not result.is_finite() or result.rad() > arb(10) ** (-digits):
raise ArithmeticError('Hi quadrature failed its absolute error target')
return result
def controls():
with precision(working_precision(140)):
tolerance = arb('1e-150')
polynomial = acb.integral(lambda position, analytic: position ** 2,
0, 3, abs_tol=tolerance, rel_tol=tolerance)
exponential = acb.integral(lambda position, analytic: (-position).exp(),
0, 1, abs_tol=tolerance, rel_tol=tolerance)
for result in (polynomial, exponential):
if not result.is_finite() or result.real.rad() > tolerance:
raise ArithmeticError('quadrature control nonfinite or too wide')
if not polynomial.real.contains(9):
raise ArithmeticError('polynomial control failed')
if not (exponential.real - 1 + arb(-1).exp()).contains(0):
raise ArithmeticError('exponential control failed')
expected = 1 / (arb(3) ** arb('2/3') * arb('2/3').gamma())
if not (arb(0).airy_ai() - expected).contains(0):
raise ArithmeticError('Airy origin gamma identity failed')
if arb(0).erf() != 0 or arb(0).erfc() != 1 or arb(0).fresnel_s() != 0:
raise ArithmeticError('exact-anchor control failed')
with mp.workdps(160):
if abs(mp.quad(lambda position: position ** 2, [0, 3]) - 9) > mp.mpf('1e-150'):
raise ArithmeticError('mpmath quadrature control failed')
return ['integral(t^2,0,3)=9', 'integral(exp(-t),0,1)=1-exp(-1)',
'Ai(0) gamma identity', 'exact erf/erfc/Fresnel origins',
'independent mpmath polynomial quadrature']
def root_bracket(identity, argument):
prefixes = {'ai-prime-zero-': 1, 'bi-prime-zero-': 3,
'ai-zero-': 0, 'bi-zero-': 2}
for prefix, component in prefixes.items():
if identity.startswith(prefix):
endpoints = as_interval(argument)
with precision(ctx.prec * 2):
lower = ball(endpoints.lower).airy()[component]
upper = ball(endpoints.upper).airy()[component]
if sign(lower) * sign(upper) != -1:
raise ArithmeticError('Airy root enclosure failed endpoint sign test')
return True
return identity.startswith('gi-critical-')
def transformed_integral(kind, identity, family, written, digits):
if family not in ('unit-quadratic', 'unit-gaussian', 'simple-integral',
'complementary-phase'):
return None
with mp.workdps(digits + 70):
argument = mp_argument(identity)
if family == 'unit-quadratic':
endpoint = int(identity.split('*')[0])
function = mp.sin if kind == 'fresnel_s' else mp.cos
value = mp.sqrt(2 / mp.pi) * mp.quad(
lambda position: function(position ** 2), [0, endpoint])
elif family == 'unit-gaussian':
coefficient = argument ** 2
integral = mp.quad(lambda position: mp.exp(-coefficient * position ** 2), [0, 1])
value = 2 * mp.sqrt(coefficient / mp.pi) * integral
if kind == 'erfc':
value = 1 - value
elif kind == 'si':
value = mp.quad(lambda position: mp.sin(argument * position) / position
if position else argument, [0, 1])
else:
value = mp.euler + mp.log(argument) + mp.quad(
lambda position: (mp.cos(argument * position) - 1) / position
if position else 0, [0, 1])
lower, upper = bounds(written)
if not mp_fraction(lower) <= value <= mp_fraction(upper):
raise ArithmeticError('transformed integral disagreement: ' + identity)
return 'independent finite-interval mpmath quadrature at %d dps' % (digits + 70)
def check(generator, numbers, verbose=True):
if set(numbers) != set(generator.entries):
raise ArithmeticError('stored and default-enumerated argument sets differ')
control_records = controls()
kind = KINDS[generator.table]
records = []
for identity, selection in generator.entries.items():
entry = numbers[identity]
written = entry['number'] if isinstance(entry, dict) else entry
digits = selection['digits']
generated = certified_value(kind, identity, digits)
assert_contained(written, generated, identity)
record = {'argument': identity, 'family': selection['family'],
'digits': digits, 'certified_interval_contained': True}
with precision(working_precision(digits + 40)):
argument = argument_for(identity, digits + 40)
value = evaluate(kind, argument, digits + 40)
if not value.is_finite():
raise ArithmeticError('nonfinite structural-check value')
record['arb_accuracy_bits'] = None if value.is_exact() else value.rel_accuracy_bits()
record['root_bracket_checked'] = root_bracket(identity, argument)
if kind in ('gi', 'hi'):
integral = quadrature_hi(argument, digits + 30)
independent = integral if kind == 'hi' else argument.airy_bi() - integral
assert_contained(written, as_interval(independent), identity + ': quadrature')
record['quadrature_interval_contained'] = True
record['quadrature_accuracy_bits'] = independent.rel_accuracy_bits()
second = scorer(argument, kind, derivative=2, digits=digits + 40)
forcing = (-1 if kind == 'gi' else 1) / arb.pi()
if not (second - argument * value - forcing).contains(0):
raise ArithmeticError('Scorer differential equation failed')
if not (scorer(argument, 'gi', digits=digits + 40)
+ scorer(argument, 'hi', digits=digits + 40)
- argument.airy_bi()).contains(0):
raise ArithmeticError('Gi+Hi=Bi failed')
if kind in ('ai', 'ai_prime', 'bi', 'bi_prime'):
airy = argument.airy()
if not (airy[0] * airy[3] - airy[1] * airy[2] - 1 / arb.pi()).contains(0):
raise ArithmeticError('Airy Wronskian failed')
if kind in ('erf', 'erfc'):
if not (argument.erf() + argument.erfc() - 1).contains(0):
raise ArithmeticError('Gaussian complement failed')
if not (argument.erf() + (-argument).erf()).contains(0):
raise ArithmeticError('Gaussian reflection failed')
if kind in ('si', 'fresnel_s', 'fresnel_c'):
if not (evaluate(kind, argument) + evaluate(kind, -argument)).contains(0):
raise ArithmeticError('odd reflection failed')
precisions = [digits + 50, digits + 100]
for decimal_digits in precisions:
with mp.workdps(decimal_digits):
independent = mp_value(kind, mp_argument(identity))
lower, upper = bounds(written)
if not mp.isfinite(independent) or not (
mp_fraction(lower) <= independent <= mp_fraction(upper)):
raise ArithmeticError('independent mpmath mismatch: ' + identity)
record['mpmath_decimal_precisions'] = precisions
record['transformed_integral'] = transformed_integral(
kind, identity, selection['family'], written, digits)
records.append(record)
if verbose:
print(generator.table, identity, 'checked', flush=True)
evidence = {'table': generator.table, 'entries': len(numbers),
'controls': control_records, 'records': records,
'versions': {package: importlib.metadata.version(package)
for package in ('numberdb', 'python-flint', 'mpmath')}}
evidence['source_sha256'] = {
name: hashlib.sha256(Path(__file__).with_name(name).read_bytes()).hexdigest()
for name in ('balanced_values.py', 'balanced_checks.py')}
return evidence
if __name__ == '__main__':
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--numbers', type=Path)
parser.add_argument('--output', type=Path, default=Path('balanced-checks-evidence.json'))
options = parser.parse_args()
data = (json.loads(options.numbers.read_text()) if options.numbers
else json.loads(Path('table.json').read_text())['Numbers'])
evidence = check(BalancedValues(), data)
options.output.write_text(json.dumps(evidence, indent=2) + '\n')