back to table · edit · history · where entries came from · files · download
6809 bytes, as of the version from 2026-10-26 15:42 (current). Recorded here, not run.
"""Generators for individually justified special-function selections.
Install numberdb, python-flint==0.6.0 and mpmath==1.3.0. Each table's
generate.py verifies by default and never writes without an explicit --fill.
Rational and square-root arguments are enclosed, never rounded before use.
"""
import argparse
import math
from contextlib import contextmanager
from fractions import Fraction
import numberdb
from flint import arb, arb_series, ctx, fmpq
PARAMETERS = {
"T14": ("s",), "T51": ("a", "b"), "T94": ("x", "s"),
"T176": ("s", "x"), "T178": ("r", "x"), "T228": ("s",),
"T246": ("s",), "T351": ("a", "x"), "T352": ("a", "x"),
"T353": ("x", "a", "b"),
}
WORKING_GUARD = 160
@contextmanager
def precision(digits):
previous = ctx.prec
ctx.prec = numberdb.bits(digits, losing=WORKING_GUARD)
try:
yield
finally:
ctx.prec = previous
def rational(value):
return arb(str(Fraction(value)))
def argument(identity):
if identity == "sqrt(2)":
return arb(2).sqrt()
if identity == "1/sqrt(2)":
return 1 / arb(2).sqrt()
return rational(identity)
def endpoint(value):
mantissa, exponent = value.man_exp()
return Fraction(int(mantissa)) * Fraction(2) ** int(exponent)
def enclosure(value, digits):
if not value.is_finite():
raise ArithmeticError("nonfinite result")
if value.rel_accuracy_bits() < numberdb.bits(digits + 10, losing=0):
raise ArithmeticError("insufficient relative accuracy")
return numberdb.RealInterval(endpoint(value.lower()), endpoint(value.upper()))
def bernoulli(index):
return Fraction(str(fmpq.bernoulli(index)))
def bernoulli_polynomial(index, value):
return sum(math.comb(index, position) * bernoulli(position) * value ** (index - position)
for position in range(index + 1))
def order_key(value):
return (abs(Fraction(value)), Fraction(value) < 0)
def selection(tid):
if tid in ("T14", "T228", "T246"):
bound = 40 if tid == "T14" else 20
indices = [Fraction(integer) for integer in range(-bound, bound + 1)
if integer != 1 and not (tid == "T246" and integer < 0 and integer % 2 == 0)]
indices.append(Fraction(1, 2))
if tid == "T14":
indices.extend([Fraction(3, 2), Fraction(5, 2)])
for index in sorted(indices, key=order_key):
yield {"s": str(index)}
elif tid == "T51":
yield {"a": "1", "b": "sqrt(2)"}
yield {"a": "1/sqrt(2)", "b": "1"}
elif tid == "T94":
for denominator in range(2, 6):
for numerator in range(1, denominator):
if math.gcd(numerator, denominator) != 1:
continue
indices = [Fraction(integer) for integer in range(-10, 11) if integer != 1]
indices.append(Fraction(1, 2))
for index in sorted(indices, key=order_key):
yield {"x": str(Fraction(numerator, denominator)), "s": str(index)}
elif tid == "T176":
for order in (2, 3, 4):
identities = ("1", "-1", "1/2", "1/4", "-1/2", "1/3", "1/9") if order == 2 else ("1", "-1", "1/2")
for identity in identities:
yield {"s": str(order), "x": identity}
elif tid == "T178":
for signature in (2, 3, 4, 6):
yield {"r": str(signature), "x": "1/2"}
elif tid in ("T351", "T352"):
for dimension in range(1, 11):
identity = str(Fraction(dimension, 2))
yield {"a": identity, "x": identity}
elif tid == "T353":
for dimension in range(1, 11):
identity = str(Fraction(dimension, 2))
yield {"x": "1/2", "a": identity, "b": identity}
else:
raise ValueError("table outside repair assignment")
def evaluate(tid, params):
if tid in ("T14", "T94"):
index = Fraction(params["s"])
if index.denominator == 1 and index <= 0:
if tid == "T94":
degree = 1 - int(index)
return -bernoulli_polynomial(degree, Fraction(params["x"])) / degree
if index == 0:
return Fraction(-1, 2)
return -bernoulli(1 - int(index)) / (1 - int(index))
return rational(index).zeta() if tid == "T14" else rational(index).zeta(rational(params["x"]))
if tid in ("T228", "T246"):
index = rational(params["s"])
derivative = arb_series([index, 1], 2).zeta()[1]
return derivative if tid == "T228" else derivative / index.zeta()
if tid == "T51":
return argument(params["a"]).agm(argument(params["b"]))
if tid == "T176":
if params["x"] == "1":
return arb(int(params["s"])).zeta()
return rational(params["x"]).polylog(int(params["s"]))
if tid == "T178":
upper = 1 / arb(int(params["r"]))
return rational(params["x"]).hypgeom_2f1(upper, 1 - upper, 1)
if tid in ("T351", "T352"):
shape, point = rational(params["a"]), rational(params["x"])
return point.gamma_upper(shape) if tid == "T351" else point.gamma_lower(shape)
if tid == "T353":
first, second, point = (Fraction(params[name]) for name in ("a", "b", "x"))
if first == second and point == Fraction(1, 2) and first.denominator == 1:
integer = first.numerator
return Fraction(math.factorial(integer - 1) ** 2, 2 * math.factorial(2 * integer - 1))
return rational(point).beta_lower(rational(first), rational(second))
raise ValueError("unknown table")
class SelectedValues(numberdb.Generator):
type = "R"
digits = 100
rigour = "proven"
files = ("generate.py", "values.py", "checks.py", "checks-evidence.json", "selection.md")
def __init__(self, tid):
if tid not in PARAMETERS:
raise ValueError("unknown table")
self.table = tid
self.parameters = PARAMETERS[tid]
def enumerate(self):
yield from selection(self.table)
def value(self, params, digits):
with precision(digits):
value = evaluate(self.table, params)
if isinstance(value, Fraction):
return value
return enclosure(value, digits)
def run(tid):
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--fill", action="store_true")
parser.add_argument("--remove-unselected", action="store_true")
arguments = parser.parse_args()
generator = SelectedValues(tid)
if arguments.fill:
print(generator.publish(removing=arguments.remove_unselected,
message="Recompute the documented intrinsic selection"))
else:
report = generator.verify(sample=None)
print(report)
raise SystemExit(0 if report.ok else 1)