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