math/special_functions/bessel_prime.hpp

0001 //  Copyright (c) 2013 Anton Bikineev
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005 //
0006 #ifndef BOOST_MATH_BESSEL_DERIVATIVES_HPP
0007 #define BOOST_MATH_BESSEL_DERIVATIVES_HPP
0008
0009 #ifdef _MSC_VER
0010 #  pragma once
0011 #endif
0012
0013 #include <boost/math/special_functions/math_fwd.hpp>
0014 #include <boost/math/special_functions/bessel.hpp>
0015 #include <boost/math/special_functions/detail/bessel_jy_derivatives_asym.hpp>
0016 #include <boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp>
0017 #include <boost/math/special_functions/detail/bessel_derivatives_linear.hpp>
0018
0019 namespace boost{ namespace math{
0020
0021 namespace detail{
0022
0023 template <class Tag, class T, class Policy>
0024 inline T cyl_bessel_j_prime_imp(T v, T x, const Policy& pol)
0025 {
0026    static const char* const function = "boost::math::cyl_bessel_j_prime<%1%>(%1%,%1%)";
0027    BOOST_MATH_STD_USING
0028    //
0029    // Prevent complex result:
0030    //
0031    if ((x < 0) && (floor(v) != v))
0032       return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function requires x >= 0", x, pol);
0033    //
0034    // Special cases for x == 0:
0035    //
0036    if (x == 0)
0037    {
0038       if (v == 1)
0039          return static_cast<T>(0.5);
0040       else if (v == -1)
0041          return static_cast<T>(-0.5);
0042       else if (floor(v) == v || v > 1)
0043          return 0;
0044       else
0045          return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function is indeterminate for this order", x, pol);
0046    }
0047    //
0048    // Special case for large x: use asymptotic expansion:
0049    //
0050    if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x))
0051       return boost::math::detail::asymptotic_bessel_j_derivative_large_x_2(v, x, pol);
0052    //
0053    // Special case for small x: use Taylor series:
0054    //
0055    if ((abs(x) < 5) || (abs(v) > x * x / 4))
0056    {
0057       bool inversed = false;
0058       if (floor(v) == v && v < 0)
0059       {
0060          v = -v;
0061          if (itrunc(v, pol) & 1)
0062             inversed = true;
0063       }
0064       T r = boost::math::detail::bessel_j_derivative_small_z_series(v, x, pol);
0065       return inversed ? T(-r) : r;
0066    }
0067    //
0068    // Special case for v == 0:
0069    //
0070    if (v == 0)
0071       return -boost::math::detail::cyl_bessel_j_imp<T>(1, x, Tag(), pol);
0072    //
0073    // Default case:
0074    //
0075    return boost::math::detail::bessel_j_derivative_linear(v, x, Tag(), pol);
0076 }
0077
0078 template <class T, class Policy>
0079 inline T sph_bessel_j_prime_imp(unsigned v, T x, const Policy& pol)
0080 {
0081    static const char* const function = "boost::math::sph_bessel_prime<%1%>(%1%,%1%)";
0082    //
0083    // Prevent complex result:
0084    //
0085    if (x < 0)
0086       return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function requires x >= 0.", x, pol);
0087    //
0088    // Special case for v == 0:
0089    //
0090    if (v == 0)
0091       return (x == 0) ? boost::math::policies::raise_overflow_error<T>(function, nullptr, pol)
0092          : static_cast<T>(-boost::math::detail::sph_bessel_j_imp<T>(1, x, pol));
0093    //
0094    // Special case for x == 0 and v > 0:
0095    //
0096    if (x == 0)
0097       return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function is indeterminate for this order", x, pol);
0098    //
0099    // Default case:
0100    //
0101    return boost::math::detail::sph_bessel_j_derivative_linear(v, x, pol);
0102 }
0103
0104 template <class T, class Policy>
0105 inline T cyl_bessel_i_prime_imp(T v, T x, const Policy& pol)
0106 {
0107    static const char* const function = "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)";
0108    BOOST_MATH_STD_USING
0109    //
0110    // Prevent complex result:
0111    //
0112    if (x < 0 && floor(v) != v)
0113       return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function requires x >= 0", x, pol);
0114    //
0115    // Special cases for x == 0:
0116    //
0117    if (x == 0)
0118    {
0119       if (v == 1 || v == -1)
0120          return static_cast<T>(0.5);
0121       else if (floor(v) == v || v > 1)
0122          return 0;
0123       else
0124          return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function is indeterminate for this order", x, pol);
0125    }
0126    //
0127    // Special case for v == 0:
0128    //
0129    if (v == 0)
0130       return boost::math::detail::cyl_bessel_i_imp<T>(1, x, pol);
0131    //
0132    // Default case:
0133    //
0134    return boost::math::detail::bessel_i_derivative_linear(v, x, pol);
0135 }
0136
0137 template <class Tag, class T, class Policy>
0138 inline T cyl_bessel_k_prime_imp(T v, T x, const Policy& pol)
0139 {
0140    //
0141    // Prevent complex and indeterminate results:
0142    //
0143    if (x <= 0)
0144       return boost::math::policies::raise_domain_error<T>("boost::math::cyl_bessel_k_prime<%1%>(%1%,%1%)", "Got x = %1%, but function requires x > 0", x, pol);
0145    //
0146    // Special case for v == 0:
0147    //
0148    if (v == 0)
0149       return -boost::math::detail::cyl_bessel_k_imp<T>(1, x, Tag(), pol);
0150    //
0151    // Default case:
0152    //
0153    return boost::math::detail::bessel_k_derivative_linear(v, x, Tag(), pol);
0154 }
0155
0156 template <class Tag, class T, class Policy>
0157 inline T cyl_neumann_prime_imp(T v, T x, const Policy& pol)
0158 {
0159    BOOST_MATH_STD_USING
0160    //
0161    // Prevent complex and indeterminate results:
0162    //
0163    if (x <= 0)
0164       return boost::math::policies::raise_domain_error<T>("boost::math::cyl_neumann_prime<%1%>(%1%,%1%)", "Got x = %1%, but function requires x > 0", x, pol);
0165    //
0166    // Special case for large x: use asymptotic expansion:
0167    //
0168    if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x))
0169       return boost::math::detail::asymptotic_bessel_y_derivative_large_x_2(v, x, pol);
0170    //
0171    // Special case for small x: use Taylor series:
0172    //
0173    if (v > 0 && floor(v) != v)
0174    {
0175       const T eps = boost::math::policies::get_epsilon<T, Policy>();
0176       if (log(eps / 2) > v * log((x * x) / (v * 4)))
0177          return boost::math::detail::bessel_y_derivative_small_z_series(v, x, pol);
0178    }
0179    //
0180    // Special case for v == 0:
0181    //
0182    if (v == 0)
0183       return -boost::math::detail::cyl_neumann_imp<T>(1, x, Tag(), pol);
0184    //
0185    // Default case:
0186    //
0187    return boost::math::detail::bessel_y_derivative_linear(v, x, Tag(), pol);
0188 }
0189
0190 template <class T, class Policy>
0191 inline T sph_neumann_prime_imp(unsigned v, T x, const Policy& pol)
0192 {
0193    //
0194    // Prevent complex and indeterminate result:
0195    //
0196    if (x <= 0)
0197       return boost::math::policies::raise_domain_error<T>("boost::math::sph_neumann_prime<%1%>(%1%,%1%)", "Got x = %1%, but function requires x > 0.", x, pol);
0198    //
0199    // Special case for v == 0:
0200    //
0201    if (v == 0)
0202       return -boost::math::detail::sph_neumann_imp<T>(1, x, pol);
0203    //
0204    // Default case:
0205    //
0206    return boost::math::detail::sph_neumann_derivative_linear(v, x, pol);
0207 }
0208
0209 } // namespace detail
0210
0211 template <class T1, class T2, class Policy>
0212 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j_prime(T1 v, T2 x, const Policy& /* pol */)
0213 {
0214    BOOST_FPU_EXCEPTION_GUARD
0215    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0216    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
0217    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0218    typedef typename policies::normalise<
0219       Policy,
0220       policies::promote_float<false>,
0221       policies::promote_double<false>,
0222       policies::discrete_quantile<>,
0223       policies::assert_undefined<> >::type forwarding_policy;
0224    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_j_prime<%1%,%1%>(%1%,%1%)");
0225 }
0226
0227 template <class T1, class T2>
0228 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_j_prime(T1 v, T2 x)
0229 {
0230    return cyl_bessel_j_prime(v, x, policies::policy<>());
0231 }
0232
0233 template <class T, class Policy>
0234 inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel_prime(unsigned v, T x, const Policy& /* pol */)
0235 {
0236    BOOST_FPU_EXCEPTION_GUARD
0237    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
0238    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0239    typedef typename policies::normalise<
0240       Policy,
0241       policies::promote_float<false>,
0242       policies::promote_double<false>,
0243       policies::discrete_quantile<>,
0244       policies::assert_undefined<> >::type forwarding_policy;
0245    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%)");
0246 }
0247
0248 template <class T>
0249 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_bessel_prime(unsigned v, T x)
0250 {
0251    return sph_bessel_prime(v, x, policies::policy<>());
0252 }
0253
0254 template <class T1, class T2, class Policy>
0255 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i_prime(T1 v, T2 x, const Policy& /* pol */)
0256 {
0257    BOOST_FPU_EXCEPTION_GUARD
0258    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0259    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0260    typedef typename policies::normalise<
0261       Policy,
0262       policies::promote_float<false>,
0263       policies::promote_double<false>,
0264       policies::discrete_quantile<>,
0265       policies::assert_undefined<> >::type forwarding_policy;
0266    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_prime_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)");
0267 }
0268
0269 template <class T1, class T2>
0270 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_i_prime(T1 v, T2 x)
0271 {
0272    return cyl_bessel_i_prime(v, x, policies::policy<>());
0273 }
0274
0275 template <class T1, class T2, class Policy>
0276 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k_prime(T1 v, T2 x, const Policy& /* pol */)
0277 {
0278    BOOST_FPU_EXCEPTION_GUARD
0279    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0280    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
0281    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0282    typedef typename policies::normalise<
0283       Policy,
0284       policies::promote_float<false>,
0285       policies::promote_double<false>,
0286       policies::discrete_quantile<>,
0287       policies::assert_undefined<> >::type forwarding_policy;
0288    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_k_prime<%1%,%1%>(%1%,%1%)");
0289 }
0290
0291 template <class T1, class T2>
0292 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_k_prime(T1 v, T2 x)
0293 {
0294    return cyl_bessel_k_prime(v, x, policies::policy<>());
0295 }
0296
0297 template <class T1, class T2, class Policy>
0298 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann_prime(T1 v, T2 x, const Policy& /* pol */)
0299 {
0300    BOOST_FPU_EXCEPTION_GUARD
0301    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0302    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
0303    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0304    typedef typename policies::normalise<
0305       Policy,
0306       policies::promote_float<false>,
0307       policies::promote_double<false>,
0308       policies::discrete_quantile<>,
0309       policies::assert_undefined<> >::type forwarding_policy;
0310    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_neumann_prime<%1%,%1%>(%1%,%1%)");
0311 }
0312
0313 template <class T1, class T2>
0314 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_neumann_prime(T1 v, T2 x)
0315 {
0316    return cyl_neumann_prime(v, x, policies::policy<>());
0317 }
0318
0319 template <class T, class Policy>
0320 inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann_prime(unsigned v, T x, const Policy& /* pol */)
0321 {
0322    BOOST_FPU_EXCEPTION_GUARD
0323    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
0324    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0325    typedef typename policies::normalise<
0326       Policy,
0327       policies::promote_float<false>,
0328       policies::promote_double<false>,
0329       policies::discrete_quantile<>,
0330       policies::assert_undefined<> >::type forwarding_policy;
0331    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%)");
0332 }
0333
0334 template <class T>
0335 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann_prime(unsigned v, T x)
0336 {
0337    return sph_neumann_prime(v, x, policies::policy<>());
0338 }
0339
0340 } // namespace math
0341 } // namespace boost
0342
0343 #endif // BOOST_MATH_BESSEL_DERIVATIVES_HPP