92 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_KN_HPP7#define BOOST_MATH_BESSEL_KN_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/math/tools/config.hpp>14#include <boost/math/tools/precision.hpp>15#include <boost/math/policies/error_handling.hpp>16#include <boost/math/special_functions/detail/bessel_k0.hpp>17#include <boost/math/special_functions/detail/bessel_k1.hpp>18#include <boost/math/special_functions/sign.hpp>19#include <boost/math/policies/error_handling.hpp>20 21// Modified Bessel function of the second kind of integer order22// K_n(z) is the dominant solution, forward recurrence always OK (though unstable)23 24namespace boost { namespace math { namespace detail{25 26template <typename T, typename Policy>27BOOST_MATH_GPU_ENABLED T bessel_kn(int n, T x, const Policy& pol)28{29 BOOST_MATH_STD_USING30 T value, current, prev;31 32 using namespace boost::math::tools;33 34 constexpr auto function = "boost::math::bessel_kn<%1%>(%1%,%1%)";35 36 if (x < 0)37 {38 return policies::raise_domain_error<T>(function, "Got x = %1%, but argument x must be non-negative, complex number result not supported.", x, pol);39 }40 if (x == 0)41 {42 return (n == 0) ? 43 policies::raise_overflow_error<T>(function, nullptr, pol) 44 : policies::raise_domain_error<T>(function, "Got x = %1%, but argument x must be positive, complex number result not supported.", x, pol);45 }46 47 if (n < 0)48 {49 n = -n; // K_{-n}(z) = K_n(z)50 }51 if (n == 0)52 {53 value = bessel_k0(x);54 }55 else if (n == 1)56 {57 value = bessel_k1(x);58 }59 else60 {61 prev = bessel_k0(x);62 current = bessel_k1(x);63 int k = 1;64 BOOST_MATH_ASSERT(k < n);65 T scale = 1;66 do67 {68 T fact = 2 * k / x;69 if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))70 {71 scale /= current;72 prev /= current;73 current = 1;74 }75 value = fact * current + prev;76 prev = current;77 current = value;78 ++k;79 }80 while(k < n);81 if (tools::max_value<T>() * scale < fabs(value))82 return ((boost::math::signbit)(scale) ? -1 : 1) * sign(value) * policies::raise_overflow_error<T>(function, nullptr, pol);83 value /= scale;84 }85 return value;86}87 88}}} // namespaces89 90#endif // BOOST_MATH_BESSEL_KN_HPP91 92