math/special_functions/bessel.hpp

0001 //  Copyright (c) 2007, 2013 John Maddock
0002 //  Copyright Christopher Kormanyos 2013.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006 //
0007 // This header just defines the function entry points, and adds dispatch
0008 // to the right implementation method.  Most of the implementation details
0009 // are in separate headers and copyright Xiaogang Zhang.
0010 //
0011 #ifndef BOOST_MATH_BESSEL_HPP
0012 #define BOOST_MATH_BESSEL_HPP
0013
0014 #ifdef _MSC_VER
0015 #  pragma once
0016 #endif
0017
0018 #include <limits>
0019 #include <boost/math/special_functions/math_fwd.hpp>
0020 #include <boost/math/special_functions/detail/bessel_jy.hpp>
0021 #include <boost/math/special_functions/detail/bessel_jn.hpp>
0022 #include <boost/math/special_functions/detail/bessel_yn.hpp>
0023 #include <boost/math/special_functions/detail/bessel_jy_zero.hpp>
0024 #include <boost/math/special_functions/detail/bessel_ik.hpp>
0025 #include <boost/math/special_functions/detail/bessel_i0.hpp>
0026 #include <boost/math/special_functions/detail/bessel_i1.hpp>
0027 #include <boost/math/special_functions/detail/bessel_kn.hpp>
0028 #include <boost/math/special_functions/detail/iconv.hpp>
0029 #include <boost/math/special_functions/sin_pi.hpp>
0030 #include <boost/math/special_functions/cos_pi.hpp>
0031 #include <boost/math/special_functions/sinc.hpp>
0032 #include <boost/math/special_functions/trunc.hpp>
0033 #include <boost/math/special_functions/round.hpp>
0034 #include <boost/math/tools/rational.hpp>
0035 #include <boost/math/tools/promotion.hpp>
0036 #include <boost/math/tools/series.hpp>
0037 #include <boost/math/tools/roots.hpp>
0038
0039 #ifdef _MSC_VER
0040 # pragma warning(push)
0041 # pragma warning(disable: 6326) // potential comparison of a constant with another constant
0042 #endif
0043
0044 namespace boost{ namespace math{
0045
0046 namespace detail{
0047
0048 template <class T, class Policy>
0049 struct sph_bessel_j_small_z_series_term
0050 {
0051    typedef T result_type;
0052
0053    sph_bessel_j_small_z_series_term(unsigned v_, T x)
0054       : N(0), v(v_)
0055    {
0056       BOOST_MATH_STD_USING
0057       mult = x / 2;
0058       if(v + 3 > max_factorial<T>::value)
0059       {
0060          term = v * log(mult) - boost::math::lgamma(v+1+T(0.5f), Policy());
0061          term = exp(term);
0062       }
0063       else
0064          term = pow(mult, T(v)) / boost::math::tgamma(v+1+T(0.5f), Policy());
0065       mult *= -mult;
0066    }
0067    T operator()()
0068    {
0069       T r = term;
0070       ++N;
0071       term *= mult / (N * T(N + v + 0.5f));
0072       return r;
0073    }
0074 private:
0075    unsigned N;
0076    unsigned v;
0077    T mult;
0078    T term;
0079 };
0080
0081 template <class T, class Policy>
0082 inline T sph_bessel_j_small_z_series(unsigned v, T x, const Policy& pol)
0083 {
0084    BOOST_MATH_STD_USING // ADL of std names
0085    sph_bessel_j_small_z_series_term<T, Policy> s(v, x);
0086    std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0087
0088    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0089
0090    policies::check_series_iterations<T>("boost::math::sph_bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
0091    return result * sqrt(constants::pi<T>() / 4);
0092 }
0093
0094 template <class T, class Policy>
0095 T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol)
0096 {
0097    BOOST_MATH_STD_USING
0098    static const char* function = "boost::math::bessel_j<%1%>(%1%,%1%)";
0099    if(x < 0)
0100    {
0101       // better have integer v:
0102       if (floor(v) == v)
0103       {
0104          // LCOV_EXCL_START
0105          // This branch is hit by multiprecision types only, and is
0106          // tested by our real_concept tests, but thee are excluded from coverage
0107          // due to time constraints.
0108          T r = cyl_bessel_j_imp(v, T(-x), t, pol);
0109          if (iround(v, pol) & 1)
0110             r = -r;
0111          return r;
0112          // LCOV_EXCL_STOP
0113       }
0114       else
0115          return policies::raise_domain_error<T>(function, "Got x = %1%, but we need x >= 0", x, pol);
0116    }
0117
0118    T result_J, y; // LCOV_EXCL_LINE
0119    bessel_jy(v, x, &result_J, &y, need_j, pol);
0120    return result_J;
0121 }
0122
0123 template <class T, class Policy>
0124 inline T cyl_bessel_j_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
0125 {
0126    BOOST_MATH_STD_USING  // ADL of std names.
0127    int ival = detail::iconv(v, pol);
0128    // If v is an integer, use the integer recursion
0129    // method, both that and Steeds method are O(v):
0130    if((0 == v - ival))
0131    {
0132       return bessel_jn(ival, x, pol);
0133    }
0134    return cyl_bessel_j_imp(v, x, bessel_no_int_tag(), pol);
0135 }
0136
0137 template <class T, class Policy>
0138 inline T cyl_bessel_j_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
0139 {
0140    BOOST_MATH_STD_USING
0141    return bessel_jn(v, x, pol);
0142 }
0143
0144 template <class T, class Policy>
0145 inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol)
0146 {
0147    BOOST_MATH_STD_USING // ADL of std names
0148    if(x < 0)
0149       return policies::raise_domain_error<T>("boost::math::sph_bessel_j<%1%>(%1%,%1%)", "Got x = %1%, but function requires x > 0.", x, pol);
0150    //
0151    // Special case, n == 0 resolves down to the sinus cardinal of x:
0152    //
0153    if(n == 0)
0154       return (boost::math::isinf)(x) ? T(0) : boost::math::sinc_pi(x, pol);
0155    //
0156    // Special case for x == 0:
0157    //
0158    if(x == 0)
0159       return 0;
0160    //
0161    // When x is small we may end up with 0/0, use series evaluation
0162    // instead, especially as it converges rapidly:
0163    //
0164    if(x < 1)
0165       return sph_bessel_j_small_z_series(n, x, pol);
0166    //
0167    // Default case is just a naive evaluation of the definition:
0168    //
0169    return sqrt(constants::pi<T>() / (2 * x))
0170       * cyl_bessel_j_imp(T(T(n)+T(0.5f)), x, bessel_no_int_tag(), pol);
0171 }
0172
0173 template <class T, class Policy>
0174 T cyl_bessel_i_imp(T v, T x, const Policy& pol)
0175 {
0176    //
0177    // This handles all the bessel I functions, note that we don't optimise
0178    // for integer v, other than the v = 0 or 1 special cases, as Millers
0179    // algorithm is at least as inefficient as the general case (the general
0180    // case has better error handling too).
0181    //
0182    BOOST_MATH_STD_USING
0183    static const char* function = "boost::math::cyl_bessel_i<%1%>(%1%,%1%)";
0184    if(x < 0)
0185    {
0186       // better have integer v:
0187       if(floor(v) == v)
0188       {
0189          T r = cyl_bessel_i_imp(v, T(-x), pol);
0190          if(iround(v, pol) & 1)
0191             r = -r;
0192          return r;
0193       }
0194       else
0195          return policies::raise_domain_error<T>(function, "Got x = %1%, but we need x >= 0", x, pol);
0196    }
0197    if(x == 0)
0198    {
0199       if(v < 0)
0200          return floor(v) == v ? static_cast<T>(0) : policies::raise_overflow_error<T>(function, nullptr, pol);
0201       return (v == 0) ? static_cast<T>(1) : static_cast<T>(0);
0202    }
0203    if(v == 0.5f)
0204    {
0205       // common special case, note try and avoid overflow in exp(x):
0206       if(x >= tools::log_max_value<T>())
0207       {
0208          T e = exp(x / 2);
0209          return e * (e / sqrt(2 * x * constants::pi<T>()));
0210       }
0211       return sqrt(2 / (x * constants::pi<T>())) * sinh(x);
0212    }
0213    if((policies::digits<T, Policy>() <= 113) && (std::numeric_limits<T>::digits <= 113) && (std::numeric_limits<T>::radix == 2))
0214    {
0215       if(v == 0)
0216       {
0217          return bessel_i0(x);
0218       }
0219       if(v == 1)
0220       {
0221          return bessel_i1(x);
0222       }
0223    }
0224    if((v > 0) && (x / v < 0.25))
0225       return bessel_i_small_z_series(v, x, pol);
0226    T result_I, result_K; // LCOV_EXCL_LINE
0227    bessel_ik(v, x, &result_I, &result_K, need_i, pol);
0228    return result_I;
0229 }
0230
0231 template <class T, class Policy>
0232 inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag& /* t */, const Policy& pol)
0233 {
0234    static const char* function = "boost::math::cyl_bessel_k<%1%>(%1%,%1%)";
0235    BOOST_MATH_STD_USING
0236    if(x < 0)
0237    {
0238       return policies::raise_domain_error<T>(function, "Got x = %1%, but we need x > 0", x, pol);
0239    }
0240    if(x == 0)
0241    {
0242       return (v == 0) ? policies::raise_overflow_error<T>(function, nullptr, pol)
0243          : policies::raise_domain_error<T>(function, "Got x = %1%, but we need x > 0", x, pol);
0244    }
0245    T result_I, result_K; // LCOV_EXCL_LINE
0246    bessel_ik(v, x, &result_I, &result_K, need_k, pol);
0247    return result_K;
0248 }
0249
0250 template <class T, class Policy>
0251 inline T cyl_bessel_k_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
0252 {
0253    BOOST_MATH_STD_USING
0254    if((floor(v) == v))
0255    {
0256       return bessel_kn(itrunc(v), x, pol);
0257    }
0258    return cyl_bessel_k_imp(v, x, bessel_no_int_tag(), pol);
0259 }
0260
0261 template <class T, class Policy>
0262 inline T cyl_bessel_k_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
0263 {
0264    return bessel_kn(v, x, pol);
0265 }
0266
0267 template <class T, class Policy>
0268 inline T cyl_neumann_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol)
0269 {
0270    static const char* function = "boost::math::cyl_neumann<%1%>(%1%,%1%)";
0271
0272    BOOST_MATH_INSTRUMENT_VARIABLE(v);
0273    BOOST_MATH_INSTRUMENT_VARIABLE(x);
0274
0275    if(x <= 0)
0276    {
0277       return (v == 0) && (x == 0) ?
0278          -policies::raise_overflow_error<T>(function, nullptr, pol) // LCOV_EXCL_LINE MP case only here, not tested in code coverage as it takes too long.
0279          : policies::raise_domain_error<T>(function, "Got x = %1%, but result is complex for x <= 0", x, pol);
0280    }
0281    T result_J, y; // LCOV_EXCL_LINE
0282    bessel_jy(v, x, &result_J, &y, need_y, pol);
0283    //
0284    // Post evaluation check for internal overflow during evaluation,
0285    // can occur when x is small and v is large, in which case the result
0286    // is -INF:
0287    //
0288    if(!(boost::math::isfinite)(y))
0289       return -policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE defensive programming?  Might be dead code now...
0290    return y;
0291 }
0292
0293 template <class T, class Policy>
0294 inline T cyl_neumann_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
0295 {
0296    BOOST_MATH_STD_USING
0297
0298    BOOST_MATH_INSTRUMENT_VARIABLE(v);
0299    BOOST_MATH_INSTRUMENT_VARIABLE(x);
0300
0301    if(floor(v) == v)
0302    {
0303       T r = bessel_yn(itrunc(v, pol), x, pol);
0304       BOOST_MATH_INSTRUMENT_VARIABLE(r);
0305       return r;
0306    }
0307    T r = cyl_neumann_imp<T>(v, x, bessel_no_int_tag(), pol);
0308    BOOST_MATH_INSTRUMENT_VARIABLE(r);
0309    return r;
0310 }
0311
0312 template <class T, class Policy>
0313 inline T cyl_neumann_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
0314 {
0315    return bessel_yn(v, x, pol);
0316 }
0317
0318 template <class T, class Policy>
0319 inline T sph_neumann_imp(unsigned v, T x, const Policy& pol)
0320 {
0321    BOOST_MATH_STD_USING // ADL of std names
0322    static const char* function = "boost::math::sph_neumann<%1%>(%1%,%1%)";
0323    //
0324    // Nothing much to do here but check for errors, and
0325    // evaluate the function's definition directly:
0326    //
0327    if(x < 0)
0328       return policies::raise_domain_error<T>(function, "Got x = %1%, but function requires x > 0.", x, pol);
0329
0330    if(x < 2 * tools::min_value<T>())
0331       return -policies::raise_overflow_error<T>(function, nullptr, pol);
0332
0333    T result = cyl_neumann_imp(T(T(v)+0.5f), x, bessel_no_int_tag(), pol);
0334    T tx = sqrt(constants::pi<T>() / (2 * x));
0335
0336    if((tx > 1) && (tools::max_value<T>() / tx < fabs(result)))
0337       return -policies::raise_overflow_error<T>(function, nullptr, pol);
0338
0339    return result * tx;
0340 }
0341
0342 template <class T, class Policy>
0343 inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol)
0344 {
0345    BOOST_MATH_STD_USING // ADL of std names, needed for floor.
0346
0347    static const char* function = "boost::math::cyl_bessel_j_zero<%1%>(%1%, int)";
0348
0349    const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
0350
0351    // Handle non-finite order.
0352    if (!(boost::math::isfinite)(v) )
0353    {
0354      return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
0355    }
0356
0357    // Handle negative rank.
0358    if(m < 0)
0359    {
0360       // Zeros of Jv(x) with negative rank are not defined and requesting one raises a domain error.
0361       return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", static_cast<T>(m), pol);
0362    }
0363
0364    // Get the absolute value of the order.
0365    const bool order_is_negative = (v < 0);
0366    const T vv((!order_is_negative) ? v : T(-v));
0367
0368    // Check if the order is very close to zero or very close to an integer.
0369    const bool order_is_zero    = (vv < half_epsilon);
0370    const bool order_is_integer = ((vv - floor(vv)) < half_epsilon);
0371
0372    if(m == 0)
0373    {
0374       if(order_is_zero)
0375       {
0376          // The zero'th zero of J0(x) is not defined and requesting it raises a domain error.
0377          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);
0378       }
0379
0380       // The zero'th zero of Jv(x) for v < 0 is not defined
0381       // unless the order is a negative integer.
0382       if(order_is_negative && (!order_is_integer))
0383       {
0384          // For non-integer, negative order, requesting the zero'th zero raises a domain error.
0385          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);
0386       }
0387
0388       // The zero'th zero does exist and its value is zero.
0389       return T(0);
0390    }
0391
0392    // Set up the initial guess for the upcoming root-finding.
0393    // If the order is a negative integer, then use the corresponding
0394    // positive integer for the order.
0395    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);
0396
0397    // Select the maximum allowed iterations from the policy.
0398    std::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
0399
0400    const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U));
0401
0402    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
0403    const T jvm =
0404       boost::math::tools::newton_raphson_iterate(
0405          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),
0406          guess_root,
0407          T(guess_root - delta_lo),
0408          T(guess_root + 0.2F),
0409          policies::digits<T, Policy>(),
0410          number_of_iterations);
0411
0412    if(number_of_iterations >= policies::get_max_root_iterations<Policy>())
0413    {
0414       return policies::raise_evaluation_error<T>(function, "Unable to locate root in a reasonable time: Current best guess is %1%", jvm, Policy()); // LCOV_EXCL_LINE
0415    }
0416
0417    return jvm;
0418 }
0419
0420 template <class T, class Policy>
0421 inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
0422 {
0423    BOOST_MATH_STD_USING // ADL of std names, needed for floor.
0424
0425    static const char* function = "boost::math::cyl_neumann_zero<%1%>(%1%, int)";
0426
0427    // Handle non-finite order.
0428    if (!(boost::math::isfinite)(v) )
0429    {
0430      return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
0431    }
0432
0433    // Handle negative rank.
0434    if(m < 0)
0435    {
0436       return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", static_cast<T>(m), pol);
0437    }
0438
0439    const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
0440
0441    // Get the absolute value of the order.
0442    const bool order_is_negative = (v < 0);
0443    const T vv((!order_is_negative) ? v : T(-v));
0444
0445    const bool order_is_integer = ((vv - floor(vv)) < half_epsilon);
0446
0447    // For negative integers, use reflection to positive integer order.
0448    if(order_is_negative && order_is_integer)
0449       return boost::math::detail::cyl_neumann_zero_imp(vv, m, pol);
0450
0451    // Check if the order is very close to a negative half-integer.
0452    const T delta_half_integer(vv - (floor(vv) + 0.5F));
0453
0454    const bool order_is_negative_half_integer =
0455       (order_is_negative && ((delta_half_integer > -half_epsilon) && (delta_half_integer < +half_epsilon)));
0456
0457    // The zero'th zero of Yv(x) for v < 0 is not defined
0458    // unless the order is a negative integer.
0459    if((m == 0) && (!order_is_negative_half_integer))
0460    {
0461       // For non-integer, negative order, requesting the zero'th zero raises a domain error.
0462       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);
0463    }
0464
0465    // For negative half-integers, use the corresponding
0466    // spherical Bessel function of positive half-integer order.
0467    if(order_is_negative_half_integer)
0468       return boost::math::detail::cyl_bessel_j_zero_imp(vv, m, pol);
0469
0470    // Set up the initial guess for the upcoming root-finding.
0471    // If the order is a negative integer, then use the corresponding
0472    // positive integer for the order.
0473    const T guess_root = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess<T, Policy>(v, m, pol);
0474
0475    // Select the maximum allowed iterations from the policy.
0476    std::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
0477
0478    const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U));
0479
0480    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
0481    const T yvm =
0482       boost::math::tools::newton_raphson_iterate(
0483          boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv_and_yv_prime<T, Policy>(v, pol),
0484          guess_root,
0485          T(guess_root - delta_lo),
0486          T(guess_root + 0.2F),
0487          policies::digits<T, Policy>(),
0488          number_of_iterations);
0489
0490    if(number_of_iterations >= policies::get_max_root_iterations<Policy>())
0491    {
0492       return policies::raise_evaluation_error<T>(function, "Unable to locate root in a reasonable time: Current best guess is %1%", yvm, Policy()); //LCOV_EXCL_LINE
0493    }
0494
0495    return yvm;
0496 }
0497
0498 } // namespace detail
0499
0500 template <class T1, class T2, class Policy>
0501 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j(T1 v, T2 x, const Policy& /* pol */)
0502 {
0503    BOOST_FPU_EXCEPTION_GUARD
0504    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0505    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
0506    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0507    typedef typename policies::normalise<
0508       Policy,
0509       policies::promote_float<false>,
0510       policies::promote_double<false>,
0511       policies::discrete_quantile<>,
0512       policies::assert_undefined<> >::type forwarding_policy;
0513    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%)");
0514 }
0515
0516 template <class T1, class T2>
0517 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_j(T1 v, T2 x)
0518 {
0519    return cyl_bessel_j(v, x, policies::policy<>());
0520 }
0521
0522 template <class T, class Policy>
0523 inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel(unsigned v, T x, const Policy& /* pol */)
0524 {
0525    BOOST_FPU_EXCEPTION_GUARD
0526    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
0527    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0528    typedef typename policies::normalise<
0529       Policy,
0530       policies::promote_float<false>,
0531       policies::promote_double<false>,
0532       policies::discrete_quantile<>,
0533       policies::assert_undefined<> >::type forwarding_policy;
0534    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%)");
0535 }
0536
0537 template <class T>
0538 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_bessel(unsigned v, T x)
0539 {
0540    return sph_bessel(v, x, policies::policy<>());
0541 }
0542
0543 template <class T1, class T2, class Policy>
0544 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i(T1 v, T2 x, const Policy& /* pol */)
0545 {
0546    BOOST_FPU_EXCEPTION_GUARD
0547    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0548    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0549    typedef typename policies::normalise<
0550       Policy,
0551       policies::promote_float<false>,
0552       policies::promote_double<false>,
0553       policies::discrete_quantile<>,
0554       policies::assert_undefined<> >::type forwarding_policy;
0555    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%)");
0556 }
0557
0558 template <class T1, class T2>
0559 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_i(T1 v, T2 x)
0560 {
0561    return cyl_bessel_i(v, x, policies::policy<>());
0562 }
0563
0564 template <class T1, class T2, class Policy>
0565 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k(T1 v, T2 x, const Policy& /* pol */)
0566 {
0567    BOOST_FPU_EXCEPTION_GUARD
0568    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0569    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag128 tag_type;
0570    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0571    typedef typename policies::normalise<
0572       Policy,
0573       policies::promote_float<false>,
0574       policies::promote_double<false>,
0575       policies::discrete_quantile<>,
0576       policies::assert_undefined<> >::type forwarding_policy;
0577    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%)");
0578 }
0579
0580 template <class T1, class T2>
0581 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_k(T1 v, T2 x)
0582 {
0583    return cyl_bessel_k(v, x, policies::policy<>());
0584 }
0585
0586 template <class T1, class T2, class Policy>
0587 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann(T1 v, T2 x, const Policy& /* pol */)
0588 {
0589    BOOST_FPU_EXCEPTION_GUARD
0590    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0591    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
0592    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0593    typedef typename policies::normalise<
0594       Policy,
0595       policies::promote_float<false>,
0596       policies::promote_double<false>,
0597       policies::discrete_quantile<>,
0598       policies::assert_undefined<> >::type forwarding_policy;
0599    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%)");
0600 }
0601
0602 template <class T1, class T2>
0603 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_neumann(T1 v, T2 x)
0604 {
0605    return cyl_neumann(v, x, policies::policy<>());
0606 }
0607
0608 template <class T, class Policy>
0609 inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann(unsigned v, T x, const Policy& /* pol */)
0610 {
0611    BOOST_FPU_EXCEPTION_GUARD
0612    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
0613    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0614    typedef typename policies::normalise<
0615       Policy,
0616       policies::promote_float<false>,
0617       policies::promote_double<false>,
0618       policies::discrete_quantile<>,
0619       policies::assert_undefined<> >::type forwarding_policy;
0620    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%)");
0621 }
0622
0623 template <class T>
0624 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann(unsigned v, T x)
0625 {
0626    return sph_neumann(v, x, policies::policy<>());
0627 }
0628
0629 template <class T, class Policy>
0630 inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_bessel_j_zero(T v, int m, const Policy& /* pol */)
0631 {
0632    BOOST_FPU_EXCEPTION_GUARD
0633    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
0634    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0635    typedef typename policies::normalise<
0636       Policy,
0637       policies::promote_float<false>,
0638       policies::promote_double<false>,
0639       policies::discrete_quantile<>,
0640       policies::assert_undefined<> >::type forwarding_policy;
0641
0642    static_assert(    false == std::numeric_limits<T>::is_specialized
0643                            || (   true  == std::numeric_limits<T>::is_specialized
0644                                && false == std::numeric_limits<T>::is_integer),
0645                            "Order must be a floating-point type.");
0646
0647    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%)");
0648 }
0649
0650 template <class T>
0651 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_bessel_j_zero(T v, int m)
0652 {
0653    static_assert(    false == std::numeric_limits<T>::is_specialized
0654                            || (   true  == std::numeric_limits<T>::is_specialized
0655                                && false == std::numeric_limits<T>::is_integer),
0656                            "Order must be a floating-point type.");
0657
0658    return cyl_bessel_j_zero<T, policies::policy<> >(v, m, policies::policy<>());
0659 }
0660
0661 template <class T, class OutputIterator, class Policy>
0662 inline OutputIterator cyl_bessel_j_zero(T v,
0663                               int start_index,
0664                               unsigned number_of_zeros,
0665                               OutputIterator out_it,
0666                               const Policy& pol)
0667 {
0668    static_assert(    false == std::numeric_limits<T>::is_specialized
0669                            || (   true  == std::numeric_limits<T>::is_specialized
0670                                && false == std::numeric_limits<T>::is_integer),
0671                            "Order must be a floating-point type.");
0672
0673    for(int i = 0; i < static_cast<int>(number_of_zeros); ++i)
0674    {
0675       *out_it = boost::math::cyl_bessel_j_zero(v, start_index + i, pol);
0676       ++out_it;
0677    }
0678    return out_it;
0679 }
0680
0681 template <class T, class OutputIterator>
0682 inline OutputIterator cyl_bessel_j_zero(T v,
0683                               int start_index,
0684                               unsigned number_of_zeros,
0685                               OutputIterator out_it)
0686 {
0687    return cyl_bessel_j_zero(v, start_index, number_of_zeros, out_it, policies::policy<>());
0688 }
0689
0690 template <class T, class Policy>
0691 inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_neumann_zero(T v, int m, const Policy& /* pol */)
0692 {
0693    BOOST_FPU_EXCEPTION_GUARD
0694    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
0695    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0696    typedef typename policies::normalise<
0697       Policy,
0698       policies::promote_float<false>,
0699       policies::promote_double<false>,
0700       policies::discrete_quantile<>,
0701       policies::assert_undefined<> >::type forwarding_policy;
0702
0703    static_assert(    false == std::numeric_limits<T>::is_specialized
0704                            || (   true  == std::numeric_limits<T>::is_specialized
0705                                && false == std::numeric_limits<T>::is_integer),
0706                            "Order must be a floating-point type.");
0707
0708    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%)");
0709 }
0710
0711 template <class T>
0712 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_neumann_zero(T v, int m)
0713 {
0714    static_assert(    false == std::numeric_limits<T>::is_specialized
0715                            || (   true  == std::numeric_limits<T>::is_specialized
0716                                && false == std::numeric_limits<T>::is_integer),
0717                            "Order must be a floating-point type.");
0718
0719    return cyl_neumann_zero<T, policies::policy<> >(v, m, policies::policy<>());
0720 }
0721
0722 template <class T, class OutputIterator, class Policy>
0723 inline OutputIterator cyl_neumann_zero(T v,
0724                              int start_index,
0725                              unsigned number_of_zeros,
0726                              OutputIterator out_it,
0727                              const Policy& pol)
0728 {
0729    static_assert(    false == std::numeric_limits<T>::is_specialized
0730                            || (   true  == std::numeric_limits<T>::is_specialized
0731                                && false == std::numeric_limits<T>::is_integer),
0732                            "Order must be a floating-point type.");
0733
0734    for(int i = 0; i < static_cast<int>(number_of_zeros); ++i)
0735    {
0736       *out_it = boost::math::cyl_neumann_zero(v, start_index + i, pol);
0737       ++out_it;
0738    }
0739    return out_it;
0740 }
0741
0742 template <class T, class OutputIterator>
0743 inline OutputIterator cyl_neumann_zero(T v,
0744                              int start_index,
0745                              unsigned number_of_zeros,
0746                              OutputIterator out_it)
0747 {
0748    return cyl_neumann_zero(v, start_index, number_of_zeros, out_it, policies::policy<>());
0749 }
0750
0751 } // namespace math
0752 } // namespace boost
0753
0754 #ifdef _MSC_VER
0755 # pragma warning(pop)
0756 #endif
0757
0758 #endif // BOOST_MATH_BESSEL_HPP
0759
0760