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
39namespace boost{ namespace math{
40
41namespace detail{
42
43template <class T, class Policy>
44struct 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 }
69private:
70 unsigned N;
71 unsigned v;
72 T mult;
73 T term;
74};
75
76template <class T, class Policy>
77inline 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>();
82#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
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
92template <class T, class Policy>
93T 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
118template <class T, class Policy>
119inline 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
132template <class T, class Policy>
133inline 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
139template <class T, class Policy>
140inline 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
170template <class T, class Policy>
171T 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 }
209 if(policies::digits<T, Policy>() <= 113)
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
227template <class T, class Policy>
228inline 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
250template <class T, class Policy>
251inline 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
261template <class T, class Policy>
262inline 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
267template <class T, class Policy>
268inline 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
295template <class T, class Policy>
296inline 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
314template <class T, class Policy>
315inline 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
320template <class T, class Policy>
321inline 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
346template <class T, class Policy>
347inline 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
425template <class T, class Policy>
426inline 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
506template <class T1, class T2, class Policy>
507inline 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
522template <class T1, class T2>
523inline 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
528template <class T, class Policy>
529inline 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
543template <class T>
544inline 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
549template <class T1, class T2, class Policy>
550inline 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
564template <class T1, class T2>
565inline 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
570template <class T1, class T2, class Policy>
571inline 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;
575 typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag128 tag_type;
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
586template <class T1, class T2>
587inline 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
592template <class T1, class T2, class Policy>
593inline 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
608template <class T1, class T2>
609inline 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
614template <class T, class Policy>
615inline 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
629template <class T>
630inline 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
635template <class T, class Policy>
636inline 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
656template <class T>
657inline 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
667template <class T, class OutputIterator, class Policy>
668inline 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
687template <class T, class OutputIterator>
688inline 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
696template <class T, class Policy>
697inline 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
717template <class T>
718inline 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
728template <class T, class OutputIterator, class Policy>
729inline 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
748template <class T, class OutputIterator>
749inline 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
763

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