]>
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_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 = static_cast<T>((n & 0x1) ? -1 : 1); // Y_{-n}(z) = (-1)^n Y_n(z) | |
49 | n = -n; | |
50 | } | |
51 | else | |
52 | { | |
53 | factor = 1; | |
54 | } | |
55 | if(x < policies::get_epsilon<T, Policy>()) | |
56 | { | |
57 | T scale = 1; | |
58 | value = bessel_yn_small_z(n, x, &scale, pol); | |
59 | if(tools::max_value<T>() * fabs(scale) < fabs(value)) | |
60 | return boost::math::sign(scale) * boost::math::sign(value) * policies::raise_overflow_error<T>(function, 0, pol); | |
61 | value /= scale; | |
62 | } | |
63 | else if(asymptotic_bessel_large_x_limit(n, x)) | |
64 | { | |
20effc67 | 65 | value = factor * asymptotic_bessel_y_large_x_2(static_cast<T>(abs(n)), x, pol); |
7c673cae FG |
66 | } |
67 | else if (n == 0) | |
68 | { | |
69 | value = bessel_y0(x, pol); | |
70 | } | |
71 | else if (n == 1) | |
72 | { | |
73 | value = factor * bessel_y1(x, pol); | |
74 | } | |
75 | else | |
76 | { | |
77 | prev = bessel_y0(x, pol); | |
78 | current = bessel_y1(x, pol); | |
79 | int k = 1; | |
80 | BOOST_ASSERT(k < n); | |
81 | policies::check_series_iterations<T>("boost::math::bessel_y_n<%1%>(%1%,%1%)", n, pol); | |
82 | T mult = 2 * k / x; | |
83 | value = mult * current - prev; | |
84 | prev = current; | |
85 | current = value; | |
86 | ++k; | |
87 | if((mult > 1) && (fabs(current) > 1)) | |
88 | { | |
89 | prev /= current; | |
90 | factor /= current; | |
91 | value /= current; | |
92 | current = 1; | |
93 | } | |
94 | while(k < n) | |
95 | { | |
96 | mult = 2 * k / x; | |
97 | value = mult * current - prev; | |
98 | prev = current; | |
99 | current = value; | |
100 | ++k; | |
101 | } | |
102 | if(fabs(tools::max_value<T>() * factor) < fabs(value)) | |
103 | return sign(value) * sign(factor) * policies::raise_overflow_error<T>(function, 0, pol); | |
104 | value /= factor; | |
105 | } | |
106 | return value; | |
107 | } | |
108 | ||
109 | }}} // namespaces | |
110 | ||
111 | #endif // BOOST_MATH_BESSEL_YN_HPP | |
112 |