back to table · edit · history · where entries came from · files · download
9794 bytes, as of the version from 2026-10-27 01:00 (current). Recorded here, not run.
"""Independent checks of every selected entry; no website writes.
python generate.py --document table.yaml
Every decimal is tested as the interval implied by its final place, not as
an exact midpoint. Finite and width guards precede enclosure comparisons.
Separate mpmath computations at two precisions supplement Arb certificates.
"""
import argparse
import json
import math
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import mpmath as mp
import numberdb
from flint import arb
from numberdb._write import to_text
from balanced_values import BalancedValues, PARAMETERS, endpoint, evaluate, flatten, precision, rational
def bounds(text):
midpoint = Fraction(text)
if "." not in text and "e" not in text.lower():
return midpoint, midpoint, Fraction(0)
unit = Fraction(10) ** Decimal(text).as_tuple().exponent
return midpoint - unit, midpoint + unit, unit
def mp_fraction(value):
return mp.mpf(value.numerator) / value.denominator
def mp_argument(identity):
if identity == "1/sqrt(2)":
return 1 / mp.sqrt(2)
if identity.startswith("sqrt(") and identity.endswith(")"):
return mp.sqrt(mp_fraction(Fraction(identity[5:-1])))
return mp_fraction(Fraction(identity))
def independent(tid, params):
values = {name: mp_argument(value) for name, value in params.items()}
if tid == "T51":
return values["a"] * mp.pi / (2 * mp.ellipk(1 - (values["b"] / values["a"]) ** 2))
if tid == "T176":
return mp.polylog(values["s"], values["x"])
if tid == "T178":
return mp.hyp2f1(1 / values["r"], 1 - 1 / values["r"], 1, values["x"])
if tid == "T250":
return mp.polygamma(int(params["n"]), values["x"])
if tid in ("T351", "T352"):
shape, point = values["a"], values["x"]
twice_shape = int(2 * shape)
if twice_shape % 2 == 0:
integer = twice_shape // 2
upper = mp.factorial(integer - 1) * mp.exp(-point) * mp.fsum(
point ** power / mp.factorial(power) for power in range(integer))
else:
upper = mp.sqrt(mp.pi) * mp.erfc(mp.sqrt(point))
current = mp.mpf("0.5")
while current < shape:
upper = current * upper + point ** current * mp.exp(-point)
current += 1
return upper if tid == "T351" else mp.gamma(shape) - upper
if tid == "T353":
return mp.betainc(values["a"], values["b"], 0, values["x"])
raise ValueError("unknown independent evaluator")
def identity(tid, params):
values = {name: rational(value) for name, value in params.items()
if "sqrt" not in value}
if tid == "T178" and params["x"] == "1/2":
upper = 1 / values["r"]
return arb.pi().sqrt() / (((upper + 1) / 2).gamma() * (1 - upper / 2).gamma())
if tid == "T176" and params == {"s": "2", "x": "1/2"}:
return arb.pi() ** 2 / 12 - arb(2).log() ** 2 / 2
if tid == "T353":
first, second, point = values["a"], values["b"], values["x"]
total = first.gamma() * second.gamma() / (first + second).gamma()
return total - (1 - point).beta_lower(second, first)
if tid in ("T351", "T352"):
total = values["a"].gamma()
complement = (values["x"].gamma_lower(values["a"]) if tid == "T351"
else values["x"].gamma_upper(values["a"]))
return total - complement
if tid == "T250" and params["x"] in ("1", "1/2"):
order = int(params["n"])
factor = 1 if params["x"] == "1" else 2 ** (order + 1) - 1
return (-1) ** (order + 1) * math.factorial(order) * factor * arb(order + 1).zeta()
return None
def contains_written(text, value, digits):
lower, upper, unit = bounds(text)
if isinstance(value, Fraction):
if not lower <= value <= upper:
raise AssertionError("exact result outside written interval")
return {"exact_result": str(value), "stored_exact": not bool(unit)}
if not value.is_finite():
raise AssertionError("nonfinite ball")
if not unit:
if value == rational(lower):
return {"stored_exact": True, "exact_ball": True}
raise AssertionError("inexact ball cannot certify an exact decimal-free value")
if value.rel_accuracy_bits() < numberdb.bits(digits + 10, losing=0):
raise AssertionError("uninformative or insufficiently precise ball")
measured_lower, measured_upper = endpoint(value.lower()), endpoint(value.upper())
if not lower < measured_lower <= measured_upper < upper:
raise AssertionError("stored interval does not contain the result enclosure")
if lower + 10 * unit < measured_lower <= measured_upper < upper + 10 * unit:
raise AssertionError("shifted-decimal control was not rejected")
return {"decimal_contains_enclosure": True,
"certified_accuracy_bits": value.rel_accuracy_bits(),
"shifted_decimal_rejected": True}
def controls():
with precision(100):
gamma_value = evaluate("T351", {"a": "1", "x": "1"})
target = arb(-1).exp()
if not gamma_value.is_finite() or gamma_value.rel_accuracy_bits() < 360:
raise AssertionError("nonzero gamma control is not informative")
if not (gamma_value - target).contains(0) or gamma_value.contains(0):
raise AssertionError("nonzero gamma control failed")
beta_value = evaluate("T353", {"x": "1/2", "a": "1/2", "b": "1/2"})
if not beta_value.is_finite() or beta_value.rel_accuracy_bits() < 360:
raise AssertionError("arcsine control is not informative")
if not (beta_value - arb.pi() / 2).contains(0):
raise AssertionError("arcsine control failed")
if evaluate("T353", {"x": "1/2", "a": "2", "b": "2"}) != Fraction(1, 12):
raise AssertionError("rational beta control failed")
if evaluate("T178", {"r": "3", "x": "0"}) != 1:
raise AssertionError("hypergeometric normalization control failed")
if evaluate("T176", {"s": "2", "x": "0"}) != 0:
raise AssertionError("polylogarithm exact-zero control failed")
return {"gamma_1_1_equals_exp_minus_1": True,
"beta_half_half_at_half_equals_pi_over_2": True,
"rational_beta_2_2_at_half_equals_1_over_12": True,
"hypergeometric_at_zero_equals_1": True,
"polylogarithm_at_zero_equals_0": True}
def check(document, manifest):
tid, digits = manifest["tid"], manifest.get("digits", 100)
rows = flatten(document)
expected = [tuple(params[name] for name in PARAMETERS[tid]) for params in manifest["parameters"]]
if len(expected) != len(set(expected)) or set(expected) != set(rows):
raise AssertionError("enumeration, duplicates or stored identities disagree")
report = {"tid": tid, "entries": len(rows), "digits": digits,
"controls": controls(), "records": [], "rigour": "proven"}
for params in manifest["parameters"]:
key = tuple(params[name] for name in PARAMETERS[tid])
text = rows[key]
lower, upper, unit = bounds(text)
wanted = max(digits, len(Decimal(text).as_tuple().digits)) if unit else digits
record = {"parameters": params, "number": text}
with precision(wanted):
measured = evaluate(tid, params)
record.update(contains_written(text, measured, wanted))
other = identity(tid, params)
if other is not None and unit:
contains_written(text, other, wanted)
record["identity_check"] = True
for working in (wanted + 50, wanted + 100):
with mp.workdps(working):
measured = independent(tid, params)
if not mp.isfinite(measured) or abs(mp.im(measured)) > mp.mpf(10) ** (-wanted - 20):
raise AssertionError("nonfinite or non-real independent value: " + str(params))
measured = mp.re(measured)
if unit:
if not mp_fraction(lower) < measured < mp_fraction(upper):
raise AssertionError("independent value outside written interval: " + str(params))
elif abs(measured - mp_fraction(lower)) > mp.mpf(10) ** (-wanted - 20) * max(1, abs(measured)):
raise AssertionError("independent exact-value mismatch: " + str(params))
record["independent_precisions"] = [wanted + 50, wanted + 100]
report["records"].append(record)
accuracies = [row["certified_accuracy_bits"] for row in report["records"] if "certified_accuracy_bits" in row]
report["minimum_certified_accuracy_bits"] = min(accuracies) if accuracies else None
report["working_guard_bits"] = 160
report["identity_checks"] = sum(row.get("identity_check", False) for row in report["records"])
report["exact_rows"] = sum(row.get("stored_exact", False) for row in report["records"])
return report
def run():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--document", type=Path, default=Path("table.yaml"))
parser.add_argument("--compute", action="store_true")
arguments = parser.parse_args()
generator = BalancedValues()
if arguments.compute:
for params in generator.enumerate():
print(json.dumps({"parameters": params,
"number": to_text(generator.value(params, generator.digits), generator.digits)}))
return
if arguments.document.suffix == ".json":
document = json.loads(arguments.document.read_text())
else:
import yaml
document = yaml.safe_load(arguments.document.read_text())
print(json.dumps(check(document, generator.manifest), indent=2))