Cylindrical bessel fuction gives wrong result if n=-1 and z>1

[Ivan Lisenkov <ivan@…>]

If the order of Bessel function is -1 and argument is more than 1.0 the algorithm goes on wrong branch, giving always 0 as result. Here a program to illustrate the problem:

#include <iostream>
#include <boost/math/special_functions/bessel.hpp>

int main()
{
	std::cout.setf (std::ios::scientific | std::ios::left);
	std::cout.precision(5);
		
	double z=0.99999;	
	std::cout<<"z="<<z<<" J_{-1}(z)="<<boost::math::cyl_bessel_j(-1, z)<<'\n';
	z=1.00007e+00;
	std::cout<<"z="<<z<<" J_{-1}(z)="<<boost::math::cyl_bessel_j(-1, z)<<'\n';

	return 0;
}

The output of the program is:

z=9.99990e-01 J_{-1}(z)=-4.40047e-01
z=1.00007e+00 J_{-1}(z)=-0.00000e+00

I have created a patch against the latest subverversion:

Index: boost/math/special_functions/detail/bessel_jn.hpp
===================================================================
--- boost/math/special_functions/detail/bessel_jn.hpp   (revision 51047)
+++ boost/math/special_functions/detail/bessel_jn.hpp   (working copy)
@@ -32,9 +32,12 @@
     {
         return bessel_j0(x);
     }
-    if (n == 1)
+    if ((n == 1) || (n == -1))
     {
-        return bessel_j1(x);
+        factor = n;  // J_{-1}(z) = -1*J_1(z)
+        value = bessel_j1(x);
+        value *=factor;
+        return value;
     }
     if (n < 0)
     {

The output after the patch is:

z=9.99990e-01 J_{-1}(z)=-4.40047e-01
z=1.00007e+00 J_{-1}(z)=-4.40073e-01

This is correspond to GNU Octave:

octave:32> besselj(-1,0.99999)
ans = -0.44005 - 0.00000i

octave:26> besselj(-1,1.00007)
ans = -0.44007 - 0.00000i

[John Maddock]

(In [51059]) Fix bug in bessel_jn for n == -1. Add new test case. Checked that the other Bessel functions do not have the same issue. Checked that real-valued -1 argument is fixed OK as well as integer argument. Fixes #2733.