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