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