compute the logarithm of Bessel functions of large argument and order

lnbessel.c

// compute the logarithm of Bessel functions of large argument and order
//
// author: Nicolas Tessore <[email protected]>
// date: 6 Dec 2020
//

#include <math.h>


// compute the Bessel function ratio J_{nu-1}(x)/J_nu(x) from a continued
// fraction
//
// reference: Rothwell, Commun. Numer. Meth. Engng 2008; 24:237-249
//
double rjnu(double nu, double x)
{
    int m;
    double r, a, c, d, delt;

    const double eps = 1e-15;
    const double tiny = 1e-30;

    r = 2*nu/x;
    c = r;
    d = 0;
    for(m = 1;; ++m)
    {
        nu += 1;
        a = (1-2*(m&1))*2*nu/x;
        c = a + 1/c;
        if(c == 0)
            c = tiny;
        d = a + d;
        if(d == 0)
            d = tiny;
        d = 1/d;
        delt = c*d;
        r = r*delt;
        if(fabs(delt - 1) < eps)
            break;
    }

    return r;
}


// compute Bessel function logarithm ln[J_{nu+n}(x)] from ln[J_nu(x)]
//
// reference: Rothwell, Commun. Numer. Meth. Engng 2008; 24:237-249
// + explicit tracking of sign change if `sgn` is not NULL
//
double lnjnupn(double nu, int n, double x, double lnjnu, double* sgn)
{
    int k;
    double r, s, c, y, t;

    r = rjnu(nu+n, x);
    k = 0;

    // Kahan summation
    s = lnjnu;
    c = 0;
    for(; n > 0; --n)
    {
        if(r < 0)
            k ^= 1;

        y = -log(fabs(r)) - c;
        t = s + y;
        c = (t - s) - y;
        s = t;

        r = 2*(nu+n-1)/x - 1/r;
    }

    if(sgn)
        *sgn = k ? -1 : +1;

    return s;
}


// test program for spherical Bessel functions -- delete

#include <stdio.h>
int main()
{
    // argument
    double x = 5.0;

    // order of start value
    double nu = 0;
    // abs log of start value: Log[Abs[SphericalBesselJ[0, 5.0`20]]]
    double lnjnu = -1.6513810824626862774;
    // track sign change
    double sgn;

    // how many times to raise the order
    int n = 10000;

    // compute the result: note +1/2 order of spherical Bessel function
    lnjnu = lnjnupn(nu+0.5, n, x, lnjnu, &sgn);

    // logjnu now contains j(nu+n, z) and sgn is the sign change from -1
    printf("jnu = %+f exp(%.16f)\n", -sgn, lnjnu);

    // the true value: Log[SphericalBesselJ[10000, 5.0`20]]
    printf("     [%+f exp(%.16f)]\n", 1.0, -72950.74713290145977);

    return 0;
}

Note that the lnjnupn function is good to compute a single value ln[J_{nu+n}(x)] from ln[J_nu(x)]. If all values ln[J_{nu+k}(x)] for k = 1, ..., n are desired, it is much more efficient to compute first all ratios rk from r and subsequently compute the partial sums forwards.