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