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

source code of include/boost/math/special_functions/bessel.hpp