back to table · edit · history · where entries came from · files · download
8810 bytes, as of the version from 2026-10-27 01:00 (current). Recorded here, not run.
"""Reproduce the full balanced selection, using exact inputs and Arb balls.
python -m pip install numberdb python-flint==0.6.0 mpmath==1.3.0
python generate.py --compute
python generate.py --document table.yaml
The selection manifest contains all retained identities as well as additions.
No stored decimal is used as a numerical input. The 160-bit working guard is
checked against the achieved relative accuracy of every nonzero ball.
This module never writes to the website.
"""
import json
import math
from contextlib import contextmanager
from fractions import Fraction
from pathlib import Path
import numberdb
from flint import acb, arb, ctx
PARAMETERS = {
"T51": ("a", "b"),
"T176": ("s", "x"),
"T178": ("r", "x"),
"T250": ("n", "x"),
"T351": ("a", "x"),
"T352": ("a", "x"),
"T353": ("x", "a", "b"),
}
GUARD_BITS = 160
METADATA_KEYS = {"equals", "comment", "param-latex", "proof", "number"}
@contextmanager
def precision(digits):
previous = ctx.prec
ctx.prec = numberdb.bits(digits, losing=GUARD_BITS)
try:
yield
finally:
ctx.prec = previous
def rational(value):
return arb(str(Fraction(value)))
def argument(identity):
if identity == "1/sqrt(2)":
return 1 / arb(2).sqrt()
if identity.startswith("sqrt(") and identity.endswith(")"):
return rational(identity[5:-1]).sqrt()
return rational(identity)
def endpoint(value):
mantissa, exponent = value.man_exp()
return Fraction(int(mantissa)) * Fraction(2) ** int(exponent)
def flatten(document, leaves=False):
result = {}
depth = len(document["Parameters"])
def descend(node, keys):
if len(keys) == depth:
if isinstance(node, dict) and "number" not in node:
return
result[tuple(keys)] = node if leaves else str(node["number"] if isinstance(node, dict) else node)
elif isinstance(node, dict):
for key, value in node.items():
if key not in METADATA_KEYS:
descend(value, keys + [key])
descend(document["Numbers"], [])
return result
def additions(tid):
if tid == "T51":
yield {"a": "1", "b": "1"}
for integer in range(2, 17):
yield {"a": "1", "b": str(integer)}
for integer in range(2, 16):
if math.isqrt(integer) ** 2 != integer:
yield {"a": "1", "b": "sqrt(%d)" % integer}
elif tid == "T176":
for order in range(2, 13):
for point in ("1", "-1", "1/2", "-1/2"):
yield {"s": str(order), "x": point}
for order in (2, 3, 4):
for point in ("1/3", "2/3", "1/4", "3/4", "-1/3", "-1/4", "-2", "-3"):
yield {"s": str(order), "x": point}
for twice_order in (1, 3, 5, 7):
for point in ("-1", "-1/2", "1/2"):
yield {"s": str(Fraction(twice_order, 2)), "x": point}
if twice_order > 1:
yield {"s": str(Fraction(twice_order, 2)), "x": "1"}
yield {"s": "2", "x": "0"}
elif tid == "T178":
for signature in (2, 3, 4, 6):
for point in ("0", "1/2", "1/4", "3/4", "1/3", "2/3", "1/9", "8/9"):
yield {"r": str(signature), "x": point}
elif tid in ("T351", "T352"):
for dimension in range(1, 21):
shape = Fraction(dimension, 2)
for point in dict.fromkeys((shape / 2, shape, 2 * shape, Fraction(1))):
yield {"a": str(shape), "x": str(point)}
for threshold in range(1, 7):
yield {"a": "1/2", "x": str(Fraction(threshold ** 2, 2))}
yield {"a": "1", "x": "0"}
elif tid == "T250":
for order in range(3, 13):
for point in ("1", "1/2", "1/3", "2/3", "1/4", "3/4"):
yield {"n": str(order), "x": point}
elif tid == "T353":
for dimension in range(1, 21):
shape = Fraction(dimension, 2)
yield {"x": "1/2", "a": str(shape), "b": str(shape)}
for dimension in range(2, 13):
shape = Fraction(dimension, 2)
other = Fraction(1, 2)
for first, second in ((shape, other), (other, shape)):
yield {"x": str(first / (first + second)), "a": str(first), "b": str(second)}
for dimension in range(1, 13):
shape = Fraction(dimension, 2)
for first, second in ((shape, shape + 1), (shape + 1, shape)):
for point in (Fraction(1, 2), first / (first + second)):
yield {"x": str(point), "a": str(first), "b": str(second)}
for point in ("0", "1/4", "3/4", "1"):
yield {"x": point, "a": "1/2", "b": "1/2"}
else:
raise ValueError("table is outside the expansion assignment")
def exact_beta(first, second, point):
return sum(Fraction((-1) ** index * math.comb(second - 1, index), first + index)
* point ** (first + index) for index in range(second))
def evaluate(tid, params):
if tid == "T51":
first, second = argument(params["a"]), argument(params["b"])
if not (first > 0 and second > 0):
raise ValueError("AGM inputs must be positive")
if params["a"] == params["b"]:
return Fraction(params["a"])
return first.agm(second)
if tid == "T176":
order, point = Fraction(params["s"]), Fraction(params["x"])
if point > 1 or (point == 1 and order <= 1):
raise ValueError("polylogarithm branch or endpoint not in the real domain")
if point == 0:
return Fraction(0)
if point == 1:
return rational(order).zeta()
if point == -1:
return -(1 - arb(2) ** rational(1 - order)) * rational(order).zeta()
return rational(point).polylog(rational(order))
if tid == "T178":
signature = int(params["r"])
point = Fraction(params["x"])
if signature not in (2, 3, 4, 6) or not 0 <= point < 1:
raise ValueError("invalid signature or argument")
if point == 0:
return Fraction(1)
upper = arb(1) / signature
return rational(point).hypgeom_2f1(upper, 1 - upper, 1, abc=True)
if tid in ("T351", "T352"):
shape, point = Fraction(params["a"]), Fraction(params["x"])
if shape <= 0 or point < 0:
raise ValueError("incomplete gamma requires positive shape and nonnegative argument")
if point == 0:
if tid == "T352":
return Fraction(0)
if shape.denominator == 1:
return Fraction(math.factorial(int(shape) - 1))
return rational(shape).gamma()
return (rational(point).gamma_upper(rational(shape)) if tid == "T351"
else rational(point).gamma_lower(rational(shape)))
if tid == "T250":
order, point = int(params["n"]), Fraction(params["x"])
if order < 1 or (point.denominator == 1 and point <= 0):
raise ValueError("polygamma order or pole invalid")
value = ((-1) ** (order + 1) * math.factorial(order)
* acb(order + 1).zeta(acb(rational(point))))
if not value.is_finite() or not value.imag.contains(0):
raise ArithmeticError("non-real polygamma enclosure")
return value.real
if tid == "T353":
first, second, point = (Fraction(params[name]) for name in ("a", "b", "x"))
if first <= 0 or second <= 0 or not 0 <= point <= 1:
raise ValueError("incomplete beta outside the defining integral domain")
if point == 0:
return Fraction(0)
if first.denominator == second.denominator == 1:
return exact_beta(int(first), int(second), point)
return rational(point).beta_lower(rational(first), rational(second))
raise ValueError("unknown table")
class BalancedValues(numberdb.Generator):
type = "R"
digits = 100
rigour = "proven"
def __init__(self, manifest=None):
self.manifest = manifest or json.loads(Path(__file__).with_name("balanced-selection.json").read_text())
self.table = self.manifest["tid"]
self.parameters = PARAMETERS[self.table]
self.digits = self.manifest.get("digits", 100)
def enumerate(self):
yield from self.manifest["parameters"]
def value(self, params, digits):
with precision(digits):
value = evaluate(self.table, params)
if isinstance(value, Fraction):
return value
if not value.is_finite() or value.rel_accuracy_bits() < numberdb.bits(digits + 10, losing=0):
raise ArithmeticError("nonfinite or insufficiently accurate enclosure")
return numberdb.RealInterval(endpoint(value.lower()), endpoint(value.upper()))