back to table · edit · history · where entries came from · files · download
16451 bytes, as of the version from 2026-10-27 01:05. Recorded here, not run.
"""Reproduce the complete attached elliptic selection in ball arithmetic.
python -m pip install numberdb python-flint mpmath
python generate.py --document table.json --output checks.json
Download selection.json beside this file. It lists exactly the numerical
entry identities, including retained rational references. The generator
never writes the website. Existing comments and referral-only branches are
preserved in table.json, not invented by this numerical implementation.
Arb evaluates rational arguments, pi, square roots and quarter periods as
balls at at least 640 bits. A symbolic K(m)/d argument denotes the complete
first-kind integral at the row's exact m divided by d, never a rounded
decimal. Its uncertainty propagates through the theta quotient. The ball
is converted to NumberDB's RealInterval through exact binary midpoint and
radius fractions; no decimal rounding is introduced in the enclosure.
Existing rational selections are retained rather than extended as grids.
New negative-parameter references use imaginary moduli i/2, i, 2i; new
third-kind slices use n=-sqrt(m) or n=m<0. Incomplete-integral additions
cover m=-1 and elementary m=1 degenerations. Jacobi sn additions are real
division points. Negative rational nomes use theta4(-q)=theta3(q).
These finite reference extents are choices, not literature claims about
every point. See DLMF 19.6, 19.7, 20.2, 22.5, 23.15 and selection.json.
The check command independently evaluates every selected number with
mpmath at two precisions. New integral values are also checked by direct
quadrature; division points are checked by inversion of the defining
integral and the half-period radical. Measurements reject nonfinite or
wide balls and test the full printed interval, not just its midpoint.
"""
import argparse
import json
from contextlib import contextmanager
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import mpmath as mp
from flint import acb, arb, ctx
import numberdb
from numberdb._compare import _claim
from numberdb._write import to_text
PARAMETERS = {"T25": ("m",), "T26": ("m",), "T59": ("n", "m"),
"T349": ("phi", "m"), "T350": ("phi", "m"),
"T424": ("m", "u"), "T425": ("m",), "T426": ("q",),
**{tid: ("r",) for tid in ("T296", "T298", "T299", "T504",
"T505", "T506", "T507")}}
@contextmanager
def working(bits=640):
previous = ctx.prec
ctx.prec = bits
try:
yield
finally:
ctx.prec = previous
def rational(text):
value = Fraction(str(text))
return arb(value.numerator) / value.denominator
def mp_rational(text):
value = Fraction(str(text))
return mp.mpf(value.numerator) / value.denominator
def finite_real(value, minimum_bits=400):
if isinstance(value, acb):
if not value.is_finite() or not value.imag.contains(0):
raise ArithmeticError("Nonfinite or nonreal complex enclosure")
value = value.real
if not value.is_finite() or value.rel_accuracy_bits() < minimum_bits:
raise ArithmeticError("Nonfinite or insufficiently accurate real enclosure")
return value
def interval(value):
def exact(part):
mantissa, exponent = part.man_exp()
return Fraction(int(mantissa)) * Fraction(2) ** int(exponent)
midpoint, radius = exact(value.mid()), exact(value.rad())
return numberdb.RealInterval(midpoint - radius, midpoint + radius)
def division(argument):
if argument.startswith("K(m)/"):
denominator = int(argument.partition("/")[2])
if denominator not in (2, 3, 4):
raise ValueError("Unvalidated division-point denominator")
return denominator
return None
def rogers_product(nome):
count = 120
product = nome ** (arb(1) / 5)
for index in range(count):
product *= (1 - nome ** (5 * index + 1)) * (1 - nome ** (5 * index + 4))
product /= (1 - nome ** (5 * index + 2)) * (1 - nome ** (5 * index + 3))
first = 5 * count + 1
bound = nome ** first * (1 + nome + nome ** 2 + nome ** 3)
bound /= (1 - nome ** 5) * (1 - nome ** first)
return product + arb(0, (abs(product) * (bound.exp() - 1)).upper())
def ball_value(tid, params):
if tid == "T26" and params["m"] == "1":
return 1
if tid == "T350" and params["m"] == "1":
if params["phi"] == "1/6":
return Fraction(1, 2)
if params["phi"] == "1/2":
return 1
if tid == "T426" and params["q"] == "0":
return 1
if tid == "T299" and params["r"] == "1":
return Fraction(1, 2)
if tid == "T25":
return finite_real(acb(rational(params["m"])).elliptic_k())
if tid == "T26":
return finite_real(acb(rational(params["m"])).elliptic_e())
if tid == "T59":
return finite_real(acb.elliptic_pi(rational(params["n"]), rational(params["m"])))
if tid in ("T349", "T350"):
function = acb.elliptic_f if tid == "T349" else acb.elliptic_e_inc
return finite_real(function(rational(params["phi"]), rational(params["m"]), pi=True))
if tid in ("T424", "T425"):
parameter = acb(rational(params["m"]))
period = parameter.elliptic_k()
tau = acb(0, 1) * (1 - parameter).elliptic_k() / period
if tid == "T425":
return finite_real((acb(0, arb.pi()) * tau).exp())
denominator = division(params["u"])
argument = period / denominator if denominator else acb(rational(params["u"]))
constants = acb(0).modular_theta(tau)
shifted = (argument / (2 * period)).modular_theta(tau)
return finite_real(constants[2] / constants[1] * shifted[0] / shifted[3])
if tid == "T426":
nome = rational(params["q"])
tau = acb(0, -abs(nome).log() / arb.pi())
return finite_real(acb(0).modular_theta(tau)[2 if nome < 0 else 3])
root = rational(params["r"]).sqrt()
tau = acb(0, root)
nome = (-arb.pi() * root).exp()
if tid == "T298":
return finite_real(rogers_product(nome))
if tid in ("T296", "T299"):
parameter = finite_real(tau.modular_lambda())
first = finite_real(acb(parameter).elliptic_k())
if tid == "T296":
return first
second = finite_real(acb(parameter).elliptic_e())
return finite_real(arb.pi() / (4 * first ** 2) + root * (1 - second / first))
if tid == "T506":
return finite_real(tau.modular_eta())
return finite_real(acb(0).modular_theta(tau)[{"T504": 2, "T505": 3, "T507": 1}[tid]])
def mp_value(tid, params, digits):
with mp.workdps(digits):
if tid in ("T25", "T26"):
return +(mp.ellipk if tid == "T25" else mp.ellipe)(mp_rational(params["m"]))
if tid == "T59":
return +mp.ellippi(mp_rational(params["n"]), mp_rational(params["m"]))
if tid in ("T349", "T350"):
return +(mp.ellipf if tid == "T349" else mp.ellipe)(
mp.pi * mp_rational(params["phi"]), mp_rational(params["m"]))
if tid == "T424":
parameter = mp_rational(params["m"])
denominator = division(params["u"])
argument = mp.ellipk(parameter) / denominator if denominator else mp_rational(params["u"])
return +mp.ellipfun("sn", argument, parameter)
if tid == "T425":
parameter = mp_rational(params["m"])
return +mp.exp(-mp.pi * mp.ellipk(1 - parameter) / mp.ellipk(parameter))
if tid == "T426":
return +mp.jtheta(4, 0, mp_rational(params["q"]))
root = mp.sqrt(mp_rational(params["r"]))
nome = mp.exp(-mp.pi * root)
if tid == "T298":
denominator = mp.mpf(1)
for index in range(200, 0, -1):
denominator = 1 + nome ** index / denominator
return +(nome ** (mp.mpf(1) / 5) / denominator)
if tid == "T506":
return +(nome ** (mp.mpf(1) / 12) * mp.qp(nome ** 2))
if tid in ("T504", "T505", "T507"):
return +mp.jtheta({"T504": 3, "T505": 4, "T507": 2}[tid], 0, nome)
parameter = (mp.jtheta(2, 0, nome) / mp.jtheta(3, 0, nome)) ** 4
first = mp.ellipk(parameter)
if tid == "T296":
return +first
return +(mp.pi / (4 * first ** 2) + root * (1 - mp.ellipe(parameter) / first))
def within(value, written):
center, radius = _claim(str(written))
if isinstance(value, (int, Fraction)):
return abs(Fraction(value) - Fraction(center)) <= Fraction(radius)
if not radius:
return value.contains(rational(str(center)))
return abs(value - rational(str(center))) < rational(str(radius))
def checks(tid, params, written, new=False):
first, second = (mp_value(tid, params, digits) for digits in (160, 210))
center, radius = _claim(str(written))
with mp.workdps(220):
if not all(mp.isfinite(value) for value in (first, second)):
raise ArithmeticError("Nonfinite independent value")
if abs(first - second) > mp.mpf("1e-145") * max(abs(second), mp.mpf("1e-50")):
raise ArithmeticError("Independent working precisions disagree")
tolerance = mp.mpf(str(radius)) if radius else mp.mpf("1e-145")
discrepancy = abs(second - mp.mpf(str(center)))
if discrepancy > tolerance:
raise ArithmeticError("Independent value outside printed interval")
evidence = {"mp_error_ulps": mp.nstr(discrepancy / tolerance, 16)}
if new:
evidence.update(structural_check(tid, params, written))
return evidence
def structural_check(tid, params, written):
value = ball_value(tid, params)
with mp.workdps(160):
target = mp_value(tid, params, 160)
if tid in ("T25", "T26", "T59", "T349", "T350"):
parameter = mp_rational(params["m"])
endpoint = mp.pi * mp_rational(params["phi"]) if "phi" in params else mp.pi / 2
if tid in ("T26", "T350"):
integrand = lambda angle: mp.sqrt(1 - parameter * mp.sin(angle) ** 2)
elif tid == "T59":
characteristic = mp_rational(params["n"])
integrand = lambda angle: 1 / ((1 - characteristic * mp.sin(angle) ** 2)
* mp.sqrt(1 - parameter * mp.sin(angle) ** 2))
else:
integrand = lambda angle: 1 / mp.sqrt(1 - parameter * mp.sin(angle) ** 2)
quadrature = mp.quad(integrand, [0, endpoint / 2, endpoint])
if abs(quadrature - target) > mp.mpf("1e-140") * max(1, abs(target)):
raise ArithmeticError("Defining-integral quadrature failed")
evidence = {"defining_integral_quadrature": True}
if tid in ("T25", "T26"):
parameter = rational(params["m"])
transformed = parameter / (parameter - 1)
expected = (acb(transformed).elliptic_k() / (1 - parameter).sqrt()
if tid == "T25" else acb(transformed).elliptic_e() * (1 - parameter).sqrt())
if not value.overlaps(finite_real(expected)):
raise ArithmeticError("Imaginary-modulus transformation failed")
evidence["imaginary_modulus_transformation"] = True
if tid == "T59":
characteristic, parameter = rational(params["n"]), rational(params["m"])
if params["n"] == params["m"]:
expected = acb(parameter).elliptic_e() / (1 - parameter)
else:
if not (characteristic ** 2 - parameter).contains(0):
raise ArithmeticError("Not a selected reduction slice")
expected = arb.pi() / (4 * (1 - characteristic)) + acb(parameter).elliptic_k() / 2
if not value.overlaps(finite_real(expected)):
raise ArithmeticError("Third-kind reduction failed")
evidence["third_kind_reduction"] = True
return evidence
if tid == "T424":
denominator = division(params["u"])
parameter = mp_rational(params["m"])
recovered = mp.quad(lambda angle: 1 / mp.sqrt(1 - parameter * mp.sin(angle) ** 2),
[0, mp.asin(target)])
if abs(recovered - mp.ellipk(parameter) / denominator) > mp.mpf("1e-140"):
raise ArithmeticError("Inverse defining integral failed")
if denominator == 2:
radical = 1 / (1 + (1 - rational(params["m"])).sqrt()).sqrt()
if not value.overlaps(finite_real(radical)):
raise ArithmeticError("Half-quarter-period radical failed")
return {"inverse_defining_integral": True, "half_quarter_radical": denominator == 2}
if tid == "T426":
nome = mp_rational(params["q"])
series = 1 + 2 * mp.fsum((-1) ** index * nome ** (index * index) for index in range(1, 120))
if abs(series - target) > mp.mpf("1e-140"):
raise ArithmeticError("Independent defining theta series failed")
return {"defining_series": True, "series_tail_below": "1e-400"}
raise ValueError("Missing structural check")
def controls():
with mp.workdps(160):
integral = mp.quad(lambda angle: mp.cos(angle), [0, mp.pi / 2])
if abs(integral - 1) > mp.mpf("1e-145"):
raise ArithmeticError("Quadrature control failed")
if ball_value("T26", {"m": "1"}) != 1 or ball_value("T426", {"q": "0"}) != 1:
raise ArithmeticError("Exact endpoint controls failed")
value = ball_value("T25", {"m": "0"})
if not value.overlaps(finite_real(arb.pi() / 2)):
raise ArithmeticError("Parameter-convention control failed")
if within(value, "3.14159"):
raise ArithmeticError("Wrong-value control accepted")
for invalid in (arb("nan"), arb(1, 2)):
try:
finite_real(invalid)
except ArithmeticError:
pass
else:
raise ArithmeticError("Nonfinite/wide-ball control accepted")
class EllipticSelection(numberdb.Generator):
type = "R"
digits = 100
rigour = "proven"
def __init__(self, path=None):
path = path or Path(__file__).with_name("selection.json")
self.selection = json.loads(Path(path).read_text())
self.table = self.selection["tid"]
self.parameters = PARAMETERS[self.table]
def enumerate(self):
for row in self.selection["rows"]:
yield row["parameters"]
def value(self, params, digits):
with working(max(640, numberdb.bits(digits, losing=160))):
value = ball_value(self.table, params)
return interval(value) if isinstance(value, arb) else value
if __name__ == "__main__":
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--selection", type=Path, default=Path(__file__).with_name("selection.json"))
parser.add_argument("--document", type=Path)
parser.add_argument("--output", type=Path, default=Path("checks.json"))
arguments = parser.parse_args()
generator = EllipticSelection(arguments.selection)
document = json.loads(arguments.document.read_text()) if arguments.document else None
results = []
with working(640):
controls()
for selected in generator.selection["rows"]:
parameters = selected["parameters"]
value = generator.value(parameters, selected.get("digits", 100))
written = to_text(value, selected.get("digits", 100))
if document is not None:
leaf = document["Numbers"]
for name in generator.parameters:
leaf = leaf[parameters[name]]
written = leaf["number"] if isinstance(leaf, dict) else leaf
if not within(ball_value(generator.table, parameters), written):
raise ArithmeticError("Existing value failed full enclosure check")
results.append({"parameters": parameters, "number": written,
**checks(generator.table, parameters, written, selected.get("new", False))})
arguments.output.write_text(json.dumps({"tid": generator.table, "controls_passed": True,
"rows": results}, indent=2) + "\n")
print(generator.table, len(results), "values checked", flush=True)