generate.py

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

2928 bytes, as of the version from 2026-10-25 16:01. Recorded here, not run.

"""Values of the Kelvin function ker_nu(x) -- numberdb.org/T493.

Run with SageMath (or a modular passagemath installation):
    sage -pip install numberdb
    sage -python generate.py                   # verify the stored draft
    sage -python generate.py --fill            # fill with NUMBERDB_API_KEY set
    sage -python checks.py                     # independent checks, no key needed
    sage -python checks.py --stored            # check the stored draft (key needed)

This stores the DLMF 10.61.2 convention, the real part of
exp(-nu*pi*i/2) K_nu(x exp(pi*i/4)), on the principal branch of K.
The arguments and orders follow the other Kelvin tables in family #205.
All input rounding, rotation and special-function errors are enclosed in
complex-ball arithmetic. Each value is returned as a real ball at 397 bits
for 100 requested decimal digits. checks.py measures every radius and checks
each written interval against independent certified quadrature.

--fill only writes entries to the existing draft. It does not publish it or
offer it for review. The key may instead be read from stdin by setting
NUMBERDB_KEY_FROM_STDIN=1. NUMBERDB_PUBLISH=1 is the wrapper's fill switch.
"""
import os
import sys

import numberdb.sage as numberdb
from sage.rings.complex_arb import ComplexBallField
from sage.rings.rational_field import QQ

WORKING_GUARD = 64


def arguments():
    return sorted({QQ(a) / b for b in range(1, 5)
                   for a in range(1, 20 * b + 1)})


def kelvin_pair(nu, x, precision):
    field = ComplexBallField(precision)
    order = QQ(nu)
    i = field(0, 1)
    phase = ((-i)**int(order) if order.denominator() == 1
             else (-i * field.pi() * field(order) / 2).exp())
    z = field(QQ(x)) * (1 + i) / field(2).sqrt()
    return phase * z.bessel_K(field(order))


class KelvinKer(numberdb.Generator):
    table = 'T493'
    parameters = ('nu', 'x')
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'check.py')

    def enumerate(self):
        for nu in range(3):
            for x in arguments():
                yield {'nu': str(nu), 'x': str(x)}

    def value(self, params, digits):
        value = kelvin_pair(params['nu'], params['x'],
                            numberdb.bits(digits, losing=WORKING_GUARD)).real()
        if not value.is_finite():
            raise ArithmeticError('nonfinite Kelvin value at %r' % params)
        return value


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    generator = KelvinKer()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill certified Kelvin ker values with independent integral checks'))
    else:
        result = generator.verify(sample=None)
        print(result)
        sys.exit(0 if result.ok else 1)