]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright (c) 2013 Anton Bikineev |
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_JY_DERIVATIVES_SERIES_HPP | |
7 | #define BOOST_MATH_BESSEL_JY_DERIVATIVES_SERIES_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | namespace boost{ namespace math{ namespace detail{ | |
14 | ||
15 | template <class T, class Policy> | |
16 | struct bessel_j_derivative_small_z_series_term | |
17 | { | |
18 | typedef T result_type; | |
19 | ||
20 | bessel_j_derivative_small_z_series_term(T v_, T x) | |
21 | : N(0), v(v_), term(1), mult(x / 2) | |
22 | { | |
23 | mult *= -mult; | |
24 | // iterate if v == 0; otherwise result of | |
25 | // first term is 0 and tools::sum_series stops | |
26 | if (v == 0) | |
27 | iterate(); | |
28 | } | |
29 | T operator()() | |
30 | { | |
31 | T r = term * (v + 2 * N); | |
32 | iterate(); | |
33 | return r; | |
34 | } | |
35 | private: | |
36 | void iterate() | |
37 | { | |
38 | ++N; | |
39 | term *= mult / (N * (N + v)); | |
40 | } | |
41 | unsigned N; | |
42 | T v; | |
43 | T term; | |
44 | T mult; | |
45 | }; | |
46 | // | |
47 | // Series evaluation for BesselJ'(v, z) as z -> 0. | |
48 | // It's derivative of http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/ | |
49 | // Converges rapidly for all z << v. | |
50 | // | |
51 | template <class T, class Policy> | |
52 | inline T bessel_j_derivative_small_z_series(T v, T x, const Policy& pol) | |
53 | { | |
54 | BOOST_MATH_STD_USING | |
55 | T prefix; | |
56 | if (v < boost::math::max_factorial<T>::value) | |
57 | { | |
58 | prefix = pow(x / 2, v - 1) / 2 / boost::math::tgamma(v + 1, pol); | |
59 | } | |
60 | else | |
61 | { | |
62 | prefix = (v - 1) * log(x / 2) - constants::ln_two<T>() - boost::math::lgamma(v + 1, pol); | |
63 | prefix = exp(prefix); | |
64 | } | |
65 | if (0 == prefix) | |
66 | return prefix; | |
67 | ||
68 | bessel_j_derivative_small_z_series_term<T, Policy> s(v, x); | |
69 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
70 | #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) | |
71 | T zero = 0; | |
72 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero); | |
73 | #else | |
74 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
75 | #endif | |
76 | boost::math::policies::check_series_iterations<T>("boost::math::bessel_j_derivative_small_z_series<%1%>(%1%,%1%)", max_iter, pol); | |
77 | return prefix * result; | |
78 | } | |
79 | ||
80 | template <class T, class Policy> | |
81 | struct bessel_y_derivative_small_z_series_term_a | |
82 | { | |
83 | typedef T result_type; | |
84 | ||
85 | bessel_y_derivative_small_z_series_term_a(T v_, T x) | |
86 | : N(0), v(v_) | |
87 | { | |
88 | mult = x / 2; | |
89 | mult *= -mult; | |
90 | term = 1; | |
91 | } | |
92 | T operator()() | |
93 | { | |
94 | T r = term * (-v + 2 * N); | |
95 | ++N; | |
96 | term *= mult / (N * (N - v)); | |
97 | return r; | |
98 | } | |
99 | private: | |
100 | unsigned N; | |
101 | T v; | |
102 | T mult; | |
103 | T term; | |
104 | }; | |
105 | ||
106 | template <class T, class Policy> | |
107 | struct bessel_y_derivative_small_z_series_term_b | |
108 | { | |
109 | typedef T result_type; | |
110 | ||
111 | bessel_y_derivative_small_z_series_term_b(T v_, T x) | |
112 | : N(0), v(v_) | |
113 | { | |
114 | mult = x / 2; | |
115 | mult *= -mult; | |
116 | term = 1; | |
117 | } | |
118 | T operator()() | |
119 | { | |
120 | T r = term * (v + 2 * N); | |
121 | ++N; | |
122 | term *= mult / (N * (N + v)); | |
123 | return r; | |
124 | } | |
125 | private: | |
126 | unsigned N; | |
127 | T v; | |
128 | T mult; | |
129 | T term; | |
130 | }; | |
131 | // | |
132 | // Series form for BesselY' as z -> 0, | |
133 | // It's derivative of http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/01/0003/ | |
134 | // This series is only useful when the second term is small compared to the first | |
135 | // otherwise we get catestrophic cancellation errors. | |
136 | // | |
137 | // Approximating tgamma(v) by v^v, and assuming |tgamma(-z)| < eps we end up requiring: | |
138 | // eps/2 * v^v(x/2)^-v > (x/2)^v or log(eps/2) > v log((x/2)^2/v) | |
139 | // | |
140 | template <class T, class Policy> | |
141 | inline T bessel_y_derivative_small_z_series(T v, T x, const Policy& pol) | |
142 | { | |
143 | BOOST_MATH_STD_USING | |
144 | static const char* function = "bessel_y_derivative_small_z_series<%1%>(%1%,%1%)"; | |
145 | T prefix; | |
146 | T gam; | |
147 | T p = log(x / 2); | |
148 | T scale = 1; | |
149 | bool need_logs = (v >= boost::math::max_factorial<T>::value) || (boost::math::tools::log_max_value<T>() / v < fabs(p)); | |
150 | if (!need_logs) | |
151 | { | |
152 | gam = boost::math::tgamma(v, pol); | |
153 | p = pow(x / 2, v + 1) * 2; | |
154 | if (boost::math::tools::max_value<T>() * p < gam) | |
155 | { | |
156 | scale /= gam; | |
157 | gam = 1; | |
158 | if (boost::math::tools::max_value<T>() * p < gam) | |
159 | { | |
160 | // This term will overflow to -INF, when combined with the series below it becomes +INF: | |
161 | return boost::math::policies::raise_overflow_error<T>(function, 0, pol); | |
162 | } | |
163 | } | |
164 | prefix = -gam / (boost::math::constants::pi<T>() * p); | |
165 | } | |
166 | else | |
167 | { | |
168 | gam = boost::math::lgamma(v, pol); | |
169 | p = (v + 1) * p + constants::ln_two<T>(); | |
170 | prefix = gam - log(boost::math::constants::pi<T>()) - p; | |
171 | if (boost::math::tools::log_max_value<T>() < prefix) | |
172 | { | |
173 | prefix -= log(boost::math::tools::max_value<T>() / 4); | |
174 | scale /= (boost::math::tools::max_value<T>() / 4); | |
175 | if (boost::math::tools::log_max_value<T>() < prefix) | |
176 | { | |
177 | return boost::math::policies::raise_overflow_error<T>(function, 0, pol); | |
178 | } | |
179 | } | |
180 | prefix = -exp(prefix); | |
181 | } | |
182 | bessel_y_derivative_small_z_series_term_a<T, Policy> s(v, x); | |
183 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
184 | #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) | |
185 | T zero = 0; | |
186 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero); | |
187 | #else | |
188 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
189 | #endif | |
190 | boost::math::policies::check_series_iterations<T>("boost::math::bessel_y_derivative_small_z_series<%1%>(%1%,%1%)", max_iter, pol); | |
191 | result *= prefix; | |
192 | ||
193 | p = pow(x / 2, v - 1) / 2; | |
194 | if (!need_logs) | |
195 | { | |
196 | prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v) * p / boost::math::constants::pi<T>(); | |
197 | } | |
198 | else | |
199 | { | |
200 | int sgn; | |
201 | prefix = boost::math::lgamma(-v, &sgn, pol) + (v - 1) * log(x / 2) - constants::ln_two<T>(); | |
202 | prefix = exp(prefix) * sgn / boost::math::constants::pi<T>(); | |
203 | } | |
204 | bessel_y_derivative_small_z_series_term_b<T, Policy> s2(v, x); | |
205 | max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
206 | #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) | |
207 | T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero); | |
208 | #else | |
209 | T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
210 | #endif | |
211 | result += scale * prefix * b; | |
212 | return result; | |
213 | } | |
214 | ||
215 | // Calculating of BesselY'(v,x) with small x (x < epsilon) and integer x using derivatives | |
216 | // of formulas in http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/ | |
217 | // seems to lose precision. Instead using linear combination of regular Bessel is preferred. | |
218 | ||
219 | }}} // namespaces | |
220 | ||
221 | #endif // BOOST_MATH_BESSEL_JY_DERIVATVIES_SERIES_HPP |