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