Chris@16: //  Copyright (c) 2007, 2013 John Maddock
Chris@16: //  Copyright Christopher Kormanyos 2013.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: //
Chris@16: // This header just defines the function entry points, and adds dispatch
Chris@16: // to the right implementation method.  Most of the implementation details
Chris@16: // are in separate headers and copyright Xiaogang Zhang.
Chris@16: //
Chris@16: #ifndef BOOST_MATH_BESSEL_HPP
Chris@16: #define BOOST_MATH_BESSEL_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #  pragma once
Chris@16: #endif
Chris@16: 
Chris@101: #include <limits>
Chris@101: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_jy.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_jn.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_yn.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_jy_zero.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_ik.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_i0.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_i1.hpp>
Chris@16: #include <boost/math/special_functions/detail/bessel_kn.hpp>
Chris@16: #include <boost/math/special_functions/detail/iconv.hpp>
Chris@16: #include <boost/math/special_functions/sin_pi.hpp>
Chris@16: #include <boost/math/special_functions/cos_pi.hpp>
Chris@16: #include <boost/math/special_functions/sinc.hpp>
Chris@16: #include <boost/math/special_functions/trunc.hpp>
Chris@16: #include <boost/math/special_functions/round.hpp>
Chris@16: #include <boost/math/tools/rational.hpp>
Chris@16: #include <boost/math/tools/promotion.hpp>
Chris@16: #include <boost/math/tools/series.hpp>
Chris@16: #include <boost/math/tools/roots.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace math{
Chris@16: 
Chris@16: namespace detail{
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: struct sph_bessel_j_small_z_series_term
Chris@16: {
Chris@16:    typedef T result_type;
Chris@16: 
Chris@16:    sph_bessel_j_small_z_series_term(unsigned v_, T x)
Chris@16:       : N(0), v(v_)
Chris@16:    {
Chris@16:       BOOST_MATH_STD_USING
Chris@16:       mult = x / 2;
Chris@16:       if(v + 3 > max_factorial<T>::value)
Chris@16:       {
Chris@16:          term = v * log(mult) - boost::math::lgamma(v+1+T(0.5f), Policy());
Chris@16:          term = exp(term);
Chris@16:       }
Chris@16:       else
Chris@16:          term = pow(mult, T(v)) / boost::math::tgamma(v+1+T(0.5f), Policy());
Chris@16:       mult *= -mult;
Chris@16:    }
Chris@16:    T operator()()
Chris@16:    {
Chris@16:       T r = term;
Chris@16:       ++N;
Chris@16:       term *= mult / (N * T(N + v + 0.5f));
Chris@16:       return r;
Chris@16:    }
Chris@16: private:
Chris@16:    unsigned N;
Chris@16:    unsigned v;
Chris@16:    T mult;
Chris@16:    T term;
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T sph_bessel_j_small_z_series(unsigned v, T x, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@16:    sph_bessel_j_small_z_series_term<T, Policy> s(v, x);
Chris@16:    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
Chris@16: #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
Chris@16:    T zero = 0;
Chris@16:    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
Chris@16: #else
Chris@16:    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
Chris@16: #endif
Chris@16:    policies::check_series_iterations<T>("boost::math::sph_bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
Chris@16:    return result * sqrt(constants::pi<T>() / 4);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    static const char* function = "boost::math::bessel_j<%1%>(%1%,%1%)";
Chris@16:    if(x < 0)
Chris@16:    {
Chris@16:       // better have integer v:
Chris@16:       if(floor(v) == v)
Chris@16:       {
Chris@16:          T r = cyl_bessel_j_imp(v, T(-x), t, pol);
Chris@16:          if(iround(v, pol) & 1)
Chris@16:             r = -r;
Chris@16:          return r;
Chris@16:       }
Chris@16:       else
Chris@16:          return policies::raise_domain_error<T>(
Chris@16:             function,
Chris@16:             "Got x = %1%, but we need x >= 0", x, pol);
Chris@16:    }
Chris@16:    
Chris@16:    T j, y;
Chris@16:    bessel_jy(v, x, &j, &y, need_j, pol);
Chris@16:    return j;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T cyl_bessel_j_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // ADL of std names.
Chris@16:    int ival = detail::iconv(v, pol);
Chris@16:    // If v is an integer, use the integer recursion
Chris@16:    // method, both that and Steeds method are O(v):
Chris@16:    if((0 == v - ival))
Chris@16:    {
Chris@16:       return bessel_jn(ival, x, pol);
Chris@16:    }
Chris@16:    return cyl_bessel_j_imp(v, x, bessel_no_int_tag(), pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T cyl_bessel_j_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    return bessel_jn(v, x, pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@16:    if(x < 0)
Chris@16:       return policies::raise_domain_error<T>(
Chris@16:          "boost::math::sph_bessel_j<%1%>(%1%,%1%)",
Chris@16:          "Got x = %1%, but function requires x > 0.", x, pol);
Chris@16:    //
Chris@16:    // Special case, n == 0 resolves down to the sinus cardinal of x:
Chris@16:    //
Chris@16:    if(n == 0)
Chris@16:       return boost::math::sinc_pi(x, pol);
Chris@16:    //
Chris@16:    // Special case for x == 0:
Chris@16:    //
Chris@16:    if(x == 0)
Chris@16:       return 0;
Chris@16:    //
Chris@16:    // When x is small we may end up with 0/0, use series evaluation
Chris@16:    // instead, especially as it converges rapidly:
Chris@16:    //
Chris@16:    if(x < 1)
Chris@16:       return sph_bessel_j_small_z_series(n, x, pol);
Chris@16:    //
Chris@16:    // Default case is just a naive evaluation of the definition:
Chris@16:    //
Chris@16:    return sqrt(constants::pi<T>() / (2 * x)) 
Chris@16:       * cyl_bessel_j_imp(T(T(n)+T(0.5f)), x, bessel_no_int_tag(), pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T cyl_bessel_i_imp(T v, T x, const Policy& pol)
Chris@16: {
Chris@16:    //
Chris@16:    // This handles all the bessel I functions, note that we don't optimise
Chris@16:    // for integer v, other than the v = 0 or 1 special cases, as Millers
Chris@16:    // algorithm is at least as inefficient as the general case (the general
Chris@16:    // case has better error handling too).
Chris@16:    //
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    if(x < 0)
Chris@16:    {
Chris@16:       // better have integer v:
Chris@16:       if(floor(v) == v)
Chris@16:       {
Chris@16:          T r = cyl_bessel_i_imp(v, T(-x), pol);
Chris@16:          if(iround(v, pol) & 1)
Chris@16:             r = -r;
Chris@16:          return r;
Chris@16:       }
Chris@16:       else
Chris@16:          return policies::raise_domain_error<T>(
Chris@16:          "boost::math::cyl_bessel_i<%1%>(%1%,%1%)",
Chris@16:             "Got x = %1%, but we need x >= 0", x, pol);
Chris@16:    }
Chris@16:    if(x == 0)
Chris@16:    {
Chris@101:       return (v == 0) ? static_cast<T>(1) : static_cast<T>(0);
Chris@16:    }
Chris@16:    if(v == 0.5f)
Chris@16:    {
Chris@16:       // common special case, note try and avoid overflow in exp(x):
Chris@16:       if(x >= tools::log_max_value<T>())
Chris@16:       {
Chris@16:          T e = exp(x / 2);
Chris@16:          return e * (e / sqrt(2 * x * constants::pi<T>()));
Chris@16:       }
Chris@16:       return sqrt(2 / (x * constants::pi<T>())) * sinh(x);
Chris@16:    }
Chris@16:    if(policies::digits<T, Policy>() <= 64)
Chris@16:    {
Chris@16:       if(v == 0)
Chris@16:       {
Chris@16:          return bessel_i0(x);
Chris@16:       }
Chris@16:       if(v == 1)
Chris@16:       {
Chris@16:          return bessel_i1(x);
Chris@16:       }
Chris@16:    }
Chris@16:    if((v > 0) && (x / v < 0.25))
Chris@16:       return bessel_i_small_z_series(v, x, pol);
Chris@16:    T I, K;
Chris@16:    bessel_ik(v, x, &I, &K, need_i, pol);
Chris@16:    return I;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag& /* t */, const Policy& pol)
Chris@16: {
Chris@16:    static const char* function = "boost::math::cyl_bessel_k<%1%>(%1%,%1%)";
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    if(x < 0)
Chris@16:    {
Chris@16:       return policies::raise_domain_error<T>(
Chris@16:          function,
Chris@16:          "Got x = %1%, but we need x > 0", x, pol);
Chris@16:    }
Chris@16:    if(x == 0)
Chris@16:    {
Chris@16:       return (v == 0) ? policies::raise_overflow_error<T>(function, 0, pol)
Chris@16:          : policies::raise_domain_error<T>(
Chris@16:          function,
Chris@16:          "Got x = %1%, but we need x > 0", x, pol);
Chris@16:    }
Chris@16:    T I, K;
Chris@16:    bessel_ik(v, x, &I, &K, need_k, pol);
Chris@16:    return K;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T cyl_bessel_k_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    if((floor(v) == v))
Chris@16:    {
Chris@16:       return bessel_kn(itrunc(v), x, pol);
Chris@16:    }
Chris@16:    return cyl_bessel_k_imp(v, x, bessel_no_int_tag(), pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T cyl_bessel_k_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
Chris@16: {
Chris@16:    return bessel_kn(v, x, pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T cyl_neumann_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol)
Chris@16: {
Chris@16:    static const char* function = "boost::math::cyl_neumann<%1%>(%1%,%1%)";
Chris@16: 
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(v);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(x);
Chris@16: 
Chris@16:    if(x <= 0)
Chris@16:    {
Chris@16:       return (v == 0) && (x == 0) ?
Chris@16:          policies::raise_overflow_error<T>(function, 0, pol)
Chris@16:          : policies::raise_domain_error<T>(
Chris@16:                function,
Chris@16:                "Got x = %1%, but result is complex for x <= 0", x, pol);
Chris@16:    }
Chris@16:    T j, y;
Chris@16:    bessel_jy(v, x, &j, &y, need_y, pol);
Chris@16:    // 
Chris@16:    // Post evaluation check for internal overflow during evaluation,
Chris@16:    // can occur when x is small and v is large, in which case the result
Chris@16:    // is -INF:
Chris@16:    //
Chris@16:    if(!(boost::math::isfinite)(y))
Chris@16:       return -policies::raise_overflow_error<T>(function, 0, pol);
Chris@16:    return y;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T cyl_neumann_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(v);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(x);
Chris@16: 
Chris@16:    if(floor(v) == v)
Chris@16:    {
Chris@16:       if(asymptotic_bessel_large_x_limit(v, x))
Chris@16:       {
Chris@16:          T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
Chris@16:          if((v < 0) && (itrunc(v, pol) & 1))
Chris@16:             r = -r;
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(r);
Chris@16:          return r;
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          T r = bessel_yn(itrunc(v, pol), x, pol);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(r);
Chris@16:          return r;
Chris@16:       }
Chris@16:    }
Chris@16:    T r = cyl_neumann_imp<T>(v, x, bessel_no_int_tag(), pol);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(r);
Chris@16:    return r;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T cyl_neumann_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(v);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(x);
Chris@16: 
Chris@16:    if(asymptotic_bessel_large_x_limit(T(v), x))
Chris@16:    {
Chris@16:       T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
Chris@16:       if((v < 0) && (v & 1))
Chris@16:          r = -r;
Chris@16:       return r;
Chris@16:    }
Chris@16:    else
Chris@16:       return bessel_yn(v, x, pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T sph_neumann_imp(unsigned v, T x, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@16:    static const char* function = "boost::math::sph_neumann<%1%>(%1%,%1%)";
Chris@16:    //
Chris@16:    // Nothing much to do here but check for errors, and
Chris@16:    // evaluate the function's definition directly:
Chris@16:    //
Chris@16:    if(x < 0)
Chris@16:       return policies::raise_domain_error<T>(
Chris@16:          function,
Chris@16:          "Got x = %1%, but function requires x > 0.", x, pol);
Chris@16: 
Chris@16:    if(x < 2 * tools::min_value<T>())
Chris@16:       return -policies::raise_overflow_error<T>(function, 0, pol);
Chris@16: 
Chris@16:    T result = cyl_neumann_imp(T(T(v)+0.5f), x, bessel_no_int_tag(), pol);
Chris@16:    T tx = sqrt(constants::pi<T>() / (2 * x));
Chris@16: 
Chris@16:    if((tx > 1) && (tools::max_value<T>() / tx < result))
Chris@16:       return -policies::raise_overflow_error<T>(function, 0, pol);
Chris@16: 
Chris@16:    return result * tx;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names, needed for floor.
Chris@16: 
Chris@16:    static const char* function = "boost::math::cyl_bessel_j_zero<%1%>(%1%, int)";
Chris@16: 
Chris@16:    const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
Chris@16: 
Chris@16:    // Handle non-finite order.
Chris@16:    if (!(boost::math::isfinite)(v) )
Chris@16:    {
Chris@16:      return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
Chris@16:    }
Chris@16: 
Chris@16:    // Handle negative rank.
Chris@16:    if(m < 0)
Chris@16:    {
Chris@16:       // Zeros of Jv(x) with negative rank are not defined and requesting one raises a domain error.
Chris@101:       return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", static_cast<T>(m), pol);
Chris@16:    }
Chris@16: 
Chris@16:    // Get the absolute value of the order.
Chris@16:    const bool order_is_negative = (v < 0);
Chris@16:    const T vv((!order_is_negative) ? v : T(-v));
Chris@16: 
Chris@16:    // Check if the order is very close to zero or very close to an integer.
Chris@16:    const bool order_is_zero    = (vv < half_epsilon);
Chris@16:    const bool order_is_integer = ((vv - floor(vv)) < half_epsilon);
Chris@16: 
Chris@16:    if(m == 0)
Chris@16:    {
Chris@16:       if(order_is_zero)
Chris@16:       {
Chris@16:          // The zero'th zero of J0(x) is not defined and requesting it raises a domain error.
Chris@101:          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);
Chris@16:       }
Chris@16: 
Chris@16:       // The zero'th zero of Jv(x) for v < 0 is not defined
Chris@16:       // unless the order is a negative integer.
Chris@16:       if(order_is_negative && (!order_is_integer))
Chris@16:       {
Chris@16:          // For non-integer, negative order, requesting the zero'th zero raises a domain error.
Chris@101:          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);
Chris@16:       }
Chris@16: 
Chris@16:       // The zero'th zero does exist and its value is zero.
Chris@16:       return T(0);
Chris@16:    }
Chris@16: 
Chris@16:    // Set up the initial guess for the upcoming root-finding.
Chris@16:    // If the order is a negative integer, then use the corresponding
Chris@16:    // positive integer for the order.
Chris@16:    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);
Chris@16: 
Chris@16:    // Select the maximum allowed iterations from the policy.
Chris@16:    boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
Chris@16: 
Chris@16:    const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U));
Chris@16: 
Chris@16:    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
Chris@16:    const T jvm =
Chris@16:       boost::math::tools::newton_raphson_iterate(
Chris@16:          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),
Chris@16:          guess_root,
Chris@16:          T(guess_root - delta_lo),
Chris@16:          T(guess_root + 0.2F),
Chris@101:          policies::digits<T, Policy>(),
Chris@16:          number_of_iterations);
Chris@16: 
Chris@16:    if(number_of_iterations >= policies::get_max_root_iterations<Policy>())
Chris@16:    {
Chris@101:       return policies::raise_evaluation_error<T>(function, "Unable to locate root in a reasonable time:"
Chris@16:          "  Current best guess is %1%", jvm, Policy());
Chris@16:    }
Chris@16: 
Chris@16:    return jvm;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names, needed for floor.
Chris@16: 
Chris@16:    static const char* function = "boost::math::cyl_neumann_zero<%1%>(%1%, int)";
Chris@16: 
Chris@16:    // Handle non-finite order.
Chris@16:    if (!(boost::math::isfinite)(v) )
Chris@16:    {
Chris@16:      return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
Chris@16:    }
Chris@16: 
Chris@16:    // Handle negative rank.
Chris@16:    if(m < 0)
Chris@16:    {
Chris@101:       return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", static_cast<T>(m), pol);
Chris@16:    }
Chris@16: 
Chris@16:    const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
Chris@16: 
Chris@16:    // Get the absolute value of the order.
Chris@16:    const bool order_is_negative = (v < 0);
Chris@16:    const T vv((!order_is_negative) ? v : T(-v));
Chris@16: 
Chris@16:    const bool order_is_integer = ((vv - floor(vv)) < half_epsilon);
Chris@16: 
Chris@16:    // For negative integers, use reflection to positive integer order.
Chris@16:    if(order_is_negative && order_is_integer)
Chris@16:       return boost::math::detail::cyl_neumann_zero_imp(vv, m, pol);
Chris@16: 
Chris@16:    // Check if the order is very close to a negative half-integer.
Chris@16:    const T delta_half_integer(vv - (floor(vv) + 0.5F));
Chris@16: 
Chris@16:    const bool order_is_negative_half_integer =
Chris@16:       (order_is_negative && ((delta_half_integer > -half_epsilon) && (delta_half_integer < +half_epsilon)));
Chris@16: 
Chris@16:    // The zero'th zero of Yv(x) for v < 0 is not defined
Chris@16:    // unless the order is a negative integer.
Chris@16:    if((m == 0) && (!order_is_negative_half_integer))
Chris@16:    {
Chris@16:       // For non-integer, negative order, requesting the zero'th zero raises a domain error.
Chris@101:       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);
Chris@16:    }
Chris@16: 
Chris@16:    // For negative half-integers, use the corresponding
Chris@16:    // spherical Bessel function of positive half-integer order.
Chris@16:    if(order_is_negative_half_integer)
Chris@16:       return boost::math::detail::cyl_bessel_j_zero_imp(vv, m, pol);
Chris@16: 
Chris@16:    // Set up the initial guess for the upcoming root-finding.
Chris@16:    // If the order is a negative integer, then use the corresponding
Chris@16:    // positive integer for the order.
Chris@16:    const T guess_root = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess<T, Policy>(v, m, pol);
Chris@16: 
Chris@16:    // Select the maximum allowed iterations from the policy.
Chris@16:    boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
Chris@16: 
Chris@16:    const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U));
Chris@16: 
Chris@16:    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
Chris@16:    const T yvm =
Chris@16:       boost::math::tools::newton_raphson_iterate(
Chris@16:          boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv_and_yv_prime<T, Policy>(v, pol),
Chris@16:          guess_root,
Chris@16:          T(guess_root - delta_lo),
Chris@16:          T(guess_root + 0.2F),
Chris@101:          policies::digits<T, Policy>(),
Chris@16:          number_of_iterations);
Chris@16: 
Chris@16:    if(number_of_iterations >= policies::get_max_root_iterations<Policy>())
Chris@16:    {
Chris@101:       return policies::raise_evaluation_error<T>(function, "Unable to locate root in a reasonable time:"
Chris@16:          "  Current best guess is %1%", yvm, Policy());
Chris@16:    }
Chris@16: 
Chris@16:    return yvm;
Chris@16: }
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j(T1 v, T2 x, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
Chris@16:    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_j(T1 v, T2 x)
Chris@16: {
Chris@16:    return cyl_bessel_j(v, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel(unsigned v, T x, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_bessel(unsigned v, T x)
Chris@16: {
Chris@16:    return sph_bessel(v, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i(T1 v, T2 x, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_i(T1 v, T2 x)
Chris@16: {
Chris@16:    return cyl_bessel_i(v, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k(T1 v, T2 x, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
Chris@16:    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_k(T1 v, T2 x)
Chris@16: {
Chris@16:    return cyl_bessel_k(v, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann(T1 v, T2 x, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
Chris@16:    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_neumann(T1 v, T2 x)
Chris@16: {
Chris@16:    return cyl_neumann(v, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann(unsigned v, T x, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann(unsigned v, T x)
Chris@16: {
Chris@16:    return sph_neumann(v, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_bessel_j_zero(T v, int m, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@101: 
Chris@101:    BOOST_STATIC_ASSERT_MSG(    false == std::numeric_limits<T>::is_specialized
Chris@101:                            || (   true  == std::numeric_limits<T>::is_specialized
Chris@101:                                && false == std::numeric_limits<T>::is_integer),
Chris@101:                            "Order must be a floating-point type.");
Chris@101: 
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_bessel_j_zero(T v, int m)
Chris@16: {
Chris@101:    BOOST_STATIC_ASSERT_MSG(    false == std::numeric_limits<T>::is_specialized
Chris@101:                            || (   true  == std::numeric_limits<T>::is_specialized
Chris@101:                                && false == std::numeric_limits<T>::is_integer),
Chris@101:                            "Order must be a floating-point type.");
Chris@101: 
Chris@16:    return cyl_bessel_j_zero<T, policies::policy<> >(v, m, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class OutputIterator, class Policy>
Chris@16: inline OutputIterator cyl_bessel_j_zero(T v,
Chris@16:                               int start_index,
Chris@16:                               unsigned number_of_zeros,
Chris@16:                               OutputIterator out_it,
Chris@16:                               const Policy& pol)
Chris@16: {
Chris@101:    BOOST_STATIC_ASSERT_MSG(    false == std::numeric_limits<T>::is_specialized
Chris@101:                            || (   true  == std::numeric_limits<T>::is_specialized
Chris@101:                                && false == std::numeric_limits<T>::is_integer),
Chris@101:                            "Order must be a floating-point type.");
Chris@101: 
Chris@101:    for(int i = 0; i < static_cast<int>(number_of_zeros); ++i)
Chris@16:    {
Chris@16:       *out_it = boost::math::cyl_bessel_j_zero(v, start_index + i, pol);
Chris@16:       ++out_it;
Chris@16:    }
Chris@16:    return out_it;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class OutputIterator>
Chris@16: inline OutputIterator cyl_bessel_j_zero(T v,
Chris@16:                               int start_index,
Chris@16:                               unsigned number_of_zeros,
Chris@16:                               OutputIterator out_it)
Chris@16: {
Chris@16:    return cyl_bessel_j_zero(v, start_index, number_of_zeros, out_it, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_neumann_zero(T v, int m, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@101: 
Chris@101:    BOOST_STATIC_ASSERT_MSG(    false == std::numeric_limits<T>::is_specialized
Chris@101:                            || (   true  == std::numeric_limits<T>::is_specialized
Chris@101:                                && false == std::numeric_limits<T>::is_integer),
Chris@101:                            "Order must be a floating-point type.");
Chris@101: 
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_neumann_zero(T v, int m)
Chris@16: {
Chris@101:    BOOST_STATIC_ASSERT_MSG(    false == std::numeric_limits<T>::is_specialized
Chris@101:                            || (   true  == std::numeric_limits<T>::is_specialized
Chris@101:                                && false == std::numeric_limits<T>::is_integer),
Chris@101:                            "Order must be a floating-point type.");
Chris@101: 
Chris@16:    return cyl_neumann_zero<T, policies::policy<> >(v, m, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class OutputIterator, class Policy>
Chris@16: inline OutputIterator cyl_neumann_zero(T v,
Chris@16:                              int start_index,
Chris@16:                              unsigned number_of_zeros,
Chris@16:                              OutputIterator out_it,
Chris@16:                              const Policy& pol)
Chris@16: {
Chris@101:    BOOST_STATIC_ASSERT_MSG(    false == std::numeric_limits<T>::is_specialized
Chris@101:                            || (   true  == std::numeric_limits<T>::is_specialized
Chris@101:                                && false == std::numeric_limits<T>::is_integer),
Chris@101:                            "Order must be a floating-point type.");
Chris@101: 
Chris@101:    for(int i = 0; i < static_cast<int>(number_of_zeros); ++i)
Chris@16:    {
Chris@16:       *out_it = boost::math::cyl_neumann_zero(v, start_index + i, pol);
Chris@16:       ++out_it;
Chris@16:    }
Chris@16:    return out_it;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class OutputIterator>
Chris@16: inline OutputIterator cyl_neumann_zero(T v,
Chris@16:                              int start_index,
Chris@16:                              unsigned number_of_zeros,
Chris@16:                              OutputIterator out_it)
Chris@16: {
Chris@16:    return cyl_neumann_zero(v, start_index, number_of_zeros, out_it, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_MATH_BESSEL_HPP
Chris@16: 
Chris@16: