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
20namespace boost { namespace math { namespace detail{
21
22template <typename T, typename Policy>
23T 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);
49 }
50 else if (n == 1)
51 {
52 value = bessel_k1(x);
53 }
54 else
55 {
56 prev = bessel_k0(x);
57 current = bessel_k1(x);
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
87

source code of include/boost/math/special_functions/detail/bessel_kn.hpp