back to table · edit · history · where entries came from · files · download
10669 bytes, as of the version from 2026-10-27 01:12 (current). Recorded here, not run.
"""Certified Kelvin/Struve extrema and Coulomb turning-point values.
Install: python -m pip install numberdb python-flint mpmath
Use the attached generate.py to verify its table without writing.
Critical arguments are defined by an isolating rational interval, not a
rounded decimal or a claimed global root ordinal. Opposite endpoint signs
give existence. Gronwall bounds and the differential equation prove that the
second derivative has a fixed nonzero sign throughout the isolating interval,
giving uniqueness and a genuine extremum. No completeness claim is made
outside these intervals. mpmath proposes a narrow bracket but never certifies
it. Arb certifies both its endpoints and the propagated argument uncertainty.
"""
import json
import math
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 acb, arb, ctx
KINDS = {'T490': 'ber', 'T491': 'bei', 'T492': 'struveh', 'T493': 'ker',
'T494': 'struvel', 'T509': 'kei', 'T510': 'coulombf', 'T511': 'coulombg'}
GUARD_BITS = 440
@contextmanager
def precision(bits):
previous = ctx.prec
ctx.prec = bits
try:
yield
finally:
ctx.prec = previous
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
def interval(lower, upper):
lower, upper = Fraction(lower), Fraction(upper)
return arb(str((lower + upper) / 2), str((upper - lower) / 2))
def struve(kind, order, argument):
order = arb(str(order))
half = arb('1/2')
scaled = acb(argument * argument / 4)
if kind == 'struveh':
scaled = -scaled
result = ((acb(argument) / 2) ** (order + 1)
* scaled.hypgeom([1], [3 * half, order + 3 * half])
/ ((3 * half).gamma() * (order + 3 * half).gamma()))
if not result.is_finite() or not result.imag.contains(0):
raise ArithmeticError('Struve result is not a finite real enclosure')
return result.real
def kelvin_pair(kind, argument):
rotation = acb(-1 if kind in ('ber', 'bei') else 1, 1) / arb(2).sqrt()
rotated = argument * rotation
if kind in ('ber', 'bei'):
value = rotated.bessel_j(0)
derivative = -rotation * rotated.bessel_j(1)
else:
value = rotated.bessel_k(0)
derivative = -rotation * rotated.bessel_k(1)
return value, derivative
def values_and_derivatives(kind, order, argument):
if kind in ('ber', 'bei', 'ker', 'kei'):
if int(order) != 0:
raise ValueError('This selection contains only order-zero Kelvin functions')
value, derivative = kelvin_pair(kind, argument)
second = acb(0, 1) * value - derivative / argument
component = 'real' if kind in ('ber', 'ker') else 'imag'
return tuple(getattr(item, component) for item in (value, derivative, second))
if kind != 'struveh':
raise ValueError('No stationary-point family implemented for ' + kind)
value = struve(kind, order, argument)
derivative = struve(kind, order - 1, argument) - order * value / argument
forcing = (argument / 2) ** (order - 1) / (
arb.pi().sqrt() * arb(str(Fraction(order) + Fraction(1, 2))).gamma())
second = forcing - derivative / argument - (1 - order * order / argument**2) * value
return value, derivative, second
def second_derivative_on_interval(kind, order, lower, upper):
domain = interval(lower, upper)
middle = arb(str((Fraction(lower) + Fraction(upper)) / 2))
half_width = arb(str((Fraction(upper) - Fraction(lower)) / 2))
left = arb(str(lower))
if kind in ('ber', 'bei', 'ker', 'kei'):
value, derivative = kelvin_pair(kind, middle)
components = (value.real, value.imag, derivative.real, derivative.imag)
norm = max(component.abs_upper() for component in components)
matrix_bound = 1 + 1 / left
variation = half_width * matrix_bound * norm * (matrix_bound * half_width).exp()
enclosed = [component + arb(0, variation.abs_upper()) for component in components]
if kind in ('ber', 'ker'):
return -enclosed[2] / domain - enclosed[1]
return -enclosed[3] / domain + enclosed[0]
value, derivative, unused = values_and_derivatives(kind, order, middle)
matrix_bound = 1 + 1 / left + order * order / left**2
forcing = (domain / 2) ** (order - 1) / (
arb.pi().sqrt() * arb(str(Fraction(order) + Fraction(1, 2))).gamma())
norm = max(value.abs_upper(), derivative.abs_upper())
variation = half_width * (matrix_bound * norm + forcing.abs_upper()) * (
matrix_bound * half_width).exp()
value += arb(0, variation.abs_upper())
derivative += arb(0, variation.abs_upper())
return forcing - derivative / domain - (1 - order * order / domain**2) * value
def mp_value(kind, order, argument):
return getattr(mp, kind)(order, argument)
def mp_derivative(kind, order, argument):
if kind in ('ber', 'bei', 'ker', 'kei'):
rotation = mp.exp((3 if kind in ('ber', 'bei') else 1) * mp.pi * 1j / 4)
bessel = mp.besselj if kind in ('ber', 'bei') else mp.besselk
derivative = -rotation * bessel(1, argument * rotation)
return derivative.real if kind in ('ber', 'ker') else derivative.imag
return mp.struveh(order - 1, argument) - order * mp.struveh(order, argument) / argument
@lru_cache(maxsize=256)
def critical_argument(kind, order, lower, upper, bits):
if bits != ctx.prec:
raise ArithmeticError('Root cache precision mismatch')
endpoint_signs = [sign(values_and_derivatives(kind, order, arb(endpoint))[1])
for endpoint in (lower, upper)]
if endpoint_signs[0] * endpoint_signs[1] != -1:
raise ArithmeticError('Root interval lacks opposite endpoint signs')
second = second_derivative_on_interval(kind, order, lower, upper)
if sign(second) == 0:
raise ArithmeticError('Root uniqueness not proved on isolating interval')
decimal_precision = math.floor(bits * math.log10(2))
with mp.workdps(70):
endpoints = [mp.mpf(str(Fraction(endpoint).numerator)) /
Fraction(endpoint).denominator for endpoint in (lower, upper)]
candidate = mp.findroot(lambda argument: mp_derivative(kind, order, argument), endpoints,
solver='anderson')
approximation = Fraction(mp.nstr(candidate, 65))
for iteration in range(5):
point = arb(str(approximation))
unused, derivative, curvature = values_and_derivatives(kind, order, point)
approximation = exact_endpoint((point - derivative / curvature).mid())
radius = Fraction(1, 10**(decimal_precision - 60))
narrow_lower = approximation - radius
narrow_upper = approximation + radius
if not Fraction(lower) < narrow_lower < narrow_upper < Fraction(upper):
raise ArithmeticError('Candidate escaped its defining interval')
narrow_signs = [sign(values_and_derivatives(kind, order, arb(str(endpoint)))[1])
for endpoint in (narrow_lower, narrow_upper)]
if narrow_signs != endpoint_signs:
raise ArithmeticError('Uncertified narrow bracket')
argument = interval(narrow_lower, narrow_upper)
return argument, {'lower': str(narrow_lower), 'upper': str(narrow_upper),
'endpoint_signs': narrow_signs,
'second_derivative_sign': sign(second),
'defining_interval': [lower, upper]}
def evaluate(kind, params, argument):
if kind in ('ber', 'bei', 'ker', 'kei'):
if argument == 0:
return {'ber': 1, 'bei': 0, 'kei': -arb.pi() / 4}[kind]
value, unused = kelvin_pair(kind, argument)
return value.real if kind in ('ber', 'ker') else value.imag
if kind in ('struveh', 'struvel'):
if argument == 0:
return 0
return struve(kind, int(params['nu']), argument)
result = getattr(acb(argument), 'coulomb_f' if kind == 'coulombf' else 'coulomb_g')(
int(params['ell']), acb(arb(params['eta'])))
if not result.is_finite() or not result.imag.contains(0):
raise ArithmeticError('Coulomb result is not a finite real enclosure')
return result.real
class IntrinsicValues(numberdb.Generator):
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'intrinsic_values.py', 'checks.py', 'selection.json',
'checks-evidence.json', 'selection.md')
def __init__(self, tid, selection_path=None):
self.table = tid
self.kind = KINDS[tid]
self.parameters = ('ell', 'eta', 'rho') if tid in ('T510', 'T511') else ('nu', 'x')
path = Path(selection_path) if selection_path else Path(__file__).with_name('selection.json')
self.selection = json.loads(path.read_text())
if self.selection['tid'] != tid:
raise ValueError('Selection belongs to another table')
def enumerate(self):
for row in self.selection['rows']:
yield row['params']
def argument(self, params):
row = next(row for row in self.selection['rows'] if row['params'] == params)
if 'bracket' in row:
return critical_argument(self.kind, int(params['nu']), *row['bracket'], ctx.prec)
if row['argument'] == 'turning':
coupling = arb(params['eta'])
order = int(params['ell'])
return coupling + (coupling**2 + order * (order + 1)).sqrt(), None
return arb(row['argument']), None
def value(self, params, digits):
with precision(math.ceil(digits * math.log2(10)) + GUARD_BITS):
argument, unused = self.argument(params)
result = evaluate(self.kind, params, argument)
if isinstance(result, int):
return result
if not result.is_finite() or result.rel_accuracy_bits() < math.ceil((digits + 15) * math.log2(10)):
raise ArithmeticError('Insufficient certified accuracy')
return as_interval(result)
def verify(tid):
generator = IntrinsicValues(tid)
report = generator.verify(sample=None)
print(report)
raise SystemExit(0 if report.ok else 1)