Chris@102: //  Copyright (c) 2013 Anton Bikineev
Chris@102: //  Use, modification and distribution are subject to the
Chris@102: //  Boost Software License, Version 1.0. (See accompanying file
Chris@102: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@102: //
Chris@102: #ifndef BOOST_MATH_BESSEL_DERIVATIVES_HPP
Chris@102: #define BOOST_MATH_BESSEL_DERIVATIVES_HPP
Chris@102: 
Chris@102: #ifdef _MSC_VER
Chris@102: #  pragma once
Chris@102: #endif
Chris@102: 
Chris@102: #include <boost/math/special_functions/math_fwd.hpp>
Chris@102: #include <boost/math/special_functions/bessel.hpp>
Chris@102: #include <boost/math/special_functions/detail/bessel_jy_derivatives_asym.hpp>
Chris@102: #include <boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp>
Chris@102: #include <boost/math/special_functions/detail/bessel_derivatives_linear.hpp>
Chris@102: 
Chris@102: namespace boost{ namespace math{
Chris@102: 
Chris@102: namespace detail{
Chris@102: 
Chris@102: template <class Tag, class T, class Policy>
Chris@102: inline T cyl_bessel_j_prime_imp(T v, T x, const Policy& pol)
Chris@102: {
Chris@102:    static const char* const function = "boost::math::cyl_bessel_j_prime<%1%>(%1%,%1%)";
Chris@102:    BOOST_MATH_STD_USING
Chris@102:    //
Chris@102:    // Prevent complex result:
Chris@102:    //
Chris@102:    if (x < 0 && floor(v) != v)
Chris@102:       return boost::math::policies::raise_domain_error<T>(
Chris@102:          function,
Chris@102:          "Got x = %1%, but function requires x >= 0", x, pol);
Chris@102:    //
Chris@102:    // Special cases for x == 0:
Chris@102:    //
Chris@102:    if (x == 0)
Chris@102:    {
Chris@102:       if (v == 1)
Chris@102:          return 0.5;
Chris@102:       else if (v == -1)
Chris@102:          return -0.5;
Chris@102:       else if (floor(v) == v || v > 1)
Chris@102:          return 0;
Chris@102:       else return boost::math::policies::raise_domain_error<T>(
Chris@102:          function,
Chris@102:          "Got x = %1%, but function is indeterminate for this order", x, pol);
Chris@102:    }
Chris@102:    //
Chris@102:    // Special case for large x: use asymptotic expansion:
Chris@102:    //
Chris@102:    if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x))
Chris@102:       return boost::math::detail::asymptotic_bessel_j_derivative_large_x_2(v, x);
Chris@102:    //
Chris@102:    // Special case for small x: use Taylor series:
Chris@102:    //
Chris@102:    if ((abs(x) < 5) || (abs(v) > x * x / 4))
Chris@102:    {
Chris@102:       bool inversed = false;
Chris@102:       if (floor(v) == v && v < 0)
Chris@102:       {
Chris@102:          v = -v;
Chris@102:          if (itrunc(v, pol) & 1)
Chris@102:             inversed = true;
Chris@102:       }
Chris@102:       T r = boost::math::detail::bessel_j_derivative_small_z_series(v, x, pol);
Chris@102:       return inversed ? T(-r) : r;
Chris@102:    }
Chris@102:    //
Chris@102:    // Special case for v == 0:
Chris@102:    //
Chris@102:    if (v == 0)
Chris@102:       return -boost::math::detail::cyl_bessel_j_imp<T>(1, x, Tag(), pol);
Chris@102:    //
Chris@102:    // Default case:
Chris@102:    //
Chris@102:    return boost::math::detail::bessel_j_derivative_linear(v, x, Tag(), pol);
Chris@102: }
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: inline T sph_bessel_j_prime_imp(unsigned v, T x, const Policy& pol)
Chris@102: {
Chris@102:    static const char* const function = "boost::math::sph_bessel_prime<%1%>(%1%,%1%)";
Chris@102:    //
Chris@102:    // Prevent complex result:
Chris@102:    //
Chris@102:    if (x < 0)
Chris@102:       return boost::math::policies::raise_domain_error<T>(
Chris@102:          function,
Chris@102:          "Got x = %1%, but function requires x >= 0.", x, pol);
Chris@102:    //
Chris@102:    // Special case for v == 0:
Chris@102:    //
Chris@102:    if (v == 0)
Chris@102:       return (x == 0) ? boost::math::policies::raise_overflow_error<T>(function, 0, pol)
Chris@102:          : static_cast<T>(-boost::math::detail::sph_bessel_j_imp<T>(1, x, pol));
Chris@102:    //
Chris@102:    // Special case for x == 0 and v > 0:
Chris@102:    //
Chris@102:    if (x == 0)
Chris@102:       return boost::math::policies::raise_domain_error<T>(
Chris@102:          function,
Chris@102:          "Got x = %1%, but function is indeterminate for this order", x, pol);
Chris@102:    //
Chris@102:    // Default case:
Chris@102:    //
Chris@102:    return boost::math::detail::sph_bessel_j_derivative_linear(v, x, pol);
Chris@102: }
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: inline T cyl_bessel_i_prime_imp(T v, T x, const Policy& pol)
Chris@102: {
Chris@102:    static const char* const function = "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)";
Chris@102:    BOOST_MATH_STD_USING
Chris@102:    //
Chris@102:    // Prevent complex result:
Chris@102:    //
Chris@102:    if (x < 0 && floor(v) != v)
Chris@102:       return boost::math::policies::raise_domain_error<T>(
Chris@102:          function,
Chris@102:          "Got x = %1%, but function requires x >= 0", x, pol);
Chris@102:    //
Chris@102:    // Special cases for x == 0:
Chris@102:    //
Chris@102:    if (x == 0)
Chris@102:    {
Chris@102:       if (v == 1 || v == -1)
Chris@102:          return 0.5;
Chris@102:       else if (floor(v) == v || v > 1)
Chris@102:          return 0;
Chris@102:       else return boost::math::policies::raise_domain_error<T>(
Chris@102:          function,
Chris@102:          "Got x = %1%, but function is indeterminate for this order", x, pol);
Chris@102:    }
Chris@102:    //
Chris@102:    // Special case for v == 0:
Chris@102:    //
Chris@102:    if (v == 0)
Chris@102:       return boost::math::detail::cyl_bessel_i_imp<T>(1, x, pol);
Chris@102:    //
Chris@102:    // Default case:
Chris@102:    //
Chris@102:    return boost::math::detail::bessel_i_derivative_linear(v, x, pol);
Chris@102: }
Chris@102: 
Chris@102: template <class Tag, class T, class Policy>
Chris@102: inline T cyl_bessel_k_prime_imp(T v, T x, const Policy& pol)
Chris@102: {
Chris@102:    //
Chris@102:    // Prevent complex and indeterminate results:
Chris@102:    //
Chris@102:    if (x <= 0)
Chris@102:       return boost::math::policies::raise_domain_error<T>(
Chris@102:          "boost::math::cyl_bessel_k_prime<%1%>(%1%,%1%)",
Chris@102:          "Got x = %1%, but function requires x > 0", x, pol);
Chris@102:    //
Chris@102:    // Special case for v == 0:
Chris@102:    //
Chris@102:    if (v == 0)
Chris@102:       return -boost::math::detail::cyl_bessel_k_imp<T>(1, x, Tag(), pol);
Chris@102:    //
Chris@102:    // Default case:
Chris@102:    //
Chris@102:    return boost::math::detail::bessel_k_derivative_linear(v, x, Tag(), pol);
Chris@102: }
Chris@102: 
Chris@102: template <class Tag, class T, class Policy>
Chris@102: inline T cyl_neumann_prime_imp(T v, T x, const Policy& pol)
Chris@102: {
Chris@102:    BOOST_MATH_STD_USING
Chris@102:    //
Chris@102:    // Prevent complex and indeterminate results:
Chris@102:    //
Chris@102:    if (x <= 0)
Chris@102:       return boost::math::policies::raise_domain_error<T>(
Chris@102:          "boost::math::cyl_neumann_prime<%1%>(%1%,%1%)",
Chris@102:          "Got x = %1%, but function requires x > 0", x, pol);
Chris@102:    //
Chris@102:    // Special case for large x: use asymptotic expansion:
Chris@102:    //
Chris@102:    if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x))
Chris@102:       return boost::math::detail::asymptotic_bessel_y_derivative_large_x_2(v, x);
Chris@102:    //
Chris@102:    // Special case for small x: use Taylor series:
Chris@102:    //
Chris@102:    if (v > 0 && floor(v) != v)
Chris@102:    {
Chris@102:       const T eps = boost::math::policies::get_epsilon<T, Policy>();
Chris@102:       if (log(eps / 2) > v * log((x * x) / (v * 4)))
Chris@102:          return boost::math::detail::bessel_y_derivative_small_z_series(v, x, pol);
Chris@102:    }
Chris@102:    //
Chris@102:    // Special case for v == 0:
Chris@102:    //
Chris@102:    if (v == 0)
Chris@102:       return -boost::math::detail::cyl_neumann_imp<T>(1, x, Tag(), pol);
Chris@102:    //
Chris@102:    // Default case:
Chris@102:    //
Chris@102:    return boost::math::detail::bessel_y_derivative_linear(v, x, Tag(), pol);
Chris@102: }
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: inline T sph_neumann_prime_imp(unsigned v, T x, const Policy& pol)
Chris@102: {
Chris@102:    //
Chris@102:    // Prevent complex and indeterminate result:
Chris@102:    //
Chris@102:    if (x <= 0)
Chris@102:       return boost::math::policies::raise_domain_error<T>(
Chris@102:          "boost::math::sph_neumann_prime<%1%>(%1%,%1%)",
Chris@102:          "Got x = %1%, but function requires x > 0.", x, pol);
Chris@102:    //
Chris@102:    // Special case for v == 0:
Chris@102:    //
Chris@102:    if (v == 0)
Chris@102:       return -boost::math::detail::sph_neumann_imp<T>(1, x, pol);
Chris@102:    //
Chris@102:    // Default case:
Chris@102:    //
Chris@102:    return boost::math::detail::sph_neumann_derivative_linear(v, x, pol);
Chris@102: }
Chris@102: 
Chris@102: } // namespace detail
Chris@102: 
Chris@102: template <class T1, class T2, class Policy>
Chris@102: inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j_prime(T1 v, T2 x, const Policy& /* pol */)
Chris@102: {
Chris@102:    BOOST_FPU_EXCEPTION_GUARD
Chris@102:    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
Chris@102:    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
Chris@102:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@102:    typedef typename policies::normalise<
Chris@102:       Policy,
Chris@102:       policies::promote_float<false>,
Chris@102:       policies::promote_double<false>,
Chris@102:       policies::discrete_quantile<>,
Chris@102:       policies::assert_undefined<> >::type forwarding_policy;
Chris@102:    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_prime_imp<tag_type, value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_j_prime<%1%,%1%>(%1%,%1%)");
Chris@102: }
Chris@102: 
Chris@102: template <class T1, class T2>
Chris@102: inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_j_prime(T1 v, T2 x)
Chris@102: {
Chris@102:    return cyl_bessel_j_prime(v, x, policies::policy<>());
Chris@102: }
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel_prime(unsigned v, T x, const Policy& /* pol */)
Chris@102: {
Chris@102:    BOOST_FPU_EXCEPTION_GUARD
Chris@102:    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
Chris@102:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@102:    typedef typename policies::normalise<
Chris@102:       Policy,
Chris@102:       policies::promote_float<false>,
Chris@102:       policies::promote_double<false>,
Chris@102:       policies::discrete_quantile<>,
Chris@102:       policies::assert_undefined<> >::type forwarding_policy;
Chris@102:    return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_bessel_j_prime_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_bessel_j_prime<%1%>(%1%,%1%)");
Chris@102: }
Chris@102: 
Chris@102: template <class T>
Chris@102: inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_bessel_prime(unsigned v, T x)
Chris@102: {
Chris@102:    return sph_bessel_prime(v, x, policies::policy<>());
Chris@102: }
Chris@102: 
Chris@102: template <class T1, class T2, class Policy>
Chris@102: inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i_prime(T1 v, T2 x, const Policy& /* pol */)
Chris@102: {
Chris@102:    BOOST_FPU_EXCEPTION_GUARD
Chris@102:    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
Chris@102:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@102:    typedef typename policies::normalise<
Chris@102:       Policy,
Chris@102:       policies::promote_float<false>,
Chris@102:       policies::promote_double<false>,
Chris@102:       policies::discrete_quantile<>,
Chris@102:       policies::assert_undefined<> >::type forwarding_policy;
Chris@102:    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_prime_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)");
Chris@102: }
Chris@102: 
Chris@102: template <class T1, class T2>
Chris@102: inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_i_prime(T1 v, T2 x)
Chris@102: {
Chris@102:    return cyl_bessel_i_prime(v, x, policies::policy<>());
Chris@102: }
Chris@102: 
Chris@102: template <class T1, class T2, class Policy>
Chris@102: inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k_prime(T1 v, T2 x, const Policy& /* pol */)
Chris@102: {
Chris@102:    BOOST_FPU_EXCEPTION_GUARD
Chris@102:    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
Chris@102:    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
Chris@102:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@102:    typedef typename policies::normalise<
Chris@102:       Policy,
Chris@102:       policies::promote_float<false>,
Chris@102:       policies::promote_double<false>,
Chris@102:       policies::discrete_quantile<>,
Chris@102:       policies::assert_undefined<> >::type forwarding_policy;
Chris@102:    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_prime_imp<tag_type, value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_k_prime<%1%,%1%>(%1%,%1%)");
Chris@102: }
Chris@102: 
Chris@102: template <class T1, class T2>
Chris@102: inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_k_prime(T1 v, T2 x)
Chris@102: {
Chris@102:    return cyl_bessel_k_prime(v, x, policies::policy<>());
Chris@102: }
Chris@102: 
Chris@102: template <class T1, class T2, class Policy>
Chris@102: inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann_prime(T1 v, T2 x, const Policy& /* pol */)
Chris@102: {
Chris@102:    BOOST_FPU_EXCEPTION_GUARD
Chris@102:    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
Chris@102:    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
Chris@102:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@102:    typedef typename policies::normalise<
Chris@102:       Policy,
Chris@102:       policies::promote_float<false>,
Chris@102:       policies::promote_double<false>,
Chris@102:       policies::discrete_quantile<>,
Chris@102:       policies::assert_undefined<> >::type forwarding_policy;
Chris@102:    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_prime_imp<tag_type, value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_neumann_prime<%1%,%1%>(%1%,%1%)");
Chris@102: }
Chris@102: 
Chris@102: template <class T1, class T2>
Chris@102: inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_neumann_prime(T1 v, T2 x)
Chris@102: {
Chris@102:    return cyl_neumann_prime(v, x, policies::policy<>());
Chris@102: }
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann_prime(unsigned v, T x, const Policy& /* pol */)
Chris@102: {
Chris@102:    BOOST_FPU_EXCEPTION_GUARD
Chris@102:    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
Chris@102:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@102:    typedef typename policies::normalise<
Chris@102:       Policy,
Chris@102:       policies::promote_float<false>,
Chris@102:       policies::promote_double<false>,
Chris@102:       policies::discrete_quantile<>,
Chris@102:       policies::assert_undefined<> >::type forwarding_policy;
Chris@102:    return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_neumann_prime_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_neumann_prime<%1%>(%1%,%1%)");
Chris@102: }
Chris@102: 
Chris@102: template <class T>
Chris@102: inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann_prime(unsigned v, T x)
Chris@102: {
Chris@102:    return sph_neumann_prime(v, x, policies::policy<>());
Chris@102: }
Chris@102: 
Chris@102: } // namespace math
Chris@102: } // namespace boost
Chris@102: 
Chris@102: #endif // BOOST_MATH_BESSEL_DERIVATIVES_HPP