checks.py

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

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

"""Check every curated value with mpmath at two independent precisions.

Arb provides the rigorous enclosures. mpmath comparisons are independent
agreement checks, not certificates. Derivative residuals and the Airy
Wronskian test the mathematical interpretation of the intrinsic arguments.
"""

import argparse
import importlib.metadata
import json
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 curated_values import CuratedValues, KINDS, argument_for, evaluate
from scorer import as_interval, precision, working_precision


def mp_argument(identity):
    prefixes = {'ai-prime-zero-': (mp.airyaizero, 1),
                'bi-prime-zero-': (mp.airybizero, 1),
                'ai-zero-': (mp.airyaizero, 0), 'bi-zero-': (mp.airybizero, 0)}
    for prefix, (function, derivative) in prefixes.items():
        if identity.startswith(prefix):
            return function(int(identity[len(prefix):]), derivative=derivative)
    if identity.startswith('sqrt('):
        return mp.sqrt(int(identity[5:-1]))
    if identity.endswith('*pi/2'):
        return int(identity[:-5]) * mp.pi / 2
    if identity.endswith('*pi'):
        return int(identity[:-3]) * mp.pi
    if identity.endswith('/sqrt(2)'):
        return int(identity[:-8]) / mp.sqrt(2)
    rational = Fraction(identity)
    return mp.mpf(rational.numerator) / rational.denominator


def mp_value(kind, argument):
    if kind in ('ai', 'ai_prime', 'bi', 'bi_prime'):
        function = mp.airyai if kind.startswith('ai') else mp.airybi
        return function(argument, derivative=int(kind.endswith('prime')))
    names = {'fresnel_s': 'fresnels', 'fresnel_c': 'fresnelc'}
    return getattr(mp, names.get(kind, kind))(argument)


def check(tid, live=False):
    generator = CuratedValues()
    generator.table = tid
    kind = KINDS[tid]
    stored = numberdb.table(tid)['Numbers'] if live else None
    numbers, records = {}, []
    with precision(working_precision(130)):
        for params in generator.enumerate():
            identity = params['x']
            entry = generator.value(params, 100)
            generated = to_text(entry['number'], 100)
            written = generated if stored is None else stored[identity]
            if isinstance(written, dict):
                written = written['number']
            exact = isinstance(entry['number'], int)
            midpoint = Fraction(Decimal(str(written)))
            unit = 0 if exact else Fraction(10) ** Decimal(str(written)).as_tuple().exponent
            lower, upper = midpoint - unit, midpoint + unit
            argument = argument_for(identity)
            value = evaluate(kind, argument)
            interval = as_interval(value)
            if not lower <= interval.lower <= interval.upper <= upper:
                raise ArithmeticError('certified enclosure outside stored digits: ' + identity)
            for decimal_digits in (150, 200):
                with mp.workdps(decimal_digits):
                    reference = mp_value(kind, mp_argument(identity))
                    lower_mp = mp.mpf(lower.numerator) / lower.denominator
                    upper_mp = mp.mpf(upper.numerator) / upper.denominator
                    if not lower_mp <= reference <= upper_mp:
                        raise ArithmeticError('independent mismatch: ' + identity)
            if 'zero-' in identity:
                airy_values = argument.airy()
                residual_index = {'ai': 1, 'ai_prime': 0, 'bi': 3, 'bi_prime': 2}[kind]
                if not airy_values[residual_index].contains(0):
                    raise ArithmeticError('wrong Airy root')
                wronskian = airy_values[0] * airy_values[3] - airy_values[1] * airy_values[2]
                if not (wronskian - 1 / arb.pi()).contains(0):
                    raise ArithmeticError('Airy Wronskian failure')
            if identity != '0' and kind in ('si', 'ci', 'fresnel_s', 'fresnel_c'):
                phase = argument if kind in ('si', 'ci') else arb.pi() * argument ** 2 / 2
                residual = phase.sin() if kind in ('si', 'fresnel_s') else phase.cos()
                if not residual.contains(0):
                    raise ArithmeticError('wrong stationary point')
            numbers[identity] = dict(entry, number=generated)
            records.append({'argument': identity, 'certified_interval_contained': True,
                            'mpmath_decimal_precisions': [150, 200],
                            'arb_accuracy_bits': None if exact else value.rel_accuracy_bits()})
        if stored is not None and set(stored) != set(numbers):
            raise ArithmeticError('default and stored arguments differ')
    return numbers, {'table': tid, 'digits': 100, 'entries': len(numbers),
                     'live': live, 'records': records,
                     'versions': {package: importlib.metadata.version(package)
                                  for package in ('python-flint', 'mpmath')}}


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('tables', nargs='+', choices=sorted(KINDS))
    parser.add_argument('--live', action='store_true')
    parser.add_argument('--output-dir', type=Path)
    arguments = parser.parse_args()
    for tid in arguments.tables:
        numbers, evidence = check(tid, arguments.live)
        if arguments.output_dir:
            directory = arguments.output_dir / tid
            directory.mkdir(parents=True, exist_ok=True)
            (directory / 'numbers.json').write_text(json.dumps(numbers, indent=2) + '\n')
            (directory / 'checks-evidence.json').write_text(json.dumps(evidence, indent=2) + '\n')
        print(tid, len(numbers), 'entries checked', flush=True)