]>
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_DERIVATIVES_HPP | |
7 | #define BOOST_MATH_BESSEL_DERIVATIVES_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | # pragma once | |
11 | #endif | |
12 | ||
13 | #include <boost/math/special_functions/math_fwd.hpp> | |
14 | #include <boost/math/special_functions/bessel.hpp> | |
15 | #include <boost/math/special_functions/detail/bessel_jy_derivatives_asym.hpp> | |
16 | #include <boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp> | |
17 | #include <boost/math/special_functions/detail/bessel_derivatives_linear.hpp> | |
18 | ||
19 | namespace boost{ namespace math{ | |
20 | ||
21 | namespace detail{ | |
22 | ||
23 | template <class Tag, class T, class Policy> | |
24 | inline T cyl_bessel_j_prime_imp(T v, T x, const Policy& pol) | |
25 | { | |
26 | static const char* const function = "boost::math::cyl_bessel_j_prime<%1%>(%1%,%1%)"; | |
27 | BOOST_MATH_STD_USING | |
28 | // | |
29 | // Prevent complex result: | |
30 | // | |
31 | if (x < 0 && floor(v) != v) | |
32 | return boost::math::policies::raise_domain_error<T>( | |
33 | function, | |
34 | "Got x = %1%, but function requires x >= 0", x, pol); | |
35 | // | |
36 | // Special cases for x == 0: | |
37 | // | |
38 | if (x == 0) | |
39 | { | |
40 | if (v == 1) | |
41 | return 0.5; | |
42 | else if (v == -1) | |
43 | return -0.5; | |
44 | else if (floor(v) == v || v > 1) | |
45 | return 0; | |
46 | else return boost::math::policies::raise_domain_error<T>( | |
47 | function, | |
48 | "Got x = %1%, but function is indeterminate for this order", x, pol); | |
49 | } | |
50 | // | |
51 | // Special case for large x: use asymptotic expansion: | |
52 | // | |
53 | if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x)) | |
54 | return boost::math::detail::asymptotic_bessel_j_derivative_large_x_2(v, x); | |
55 | // | |
56 | // Special case for small x: use Taylor series: | |
57 | // | |
58 | if ((abs(x) < 5) || (abs(v) > x * x / 4)) | |
59 | { | |
60 | bool inversed = false; | |
61 | if (floor(v) == v && v < 0) | |
62 | { | |
63 | v = -v; | |
64 | if (itrunc(v, pol) & 1) | |
65 | inversed = true; | |
66 | } | |
67 | T r = boost::math::detail::bessel_j_derivative_small_z_series(v, x, pol); | |
68 | return inversed ? T(-r) : r; | |
69 | } | |
70 | // | |
71 | // Special case for v == 0: | |
72 | // | |
73 | if (v == 0) | |
74 | return -boost::math::detail::cyl_bessel_j_imp<T>(1, x, Tag(), pol); | |
75 | // | |
76 | // Default case: | |
77 | // | |
78 | return boost::math::detail::bessel_j_derivative_linear(v, x, Tag(), pol); | |
79 | } | |
80 | ||
81 | template <class T, class Policy> | |
82 | inline T sph_bessel_j_prime_imp(unsigned v, T x, const Policy& pol) | |
83 | { | |
84 | static const char* const function = "boost::math::sph_bessel_prime<%1%>(%1%,%1%)"; | |
85 | // | |
86 | // Prevent complex result: | |
87 | // | |
88 | if (x < 0) | |
89 | return boost::math::policies::raise_domain_error<T>( | |
90 | function, | |
91 | "Got x = %1%, but function requires x >= 0.", x, pol); | |
92 | // | |
93 | // Special case for v == 0: | |
94 | // | |
95 | if (v == 0) | |
96 | return (x == 0) ? boost::math::policies::raise_overflow_error<T>(function, 0, pol) | |
97 | : static_cast<T>(-boost::math::detail::sph_bessel_j_imp<T>(1, x, pol)); | |
98 | // | |
99 | // Special case for x == 0 and v > 0: | |
100 | // | |
101 | if (x == 0) | |
102 | return boost::math::policies::raise_domain_error<T>( | |
103 | function, | |
104 | "Got x = %1%, but function is indeterminate for this order", x, pol); | |
105 | // | |
106 | // Default case: | |
107 | // | |
108 | return boost::math::detail::sph_bessel_j_derivative_linear(v, x, pol); | |
109 | } | |
110 | ||
111 | template <class T, class Policy> | |
112 | inline T cyl_bessel_i_prime_imp(T v, T x, const Policy& pol) | |
113 | { | |
114 | static const char* const function = "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)"; | |
115 | BOOST_MATH_STD_USING | |
116 | // | |
117 | // Prevent complex result: | |
118 | // | |
119 | if (x < 0 && floor(v) != v) | |
120 | return boost::math::policies::raise_domain_error<T>( | |
121 | function, | |
122 | "Got x = %1%, but function requires x >= 0", x, pol); | |
123 | // | |
124 | // Special cases for x == 0: | |
125 | // | |
126 | if (x == 0) | |
127 | { | |
128 | if (v == 1 || v == -1) | |
129 | return 0.5; | |
130 | else if (floor(v) == v || v > 1) | |
131 | return 0; | |
132 | else return boost::math::policies::raise_domain_error<T>( | |
133 | function, | |
134 | "Got x = %1%, but function is indeterminate for this order", x, pol); | |
135 | } | |
136 | // | |
137 | // Special case for v == 0: | |
138 | // | |
139 | if (v == 0) | |
140 | return boost::math::detail::cyl_bessel_i_imp<T>(1, x, pol); | |
141 | // | |
142 | // Default case: | |
143 | // | |
144 | return boost::math::detail::bessel_i_derivative_linear(v, x, pol); | |
145 | } | |
146 | ||
147 | template <class Tag, class T, class Policy> | |
148 | inline T cyl_bessel_k_prime_imp(T v, T x, const Policy& pol) | |
149 | { | |
150 | // | |
151 | // Prevent complex and indeterminate results: | |
152 | // | |
153 | if (x <= 0) | |
154 | return boost::math::policies::raise_domain_error<T>( | |
155 | "boost::math::cyl_bessel_k_prime<%1%>(%1%,%1%)", | |
156 | "Got x = %1%, but function requires x > 0", x, pol); | |
157 | // | |
158 | // Special case for v == 0: | |
159 | // | |
160 | if (v == 0) | |
161 | return -boost::math::detail::cyl_bessel_k_imp<T>(1, x, Tag(), pol); | |
162 | // | |
163 | // Default case: | |
164 | // | |
165 | return boost::math::detail::bessel_k_derivative_linear(v, x, Tag(), pol); | |
166 | } | |
167 | ||
168 | template <class Tag, class T, class Policy> | |
169 | inline T cyl_neumann_prime_imp(T v, T x, const Policy& pol) | |
170 | { | |
171 | BOOST_MATH_STD_USING | |
172 | // | |
173 | // Prevent complex and indeterminate results: | |
174 | // | |
175 | if (x <= 0) | |
176 | return boost::math::policies::raise_domain_error<T>( | |
177 | "boost::math::cyl_neumann_prime<%1%>(%1%,%1%)", | |
178 | "Got x = %1%, but function requires x > 0", x, pol); | |
179 | // | |
180 | // Special case for large x: use asymptotic expansion: | |
181 | // | |
182 | if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x)) | |
183 | return boost::math::detail::asymptotic_bessel_y_derivative_large_x_2(v, x); | |
184 | // | |
185 | // Special case for small x: use Taylor series: | |
186 | // | |
187 | if (v > 0 && floor(v) != v) | |
188 | { | |
189 | const T eps = boost::math::policies::get_epsilon<T, Policy>(); | |
190 | if (log(eps / 2) > v * log((x * x) / (v * 4))) | |
191 | return boost::math::detail::bessel_y_derivative_small_z_series(v, x, pol); | |
192 | } | |
193 | // | |
194 | // Special case for v == 0: | |
195 | // | |
196 | if (v == 0) | |
197 | return -boost::math::detail::cyl_neumann_imp<T>(1, x, Tag(), pol); | |
198 | // | |
199 | // Default case: | |
200 | // | |
201 | return boost::math::detail::bessel_y_derivative_linear(v, x, Tag(), pol); | |
202 | } | |
203 | ||
204 | template <class T, class Policy> | |
205 | inline T sph_neumann_prime_imp(unsigned v, T x, const Policy& pol) | |
206 | { | |
207 | // | |
208 | // Prevent complex and indeterminate result: | |
209 | // | |
210 | if (x <= 0) | |
211 | return boost::math::policies::raise_domain_error<T>( | |
212 | "boost::math::sph_neumann_prime<%1%>(%1%,%1%)", | |
213 | "Got x = %1%, but function requires x > 0.", x, pol); | |
214 | // | |
215 | // Special case for v == 0: | |
216 | // | |
217 | if (v == 0) | |
218 | return -boost::math::detail::sph_neumann_imp<T>(1, x, pol); | |
219 | // | |
220 | // Default case: | |
221 | // | |
222 | return boost::math::detail::sph_neumann_derivative_linear(v, x, pol); | |
223 | } | |
224 | ||
225 | } // namespace detail | |
226 | ||
227 | template <class T1, class T2, class Policy> | |
228 | inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j_prime(T1 v, T2 x, const Policy& /* pol */) | |
229 | { | |
230 | BOOST_FPU_EXCEPTION_GUARD | |
231 | typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; | |
232 | typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type; | |
233 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
234 | typedef typename policies::normalise< | |
235 | Policy, | |
236 | policies::promote_float<false>, | |
237 | policies::promote_double<false>, | |
238 | policies::discrete_quantile<>, | |
239 | policies::assert_undefined<> >::type forwarding_policy; | |
240 | return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_j_prime<%1%,%1%>(%1%,%1%)"); | |
241 | } | |
242 | ||
243 | template <class T1, class T2> | |
244 | inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_j_prime(T1 v, T2 x) | |
245 | { | |
246 | return cyl_bessel_j_prime(v, x, policies::policy<>()); | |
247 | } | |
248 | ||
249 | template <class T, class Policy> | |
250 | inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel_prime(unsigned v, T x, const Policy& /* pol */) | |
251 | { | |
252 | BOOST_FPU_EXCEPTION_GUARD | |
253 | typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type; | |
254 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
255 | typedef typename policies::normalise< | |
256 | Policy, | |
257 | policies::promote_float<false>, | |
258 | policies::promote_double<false>, | |
259 | policies::discrete_quantile<>, | |
260 | policies::assert_undefined<> >::type forwarding_policy; | |
261 | return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_bessel_j_prime_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_bessel_j_prime<%1%>(%1%,%1%)"); | |
262 | } | |
263 | ||
264 | template <class T> | |
265 | inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_bessel_prime(unsigned v, T x) | |
266 | { | |
267 | return sph_bessel_prime(v, x, policies::policy<>()); | |
268 | } | |
269 | ||
270 | template <class T1, class T2, class Policy> | |
271 | inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i_prime(T1 v, T2 x, const Policy& /* pol */) | |
272 | { | |
273 | BOOST_FPU_EXCEPTION_GUARD | |
274 | typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; | |
275 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
276 | typedef typename policies::normalise< | |
277 | Policy, | |
278 | policies::promote_float<false>, | |
279 | policies::promote_double<false>, | |
280 | policies::discrete_quantile<>, | |
281 | policies::assert_undefined<> >::type forwarding_policy; | |
282 | return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_prime_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)"); | |
283 | } | |
284 | ||
285 | template <class T1, class T2> | |
286 | inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_i_prime(T1 v, T2 x) | |
287 | { | |
288 | return cyl_bessel_i_prime(v, x, policies::policy<>()); | |
289 | } | |
290 | ||
291 | template <class T1, class T2, class Policy> | |
292 | inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k_prime(T1 v, T2 x, const Policy& /* pol */) | |
293 | { | |
294 | BOOST_FPU_EXCEPTION_GUARD | |
295 | typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; | |
296 | typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type; | |
297 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
298 | typedef typename policies::normalise< | |
299 | Policy, | |
300 | policies::promote_float<false>, | |
301 | policies::promote_double<false>, | |
302 | policies::discrete_quantile<>, | |
303 | policies::assert_undefined<> >::type forwarding_policy; | |
304 | return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_k_prime<%1%,%1%>(%1%,%1%)"); | |
305 | } | |
306 | ||
307 | template <class T1, class T2> | |
308 | inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_k_prime(T1 v, T2 x) | |
309 | { | |
310 | return cyl_bessel_k_prime(v, x, policies::policy<>()); | |
311 | } | |
312 | ||
313 | template <class T1, class T2, class Policy> | |
314 | inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann_prime(T1 v, T2 x, const Policy& /* pol */) | |
315 | { | |
316 | BOOST_FPU_EXCEPTION_GUARD | |
317 | typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; | |
318 | typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type; | |
319 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
320 | typedef typename policies::normalise< | |
321 | Policy, | |
322 | policies::promote_float<false>, | |
323 | policies::promote_double<false>, | |
324 | policies::discrete_quantile<>, | |
325 | policies::assert_undefined<> >::type forwarding_policy; | |
326 | return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_neumann_prime<%1%,%1%>(%1%,%1%)"); | |
327 | } | |
328 | ||
329 | template <class T1, class T2> | |
330 | inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_neumann_prime(T1 v, T2 x) | |
331 | { | |
332 | return cyl_neumann_prime(v, x, policies::policy<>()); | |
333 | } | |
334 | ||
335 | template <class T, class Policy> | |
336 | inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann_prime(unsigned v, T x, const Policy& /* pol */) | |
337 | { | |
338 | BOOST_FPU_EXCEPTION_GUARD | |
339 | typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type; | |
340 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
341 | typedef typename policies::normalise< | |
342 | Policy, | |
343 | policies::promote_float<false>, | |
344 | policies::promote_double<false>, | |
345 | policies::discrete_quantile<>, | |
346 | policies::assert_undefined<> >::type forwarding_policy; | |
347 | return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_neumann_prime_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_neumann_prime<%1%>(%1%,%1%)"); | |
348 | } | |
349 | ||
350 | template <class T> | |
351 | inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann_prime(unsigned v, T x) | |
352 | { | |
353 | return sph_neumann_prime(v, x, policies::policy<>()); | |
354 | } | |
355 | ||
356 | } // namespace math | |
357 | } // namespace boost | |
358 | ||
359 | #endif // BOOST_MATH_BESSEL_DERIVATIVES_HPP |