ceph/src/boost/boost/math/special_functions/bessel.hpp

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