balanced_generate.py

back to table · edit · history · where entries came from · files · download

9743 bytes, as of the version from 2026-10-27 01:13 (current). Recorded here, not run.

"""Recompute this table's entire attached selection, including retained entries.

Install: python -m pip install numberdb python-flint mpmath
For value-table bundles, run the attached forwarding entry point:
    $ python generate.py
Or run this standalone file directly:
    $ python balanced_generate.py
Online verification: python balanced_generate.py --verify

The default is an offline comparison with balanced-selection.json. It never
writes to the website. Every retained entry's value and annotations come from
the identified website revision, not this expansion run. All entries are
recomputed; retained numbers are not returned as if this code produced them.
The old entry point is preserved as pre-balanced-generator.py; its supporting
attachments remain unchanged. For isolated critical arguments this
module delegates certification to the current attached intrinsic_values.py,
whose digest is recorded in balanced-selection.json. Its rational bracket is
recomputed and its uncertainty propagated, never replaced by a decimal root.

Ordinary additions use exact rational inputs or Arb-enclosed algebraic and
hyperbolic inputs. The 440-bit guard is inherited from the critical-argument
certification budget, and every result must carry at least 15 surplus decimal
digits. Independent checks and measured accuracies are in balanced-checks.json.
In root-table bundles this module only supplies arithmetic helpers to
balanced_root_checks.py; the preserved Sage generator remains the entry point.
"""

import argparse
import hashlib
import importlib.util
import json
import math
from contextlib import contextmanager
from decimal import Decimal
from fractions import Fraction
from pathlib import Path

import numberdb
from flint import acb, arb, ctx


@contextmanager
def precision(bits):
    previous = ctx.prec
    ctx.prec = bits
    try:
        yield
    finally:
        ctx.prec = previous


def endpoint(value):
    mantissa, exponent = value.man_exp()
    return Fraction(int(mantissa)) * Fraction(2) ** int(exponent)


def as_interval(value):
    if not value.is_finite():
        raise ArithmeticError("Nonfinite enclosure")
    if value.is_exact():
        exact = endpoint(value)
        return numberdb.RealInterval(exact, exact)
    with precision(max(ctx.prec, int(value.rel_accuracy_bits()) + 64)):
        return numberdb.RealInterval(endpoint(value.lower()), endpoint(value.upper()))


def written_interval(text):
    text = str(text).rstrip("?")
    center = Fraction(text)
    unit = Fraction(10) ** Decimal(text).as_tuple().exponent if "." in text or "e" in text.lower() else Fraction(0)
    return center - unit, center + unit


def contains(text, value, required_bits=350):
    lower, upper = written_interval(text)
    if isinstance(value, (int, Fraction)):
        return lower <= value <= upper
    if not value.is_finite() or value.rel_accuracy_bits() < required_bits:
        return False
    interval = as_interval(value)
    return lower <= interval.lower <= interval.upper <= upper


def real_part(value):
    if not value.is_finite() or not value.imag.contains(0):
        raise ArithmeticError("Expected a finite real enclosure")
    return value.real


def argument_ball(row):
    descriptor = row["argument"]
    kind = descriptor["type"]
    if kind == "rational":
        return arb(descriptor["value"])
    if kind == "airy":
        value = arb(descriptor["value"])
        return 2 * value * value.sqrt() / 3
    if kind == "scaled_sqrt":
        return arb(descriptor["factor"]) * arb(descriptor["radicand"]).sqrt()
    if kind == "cosh":
        return arb(descriptor["value"]).cosh()
    if kind == "turning":
        order = int(row["params"]["ell"])
        coupling = arb(row["params"]["eta"])
        return coupling + (coupling**2 + order * (order + 1)).sqrt()
    raise ValueError("Unsupported argument descriptor: " + kind)


