boost/math/special_functions/detail/bessel_jn.hpp

   1 //  Copyright (c) 2006 Xiaogang Zhang
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_BESSEL_JN_HPP
   7 #define BOOST_MATH_BESSEL_JN_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/special_functions/detail/bessel_j0.hpp>
  14 #include <boost/math/special_functions/detail/bessel_j1.hpp>
  15 #include <boost/math/special_functions/detail/bessel_jy.hpp>
  16 #include <boost/math/special_functions/detail/bessel_jy_asym.hpp>
  17 #include <boost/math/special_functions/detail/bessel_jy_series.hpp>
  18
  19 // Bessel function of the first kind of integer order
  20 // J_n(z) is the minimal solution
  21 // n < abs(z), forward recurrence stable and usable
  22 // n >= abs(z), forward recurrence unstable, use Miller's algorithm
  23
  24 namespace boost { namespace math { namespace detail{
  25
  26 template <typename T, typename Policy>
  27 T bessel_jn(int n, T x, const Policy& pol)
  28 {
  29     T value(0), factor, current, prev, next;
  30
  31     BOOST_MATH_STD_USING
  32
  33     //
  34     // Reflection has to come first:
  35     //
  36     if (n < 0)
  37     {
  38         factor = (n & 0x1) ? -1 : 1;  // J_{-n}(z) = (-1)^n J_n(z)
  39         n = -n;
  40     }
  41     else
  42     {
  43         factor = 1;
  44     }
  45     if(x < 0)
  46     {
  47         factor *= (n & 0x1) ? -1 : 1;  // J_{n}(-z) = (-1)^n J_n(z)
  48         x = -x;
  49     }
  50     //
  51     // Special cases:
  52     //
  53     if (n == 0)
  54     {
  55         return factor * bessel_j0(x);
  56     }
  57     if (n == 1)
  58     {
  59         return factor * bessel_j1(x);
  60     }
  61
  62     if (x == 0)                             // n >= 2
  63     {
  64         return static_cast<T>(0);
  65     }
  66
  67     if(asymptotic_bessel_large_x_limit(T(n), x))
  68       return factor * asymptotic_bessel_j_large_x_2<T>(n, x);
  69
  70     BOOST_ASSERT(n > 1);
  71     T scale = 1;
  72     if (n < abs(x))                         // forward recurrence
  73     {
  74         prev = bessel_j0(x);
  75         current = bessel_j1(x);
  76         policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
  77         for (int k = 1; k < n; k++)
  78         {
  79             T fact = 2 * k / x;
  80             //
  81             // rescale if we would overflow or underflow:
  82             //
  83             if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))
  84             {
  85                scale /= current;
  86                prev /= current;
  87                current = 1;
  88             }
  89             value = fact * current - prev;
  90             prev = current;
  91             current = value;
  92         }
  93     }
  94     else if((x < 1) || (n > x * x / 4) || (x < 5))
  95     {
  96        return factor * bessel_j_small_z_series(T(n), x, pol);
  97     }
  98     else                                    // backward recurrence
  99     {
 100         T fn; int s;                        // fn = J_(n+1) / J_n
 101         // |x| <= n, fast convergence for continued fraction CF1
 102         boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol);
 103         prev = fn;
 104         current = 1;
 105         // Check recursion won't go on too far:
 106         policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
 107         for (int k = n; k > 0; k--)
 108         {
 109             T fact = 2 * k / x;
 110             if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))
 111             {
 112                prev /= current;
 113                scale /= current;
 114                current = 1;
 115             }
 116             next = fact * current - prev;
 117             prev = current;
 118             current = next;
 119         }
 120         value = bessel_j0(x) / current;       // normalization
 121         scale = 1 / scale;
 122     }
 123     value *= factor;
 124
 125     if(tools::max_value<T>() * scale < fabs(value))
 126        return policies::raise_overflow_error<T>("boost::math::bessel_jn<%1%>(%1%,%1%)", 0, pol);
 127
 128     return value / scale;
 129 }
 130
 131 }}} // namespaces
 132
 133 #endif // BOOST_MATH_BESSEL_JN_HPP
 134