back to table · edit · history · where entries came from · files · download
11391 bytes, as of the version from 2026-10-27 01:05 (current). Recorded here, not run.
"""Reproduce balanced Airy, Scorer, integral and Gaussian reference values.
Install with ``pip install numberdb python-flint mpmath``. Download the table
bundle, then run ``python generate.py`` to verify its entire attached selection
without networking or writing. ``python generate.py --independent`` adds the
independent checks. Website comparison is a separate fresh-readback step.
The selection manifest records all retained identities and their metadata.
No older repository module is needed. All arguments, including roots, remain
balls. Scorer values use the entire inhomogeneous Airy series with geometric
bounds on its three residue-class tails. The initial guard is 640 bits;
accuracy is measured and insufficient enclosures are rejected or escalated.
Gi root ordering uses an interval count on rational cells in [-16,1]. For a
cell midpoint m and half-width h, Gronwall bounds |Gi|+|Gi'| by
(|Gi(m)|+|Gi'(m)|+h/pi)*exp(16*h). This encloses Gi, Gi' and hence
Gi''=x*Gi-1/pi on the cell. Each root is then independently bracketed.
The certificate follows the previously attached scorer.py, not a decimal
approximation to a root. Original sources survive in revision history.
"""
import json
import math
import re
from contextlib import contextmanager
from fractions import Fraction
from functools import lru_cache
from pathlib import Path
import mpmath as mp
import numberdb
from flint import arb, ctx
KINDS = {
'T514': 'gi', 'T515': 'hi', 'T198': 'ai', 'T202': 'bi',
'T201': 'ai_prime', 'T203': 'bi_prime', 'T192': 'si', 'T193': 'ci',
'T207': 'fresnel_s', 'T208': 'fresnel_c', 'T197': 'erf', 'T204': 'erfc',
}
GUARD_BITS = 640
@contextmanager
def precision(bits):
previous = ctx.prec
ctx.prec = bits
try:
yield
finally:
ctx.prec = previous
def working_precision(digits):
if digits < 1:
raise ValueError('digits must be positive')
return math.ceil(digits * math.log2(10)) + GUARD_BITS
def ball(value):
return value if isinstance(value, arb) else arb(str(value))
def widen(value, error):
return value + arb(0, ball(error).abs_upper())
def exact_endpoint(value):
mantissa, exponent = value.man_exp()
return Fraction(int(mantissa)) * Fraction(2) ** int(exponent)
def as_interval(value):
if not value.is_finite():
raise ArithmeticError('nonfinite enclosure')
return numberdb.RealInterval(exact_endpoint(value.lower()),
exact_endpoint(value.upper()))
def sign(value):
if not value.is_finite():
raise ArithmeticError('nonfinite sign test')
return 1 if value > 0 else -1 if value < 0 else 0
@lru_cache(maxsize=16)
def scorer_initial(kind, bits):
if bits != ctx.prec or kind not in ('gi', 'hi'):
raise ValueError('invalid coefficient context')
multiplier = 1 if kind == 'gi' else 2
forcing = -1 if kind == 'gi' else 1
return (multiplier / (arb(3) ** arb('7/6') * arb('2/3').gamma()),
multiplier / (arb(3) ** arb('5/6') * arb('1/3').gamma()),
forcing / (2 * arb.pi()))
def scorer(argument, kind='gi', derivative=0, digits=100):
if derivative not in (0, 1, 2):
raise ValueError('supported derivative orders are zero through two')
argument = ball(argument)
radius = argument.abs_upper()
tolerance = arb(10) ** (-digits - 80)
answer = arb(0)
for residue, coefficient in enumerate(scorer_initial(kind, ctx.prec)):
index = residue
if index < derivative:
coefficient /= (index + 3) * (index + 2)
index += 3
factor = math.factorial(index) // math.factorial(index - derivative)
term = coefficient * factor * argument ** (index - derivative)
subtotal = arb(0)
while index < 20000:
subtotal += term
denominator = ((index + 3) * (index + 2),
index * (index + 2), index * (index - 1))[derivative]
ratio = radius ** 3 / denominator
if ratio < 1:
tail = abs(term) * ratio / (1 - ratio)
if tail < tolerance:
answer += widen(subtotal, tail)
break
term *= argument ** 3 / denominator
index += 3
else:
raise ArithmeticError('Taylor tail did not converge')
if not answer.is_finite():
raise ArithmeticError('nonfinite Taylor result')
return answer
@lru_cache(maxsize=4096)
def scorer_point(argument, bits):
if bits != ctx.prec:
raise ArithmeticError('point precision mismatch')
return scorer(argument, digits=40), scorer(argument, derivative=1, digits=40)
def local_enclosures(left, right, digits=40):
midpoint = (left + right) / 2
half_width = ball((right - left) / 2)
if max(abs(left), abs(right)) > 16:
raise ValueError('root-count certificate covers only |x| <= 16')
if digits == 40:
function, velocity = scorer_point(midpoint, ctx.prec)
else:
function = scorer(midpoint, digits=digits)
velocity = scorer(midpoint, derivative=1, digits=digits)
bound = (abs(function) + abs(velocity) + half_width / arb.pi())
bound *= (16 * half_width).exp()
function_box = widen(function, half_width * bound)
velocity_box = widen(velocity, half_width * (
ball(max(abs(left), abs(right))) * bound + 1 / arb.pi()))
acceleration = widen(ball(midpoint), half_width) * function_box - 1 / arb.pi()
return velocity_box, acceleration
def classify(left, right, depth=0):
velocity, acceleration = local_enclosures(left, right)
if sign(velocity):
return []
if sign(acceleration):
left_sign = sign(scorer_point(left, ctx.prec)[1])
right_sign = sign(scorer_point(right, ctx.prec)[1])
if left_sign and right_sign:
return [(left, right)] if left_sign != right_sign else []
if depth >= 24:
raise ArithmeticError('root count could not resolve a cell')
midpoint = (left + right) / 2
return classify(left, midpoint, depth + 1) + classify(midpoint, right, depth + 1)
@lru_cache(maxsize=1)
def critical_cells():
with precision(working_precision(50)):
positive = []
for index in range(64):
positive.extend(classify(Fraction(index, 64), Fraction(index + 1, 64)))
negative = []
for index in range(1024):
negative.extend(reversed(classify(-Fraction(index + 1, 64),
-Fraction(index, 64))))
if len(negative) == 12:
break
if len(positive) != 1 or len(negative) != 12:
raise ArithmeticError('unexpected number of critical points')
return positive + negative
def mp_fraction(value):
return mp.mpf(value.numerator) / value.denominator
def critical_point(index, digits):
if not 0 <= index <= 12:
raise ValueError('certified Gi critical indices are zero through twelve')
left, right = critical_cells()[index]
places = digits + 80
with mp.workdps(digits + 120):
candidate = mp.findroot(lambda argument: mp.diff(mp.scorergi, argument),
(mp_fraction(left), mp_fraction(right)))
center = Fraction(int(mp.nint(candidate * 10 ** places)), 10 ** places)
epsilon = Fraction(1, 10 ** places)
lower, upper = center - epsilon, center + epsilon
if not left < lower < upper < right:
raise ArithmeticError('candidate escaped its counted cell')
if sign(scorer(lower, derivative=1, digits=digits + 20)) * sign(
scorer(upper, derivative=1, digits=digits + 20)) != -1:
raise ArithmeticError('critical point is not bracketed')
if not sign(local_enclosures(lower, upper, digits=digits + 20)[1]):
raise ArithmeticError('critical point is not isolated')
return widen(ball(center), ball(epsilon))
def argument_for(identity, digits=100):
prefixes = {'ai-prime-zero-': (arb.airy_ai_zero, 1),
'bi-prime-zero-': (arb.airy_bi_zero, 1),
'ai-zero-': (arb.airy_ai_zero, 0),
'bi-zero-': (arb.airy_bi_zero, 0)}
for prefix, (function, derivative) in prefixes.items():
if identity.startswith(prefix):
index = int(identity[len(prefix):])
if index < 1:
raise ValueError('Airy indices start at one')
return function(index, derivative=derivative)
if identity.startswith('gi-critical-'):
return critical_point(int(identity[len('gi-critical-'):]), digits)
match = re.fullmatch(r'(-?\d+)/cbrt\(3\)', identity)
if match:
return arb(int(match[1])) / arb(3).root(3)
match = re.fullmatch(r'(\d+)\*sqrt\(2/pi\)', identity)
if match:
return arb(int(match[1])) * (2 / arb.pi()).sqrt()
if identity.startswith('sqrt(') and identity.endswith(')'):
return ball(Fraction(identity[5:-1])).sqrt()
if identity.endswith('*pi/2'):
return arb(int(identity[:-5])) * arb.pi() / 2
if identity.endswith('*pi'):
return arb(int(identity[:-3])) * arb.pi()
if identity.endswith('/sqrt(2)'):
return arb(int(identity[:-8])) / arb(2).sqrt()
return ball(Fraction(identity))
def evaluate(kind, argument, digits=100):
if kind in ('gi', 'hi'):
return scorer(argument, kind, digits=digits)
airy_indices = {'ai': 0, 'ai_prime': 1, 'bi': 2, 'bi_prime': 3}
if kind in airy_indices:
return argument.airy()[airy_indices[kind]]
return getattr(argument, kind)()
def certified_value(kind, identity, digits):
for attempt in range(4):
with precision(working_precision(digits) * 2 ** attempt):
argument = argument_for(identity, digits)
result = evaluate(kind, argument, digits)
if not result.is_finite():
raise ArithmeticError('nonfinite result')
if result.is_exact():
if result == 0:
return 0
if result == 1:
return 1
raise ArithmeticError('unexpected exact value')
if result.rel_accuracy_bits() >= math.ceil((digits + 10) * math.log2(10)):
return as_interval(result)
raise ArithmeticError('insufficient certified accuracy after escalation')
class BalancedValues(numberdb.Generator):
parameters = ('x',)
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'balanced_values.py', 'balanced_checks.py',
'balanced-selection.json', 'balanced-selection.md',
'balanced-numbers.json', 'balanced-checks-evidence.json')
def __init__(self, manifest=None):
path = Path(manifest) if manifest else Path(__file__).with_name('balanced-selection.json')
self.selection = json.loads(path.read_text())
self.table = self.selection['tid']
self.entries = self.selection['entries']
def enumerate(self):
for identity in self.entries:
yield {'x': identity}
def digits_for(self, params):
return self.entries[str(params['x'])]['digits']
def value(self, params, digits):
identity = str(params['x'])
metadata = dict(self.entries[identity]['metadata'])
metadata['number'] = certified_value(KINDS[self.table], identity, digits)
return metadata