]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright (c) 2007, 2013 John Maddock |
2 | // Copyright Christopher Kormanyos 2013. | |
3 | // Use, modification and distribution are subject to the | |
4 | // Boost Software License, Version 1.0. (See accompanying file | |
5 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | // | |
7 | // This header just defines the function entry points, and adds dispatch | |
8 | // to the right implementation method. Most of the implementation details | |
9 | // are in separate headers and copyright Xiaogang Zhang. | |
10 | // | |
11 | #ifndef BOOST_MATH_BESSEL_HPP | |
12 | #define BOOST_MATH_BESSEL_HPP | |
13 | ||
14 | #ifdef _MSC_VER | |
15 | # pragma once | |
16 | #endif | |
17 | ||
18 | #include <limits> | |
19 | #include <boost/math/special_functions/math_fwd.hpp> | |
20 | #include <boost/math/special_functions/detail/bessel_jy.hpp> | |
21 | #include <boost/math/special_functions/detail/bessel_jn.hpp> | |
22 | #include <boost/math/special_functions/detail/bessel_yn.hpp> | |
23 | #include <boost/math/special_functions/detail/bessel_jy_zero.hpp> | |
24 | #include <boost/math/special_functions/detail/bessel_ik.hpp> | |
25 | #include <boost/math/special_functions/detail/bessel_i0.hpp> | |
26 | #include <boost/math/special_functions/detail/bessel_i1.hpp> | |
27 | #include <boost/math/special_functions/detail/bessel_kn.hpp> | |
28 | #include <boost/math/special_functions/detail/iconv.hpp> | |
29 | #include <boost/math/special_functions/sin_pi.hpp> | |
30 | #include <boost/math/special_functions/cos_pi.hpp> | |
31 | #include <boost/math/special_functions/sinc.hpp> | |
32 | #include <boost/math/special_functions/trunc.hpp> | |
33 | #include <boost/math/special_functions/round.hpp> | |
34 | #include <boost/math/tools/rational.hpp> | |
35 | #include <boost/math/tools/promotion.hpp> | |
36 | #include <boost/math/tools/series.hpp> | |
37 | #include <boost/math/tools/roots.hpp> | |
38 | ||
39 | namespace boost{ namespace math{ | |
40 | ||
41 | namespace detail{ | |
42 | ||
43 | template <class T, class Policy> | |
44 | struct sph_bessel_j_small_z_series_term | |
45 | { | |
46 | typedef T result_type; | |
47 | ||
48 | sph_bessel_j_small_z_series_term(unsigned v_, T x) | |
49 | : N(0), v(v_) | |
50 | { | |
51 | BOOST_MATH_STD_USING | |
52 | mult = x / 2; | |
53 | if(v + 3 > max_factorial<T>::value) | |
54 | { | |
55 | term = v * log(mult) - boost::math::lgamma(v+1+T(0.5f), Policy()); | |
56 | term = exp(term); | |
57 | } | |
58 | else | |
59 | term = pow(mult, T(v)) / boost::math::tgamma(v+1+T(0.5f), Policy()); | |
60 | mult *= -mult; | |
61 | } | |
62 | T operator()() | |
63 | { | |
64 | T r = term; | |
65 | ++N; | |
66 | term *= mult / (N * T(N + v + 0.5f)); | |
67 | return r; | |
68 | } | |
69 | private: | |
70 | unsigned N; | |
71 | unsigned v; | |
72 | T mult; | |
73 | T term; | |
74 | }; | |
75 | ||
76 | template <class T, class Policy> | |
77 | inline T sph_bessel_j_small_z_series(unsigned v, T x, const Policy& pol) | |
78 | { | |
79 | BOOST_MATH_STD_USING // ADL of std names | |
80 | sph_bessel_j_small_z_series_term<T, Policy> s(v, x); | |
81 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
20effc67 | 82 | #if BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582)) |
7c673cae FG |
83 | T zero = 0; |
84 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero); | |
85 | #else | |
86 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
87 | #endif | |
88 | policies::check_series_iterations<T>("boost::math::sph_bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol); | |
89 | return result * sqrt(constants::pi<T>() / 4); | |
90 | } | |
91 | ||
92 | template <class T, class Policy> | |
93 | T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol) | |
94 | { | |
95 | BOOST_MATH_STD_USING | |
96 | static const char* function = "boost::math::bessel_j<%1%>(%1%,%1%)"; | |
97 | if(x < 0) | |
98 | { | |
99 | // better have integer v: | |
100 | if(floor(v) == v) | |
101 | { | |
102 | T r = cyl_bessel_j_imp(v, T(-x), t, pol); | |
103 | if(iround(v, pol) & 1) | |
104 | r = -r; | |
105 | return r; | |
106 | } | |
107 | else | |
108 | return policies::raise_domain_error<T>( | |
109 | function, | |
110 | "Got x = %1%, but we need x >= 0", x, pol); | |
111 | } | |
112 | ||
113 | T j, y; | |
114 | bessel_jy(v, x, &j, &y, need_j, pol); | |
115 | return j; | |
116 | } | |
117 | ||
118 | template <class T, class Policy> | |
119 | inline T cyl_bessel_j_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol) | |
120 | { | |
121 | BOOST_MATH_STD_USING // ADL of std names. | |
122 | int ival = detail::iconv(v, pol); | |
123 | // If v is an integer, use the integer recursion | |
124 | // method, both that and Steeds method are O(v): | |
125 | if((0 == v - ival)) | |
126 | { | |
127 | return bessel_jn(ival, x, pol); | |
128 | } | |
129 | return cyl_bessel_j_imp(v, x, bessel_no_int_tag(), pol); | |
130 | } | |
131 | ||
132 | template <class T, class Policy> | |
133 | inline T cyl_bessel_j_imp(int v, T x, const bessel_int_tag&, const Policy& pol) | |
134 | { | |
135 | BOOST_MATH_STD_USING | |
136 | return bessel_jn(v, x, pol); | |
137 | } | |
138 | ||
139 | template <class T, class Policy> | |
140 | inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol) | |
141 | { | |
142 | BOOST_MATH_STD_USING // ADL of std names | |
143 | if(x < 0) | |
144 | return policies::raise_domain_error<T>( | |
145 | "boost::math::sph_bessel_j<%1%>(%1%,%1%)", | |
146 | "Got x = %1%, but function requires x > 0.", x, pol); | |
147 | // | |
148 | // Special case, n == 0 resolves down to the sinus cardinal of x: | |
149 | // | |
150 | if(n == 0) | |
151 | return boost::math::sinc_pi(x, pol); | |
152 | // | |
153 | // Special case for x == 0: | |
154 | // | |
155 | if(x == 0) | |
156 | return 0; | |
157 | // | |
158 | // When x is small we may end up with 0/0, use series evaluation | |
159 | // instead, especially as it converges rapidly: | |
160 | // | |
161 | if(x < 1) | |
162 | return sph_bessel_j_small_z_series(n, x, pol); | |
163 | // | |
164 | // Default case is just a naive evaluation of the definition: | |
165 | // | |
166 | return sqrt(constants::pi<T>() / (2 * x)) | |
167 | * cyl_bessel_j_imp(T(T(n)+T(0.5f)), x, bessel_no_int_tag(), pol); | |
168 | } | |
169 | ||
170 | template <class T, class Policy> | |
171 | T cyl_bessel_i_imp(T v, T x, const Policy& pol) | |
172 | { | |
173 | // | |
174 | // This handles all the bessel I functions, note that we don't optimise | |
175 | // for integer v, other than the v = 0 or 1 special cases, as Millers | |
176 | // algorithm is at least as inefficient as the general case (the general | |
177 | // case has better error handling too). | |
178 | // | |
179 | BOOST_MATH_STD_USING | |
180 | if(x < 0) | |
181 | { | |
182 | // better have integer v: | |
183 | if(floor(v) == v) | |
184 | { | |
185 | T r = cyl_bessel_i_imp(v, T(-x), pol); | |
186 | if(iround(v, pol) & 1) | |
187 | r = -r; | |
188 | return r; | |
189 | } | |
190 | else | |
191 | return policies::raise_domain_error<T>( | |
192 | "boost::math::cyl_bessel_i<%1%>(%1%,%1%)", | |
193 | "Got x = %1%, but we need x >= 0", x, pol); | |
194 | } | |
195 | if(x == 0) | |
196 | { | |
197 | return (v == 0) ? static_cast<T>(1) : static_cast<T>(0); | |
198 | } | |
199 | if(v == 0.5f) | |
200 | { | |
201 | // common special case, note try and avoid overflow in exp(x): | |
202 | if(x >= tools::log_max_value<T>()) | |
203 | { | |
204 | T e = exp(x / 2); | |
205 | return e * (e / sqrt(2 * x * constants::pi<T>())); | |
206 | } | |
207 | return sqrt(2 / (x * constants::pi<T>())) * sinh(x); | |
208 | } | |
20effc67 | 209 | if((policies::digits<T, Policy>() <= 113) && (std::numeric_limits<T>::digits <= 113) && (std::numeric_limits<T>::radix == 2)) |
7c673cae FG |
210 | { |
211 | if(v == 0) | |
212 | { | |
213 | return bessel_i0(x); | |
214 | } | |
215 | if(v == 1) | |
216 | { | |
217 | return bessel_i1(x); | |
218 | } | |
219 | } | |
220 | if((v > 0) && (x / v < 0.25)) | |
221 | return bessel_i_small_z_series(v, x, pol); | |
222 | T I, K; | |
223 | bessel_ik(v, x, &I, &K, need_i, pol); | |
224 | return I; | |
225 | } | |
226 | ||
227 | template <class T, class Policy> | |
228 | inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag& /* t */, const Policy& pol) | |
229 | { | |
230 | static const char* function = "boost::math::cyl_bessel_k<%1%>(%1%,%1%)"; | |
231 | BOOST_MATH_STD_USING | |
232 | if(x < 0) | |
233 | { | |
234 | return policies::raise_domain_error<T>( | |
235 | function, | |
236 | "Got x = %1%, but we need x > 0", x, pol); | |
237 | } | |
238 | if(x == 0) | |
239 | { | |
240 | return (v == 0) ? policies::raise_overflow_error<T>(function, 0, pol) | |
241 | : policies::raise_domain_error<T>( | |
242 | function, | |
243 | "Got x = %1%, but we need x > 0", x, pol); | |
244 | } | |
245 | T I, K; | |
246 | bessel_ik(v, x, &I, &K, need_k, pol); | |
247 | return K; | |
248 | } | |
249 | ||
250 | template <class T, class Policy> | |
251 | inline T cyl_bessel_k_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol) | |
252 | { | |
253 | BOOST_MATH_STD_USING | |
254 | if((floor(v) == v)) | |
255 | { | |
256 | return bessel_kn(itrunc(v), x, pol); | |
257 | } | |
258 | return cyl_bessel_k_imp(v, x, bessel_no_int_tag(), pol); | |
259 | } | |
260 | ||
261 | template <class T, class Policy> | |
262 | inline T cyl_bessel_k_imp(int v, T x, const bessel_int_tag&, const Policy& pol) | |
263 | { | |
264 | return bessel_kn(v, x, pol); | |
265 | } | |
266 | ||
267 | template <class T, class Policy> | |
268 | inline T cyl_neumann_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol) | |
269 | { | |
270 | static const char* function = "boost::math::cyl_neumann<%1%>(%1%,%1%)"; | |
271 | ||
272 | BOOST_MATH_INSTRUMENT_VARIABLE(v); | |
273 | BOOST_MATH_INSTRUMENT_VARIABLE(x); | |
274 | ||
275 | if(x <= 0) | |
276 | { | |
277 | return (v == 0) && (x == 0) ? | |
278 | policies::raise_overflow_error<T>(function, 0, pol) | |
279 | : policies::raise_domain_error<T>( | |
280 | function, | |
281 | "Got x = %1%, but result is complex for x <= 0", x, pol); | |
282 | } | |
283 | T j, y; | |
284 | bessel_jy(v, x, &j, &y, need_y, pol); | |
285 | // | |
286 | // Post evaluation check for internal overflow during evaluation, | |
287 | // can occur when x is small and v is large, in which case the result | |
288 | // is -INF: | |
289 | // | |
290 | if(!(boost::math::isfinite)(y)) | |
291 | return -policies::raise_overflow_error<T>(function, 0, pol); | |
292 | return y; | |
293 | } | |
294 | ||
295 | template <class T, class Policy> | |
296 | inline T cyl_neumann_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol) | |
297 | { | |
298 | BOOST_MATH_STD_USING | |
299 | ||
300 | BOOST_MATH_INSTRUMENT_VARIABLE(v); | |
301 | BOOST_MATH_INSTRUMENT_VARIABLE(x); | |
302 | ||
303 | if(floor(v) == v) | |
304 | { | |
305 | T r = bessel_yn(itrunc(v, pol), x, pol); | |
306 | BOOST_MATH_INSTRUMENT_VARIABLE(r); | |
307 | return r; | |
308 | } | |
309 | T r = cyl_neumann_imp<T>(v, x, bessel_no_int_tag(), pol); | |
310 | BOOST_MATH_INSTRUMENT_VARIABLE(r); | |
311 | return r; | |
312 | } | |
313 | ||
314 | template <class T, class Policy> | |
315 | inline T cyl_neumann_imp(int v, T x, const bessel_int_tag&, const Policy& pol) | |
316 | { | |
317 | return bessel_yn(v, x, pol); | |
318 | } | |
319 | ||
320 | template <class T, class Policy> | |
321 | inline T sph_neumann_imp(unsigned v, T x, const Policy& pol) | |
322 | { | |
323 | BOOST_MATH_STD_USING // ADL of std names | |
324 | static const char* function = "boost::math::sph_neumann<%1%>(%1%,%1%)"; | |
325 | // | |
326 | // Nothing much to do here but check for errors, and | |
327 | // evaluate the function's definition directly: | |
328 | // | |
329 | if(x < 0) | |
330 | return policies::raise_domain_error<T>( | |
331 | function, | |
332 | "Got x = %1%, but function requires x > 0.", x, pol); | |
333 | ||
334 | if(x < 2 * tools::min_value<T>()) | |
335 | return -policies::raise_overflow_error<T>(function, 0, pol); | |
336 | ||
337 | T result = cyl_neumann_imp(T(T(v)+0.5f), x, bessel_no_int_tag(), pol); | |
338 | T tx = sqrt(constants::pi<T>() / (2 * x)); | |
339 | ||
340 | if((tx > 1) && (tools::max_value<T>() / tx < result)) | |
341 | return -policies::raise_overflow_error<T>(function, 0, pol); | |
342 | ||
343 | return result * tx; | |
344 | } | |
345 | ||
346 | template <class T, class Policy> | |
347 | inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol) | |
348 | { | |
349 | BOOST_MATH_STD_USING // ADL of std names, needed for floor. | |
350 | ||
351 | static const char* function = "boost::math::cyl_bessel_j_zero<%1%>(%1%, int)"; | |
352 | ||
353 | const T half_epsilon(boost::math::tools::epsilon<T>() / 2U); | |
354 | ||
355 | // Handle non-finite order. | |
356 | if (!(boost::math::isfinite)(v) ) | |
357 | { | |
358 | return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol); | |
359 | } | |
360 | ||
361 | // Handle negative rank. | |
362 | if(m < 0) | |
363 | { | |
364 | // Zeros of Jv(x) with negative rank are not defined and requesting one raises a domain error. | |
365 | return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", static_cast<T>(m), pol); | |
366 | } | |
367 | ||
368 | // Get the absolute value of the order. | |
369 | const bool order_is_negative = (v < 0); | |
370 | const T vv((!order_is_negative) ? v : T(-v)); | |
371 | ||
372 | // Check if the order is very close to zero or very close to an integer. | |
373 | const bool order_is_zero = (vv < half_epsilon); | |
374 | const bool order_is_integer = ((vv - floor(vv)) < half_epsilon); | |
375 | ||
376 | if(m == 0) | |
377 | { | |
378 | if(order_is_zero) | |
379 | { | |
380 | // The zero'th zero of J0(x) is not defined and requesting it raises a domain error. | |
381 | return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of J0, but the rank must be > 0 !", static_cast<T>(m), pol); | |
382 | } | |
383 | ||
384 | // The zero'th zero of Jv(x) for v < 0 is not defined | |
385 | // unless the order is a negative integer. | |
386 | if(order_is_negative && (!order_is_integer)) | |
387 | { | |
388 | // For non-integer, negative order, requesting the zero'th zero raises a domain error. | |
389 | return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Jv for negative, non-integer order, but the rank must be > 0 !", static_cast<T>(m), pol); | |
390 | } | |
391 | ||
392 | // The zero'th zero does exist and its value is zero. | |
393 | return T(0); | |
394 | } | |
395 | ||
396 | // Set up the initial guess for the upcoming root-finding. | |
397 | // If the order is a negative integer, then use the corresponding | |
398 | // positive integer for the order. | |
399 | const T guess_root = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess<T, Policy>((order_is_integer ? vv : v), m, pol); | |
400 | ||
401 | // Select the maximum allowed iterations from the policy. | |
402 | boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>(); | |
403 | ||
404 | const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U)); | |
405 | ||
406 | // Perform the root-finding using Newton-Raphson iteration from Boost.Math. | |
407 | const T jvm = | |
408 | boost::math::tools::newton_raphson_iterate( | |
409 | boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object_jv_and_jv_prime<T, Policy>((order_is_integer ? vv : v), order_is_zero, pol), | |
410 | guess_root, | |
411 | T(guess_root - delta_lo), | |
412 | T(guess_root + 0.2F), | |
413 | policies::digits<T, Policy>(), | |
414 | number_of_iterations); | |
415 | ||
416 | if(number_of_iterations >= policies::get_max_root_iterations<Policy>()) | |
417 | { | |
418 | return policies::raise_evaluation_error<T>(function, "Unable to locate root in a reasonable time:" | |
419 | " Current best guess is %1%", jvm, Policy()); | |
420 | } | |
421 | ||
422 | return jvm; | |
423 | } | |
424 | ||
425 | template <class T, class Policy> | |
426 | inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol) | |
427 | { | |
428 | BOOST_MATH_STD_USING // ADL of std names, needed for floor. | |
429 | ||
430 | static const char* function = "boost::math::cyl_neumann_zero<%1%>(%1%, int)"; | |
431 | ||
432 | // Handle non-finite order. | |
433 | if (!(boost::math::isfinite)(v) ) | |
434 | { | |
435 | return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol); | |
436 | } | |
437 | ||
438 | // Handle negative rank. | |
439 | if(m < 0) | |
440 | { | |
441 | return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", static_cast<T>(m), pol); | |
442 | } | |
443 | ||
444 | const T half_epsilon(boost::math::tools::epsilon<T>() / 2U); | |
445 | ||
446 | // Get the absolute value of the order. | |
447 | const bool order_is_negative = (v < 0); | |
448 | const T vv((!order_is_negative) ? v : T(-v)); | |
449 | ||
450 | const bool order_is_integer = ((vv - floor(vv)) < half_epsilon); | |
451 | ||
452 | // For negative integers, use reflection to positive integer order. | |
453 | if(order_is_negative && order_is_integer) | |
454 | return boost::math::detail::cyl_neumann_zero_imp(vv, m, pol); | |
455 | ||
456 | // Check if the order is very close to a negative half-integer. | |
457 | const T delta_half_integer(vv - (floor(vv) + 0.5F)); | |
458 | ||
459 | const bool order_is_negative_half_integer = | |
460 | (order_is_negative && ((delta_half_integer > -half_epsilon) && (delta_half_integer < +half_epsilon))); | |
461 | ||
462 | // The zero'th zero of Yv(x) for v < 0 is not defined | |
463 | // unless the order is a negative integer. | |
464 | if((m == 0) && (!order_is_negative_half_integer)) | |
465 | { | |
466 | // For non-integer, negative order, requesting the zero'th zero raises a domain error. | |
467 | return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Yv for negative, non-half-integer order, but the rank must be > 0 !", static_cast<T>(m), pol); | |
468 | } | |
469 | ||
470 | // For negative half-integers, use the corresponding | |
471 | // spherical Bessel function of positive half-integer order. | |
472 | if(order_is_negative_half_integer) | |
473 | return boost::math::detail::cyl_bessel_j_zero_imp(vv, m, pol); | |
474 | ||
475 | // Set up the initial guess for the upcoming root-finding. | |
476 | // If the order is a negative integer, then use the corresponding | |
477 | // positive integer for the order. | |
478 | const T guess_root = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess<T, Policy>(v, m, pol); | |
479 | ||
480 | // Select the maximum allowed iterations from the policy. | |
481 | boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>(); | |
482 | ||
483 | const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U)); | |
484 | ||
485 | // Perform the root-finding using Newton-Raphson iteration from Boost.Math. | |
486 | const T yvm = | |
487 | boost::math::tools::newton_raphson_iterate( | |
488 | boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv_and_yv_prime<T, Policy>(v, pol), | |
489 | guess_root, | |
490 | T(guess_root - delta_lo), | |
491 | T(guess_root + 0.2F), | |
492 | policies::digits<T, Policy>(), | |
493 | number_of_iterations); | |
494 | ||
495 | if(number_of_iterations >= policies::get_max_root_iterations<Policy>()) | |
496 | { | |
497 | return policies::raise_evaluation_error<T>(function, "Unable to locate root in a reasonable time:" | |
498 | " Current best guess is %1%", yvm, Policy()); | |
499 | } | |
500 | ||
501 | return yvm; | |
502 | } | |
503 | ||
504 | } // namespace detail | |
505 | ||
506 | template <class T1, class T2, class Policy> | |
507 | inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j(T1 v, T2 x, const Policy& /* pol */) | |
508 | { | |
509 | BOOST_FPU_EXCEPTION_GUARD | |
510 | typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; | |
511 | typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type; | |
512 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
513 | typedef typename policies::normalise< | |
514 | Policy, | |
515 | policies::promote_float<false>, | |
516 | policies::promote_double<false>, | |
517 | policies::discrete_quantile<>, | |
518 | policies::assert_undefined<> >::type forwarding_policy; | |
519 | return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_imp<value_type>(v, static_cast<value_type>(x), tag_type(), forwarding_policy()), "boost::math::cyl_bessel_j<%1%>(%1%,%1%)"); | |
520 | } | |
521 | ||
522 | template <class T1, class T2> | |
523 | inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_j(T1 v, T2 x) | |
524 | { | |
525 | return cyl_bessel_j(v, x, policies::policy<>()); | |
526 | } | |
527 | ||
528 | template <class T, class Policy> | |
529 | inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel(unsigned v, T x, const Policy& /* pol */) | |
530 | { | |
531 | BOOST_FPU_EXCEPTION_GUARD | |
532 | typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type; | |
533 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
534 | typedef typename policies::normalise< | |
535 | Policy, | |
536 | policies::promote_float<false>, | |
537 | policies::promote_double<false>, | |
538 | policies::discrete_quantile<>, | |
539 | policies::assert_undefined<> >::type forwarding_policy; | |
540 | return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_bessel_j_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_bessel<%1%>(%1%,%1%)"); | |
541 | } | |
542 | ||
543 | template <class T> | |
544 | inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_bessel(unsigned v, T x) | |
545 | { | |
546 | return sph_bessel(v, x, policies::policy<>()); | |
547 | } | |
548 | ||
549 | template <class T1, class T2, class Policy> | |
550 | inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i(T1 v, T2 x, const Policy& /* pol */) | |
551 | { | |
552 | BOOST_FPU_EXCEPTION_GUARD | |
553 | typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; | |
554 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
555 | typedef typename policies::normalise< | |
556 | Policy, | |
557 | policies::promote_float<false>, | |
558 | policies::promote_double<false>, | |
559 | policies::discrete_quantile<>, | |
560 | policies::assert_undefined<> >::type forwarding_policy; | |
561 | return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i<%1%>(%1%,%1%)"); | |
562 | } | |
563 | ||
564 | template <class T1, class T2> | |
565 | inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_i(T1 v, T2 x) | |
566 | { | |
567 | return cyl_bessel_i(v, x, policies::policy<>()); | |
568 | } | |
569 | ||
570 | template <class T1, class T2, class Policy> | |
571 | inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k(T1 v, T2 x, const Policy& /* pol */) | |
572 | { | |
573 | BOOST_FPU_EXCEPTION_GUARD | |
574 | typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; | |
b32b8144 | 575 | typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag128 tag_type; |
7c673cae FG |
576 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
577 | typedef typename policies::normalise< | |
578 | Policy, | |
579 | policies::promote_float<false>, | |
580 | policies::promote_double<false>, | |
581 | policies::discrete_quantile<>, | |
582 | policies::assert_undefined<> >::type forwarding_policy; | |
583 | return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_imp<value_type>(v, static_cast<value_type>(x), tag_type(), forwarding_policy()), "boost::math::cyl_bessel_k<%1%>(%1%,%1%)"); | |
584 | } | |
585 | ||
586 | template <class T1, class T2> | |
587 | inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_k(T1 v, T2 x) | |
588 | { | |
589 | return cyl_bessel_k(v, x, policies::policy<>()); | |
590 | } | |
591 | ||
592 | template <class T1, class T2, class Policy> | |
593 | inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann(T1 v, T2 x, const Policy& /* pol */) | |
594 | { | |
595 | BOOST_FPU_EXCEPTION_GUARD | |
596 | typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; | |
597 | typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type; | |
598 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
599 | typedef typename policies::normalise< | |
600 | Policy, | |
601 | policies::promote_float<false>, | |
602 | policies::promote_double<false>, | |
603 | policies::discrete_quantile<>, | |
604 | policies::assert_undefined<> >::type forwarding_policy; | |
605 | return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_imp<value_type>(v, static_cast<value_type>(x), tag_type(), forwarding_policy()), "boost::math::cyl_neumann<%1%>(%1%,%1%)"); | |
606 | } | |
607 | ||
608 | template <class T1, class T2> | |
609 | inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_neumann(T1 v, T2 x) | |
610 | { | |
611 | return cyl_neumann(v, x, policies::policy<>()); | |
612 | } | |
613 | ||
614 | template <class T, class Policy> | |
615 | inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann(unsigned v, T x, const Policy& /* pol */) | |
616 | { | |
617 | BOOST_FPU_EXCEPTION_GUARD | |
618 | typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type; | |
619 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
620 | typedef typename policies::normalise< | |
621 | Policy, | |
622 | policies::promote_float<false>, | |
623 | policies::promote_double<false>, | |
624 | policies::discrete_quantile<>, | |
625 | policies::assert_undefined<> >::type forwarding_policy; | |
626 | return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_neumann_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_neumann<%1%>(%1%,%1%)"); | |
627 | } | |
628 | ||
629 | template <class T> | |
630 | inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann(unsigned v, T x) | |
631 | { | |
632 | return sph_neumann(v, x, policies::policy<>()); | |
633 | } | |
634 | ||
635 | template <class T, class Policy> | |
636 | inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_bessel_j_zero(T v, int m, const Policy& /* pol */) | |
637 | { | |
638 | BOOST_FPU_EXCEPTION_GUARD | |
639 | typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type; | |
640 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
641 | typedef typename policies::normalise< | |
642 | Policy, | |
643 | policies::promote_float<false>, | |
644 | policies::promote_double<false>, | |
645 | policies::discrete_quantile<>, | |
646 | policies::assert_undefined<> >::type forwarding_policy; | |
647 | ||
648 | BOOST_STATIC_ASSERT_MSG( false == std::numeric_limits<T>::is_specialized | |
649 | || ( true == std::numeric_limits<T>::is_specialized | |
650 | && false == std::numeric_limits<T>::is_integer), | |
651 | "Order must be a floating-point type."); | |
652 | ||
653 | return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_zero_imp<value_type>(v, m, forwarding_policy()), "boost::math::cyl_bessel_j_zero<%1%>(%1%,%1%)"); | |
654 | } | |
655 | ||
656 | template <class T> | |
657 | inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_bessel_j_zero(T v, int m) | |
658 | { | |
659 | BOOST_STATIC_ASSERT_MSG( false == std::numeric_limits<T>::is_specialized | |
660 | || ( true == std::numeric_limits<T>::is_specialized | |
661 | && false == std::numeric_limits<T>::is_integer), | |
662 | "Order must be a floating-point type."); | |
663 | ||
664 | return cyl_bessel_j_zero<T, policies::policy<> >(v, m, policies::policy<>()); | |
665 | } | |
666 | ||
667 | template <class T, class OutputIterator, class Policy> | |
668 | inline OutputIterator cyl_bessel_j_zero(T v, | |
669 | int start_index, | |
670 | unsigned number_of_zeros, | |
671 | OutputIterator out_it, | |
672 | const Policy& pol) | |
673 | { | |
674 | BOOST_STATIC_ASSERT_MSG( false == std::numeric_limits<T>::is_specialized | |
675 | || ( true == std::numeric_limits<T>::is_specialized | |
676 | && false == std::numeric_limits<T>::is_integer), | |
677 | "Order must be a floating-point type."); | |
678 | ||
679 | for(int i = 0; i < static_cast<int>(number_of_zeros); ++i) | |
680 | { | |
681 | *out_it = boost::math::cyl_bessel_j_zero(v, start_index + i, pol); | |
682 | ++out_it; | |
683 | } | |
684 | return out_it; | |
685 | } | |
686 | ||
687 | template <class T, class OutputIterator> | |
688 | inline OutputIterator cyl_bessel_j_zero(T v, | |
689 | int start_index, | |
690 | unsigned number_of_zeros, | |
691 | OutputIterator out_it) | |
692 | { | |
693 | return cyl_bessel_j_zero(v, start_index, number_of_zeros, out_it, policies::policy<>()); | |
694 | } | |
695 | ||
696 | template <class T, class Policy> | |
697 | inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_neumann_zero(T v, int m, const Policy& /* pol */) | |
698 | { | |
699 | BOOST_FPU_EXCEPTION_GUARD | |
700 | typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type; | |
701 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
702 | typedef typename policies::normalise< | |
703 | Policy, | |
704 | policies::promote_float<false>, | |
705 | policies::promote_double<false>, | |
706 | policies::discrete_quantile<>, | |
707 | policies::assert_undefined<> >::type forwarding_policy; | |
708 | ||
709 | BOOST_STATIC_ASSERT_MSG( false == std::numeric_limits<T>::is_specialized | |
710 | || ( true == std::numeric_limits<T>::is_specialized | |
711 | && false == std::numeric_limits<T>::is_integer), | |
712 | "Order must be a floating-point type."); | |
713 | ||
714 | return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_zero_imp<value_type>(v, m, forwarding_policy()), "boost::math::cyl_neumann_zero<%1%>(%1%,%1%)"); | |
715 | } | |
716 | ||
717 | template <class T> | |
718 | inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_neumann_zero(T v, int m) | |
719 | { | |
720 | BOOST_STATIC_ASSERT_MSG( false == std::numeric_limits<T>::is_specialized | |
721 | || ( true == std::numeric_limits<T>::is_specialized | |
722 | && false == std::numeric_limits<T>::is_integer), | |
723 | "Order must be a floating-point type."); | |
724 | ||
725 | return cyl_neumann_zero<T, policies::policy<> >(v, m, policies::policy<>()); | |
726 | } | |
727 | ||
728 | template <class T, class OutputIterator, class Policy> | |
729 | inline OutputIterator cyl_neumann_zero(T v, | |
730 | int start_index, | |
731 | unsigned number_of_zeros, | |
732 | OutputIterator out_it, | |
733 | const Policy& pol) | |
734 | { | |
735 | BOOST_STATIC_ASSERT_MSG( false == std::numeric_limits<T>::is_specialized | |
736 | || ( true == std::numeric_limits<T>::is_specialized | |
737 | && false == std::numeric_limits<T>::is_integer), | |
738 | "Order must be a floating-point type."); | |
739 | ||
740 | for(int i = 0; i < static_cast<int>(number_of_zeros); ++i) | |
741 | { | |
742 | *out_it = boost::math::cyl_neumann_zero(v, start_index + i, pol); | |
743 | ++out_it; | |
744 | } | |
745 | return out_it; | |
746 | } | |
747 | ||
748 | template <class T, class OutputIterator> | |
749 | inline OutputIterator cyl_neumann_zero(T v, | |
750 | int start_index, | |
751 | unsigned number_of_zeros, | |
752 | OutputIterator out_it) | |
753 | { | |
754 | return cyl_neumann_zero(v, start_index, number_of_zeros, out_it, policies::policy<>()); | |
755 | } | |
756 | ||
757 | } // namespace math | |
758 | } // namespace boost | |
759 | ||
760 | #endif // BOOST_MATH_BESSEL_HPP | |
761 | ||
762 |