back to table · edit · history · where entries came from · files · download
14289 bytes, as of the version from 2026-10-27 01:12 (current). Recorded here, not run.
"""Independent checks: python balanced_checks.py balanced-selection.json.
Every retained value is recomputed in narrow Arb arithmetic. Every addition
is also evaluated independently with mpmath at 160 and 210 decimal digits.
Elementary, Airy, elliptic, integral and Wronskian checks test normalization.
All comparisons check finite values and meaningful widths first.
"""
import argparse
import hashlib
import json
from fractions import Fraction
from pathlib import Path
import mpmath as mp
from flint import acb, arb
from balanced_generate import (
BalancedValues, as_interval, contains, evaluate, precision, written_interval,
)
def mp_rational(text):
value = Fraction(text)
return mp.mpf(value.numerator) / value.denominator
def mp_argument(row, generator):
descriptor = row["argument"]
kind = descriptor["type"]
if kind == "rational":
return mp_rational(descriptor["value"])
if kind == "airy":
value = mp_rational(descriptor["value"])
return 2 * value * mp.sqrt(value) / 3
if kind == "scaled_sqrt":
return mp_rational(descriptor["factor"]) * mp.sqrt(mp_rational(descriptor["radicand"]))
if kind == "cosh":
return mp.cosh(mp_rational(descriptor["value"]))
if kind == "turning":
order = int(row["params"]["ell"])
coupling = mp_rational(row["params"]["eta"])
return coupling + mp.sqrt(coupling**2 + order * (order + 1))
with precision(1000):
argument, unused = generator.argument(row)
interval = as_interval(argument)
return mp_rational((interval.lower + interval.upper) / 2)
def independent_value(kind, params, argument):
if kind.startswith("coulomb"):
return getattr(mp, kind)(int(params["ell"]), mp_rational(params["eta"]), argument)
order = mp_rational(params["nu"])
if kind == "legenq":
result = mp.legenq(order, 0, argument, type=3)
if abs(result.imag) > mp.mpf("1e-140"):
raise ArithmeticError("Legendre type-3 result is not real")
return result.real
if argument == 0 and kind == "kei" and order == 0:
return -mp.pi / 4
return getattr(mp, kind)(order, argument)
def assert_residual(value, name):
if not value.is_finite() or not value.contains(0) or value.rad() >= arb("1e-160"):
raise ArithmeticError("Invalid structural check: " + name)
def half_bessel(kind, order, argument):
degree = int(Fraction(order) - Fraction(1, 2))
if kind == "besselj":
previous = argument.sin() / argument
current = argument.sin() / argument**2 - argument.cos() / argument
else:
previous = -argument.cos() / argument
current = -argument.cos() / argument**2 - argument.sin() / argument
if degree == 0:
current = previous
else:
for index in range(1, degree):
previous, current = current, (2 * index + 1) * current / argument - previous
return (2 * argument / arb.pi()).sqrt() * current
def elliptic_q(degree, argument):
parameter = 2 / (argument + 1)
first = (2 / (argument + 1)).sqrt() * acb(parameter).elliptic_k().real
second = argument * first - (2 * (argument + 1)).sqrt() * acb(parameter).elliptic_e().real
target = Fraction(degree)
if target == Fraction(-1, 2):
return first
order = Fraction(1, 2)
while order < target:
coefficient = arb(str(order))
first, second = second, ((2 * coefficient + 1) * argument * second - coefficient * first) / (coefficient + 1)
order += 1
return second
def structural_checks(generator, row, argument, value, text):
kind = generator.kind
params = row["params"]
checks = []
if argument == 0:
return ["exact_or_continuous_origin"]
order = Fraction(params.get("nu", "0"))
if kind in ("besselj", "bessely"):
first = argument.bessel_j(arb(str(order)))
second = argument.bessel_y(arb(str(order)))
first_next = argument.bessel_j(arb(str(order + 1)))
second_next = argument.bessel_y(arb(str(order + 1)))
assert_residual(first_next * second - first * second_next - 2 / (arb.pi() * argument), "Bessel Wronskian")
checks.append("bessel_wronskian")
if order.denominator == 2:
if not contains(text, half_bessel(kind, order, argument)):
raise ArithmeticError("Elementary spherical recurrence failed")
checks.append("elementary_half_order_recurrence")
if kind in ("besseli", "besselk") and row["argument"]["type"] == "airy":
airy_argument = arb(row["argument"]["value"])
airy_ai, airy_prime, airy_bi, airy_bi_prime = airy_argument.airy()
if kind == "besselk":
expected = (arb.pi() * (3 / airy_argument).sqrt() * airy_ai if order == Fraction(1, 3)
else -arb.pi() * arb(3).sqrt() * airy_prime / airy_argument)
else:
expected = ((arb(3).sqrt() * airy_bi - 3 * airy_ai) / (2 * airy_argument.sqrt())
if order == Fraction(1, 3) else
(arb(3).sqrt() * airy_bi_prime + 3 * airy_prime) / (2 * airy_argument))
if not contains(text, expected):
raise ArithmeticError("Independent Airy reduction failed")
checks.append("independent_airy_function_reduction")
if kind in ("ber", "bei", "ker", "kei") and order == Fraction(1, 2):
rotated = argument * acb(-1 if kind in ("ber", "bei") else 1, 1) / arb(2).sqrt()
if kind in ("ber", "bei"):
elementary = (2 / (arb.pi() * rotated)).sqrt() * rotated.sin()
else:
elementary = acb(0, -arb.pi() / 4).exp() * (arb.pi() / (2 * rotated)).sqrt() * (-rotated).exp()
component = elementary.real if kind in ("ber", "ker") else elementary.imag
if not contains(text, component):
raise ArithmeticError("Kelvin half-order normalization failed")
checks.append("elementary_half_order_rotation_and_phase")
if kind in ("struveh", "struvel") and order == Fraction(1, 2):
expected = (2 / (arb.pi() * argument)).sqrt() * (
1 - argument.cos() if kind == "struveh" else argument.cosh() - 1)
if not contains(text, expected):
raise ArithmeticError("Struve half-order formula failed")
checks.append("elementary_half_order_formula")
if kind == "legenq" and order.denominator == 2:
if not contains(text, elliptic_q(order, argument)):
raise ArithmeticError("Legendre elliptic-integral reduction failed")
checks.append("elliptic_integral_seed_and_degree_recurrence")
if kind.startswith("coulomb"):
angular = int(params["ell"])
coupling = arb(params["eta"])
regular = acb(argument).coulomb_f(angular, coupling).real
irregular = acb(argument).coulomb_g(angular, coupling).real
factor = (1 + coupling**2 / (angular + 1)**2).sqrt()
shifted = (angular + 1) / argument + coupling / (angular + 1)
regular_prime = shifted * regular - factor * acb(argument).coulomb_f(angular + 1, coupling).real
irregular_prime = shifted * irregular - factor * acb(argument).coulomb_g(angular + 1, coupling).real
assert_residual(irregular * regular_prime - regular * irregular_prime - 1, "Coulomb Wronskian")
checks.append("coulomb_wronskian_equals_one")
if coupling == 0:
expected = half_bessel("besselj" if kind == "coulombf" else "bessely",
Fraction(2 * angular + 1, 2), argument)
expected *= (arb.pi() * argument / 2).sqrt()
if kind == "coulombg":
expected = -expected
if not contains(text, expected):
raise ArithmeticError("Free-wave elementary recurrence failed")
checks.append("free_wave_elementary_recurrence")
if row["argument"]["type"] == "turning":
assert_residual(1 - 2 * coupling / argument - angular * (angular + 1) / argument**2, "turning radius")
checks.append("turning_coefficient_zero")
return checks
def integral_check(kind, params, argument):
order = mp_rational(params["nu"])
if kind in ("struveh", "struvel"):
trig = mp.sin if kind == "struveh" else mp.sinh
factor = 2 * (argument / 2)**order / (mp.sqrt(mp.pi) * mp.gamma(order + mp.mpf("0.5")))
return factor * mp.quad(lambda angle: trig(argument * mp.cos(angle)) * mp.sin(angle)**(2 * order),
[0, mp.pi / 4, mp.pi / 2])
if kind in ("ber", "bei"):
rotated = argument * mp.exp(3j * mp.pi / 4)
value = mp.quad(lambda angle: mp.cos(order * angle - rotated * mp.sin(angle)),
[0, mp.pi / 2, mp.pi]) / mp.pi
return value.real if kind == "ber" else value.imag
rotated = argument * mp.exp(1j * mp.pi / 4)
value = mp.exp(-1j * mp.pi * order / 2) * mp.quad(
lambda variable: mp.exp(-rotated * mp.cosh(variable)) * mp.cosh(order * variable),
[0, 1, 2, 4, 6, 8])
return value.real if kind == "ker" else value.imag
def controls():
with precision(1000), mp.workdps(150):
pi_text = "3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068"
if not contains(pi_text, arb.pi()):
raise AssertionError("positive decimal control")
if contains(pi_text[:-3] + "168", arb.pi()):
raise AssertionError("displaced value accepted")
if contains(pi_text, arb("nan")) or contains(pi_text, arb(3, 10)):
raise AssertionError("vacuous enclosure accepted")
if abs(mp.quad(lambda variable: variable, [0, 1]) - mp.mpf("0.5")) > mp.mpf("1e-140"):
raise AssertionError("nonzero integration control")
assert_residual(evaluate("legenq", {"nu": "0"}, arb(2)) - arb(3).log() / 2, "Q0 logarithm")
assert_residual(evaluate("coulombf", {"ell": "0", "eta": "0"}, arb(1)) - arb(1).sin(), "free s-wave F")
assert_residual(evaluate("coulombg", {"ell": "0", "eta": "0"}, arb(1)) - arb(1).cos(), "free s-wave G")
return ["pi", "displaced_value_rejected", "nan_rejected", "wide_ball_rejected",
"integral_t_equals_half", "Q0_logarithm", "F0_free_sine", "G0_free_cosine"]
def check(selection_path):
generator = BalancedValues(selection_path)
records = []
control_results = controls()
for row in generator.selection["rows"]:
entry = row["entry"]
text = entry["number"] if isinstance(entry, dict) else entry
digits = max(100, int(row.get("digits", 100)))
value, certificate = generator.ball(row, digits, bits=max(1000, 4 * digits + 500))
if not contains(text, value):
raise ArithmeticError("Written interval does not enclose recomputation: " + str(row["params"]))
record = {"params": row["params"], "origin": row["origin"],
"stored_interval_contains_ball": True,
"accuracy_bits": None if isinstance(value, int) else int(value.rel_accuracy_bits())}
if certificate:
record["root_certificate"] = certificate
if row["origin"] == "balanced expansion":
lower, upper = written_interval(text)
for digits in (160, 210):
with mp.workdps(digits):
argument = mp_argument(row, generator)
independent = independent_value(generator.kind, row["params"], argument)
if not mp.isfinite(independent) or not mp_rational(lower) <= independent <= mp_rational(upper):
raise ArithmeticError("Independent mpmath mismatch: " + str(row["params"]))
record["mpmath_%d_digits" % digits] = True
with precision(1000):
argument, unused = generator.argument(row)
record["structural_checks"] = structural_checks(generator, row, argument, value, text)
kind = generator.kind
check_integral = kind in ("struveh", "struvel") or (
kind in ("ber", "bei", "ker", "kei") and row["params"]["nu"] in ("0", "1")
and row["params"]["x"] in ("1", "2"))
if check_integral and row["params"].get("x") != "0":
with mp.workdps(150):
argument = mp_argument(row, generator)
integral = integral_check(kind, row["params"], argument)
independent = independent_value(kind, row["params"], argument)
relative = abs((integral - independent) / independent)
if not mp.isfinite(relative) or relative >= mp.mpf("1e-115"):
raise ArithmeticError("Independent integral mismatch: " + str(row["params"]))
record["integral_relative_error"] = mp.nstr(relative, 12)
records.append(record)
if len(records) % 50 == 0:
print(generator.table, len(records), "checked", flush=True)
return {"tid": generator.table, "all_passed": True, "controls": control_results,
"source_sha256": {name: hashlib.sha256(Path(__file__).with_name(name).read_bytes()).hexdigest()
for name in ("balanced_checks.py", "balanced_generate.py")},
"selection_sha256": hashlib.sha256(Path(selection_path).read_bytes()).hexdigest(),
"count": len(records), "independently_checked_additions": sum(
row["origin"] == "balanced expansion" for row in records),
"minimum_accuracy_bits": min(row["accuracy_bits"] for row in records if row["accuracy_bits"] is not None),
"records": records,
"limitations": "Retained values are fully Arb-recomputed; new values additionally have independent mpmath and family-identity checks. Existing manual numbers and annotations are not claimed to originate in this run."}
if __name__ == "__main__":
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("selection", type=Path)
options = parser.parse_args()
evidence = check(options.selection)
options.selection.with_name("balanced-checks.json").write_text(json.dumps(evidence, indent=2) + "\n")
print(evidence["tid"], "passed", evidence["count"], flush=True)