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             T fact = 2 * k / x;
0079             //
0080             // rescale if we would overflow or underflow:
0081             //
0082             if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))
0083             {
0084                scale /= current;
0085                prev /= current;
0086                current = 1;
0087             }
0088             value = fact * current - prev;
0089             prev = current;
0090             current = value;
0091         }
0092     }
0093     else if((x < 1) || (n > x * x / 4) || (x < 5))
0094     {
0095        return factor * bessel_j_small_z_series(T(n), x, pol);
0096     }
0097     else                                    // backward recurrence
0098     {
0099         T fn; int s;                        // fn = J_(n+1) / J_n
0100         // |x| <= n, fast convergence for continued fraction CF1
0101         boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol);
0102         prev = fn;
0103         current = 1;
0104         // Check recursion won't go on too far:
0105         policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
0106         for (int k = n; k > 0; k--)
0107         {
0108             T fact = 2 * k / x;
0109             if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))
0110             {
0111                prev /= current;
0112                scale /= current;
0113                current = 1;
0114             }
0115             next = fact * current - prev;
0116             prev = current;
0117             current = next;
0118         }
0119         value = bessel_j0(x) / current;       // normalization
0120         scale = 1 / scale;
0121     }
0122     value *= factor;
0123
0124     if(tools::max_value<T>() * scale < fabs(value))
0125        return policies::raise_overflow_error<T>("boost::math::bessel_jn<%1%>(%1%,%1%)", nullptr, pol);
0126
0127     return value / scale;
0128 }
0129
0130 }}} // namespaces
0131
0132 #endif // BOOST_MATH_BESSEL_JN_HPP
0133