128 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_BESSEL_JN_HPP7#define BOOST_MATH_BESSEL_JN_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/math/tools/config.hpp>14#include <boost/math/tools/assert.hpp>15#include <boost/math/policies/error_handling.hpp>16#include <boost/math/special_functions/gamma.hpp>17#include <boost/math/special_functions/detail/bessel_j0.hpp>18#include <boost/math/special_functions/detail/bessel_j1.hpp>19#include <boost/math/special_functions/detail/bessel_jy.hpp>20#include <boost/math/special_functions/detail/bessel_jy_asym.hpp>21#include <boost/math/special_functions/detail/bessel_jy_series.hpp>22 23// Bessel function of the first kind of integer order24// J_n(z) is the minimal solution25// n < abs(z), forward recurrence stable and usable26// n >= abs(z), forward recurrence unstable, use Miller's algorithm27 28namespace boost { namespace math { namespace detail{29 30template <typename T, typename Policy>31BOOST_MATH_GPU_ENABLED T bessel_jn(int n, T x, const Policy& pol)32{33 T value(0), factor, current, prev, next;34 35 BOOST_MATH_STD_USING36 37 //38 // Reflection has to come first:39 //40 if (n < 0)41 {42 factor = static_cast<T>((n & 0x1) ? -1 : 1); // J_{-n}(z) = (-1)^n J_n(z)43 n = -n;44 }45 else46 {47 factor = 1;48 }49 if(x < 0)50 {51 factor *= (n & 0x1) ? -1 : 1; // J_{n}(-z) = (-1)^n J_n(z)52 x = -x;53 }54 //55 // Special cases:56 //57 if(asymptotic_bessel_large_x_limit(T(n), x))58 return factor * asymptotic_bessel_j_large_x_2<T>(T(n), x, pol);59 if (n == 0)60 {61 return factor * bessel_j0(x);62 }63 if (n == 1)64 {65 return factor * bessel_j1(x);66 }67 68 if (x == 0) // n >= 269 {70 return static_cast<T>(0);71 }72 73 BOOST_MATH_ASSERT(n > 1);74 T scale = 1;75 if (n < abs(x)) // forward recurrence76 {77 prev = bessel_j0(x);78 current = bessel_j1(x);79 policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", static_cast<unsigned>(n), pol);80 for (int k = 1; k < n; k++)81 {82 value = (2 * k * current / x) - prev;83 prev = current;84 current = value;85 }86 }87 else if((x < 1) || (n > x * x / 4) || (x < 5))88 {89 return factor * bessel_j_small_z_series(T(n), x, pol);90 }91 else // backward recurrence92 {93 T fn; int s; // fn = J_(n+1) / J_n94 // |x| <= n, fast convergence for continued fraction CF195 boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol);96 prev = fn;97 current = 1;98 // Check recursion won't go on too far:99 policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", static_cast<unsigned>(n), pol);100 for (int k = n; k > 0; k--)101 {102 T fact = 2 * k / x;103 if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))104 {105 prev /= current;106 scale /= current;107 current = 1;108 }109 next = fact * current - prev;110 prev = current;111 current = next;112 }113 value = bessel_j0(x) / current; // normalization114 scale = 1 / scale;115 }116 value *= factor;117 118 if(tools::max_value<T>() * scale < fabs(value))119 return policies::raise_overflow_error<T>("boost::math::bessel_jn<%1%>(%1%,%1%)", nullptr, pol); // LCOV_EXCL_LINE we should never get here!120 121 return value / scale;122}123 124}}} // namespaces125 126#endif // BOOST_MATH_BESSEL_JN_HPP127 128