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