back to table · edit · history · where entries came from · files · download
11085 bytes, as of the version from 2026-10-27 01:02 (current). Recorded here, not run.
"""Reproduce the balanced remaining-function selections without API writes.
pip install numberdb python-flint mpmath
python generate.py
Download balanced_values.py and selection.json beside generate.py. The latter
records the complete current selection, including retained identifiers and row
annotations. No previous checkout or private helper is required.
All nonexact computations use Arb, with 400 guard bits beyond the requested
decimal precision. Every result must retain at least digits+12 decimal digits.
Exact rational arguments and certified root brackets propagate through the
evaluation. Special-function formulas follow DLMF 4.13, 6.2, 7.2 and 7.5.
The delay method extends the algorithm in the previously attached
remaining_values.py. A centered degree-N polynomial and a separate uniform
error describe each unit interval. Dividing P(y) by a+y gives coefficients
d_j=(p_j-d_(j-1))/a. The omitted derivative series on |y|<=1/2 is bounded by
|d_N| (1/2)^N/(1-1/(2a)); the inherited function error contributes E/(a-1/2).
Integration has length at most one, and the boundary error is added separately.
Ball coefficients enclose rounding error. The error is not differentiated as
though it were a constant coefficient. Degree 320 suffices for 100 digits in
this bounded selection; the independent check increases degree and precision.
Existing stored strings are preserved by the preparation process, not claimed
to have been newly authored by this generator. Verification recomputes them.
"""
import json
import math
import re
from contextlib import contextmanager
from fractions import Fraction
from functools import lru_cache
from pathlib import Path
import numberdb
from flint import arb, ctx
GUARD_BITS = 400
PARAMETERS = {
'T188': ('normalisation', 'x'), 'T191': ('x',), 'T194': ('x',),
'T195': ('x',), 'T200': ('branch', 'x'), 'T205': ('x',),
'T206': ('x',), 'T422': ('s',), 'T423': ('u',),
'T433': ('u',), 'T434': ('u',),
}
@contextmanager
def precision(bits):
previous = ctx.prec
ctx.prec = bits
try:
yield
finally:
ctx.prec = previous
def as_ball(value):
return value if isinstance(value, arb) else arb(str(value))
def 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(endpoint(value.lower()), endpoint(value.upper()))
def widen(value, error):
return value + arb(0, as_ball(error).abs_upper())
def dawson(argument):
return arb.pi().sqrt() * (-(argument * argument)).exp() * argument.erfi() / 2
def bracket_root(function, left, right, digits):
left, right = Fraction(left), Fraction(right)
left_value = function(as_ball(left))
sign = 1 if left_value > 0 else -1
if not left_value * sign > 0 or not function(as_ball(right)) * sign < 0:
raise ArithmeticError('root endpoints do not have strict opposite signs')
while right - left > Fraction(1, 10 ** (digits + 25)):
middle = (left + right) / 2
value = function(as_ball(middle)) * sign
if value > 0:
left = middle
elif value < 0:
right = middle
else:
raise ArithmeticError('undecidable root sign')
root = widen(as_ball((left + right) / 2), as_ball((right - left) / 2))
if not function(root).is_finite() or not function(root).contains(0):
raise ArithmeticError('root enclosure residual failed')
return root, left, right
@lru_cache(maxsize=16)
def dawson_root(digits, bits):
with precision(bits):
return bracket_root(lambda argument: 1 - 2 * argument * dawson(argument),
Fraction(1, 2), 1, digits)
def argument(identity, digits=100):
if identity == 'dawson-critical-1':
return dawson_root(digits, ctx.prec)[0]
if identity == 'e':
return arb(1).exp()
if identity == '-1/e':
return -arb(-1).exp()
match = re.fullmatch(r'(log|exp|sqrt)\((-?\d+)\)', identity)
if match:
operation, integer = match.groups()
return getattr(arb(integer), operation)()
match = re.fullmatch(r'-(\d+(?:/\d+)?)\*exp\(-\1\)', identity)
if match:
coefficient = as_ball(Fraction(match.group(1)))
return -coefficient * (-coefficient).exp()
match = re.fullmatch(r'-exp\(-1/(\d+)\)/(\d+)', identity)
if match and match.group(1) == match.group(2):
denominator = int(match.group(1))
return -(arb(-1) / denominator).exp() / denominator
return as_ball(Fraction(identity))
def polynomial(coefficients, offset):
value = arb(0)
for coefficient in reversed(coefficients):
value = value * offset + coefficient
return value
class Model:
def __init__(self, coefficients, error=0):
self.coefficients = coefficients
self.error = as_ball(error)
def evaluate(self, offset):
return widen(polynomial(self.coefficients, as_ball(offset)), self.error)
def constant_model(value, degree):
return Model([as_ball(value)] + [arb(0) for index in range(degree)])
def integrate_step(previous, boundary, midpoint, sign, degree):
half = arb(1) / 2
midpoint = as_ball(midpoint)
quotient = []
prior = arb(0)
for coefficient in previous.coefficients:
prior = (coefficient - prior) / midpoint
quotient.append(prior)
omitted = abs(quotient[degree]) * half ** degree / (1 - half / midpoint)
inherited = previous.error / (midpoint - half)
current = [arb(0)] + [sign * quotient[index - 1] / index
for index in range(1, degree + 1)]
current[0] = polynomial(boundary.coefficients, half) - polynomial(current, -half)
return Model(current, boundary.error + omitted + inherited)
@lru_cache(maxsize=16)
def delay_models(kind, degree, bits):
with precision(bits):
if kind == 'rho':
models = [constant_model(1, degree)]
for interval in range(1, 10):
models.append(integrate_step(models[-1], models[-1],
Fraction(2 * interval + 1, 2), -1, degree))
return models
if kind == 'buchstab':
models = [None, constant_model(1, degree)]
for interval in range(2, 8):
models.append(integrate_step(models[-1], models[-1],
Fraction(2 * interval - 1, 2), 1, degree))
return models
if kind == 'sieve':
upper, lower = [], []
for interval in range(10):
midpoint = Fraction(2 * interval - 1, 2)
upper.append(constant_model(2 * arb.const_euler().exp(), degree)
if interval <= 2 else integrate_step(
lower[interval - 1], upper[interval - 1], midpoint, 1, degree))
lower.append(constant_model(0, degree) if interval <= 1 else integrate_step(
upper[interval - 1], lower[interval - 1], midpoint, 1, degree))
return upper, lower
raise ValueError('unknown delay equation')
def delay_value(tid, identity, degree):
point = Fraction(identity)
interval = (point.numerator - 1) // point.denominator
offset = point - interval - Fraction(1, 2)
if tid == 'T433':
if point <= 1:
return 1
return delay_models('rho', degree, ctx.prec)[interval].evaluate(offset)
if tid == 'T434':
if point <= 2:
return 1 / point
return delay_models('buchstab', degree, ctx.prec)[interval].evaluate(offset) / as_ball(point)
return delay_models('sieve', degree, ctx.prec)[0][interval].evaluate(offset) / as_ball(point)
def compute(tid, params, digits=100, extra_bits=0, degree=None):
bits = numberdb.bits(digits, losing=GUARD_BITS + extra_bits)
degree = degree or max(320, 3 * digits + 20)
identity = str(params[PARAMETERS[tid][-1]])
exact_lambert = {('0', '0'): 0, ('0', 'e'): 1,
('0', '-1/e'): -1, ('-1', '-1/e'): -1}
if tid == 'T200' and (params['branch'], identity) in exact_lambert:
return exact_lambert[params['branch'], identity]
if tid in ('T194', 'T205', 'T206') and identity == '0':
return 0
if tid == 'T423' and identity == '1':
return 0
with precision(bits):
point = argument(identity, digits)
if tid == 'T188':
value = point.ei() if params['normalisation'] == 'Ei' else -(-point).ei()
elif tid == 'T191':
value = point.log().ei()
elif tid == 'T194':
value = point.shi()
elif tid == 'T195':
value = point.chi()
elif tid == 'T205':
value = point.erfi()
elif tid == 'T206':
value = dawson(point)
elif tid == 'T200':
branch = int(params['branch'])
value = point.lambertw(branch)
if branch not in (0, -1) or not (value > -1 if branch == 0 else value < -1):
raise ArithmeticError('wrong real Lambert branch')
residual = value * value.exp() - point
if not residual.is_finite() or not residual.contains(0):
raise ArithmeticError('Lambert defining equation failed')
if not residual.abs_upper() < arb(10) ** (-(digits + 10)):
raise ArithmeticError('uninformative Lambert residual')
elif tid == 'T423':
value = bracket_root(lambda root: root.exp() - 1 - point * root,
Fraction(1, 128), 8, digits)[0]
else:
value = delay_value(tid, identity, degree)
if isinstance(value, arb):
required = math.ceil((digits + 12) * math.log2(10))
if not value.is_finite() or value.rel_accuracy_bits() < required:
raise ArithmeticError('insufficient finite relative accuracy: ' + tid + ' ' + identity)
return value
class BalancedValues(numberdb.Generator):
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'balanced_values.py', 'selection.json',
'balanced_checks.py', 'coverage-checks.json', 'coverage-rationale.json')
def __init__(self, tid, selection=None):
self.table = tid
self.parameters = PARAMETERS[tid]
self.selection = selection or json.loads(Path(__file__).with_name('selection.json').read_text())
def enumerate(self):
for row in self.selection['rows']:
yield row['params']
def value(self, params, digits):
record = next(row for row in self.selection['rows'] if row['params'] == params)
value = compute(self.table, params, digits)
with precision(numberdb.bits(digits, losing=GUARD_BITS)):
result = dict(record.get('metadata', {}))
result['number'] = as_interval(value) if isinstance(value, arb) else value
return result