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_YN_HPP 7 #define BOOST_MATH_BESSEL_YN_HPP 8 9 #ifdef _MSC_VER 10 #pragma once 11 #endif 12 13 #include <boost/math/special_functions/detail/bessel_y0.hpp> 14 #include <boost/math/special_functions/detail/bessel_y1.hpp> 15 #include <boost/math/special_functions/detail/bessel_jy_series.hpp> 16 #include <boost/math/policies/error_handling.hpp> 17 18 // Bessel function of the second kind of integer order 19 // Y_n(z) is the dominant solution, forward recurrence always OK (though unstable) 20 21 namespace boost { namespace math { namespace detail{ 22 23 template <typename T, typename Policy> 24 T bessel_yn(int n, T x, const Policy& pol) 25 { 26 BOOST_MATH_STD_USING 27 T value, factor, current, prev; 28 29 using namespace boost::math::tools; 30 31 static const char* function = "boost::math::bessel_yn<%1%>(%1%,%1%)"; 32 33 if ((x == 0) && (n == 0)) 34 { 35 return -policies::raise_overflow_error<T>(function, 0, pol); 36 } 37 if (x <= 0) 38 { 39 return policies::raise_domain_error<T>(function, 40 "Got x = %1%, but x must be > 0, complex result not supported.", x, pol); 41 } 42 43 // 44 // Reflection comes first: 45 // 46 if (n < 0) 47 { 48 factor = (n & 0x1) ? -1 : 1; // Y_{-n}(z) = (-1)^n Y_n(z) 49 n = -n; 50 } 51 else 52 { 53 factor = 1; 54 } 55 56 if(x < policies::get_epsilon<T, Policy>()) 57 { 58 T scale = 1; 59 value = bessel_yn_small_z(n, x, &scale, pol); 60 if(tools::max_value<T>() * fabs(scale) < fabs(value)) 61 return boost::math::sign(scale) * boost::math::sign(value) * policies::raise_overflow_error<T>(function, 0, pol); 62 value /= scale; 63 } 64 else if (n == 0) 65 { 66 value = bessel_y0(x, pol); 67 } 68 else if (n == 1) 69 { 70 value = factor * bessel_y1(x, pol); 71 } 72 else 73 { 74 prev = bessel_y0(x, pol); 75 current = bessel_y1(x, pol); 76 int k = 1; 77 BOOST_ASSERT(k < n); 78 policies::check_series_iterations<T>("boost::math::bessel_y_n<%1%>(%1%,%1%)", n, pol); 79 do 80 { 81 T fact = 2 * k / x; 82 if((fact > 1) && ((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))) 83 { 84 prev /= current; 85 factor /= current; 86 current = 1; 87 } 88 value = fact * current - prev; 89 prev = current; 90 current = value; 91 ++k; 92 } 93 while(k < n); 94 if(fabs(tools::max_value<T>() * factor) < fabs(value)) 95 return sign(value) * sign(value) * policies::raise_overflow_error<T>(function, 0, pol); 96 value /= factor; 97 } 98 return value; 99 } 100 101 }}} // namespaces 102 103 #endif // BOOST_MATH_BESSEL_YN_HPP 104 105