Chris@16: //  Copyright (c) 2006 Xiaogang Zhang
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_BESSEL_JN_HPP
Chris@16: #define BOOST_MATH_BESSEL_JN_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/math/special_functions/detail/bessel_j0.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_j1.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_jy.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_jy_asym.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_jy_series.hpp>
Chris@16: 
Chris@16: // Bessel function of the first kind of integer order
Chris@16: // J_n(z) is the minimal solution
Chris@16: // n < abs(z), forward recurrence stable and usable
Chris@16: // n >= abs(z), forward recurrence unstable, use Miller's algorithm
Chris@16: 
Chris@16: namespace boost { namespace math { namespace detail{
Chris@16: 
Chris@16: template <typename T, typename Policy>
Chris@16: T bessel_jn(int n, T x, const Policy& pol)
Chris@16: {
Chris@16:     T value(0), factor, current, prev, next;
Chris@16: 
Chris@16:     BOOST_MATH_STD_USING
Chris@16: 
Chris@16:     //
Chris@16:     // Reflection has to come first:
Chris@16:     //
Chris@16:     if (n < 0)
Chris@16:     {
Chris@16:         factor = (n & 0x1) ? -1 : 1;  // J_{-n}(z) = (-1)^n J_n(z)
Chris@16:         n = -n;
Chris@16:     }
Chris@16:     else
Chris@16:     {
Chris@16:         factor = 1;
Chris@16:     }
Chris@16:     if(x < 0)
Chris@16:     {
Chris@16:         factor *= (n & 0x1) ? -1 : 1;  // J_{n}(-z) = (-1)^n J_n(z)
Chris@16:         x = -x;
Chris@16:     }
Chris@16:     //
Chris@16:     // Special cases:
Chris@16:     //
Chris@16:     if (n == 0)
Chris@16:     {
Chris@16:         return factor * bessel_j0(x);
Chris@16:     }
Chris@16:     if (n == 1)
Chris@16:     {
Chris@16:         return factor * bessel_j1(x);
Chris@16:     }
Chris@16: 
Chris@16:     if (x == 0)                             // n >= 2
Chris@16:     {
Chris@16:         return static_cast<T>(0);
Chris@16:     }
Chris@16: 
Chris@16:     if(asymptotic_bessel_large_x_limit(T(n), x))
Chris@16:       return factor * asymptotic_bessel_j_large_x_2<T>(n, x);
Chris@16: 
Chris@16:     BOOST_ASSERT(n > 1);
Chris@16:     T scale = 1;
Chris@16:     if (n < abs(x))                         // forward recurrence
Chris@16:     {
Chris@16:         prev = bessel_j0(x);
Chris@16:         current = bessel_j1(x);
Chris@16:         policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
Chris@16:         for (int k = 1; k < n; k++)
Chris@16:         {
Chris@16:             T fact = 2 * k / x;
Chris@16:             //
Chris@16:             // rescale if we would overflow or underflow:
Chris@16:             //
Chris@16:             if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))
Chris@16:             {
Chris@16:                scale /= current;
Chris@16:                prev /= current;
Chris@16:                current = 1;
Chris@16:             }
Chris@16:             value = fact * current - prev;
Chris@16:             prev = current;
Chris@16:             current = value;
Chris@16:         }
Chris@16:     }
Chris@16:     else if((x < 1) || (n > x * x / 4) || (x < 5))
Chris@16:     {
Chris@16:        return factor * bessel_j_small_z_series(T(n), x, pol);
Chris@16:     }
Chris@16:     else                                    // backward recurrence
Chris@16:     {
Chris@16:         T fn; int s;                        // fn = J_(n+1) / J_n
Chris@16:         // |x| <= n, fast convergence for continued fraction CF1
Chris@16:         boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol);
Chris@16:         prev = fn;
Chris@16:         current = 1;
Chris@16:         // Check recursion won't go on too far:
Chris@16:         policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
Chris@16:         for (int k = n; k > 0; k--)
Chris@16:         {
Chris@16:             T fact = 2 * k / x;
Chris@16:             if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))
Chris@16:             {
Chris@16:                prev /= current;
Chris@16:                scale /= current;
Chris@16:                current = 1;
Chris@16:             }
Chris@16:             next = fact * current - prev;
Chris@16:             prev = current;
Chris@16:             current = next;
Chris@16:         }
Chris@16:         value = bessel_j0(x) / current;       // normalization
Chris@16:         scale = 1 / scale;
Chris@16:     }
Chris@16:     value *= factor;
Chris@16: 
Chris@16:     if(tools::max_value<T>() * scale < fabs(value))
Chris@16:        return policies::raise_overflow_error<T>("boost::math::bessel_jn<%1%>(%1%,%1%)", 0, pol);
Chris@16: 
Chris@16:     return value / scale;
Chris@16: }
Chris@16: 
Chris@16: }}} // namespaces
Chris@16: 
Chris@16: #endif // BOOST_MATH_BESSEL_JN_HPP
Chris@16: