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