ceph/src/boost/libs/math/include/boost/math/special_functions/bessel_prime.hpp

   1 //  Copyright (c) 2013 Anton Bikineev
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5 //
   6 #ifndef BOOST_MATH_BESSEL_DERIVATIVES_HPP
   7 #define BOOST_MATH_BESSEL_DERIVATIVES_HPP
   8
   9 #ifdef _MSC_VER
  10 #  pragma once
  11 #endif
  12
  13 #include <boost/math/special_functions/math_fwd.hpp>
  14 #include <boost/math/special_functions/bessel.hpp>
  15 #include <boost/math/special_functions/detail/bessel_jy_derivatives_asym.hpp>
  16 #include <boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp>
  17 #include <boost/math/special_functions/detail/bessel_derivatives_linear.hpp>
  18
  19 namespace boost{ namespace math{
  20
  21 namespace detail{
  22
  23 template <class Tag, class T, class Policy>
  24 inline T cyl_bessel_j_prime_imp(T v, T x, const Policy& pol)
  25 {
  26    static const char* const function = "boost::math::cyl_bessel_j_prime<%1%>(%1%,%1%)";
  27    BOOST_MATH_STD_USING
  28    //
  29    // Prevent complex result:
  30    //
  31    if (x < 0 && floor(v) != v)
  32       return boost::math::policies::raise_domain_error<T>(
  33          function,
  34          "Got x = %1%, but function requires x >= 0", x, pol);
  35    //
  36    // Special cases for x == 0:
  37    //
  38    if (x == 0)
  39    {
  40       if (v == 1)
  41          return 0.5;
  42       else if (v == -1)
  43          return -0.5;
  44       else if (floor(v) == v || v > 1)
  45          return 0;
  46       else return boost::math::policies::raise_domain_error<T>(
  47          function,
  48          "Got x = %1%, but function is indeterminate for this order", x, pol);
  49    }
  50    //
  51    // Special case for large x: use asymptotic expansion:
  52    //
  53    if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x))
  54       return boost::math::detail::asymptotic_bessel_j_derivative_large_x_2(v, x);
  55    //
  56    // Special case for small x: use Taylor series:
  57    //
  58    if ((abs(x) < 5) || (abs(v) > x * x / 4))
  59    {
  60       bool inversed = false;
  61       if (floor(v) == v && v < 0)
  62       {
  63          v = -v;
  64          if (itrunc(v, pol) & 1)
  65             inversed = true;
  66       }
  67       T r = boost::math::detail::bessel_j_derivative_small_z_series(v, x, pol);
  68       return inversed ? T(-r) : r;
  69    }
  70    //
  71    // Special case for v == 0:
  72    //
  73    if (v == 0)
  74       return -boost::math::detail::cyl_bessel_j_imp<T>(1, x, Tag(), pol);
  75    //
  76    // Default case:
  77    //
  78    return boost::math::detail::bessel_j_derivative_linear(v, x, Tag(), pol);
  79 }
  80
  81 template <class T, class Policy>
  82 inline T sph_bessel_j_prime_imp(unsigned v, T x, const Policy& pol)
  83 {
  84    static const char* const function = "boost::math::sph_bessel_prime<%1%>(%1%,%1%)";
  85    //
  86    // Prevent complex result:
  87    //
  88    if (x < 0)
  89       return boost::math::policies::raise_domain_error<T>(
  90          function,
  91          "Got x = %1%, but function requires x >= 0.", x, pol);
  92    //
  93    // Special case for v == 0:
  94    //
  95    if (v == 0)
  96       return (x == 0) ? boost::math::policies::raise_overflow_error<T>(function, 0, pol)
  97          : static_cast<T>(-boost::math::detail::sph_bessel_j_imp<T>(1, x, pol));
  98    //
  99    // Special case for x == 0 and v > 0:
 100    //
 101    if (x == 0)
 102       return boost::math::policies::raise_domain_error<T>(
 103          function,
 104          "Got x = %1%, but function is indeterminate for this order", x, pol);
 105    //
 106    // Default case:
 107    //
 108    return boost::math::detail::sph_bessel_j_derivative_linear(v, x, pol);
 109 }
 110
 111 template <class T, class Policy>
 112 inline T cyl_bessel_i_prime_imp(T v, T x, const Policy& pol)
 113 {
 114    static const char* const function = "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)";
 115    BOOST_MATH_STD_USING
 116    //
 117    // Prevent complex result:
 118    //
 119    if (x < 0 && floor(v) != v)
 120       return boost::math::policies::raise_domain_error<T>(
 121          function,
 122          "Got x = %1%, but function requires x >= 0", x, pol);
 123    //
 124    // Special cases for x == 0:
 125    //
 126    if (x == 0)
 127    {
 128       if (v == 1 || v == -1)
 129          return 0.5;
 130       else if (floor(v) == v || v > 1)
 131          return 0;
 132       else return boost::math::policies::raise_domain_error<T>(
 133          function,
 134          "Got x = %1%, but function is indeterminate for this order", x, pol);
 135    }
 136    //
 137    // Special case for v == 0:
 138    //
 139    if (v == 0)
 140       return boost::math::detail::cyl_bessel_i_imp<T>(1, x, pol);
 141    //
 142    // Default case:
 143    //
 144    return boost::math::detail::bessel_i_derivative_linear(v, x, pol);
 145 }
 146
 147 template <class Tag, class T, class Policy>
 148 inline T cyl_bessel_k_prime_imp(T v, T x, const Policy& pol)
 149 {
 150    //
 151    // Prevent complex and indeterminate results:
 152    //
 153    if (x <= 0)
 154       return boost::math::policies::raise_domain_error<T>(
 155          "boost::math::cyl_bessel_k_prime<%1%>(%1%,%1%)",
 156          "Got x = %1%, but function requires x > 0", x, pol);
 157    //
 158    // Special case for v == 0:
 159    //
 160    if (v == 0)
 161       return -boost::math::detail::cyl_bessel_k_imp<T>(1, x, Tag(), pol);
 162    //
 163    // Default case:
 164    //
 165    return boost::math::detail::bessel_k_derivative_linear(v, x, Tag(), pol);
 166 }
 167
 168 template <class Tag, class T, class Policy>
 169 inline T cyl_neumann_prime_imp(T v, T x, const Policy& pol)
 170 {
 171    BOOST_MATH_STD_USING
 172    //
 173    // Prevent complex and indeterminate results:
 174    //
 175    if (x <= 0)
 176       return boost::math::policies::raise_domain_error<T>(
 177          "boost::math::cyl_neumann_prime<%1%>(%1%,%1%)",
 178          "Got x = %1%, but function requires x > 0", x, pol);
 179    //
 180    // Special case for large x: use asymptotic expansion:
 181    //
 182    if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x))
 183       return boost::math::detail::asymptotic_bessel_y_derivative_large_x_2(v, x);
 184    //
 185    // Special case for small x: use Taylor series:
 186    //
 187    if (v > 0 && floor(v) != v)
 188    {
 189       const T eps = boost::math::policies::get_epsilon<T, Policy>();
 190       if (log(eps / 2) > v * log((x * x) / (v * 4)))
 191          return boost::math::detail::bessel_y_derivative_small_z_series(v, x, pol);
 192    }
 193    //
 194    // Special case for v == 0:
 195    //
 196    if (v == 0)
 197       return -boost::math::detail::cyl_neumann_imp<T>(1, x, Tag(), pol);
 198    //
 199    // Default case:
 200    //
 201    return boost::math::detail::bessel_y_derivative_linear(v, x, Tag(), pol);
 202 }
 203
 204 template <class T, class Policy>
 205 inline T sph_neumann_prime_imp(unsigned v, T x, const Policy& pol)
 206 {
 207    //
 208    // Prevent complex and indeterminate result:
 209    //
 210    if (x <= 0)
 211       return boost::math::policies::raise_domain_error<T>(
 212          "boost::math::sph_neumann_prime<%1%>(%1%,%1%)",
 213          "Got x = %1%, but function requires x > 0.", x, pol);
 214    //
 215    // Special case for v == 0:
 216    //
 217    if (v == 0)
 218       return -boost::math::detail::sph_neumann_imp<T>(1, x, pol);
 219    //
 220    // Default case:
 221    //
 222    return boost::math::detail::sph_neumann_derivative_linear(v, x, pol);
 223 }
 224
 225 } // namespace detail
 226
 227 template <class T1, class T2, class Policy>
 228 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j_prime(T1 v, T2 x, const Policy& /* pol */)
 229 {
 230    BOOST_FPU_EXCEPTION_GUARD
 231    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
 232    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
 233    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 234    typedef typename policies::normalise<
 235       Policy,
 236       policies::promote_float<false>,
 237       policies::promote_double<false>,
 238       policies::discrete_quantile<>,
 239       policies::assert_undefined<> >::type forwarding_policy;
 240    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%)");
 241 }
 242
 243 template <class T1, class T2>
 244 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_j_prime(T1 v, T2 x)
 245 {
 246    return cyl_bessel_j_prime(v, x, policies::policy<>());
 247 }
 248
 249 template <class T, class Policy>
 250 inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel_prime(unsigned v, T x, const Policy& /* pol */)
 251 {
 252    BOOST_FPU_EXCEPTION_GUARD
 253    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
 254    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 255    typedef typename policies::normalise<
 256       Policy,
 257       policies::promote_float<false>,
 258       policies::promote_double<false>,
 259       policies::discrete_quantile<>,
 260       policies::assert_undefined<> >::type forwarding_policy;
 261    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%)");
 262 }
 263
 264 template <class T>
 265 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_bessel_prime(unsigned v, T x)
 266 {
 267    return sph_bessel_prime(v, x, policies::policy<>());
 268 }
 269
 270 template <class T1, class T2, class Policy>
 271 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i_prime(T1 v, T2 x, const Policy& /* pol */)
 272 {
 273    BOOST_FPU_EXCEPTION_GUARD
 274    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
 275    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 276    typedef typename policies::normalise<
 277       Policy,
 278       policies::promote_float<false>,
 279       policies::promote_double<false>,
 280       policies::discrete_quantile<>,
 281       policies::assert_undefined<> >::type forwarding_policy;
 282    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%)");
 283 }
 284
 285 template <class T1, class T2>
 286 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_i_prime(T1 v, T2 x)
 287 {
 288    return cyl_bessel_i_prime(v, x, policies::policy<>());
 289 }
 290
 291 template <class T1, class T2, class Policy>
 292 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k_prime(T1 v, T2 x, const Policy& /* pol */)
 293 {
 294    BOOST_FPU_EXCEPTION_GUARD
 295    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
 296    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
 297    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 298    typedef typename policies::normalise<
 299       Policy,
 300       policies::promote_float<false>,
 301       policies::promote_double<false>,
 302       policies::discrete_quantile<>,
 303       policies::assert_undefined<> >::type forwarding_policy;
 304    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%)");
 305 }
 306
 307 template <class T1, class T2>
 308 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_k_prime(T1 v, T2 x)
 309 {
 310    return cyl_bessel_k_prime(v, x, policies::policy<>());
 311 }
 312
 313 template <class T1, class T2, class Policy>
 314 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann_prime(T1 v, T2 x, const Policy& /* pol */)
 315 {
 316    BOOST_FPU_EXCEPTION_GUARD
 317    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
 318    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
 319    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 320    typedef typename policies::normalise<
 321       Policy,
 322       policies::promote_float<false>,
 323       policies::promote_double<false>,
 324       policies::discrete_quantile<>,
 325       policies::assert_undefined<> >::type forwarding_policy;
 326    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%)");
 327 }
 328
 329 template <class T1, class T2>
 330 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_neumann_prime(T1 v, T2 x)
 331 {
 332    return cyl_neumann_prime(v, x, policies::policy<>());
 333 }
 334
 335 template <class T, class Policy>
 336 inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann_prime(unsigned v, T x, const Policy& /* pol */)
 337 {
 338    BOOST_FPU_EXCEPTION_GUARD
 339    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
 340    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 341    typedef typename policies::normalise<
 342       Policy,
 343       policies::promote_float<false>,
 344       policies::promote_double<false>,
 345       policies::discrete_quantile<>,
 346       policies::assert_undefined<> >::type forwarding_policy;
 347    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%)");
 348 }
 349
 350 template <class T>
 351 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann_prime(unsigned v, T x)
 352 {
 353    return sph_neumann_prime(v, x, policies::policy<>());
 354 }
 355
 356 } // namespace math
 357 } // namespace boost
 358
 359 #endif // BOOST_MATH_BESSEL_DERIVATIVES_HPP