special_functions/detail/bessel_kn.hpp

0001 //  Copyright (c) 2006 Xiaogang Zhang
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_BESSEL_KN_HPP
0007 #define BOOST_MATH_BESSEL_KN_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/config.hpp>
0014 #include <boost/math/tools/precision.hpp>
0015 #include <boost/math/policies/error_handling.hpp>
0016 #include <boost/math/special_functions/detail/bessel_k0.hpp>
0017 #include <boost/math/special_functions/detail/bessel_k1.hpp>
0018 #include <boost/math/special_functions/sign.hpp>
0019 #include <boost/math/policies/error_handling.hpp>
0020
0021 // Modified Bessel function of the second kind of integer order
0022 // K_n(z) is the dominant solution, forward recurrence always OK (though unstable)
0023
0024 namespace boost { namespace math { namespace detail{
0025
0026 template <typename T, typename Policy>
0027 BOOST_MATH_GPU_ENABLED T bessel_kn(int n, T x, const Policy& pol)
0028 {
0029     BOOST_MATH_STD_USING
0030     T value, current, prev;
0031
0032     using namespace boost::math::tools;
0033
0034     constexpr auto function = "boost::math::bessel_kn<%1%>(%1%,%1%)";
0035
0036     if (x < 0)
0037     {
0038        return policies::raise_domain_error<T>(function, "Got x = %1%, but argument x must be non-negative, complex number result not supported.", x, pol);
0039     }
0040     if (x == 0)
0041     {
0042        return (n == 0) ?
0043           policies::raise_overflow_error<T>(function, nullptr, pol)
0044           : policies::raise_domain_error<T>(function, "Got x = %1%, but argument x must be positive, complex number result not supported.", x, pol);
0045     }
0046
0047     if (n < 0)
0048     {
0049         n = -n;                             // K_{-n}(z) = K_n(z)
0050     }
0051     if (n == 0)
0052     {
0053         value = bessel_k0(x);
0054     }
0055     else if (n == 1)
0056     {
0057         value = bessel_k1(x);
0058     }
0059     else
0060     {
0061        prev = bessel_k0(x);
0062        current = bessel_k1(x);
0063        int k = 1;
0064        BOOST_MATH_ASSERT(k < n);
0065        T scale = 1;
0066        do
0067        {
0068            T fact = 2 * k / x;
0069            if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
0070            {
0071               scale /= current;
0072               prev /= current;
0073               current = 1;
0074            }
0075            value = fact * current + prev;
0076            prev = current;
0077            current = value;
0078            ++k;
0079        }
0080        while(k < n);
0081        if (tools::max_value<T>() * scale < fabs(value))
0082           return ((boost::math::signbit)(scale) ? -1 : 1) * sign(value) * policies::raise_overflow_error<T>(function, nullptr, pol);
0083        value /= scale;
0084     }
0085     return value;
0086 }
0087
0088 }}} // namespaces
0089
0090 #endif // BOOST_MATH_BESSEL_KN_HPP
0091