[Boost-bugs] [Boost C++ Libraries] #2733: Cylindrical bessel fuction gives wrong result if n=-1 and z>1

Subject: [Boost-bugs] [Boost C++ Libraries] #2733: Cylindrical bessel fuction gives wrong result if n=-1 and z>1
From: Boost C++ Libraries (noreply_at_[hidden])
Date: 2009-02-06 08:48:24

#2733: Cylindrical bessel fuction gives wrong result if n=-1 and z>1
--------------------------------------------+-------------------------------
 Reporter:  Ivan Lisenkov <ivan_at_[hidden]>  |       Owner:         
     Type:  Patches                         |      Status:  new    
Milestone:  Boost 1.38.0                    |   Component:  None   
  Version:  Boost 1.37.0                    |    Severity:  Problem
 Keywords:  bessel                          |  
--------------------------------------------+-------------------------------
 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
 }}}

-- 
Ticket URL: <https://svn.boost.org/trac/boost/ticket/2733>
Boost C++ Libraries <http://www.boost.org/>
Boost provides free peer-reviewed portable C++ source libraries.