summaryrefslogtreecommitdiff log msg author committer range
blob: e3a5023c63f217aae874846622d36dd7b292ff36 (plain)
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 ``` ``````// Copyright (c) 2006 Xiaogang Zhang // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_MATH_BESSEL_KN_HPP #define BOOST_MATH_BESSEL_KN_HPP #ifdef _MSC_VER #pragma once #endif #include #include #include // Modified Bessel function of the second kind of integer order // K_n(z) is the dominant solution, forward recurrence always OK (though unstable) namespace boost { namespace math { namespace detail{ template T bessel_kn(int n, T x, const Policy& pol) { BOOST_MATH_STD_USING T value, current, prev; using namespace boost::math::tools; static const char* function = "boost::math::bessel_kn<%1%>(%1%,%1%)"; if (x < 0) { return policies::raise_domain_error(function, "Got x = %1%, but argument x must be non-negative, complex number result not supported.", x, pol); } if (x == 0) { return policies::raise_overflow_error(function, 0, pol); } if (n < 0) { n = -n; // K_{-n}(z) = K_n(z) } if (n == 0) { value = bessel_k0(x, pol); } else if (n == 1) { value = bessel_k1(x, pol); } else { prev = bessel_k0(x, pol); current = bessel_k1(x, pol); int k = 1; BOOST_ASSERT(k < n); T scale = 1; do { T fact = 2 * k / x; if((tools::max_value() - fabs(prev)) / fact < fabs(current)) { scale /= current; prev /= current; current = 1; } value = fact * current + prev; prev = current; current = value; ++k; } while(k < n); if(tools::max_value() * scale < fabs(value)) return sign(scale) * sign(value) * policies::raise_overflow_error(function, 0, pol); value /= scale; } return value; } }}} // namespaces #endif // BOOST_MATH_BESSEL_KN_HPP ``````