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