generate.py

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

2911 bytes, as of the version from 2026-10-27 02:59. Recorded here, not run.

"""Rayleigh sums of Bessel derivatives -- numberdb.org/T535.

Run with SageMath or passagemath and the numberdb package installed:
    sage -pip install numberdb          # once
    sage -python generate.py             # verify all stored entries
    sage -python checks.py               # independent local checks
    sage -python generate.py --fill      # write entries to this draft

The --fill action needs NUMBERDB_API_KEY. It writes values and attachments;
it does not offer, review or publish the website draft. Default enumeration:
nu=0,1,2,3, and n=1,...,12, in that order. See argument-selection.md.

For nu>0 the constant-term-one entire function of t is
  F_nu(t) = 2^nu Gamma(nu)/t^((nu-1)/2) J'_nu(sqrt(t)).
Its coefficient c_m is (-1)^m (nu+2m)/(nu 4^m m! (nu+1)_m).
Its product is product_k (1-t/(j'_(nu,k))^2), DLMF 10.21.16.
For nu=0 use J'_0=-J_1 and F_0(t)=2 J_1(sqrt(t))/sqrt(t).
The branch-free series defines these normalized functions also at t=0.

Comparing coefficients in -t F'_nu=F_nu sum_n tau_n t^n gives
  tau_n = -n c_n - sum_(k=1)^(n-1) c_(n-k) tau_k.
All arithmetic is exact QQ; no root truncation or decimal rounding is used.
"""

import sys
from functools import lru_cache
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ


@lru_cache(maxsize=None)
def rayleigh_derivative(nu, n):
    nu = QQ(nu)
    if nu < 0 or n < 1 or int(n) != n:
        raise ValueError('Require rational nu >= 0 and integer n >= 1')
    n = int(n)
    order = QQ(1) if nu == 0 else nu
    coefficients = [QQ(1)]
    base = QQ(1)
    for m in range(1, n + 1):
        base = -base / (QQ(4) * QQ(m) * (order + QQ(m)))
        coefficients.append(base if nu == 0 else base * (nu + 2*m) / nu)
    moments = [QQ(0)]
    for m in range(1, n + 1):
        moments.append(-QQ(m)*coefficients[m] - sum(
            (coefficients[m-k]*moments[k] for k in range(1, m)), QQ(0)))
    return moments[n]


class BesselDerivativeRayleigh(numberdb.Generator):
    table = 'T535'
    parameters = ('nu', 'n')
    type = 'Q'
    rigour = 'exact'
    digits = 100
    files = ('generate.py', 'checks.py', 'check.py', 'argument-selection.md',
             'numerical-evidence.md', 'skill-evidence.json')

    def enumerate(self):
        for nu in range(4):
            for n in range(1, 13):
                yield {'nu': str(nu), 'n': str(n)}

    def value(self, params, digits):
        order = Fraction(str(params['nu']))
        nu = QQ(order.numerator) / QQ(order.denominator)
        return rayleigh_derivative(nu, int(params['n']))


if __name__ == '__main__':
    generator = BesselDerivativeRayleigh()
    if '--fill' in sys.argv:
        print(generator.publish(message='Fill exact derivative Rayleigh sums after independent checks'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)