In besselfrac.h it seems like the cache tests in jnu(), ynu(), inu(), and knu() can, in some situations, cause jo, yo, io, or ko to be returned without calling the method where they are calculated.

I believe this code illustrates the issue

#include <iostream>
#include "Numerics/nr3.h"
#include "Numerics/besselfrac.h"

int main (int argc, char * const argv[]) {
Bessel bes1;
std::cout << "J1(0.75): " << bes1.jnu(1.0,0.75) << endl;
std::cout << "I1(0.75): " << bes1.inu(1.0,0.75) << endl;
bes1.besselik(1.0,0.75);
std::cout << "I1(0.75): " << bes1.inu(1.0,0.75) << endl << endl;

Bessel bes2;
std::cout << "I1(0.75): " << bes2.inu(1.0,0.75) << endl;
std::cout << "J1(0.75): " << bes2.jnu(1.0,0.75) << endl;
bes2.besseljy(1.0,0.75);
std::cout << "J1(0.75): " << bes2.jnu(1.0,0.75) << endl;

return 0;
}


The output on my system is:
J1(0.75): 0.349244
I1(0.75): -3.2233e-232
I1(0.75): 0.401992

I1(0.75): 0.401992
J1(0.75): 0
J1(0.75): 0.349244


I don't think besseljy() and besselik() should share xo and nuo. sphbesj() and sphbesy() can share besseljy's xo. It looks like airy needs to take the value of x into account in order to choose which version of xo and nuo to test.

-Steve Nicholson

Hi Steve,

Yes, you're absolutely correct that this is a bug. We'll post a fix under Official Bug Reports. Thanks!

Saul Teukolsky