special_functions/detail/bessel_jn.hpp

0001 //  Copyright (c) 2006 Xiaogang Zhang
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_BESSEL_JN_HPP
0007 #define BOOST_MATH_BESSEL_JN_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/special_functions/detail/bessel_j0.hpp>
0014 #include <boost/math/special_functions/detail/bessel_j1.hpp>
0015 #include <boost/math/special_functions/detail/bessel_jy.hpp>
0016 #include <boost/math/special_functions/detail/bessel_jy_asym.hpp>
0017 #include <boost/math/special_functions/detail/bessel_jy_series.hpp>
0018
0019 // Bessel function of the first kind of integer order
0020 // J_n(z) is the minimal solution
0021 // n < abs(z), forward recurrence stable and usable
0022 // n >= abs(z), forward recurrence unstable, use Miller's algorithm
0023
0024 namespace boost { namespace math { namespace detail{
0025
0026 template <typename T, typename Policy>
0027 T bessel_jn(int n, T x, const Policy& pol)
0028 {
0029     T value(0), factor, current, prev, next;
0030
0031     BOOST_MATH_STD_USING
0032
0033     //
0034     // Reflection has to come first:
0035     //
0036     if (n < 0)
0037     {
0038         factor = static_cast<T>((n & 0x1) ? -1 : 1);  // J_{-n}(z) = (-1)^n J_n(z)
0039         n = -n;
0040     }
0041     else
0042     {
0043         factor = 1;
0044     }
0045     if(x < 0)
0046     {
0047         factor *= (n & 0x1) ? -1 : 1;  // J_{n}(-z) = (-1)^n J_n(z)
0048         x = -x;
0049     }
0050     //
0051     // Special cases:
0052     //
0053     if(asymptotic_bessel_large_x_limit(T(n), x))
0054        return factor * asymptotic_bessel_j_large_x_2<T>(T(n), x, pol);
0055     if (n == 0)
0056     {
0057         return factor * bessel_j0(x);
0058     }
0059     if (n == 1)
0060     {
0061         return factor * bessel_j1(x);
0062     }
0063
0064     if (x == 0)                             // n >= 2
0065     {
0066         return static_cast<T>(0);
0067     }
0068
0069     BOOST_MATH_ASSERT(n > 1);
0070     T scale = 1;
0071     if (n < abs(x))                         // forward recurrence
0072     {
0073         prev = bessel_j0(x);
0074         current = bessel_j1(x);
0075         policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
0076         for (int k = 1; k < n; k++)
0077         {
0078             value = (2 * k * current / x) - prev;
0079             prev = current;
0080             current = value;
0081         }
0082     }
0083     else if((x < 1) || (n > x * x / 4) || (x < 5))
0084     {
0085        return factor * bessel_j_small_z_series(T(n), x, pol);
0086     }
0087     else                                    // backward recurrence
0088     {
0089         T fn; int s;                        // fn = J_(n+1) / J_n
0090         // |x| <= n, fast convergence for continued fraction CF1
0091         boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol);
0092         prev = fn;
0093         current = 1;
0094         // Check recursion won't go on too far:
0095         policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
0096         for (int k = n; k > 0; k--)
0097         {
0098             T fact = 2 * k / x;
0099             if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))
0100             {
0101                prev /= current;
0102                scale /= current;
0103                current = 1;
0104             }
0105             next = fact * current - prev;
0106             prev = current;
0107             current = next;
0108         }
0109         value = bessel_j0(x) / current;       // normalization
0110         scale = 1 / scale;
0111     }
0112     value *= factor;
0113
0114     if(tools::max_value<T>() * scale < fabs(value))
0115        return policies::raise_overflow_error<T>("boost::math::bessel_jn<%1%>(%1%,%1%)", nullptr, pol); // LCOV_EXCL_LINE we should never get here!
0116
0117     return value / scale;
0118 }
0119
0120 }}} // namespaces
0121
0122 #endif // BOOST_MATH_BESSEL_JN_HPP
0123