def evaluate(kind, params, argument):
    if kind.startswith("coulomb"):
        method = "coulomb_f" if kind == "coulombf" else "coulomb_g"
        return real_part(getattr(acb(argument), method)(int(params["ell"]), arb(params["eta"])))
    order = arb(params["nu"])
    if kind.startswith("bessel"):
        if argument == 0 and kind in ("besselj", "besseli"):
            return 1 if order == 0 else 0
        return getattr(argument, "bessel_" + kind[-1])(order)
    if kind in ("ber", "bei", "ker", "kei"):
        if argument == 0 and order == 0:
            if kind == "ber":
                return 1
            if kind == "bei":
                return 0
            if kind == "kei":
                return -arb.pi() / 4
            raise ValueError("ker has no finite origin value")
        if kind in ("ber", "bei"):
            value = (argument * acb(-1, 1) / arb(2).sqrt()).bessel_j(order)
        else:
            phase = acb(0, -arb.pi() * order / 2).exp()
            value = phase * (argument * acb(1, 1) / arb(2).sqrt()).bessel_k(order)
        return value.real if kind in ("ber", "ker") else value.imag
    if kind in ("struveh", "struvel"):
        if argument == 0 and order >= 0:
            return 0
        half = arb("1/2")
        scaled = acb(argument**2 / 4)
        if kind == "struveh":
            scaled = -scaled
        value = ((acb(argument) / 2) ** (order + 1)
                 * scaled.hypgeom([1], [3 * half, order + 3 * half])
                 / ((3 * half).gamma() * (order + 3 * half).gamma()))
        return real_part(value)
    if kind == "legenq":
        half = arb("1/2")
        return (arb.pi().sqrt() * (order + 1).gamma()
                * (1 / argument**2).hypgeom_2f1((order + 2) / 2, (order + 1) / 2, order + 3 * half)
                / (arb(2) ** (order + 1) * (order + 3 * half).gamma() * argument ** (order + 1)))
    raise ValueError("Unsupported function: " + kind)


class BalancedValues(numberdb.Generator):
    type = "R"
    digits = 100
    rigour = "proven"

    def __init__(self, selection_path=None):
        self.path = Path(selection_path) if selection_path else Path(__file__).with_name("balanced-selection.json")
        self.selection = json.loads(self.path.read_text())
        self.table = self.selection["tid"]
        self.kind = self.selection["kind"]
        self.parameters = tuple(self.selection["parameters"])
        self.rows = {tuple(row["params"][key] for key in self.parameters): row for row in self.selection["rows"]}
        if len(self.rows) != len(self.selection["rows"]):
            raise ValueError("Duplicate parameter tuple")
        self.legacy = None
        if any(row["argument"]["type"] == "critical" for row in self.selection["rows"]):
            source = self.path.with_name("intrinsic_values.py")
            expected = self.selection["critical_source_sha256"]
            if hashlib.sha256(source.read_bytes()).hexdigest() != expected:
                raise ValueError("Critical-argument source changed")
            specification = importlib.util.spec_from_file_location("retained_intrinsic_values", source)
            self.legacy = importlib.util.module_from_spec(specification)
            specification.loader.exec_module(self.legacy)

    def enumerate(self):
        for row in self.selection["rows"]:
            yield row["params"]

    def argument(self, row):
        if row["argument"]["type"] == "critical":
            if self.legacy is None:
                raise ValueError("Missing retained critical-point certificate implementation")
            order = Fraction(row["params"]["nu"])
            if self.kind not in ("besselj", "bessely"):
                if order.denominator != 1:
                    raise ValueError("Retained non-Bessel critical source requires integer order")
                order = order.numerator
            return self.legacy.critical_argument(
                self.kind, order, *row["argument"]["bracket"], ctx.prec)
        return argument_ball(row), None

    def ball(self, row, digits=100, bits=None):
        with precision(bits or math.ceil(digits * math.log2(10)) + 440):
            argument, certificate = self.argument(row)
            value = evaluate(self.kind, row["params"], argument)
            if not isinstance(value, int):
                if not value.is_finite() or value.rel_accuracy_bits() < math.ceil((digits + 15) * math.log2(10)):
                    raise ArithmeticError("Insufficient certified accuracy: " + str(row["params"]))
            return value, certificate

    def value(self, params, digits):
        row = self.rows[tuple(str(params[key]) for key in self.parameters)]
        value, unused = self.ball(row, digits)
        result = value if isinstance(value, int) else as_interval(value)
        entry = row.get("entry")
        if isinstance(entry, dict):
            return dict(entry, number=result)
        return result

    def check_local(self):
        for row in self.selection["rows"]:
            entry = row["entry"]
            text = entry["number"] if isinstance(entry, dict) else entry
            value, unused = self.ball(row, max(100, int(row.get("digits", 100))), bits=1000)
            if not contains(text, value):
                raise ArithmeticError("Stored interval mismatch: " + str(row["params"]))
        print(self.table, len(self.rows), "entries recomputed from the attached selection")


if __name__ == "__main__":
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--verify", action="store_true")
    parser.add_argument("--selection", type=Path)
    options = parser.parse_args()
    generator = BalancedValues(options.selection)
    if options.verify:
        report = generator.verify(sample=None)
        print(report)
        raise SystemExit(0 if report.ok else 1)
    generator.check_local()