Cephes is an old and respected software library full of numerical approximations for eclectic functions, which is incorporated in more modern libraries such as Scipy. According to the comments in the code, the function j1
Returns Bessel function of order one of the argument. The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 24 term Chebyshev expansion is used. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5.
Looking at the actual code (reproduced below) that doesn't sound right:
- The code is clearly splitting on 5 rather than 8
- The lists of constants involved in fitting that smaller interval isn't 24 long: the RP and RQ arrays are 4 and 8 elements long respectively
All that together, my natural conclusion is that the comment has drifted out of sync with the actual code, and it's actually doing something different. I'm trying to figure out what the comment should say, and whether it is just faulty on some numerical particulars (including the precision notes?) or whether there was a fundamental change of approach. My original end goal was to work out how to do some similar approximation for jinc directly, but at this point I'm largely caught up in the puzzle.
Relevant C code
double j1(x)
double x;
{
double w, z, p, q, xn;
w = x;
if (x < 0)
return -j1(-x);
if (w <= 5.0) {
z = x * x;
w = polevl(z, RP, 3) / p1evl(z, RQ, 8);
w = w * x * (z - Z1) * (z - Z2);
return (w);
}
w = 5.0 / x;
z = w * w;
p = polevl(z, PP, 6) / polevl(z, PQ, 6);
q = polevl(z, QP, 7) / p1evl(z, QQ, 7);
xn = x - THPIO4;
p = p * cos(xn) - w * q * sin(xn);
return (p * SQ2OPI / sqrt(x));
}
polevl and p1evl are both polynomial evaluation functions found here, the latter assuming that the first term is 1.
For completeness, here are the constant lists.
static double RP[4] = {
-8.99971225705559398224E8,
4.52228297998194034323E11,
-7.27494245221818276015E13,
3.68295732863852883286E15,
};
static double RQ[8] = {
/* 1.00000000000000000000E0, */
6.20836478118054335476E2,
2.56987256757748830383E5,
8.35146791431949253037E7,
2.21511595479792499675E10,
4.74914122079991414898E12,
7.84369607876235854894E14,
8.95222336184627338078E16,
5.32278620332680085395E18,
};
static double PP[7] = {
7.62125616208173112003E-4,
7.31397056940917570436E-2,
1.12719608129684925192E0,
5.11207951146807644818E0,
8.42404590141772420927E0,
5.21451598682361504063E0,
1.00000000000000000254E0,
};
static double PQ[7] = {
5.71323128072548699714E-4,
6.88455908754495404082E-2,
1.10514232634061696926E0,
5.07386386128601488557E0,
8.39985554327604159757E0,
5.20982848682361821619E0,
9.99999999999999997461E-1,
};
static double QP[8] = {
5.10862594750176621635E-2,
4.98213872951233449420E0,
7.58238284132545283818E1,
3.66779609360150777800E2,
7.10856304998926107277E2,
5.97489612400613639965E2,
2.11688757100572135698E2,
2.52070205858023719784E1,
};
static double QQ[7] = {
/* 1.00000000000000000000E0, */
7.42373277035675149943E1,
1.05644886038262816351E3,
4.98641058337653607651E3,
9.56231892404756170795E3,
7.99704160447350683650E3,
2.82619278517639096600E3,
3.36093607810698293419E2,
};