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