ceph/src/boost/libs/math/include/boost/math/special_functions/detail/bessel_jy_derivatives_asym.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 //
   7 // This is a partial header, do not include on it's own!!!
   8 //
   9 // Contains asymptotic expansions for derivatives of Bessel J(v,x) and Y(v,x)
  10 // functions, as x -> INF.
  11 #ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP
  12 #define BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP
  13
  14 #ifdef _MSC_VER
  15 #pragma once
  16 #endif
  17
  18 namespace boost{ namespace math{ namespace detail{
  19
  20 template <class T>
  21 inline T asymptotic_bessel_derivative_amplitude(T v, T x)
  22 {
  23    // Calculate the amplitude for J'(v,x) and I'(v,x)
  24    // for large x: see A&S 9.2.30.
  25    BOOST_MATH_STD_USING
  26    T s = 1;
  27    const T mu = 4 * v * v;
  28    T txq = 2 * x;
  29    txq *= txq;
  30
  31    s -= (mu - 3) / (2 * txq);
  32    s -= ((mu - 1) * (mu - 45)) / (txq * txq * 8);
  33
  34    return sqrt(s * 2 / (boost::math::constants::pi<T>() * x));
  35 }
  36
  37 template <class T>
  38 inline T asymptotic_bessel_derivative_phase_mx(T v, T x)
  39 {
  40    // Calculate the phase of J'(v, x) and Y'(v, x) for large x.
  41    // See A&S 9.2.31.
  42    // Note that the result returned is the phase less (x - PI(v/2 - 1/4))
  43    // which we'll factor in later when we calculate the sines/cosines of the result:
  44    const T mu = 4 * v * v;
  45    const T mu2 = mu * mu;
  46    const T mu3 = mu2 * mu;
  47    T denom = 4 * x;
  48    T denom_mult = denom * denom;
  49
  50    T s = 0;
  51    s += (mu + 3) / (2 * denom);
  52    denom *= denom_mult;
  53    s += (mu2 + (46 * mu) - 63) / (6 * denom);
  54    denom *= denom_mult;
  55    s += (mu3 + (185 * mu2) - (2053 * mu) + 1899) / (5 * denom);
  56    return s;
  57 }
  58
  59 template <class T>
  60 inline T asymptotic_bessel_y_derivative_large_x_2(T v, T x)
  61 {
  62    // See A&S 9.2.20.
  63    BOOST_MATH_STD_USING
  64    // Get the phase and amplitude:
  65    const T ampl = asymptotic_bessel_derivative_amplitude(v, x);
  66    const T phase = asymptotic_bessel_derivative_phase_mx(v, x);
  67    BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
  68    BOOST_MATH_INSTRUMENT_VARIABLE(phase);
  69    //
  70    // Calculate the sine of the phase, using
  71    // sine/cosine addition rules to factor in
  72    // the x - PI(v/2 - 1/4) term not added to the
  73    // phase when we calculated it.
  74    //
  75    const T cx = cos(x);
  76    const T sx = sin(x);
  77    const T vd2shifted = (v / 2) - 0.25f;
  78    const T ci = cos_pi(vd2shifted);
  79    const T si = sin_pi(vd2shifted);
  80    const T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si);
  81    BOOST_MATH_INSTRUMENT_CODE(sin(phase));
  82    BOOST_MATH_INSTRUMENT_CODE(cos(x));
  83    BOOST_MATH_INSTRUMENT_CODE(cos(phase));
  84    BOOST_MATH_INSTRUMENT_CODE(sin(x));
  85    return sin_phase * ampl;
  86 }
  87
  88 template <class T>
  89 inline T asymptotic_bessel_j_derivative_large_x_2(T v, T x)
  90 {
  91    // See A&S 9.2.20.
  92    BOOST_MATH_STD_USING
  93    // Get the phase and amplitude:
  94    const T ampl = asymptotic_bessel_derivative_amplitude(v, x);
  95    const T phase = asymptotic_bessel_derivative_phase_mx(v, x);
  96    BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
  97    BOOST_MATH_INSTRUMENT_VARIABLE(phase);
  98    //
  99    // Calculate the sine of the phase, using
 100    // sine/cosine addition rules to factor in
 101    // the x - PI(v/2 - 1/4) term not added to the
 102    // phase when we calculated it.
 103    //
 104    BOOST_MATH_INSTRUMENT_CODE(cos(phase));
 105    BOOST_MATH_INSTRUMENT_CODE(cos(x));
 106    BOOST_MATH_INSTRUMENT_CODE(sin(phase));
 107    BOOST_MATH_INSTRUMENT_CODE(sin(x));
 108    const T cx = cos(x);
 109    const T sx = sin(x);
 110    const T vd2shifted = (v / 2) - 0.25f;
 111    const T ci = cos_pi(vd2shifted);
 112    const T si = sin_pi(vd2shifted);
 113    const T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si);
 114    BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase);
 115    return sin_phase * ampl;
 116 }
 117
 118 template <class T>
 119 inline bool asymptotic_bessel_derivative_large_x_limit(const T& v, const T& x)
 120 {
 121    BOOST_MATH_STD_USING
 122    //
 123    // This function is the copy of math::asymptotic_bessel_large_x_limit
 124    // It means that we use the same rules for determining how x is large
 125    // compared to v.
 126    //
 127    // Determines if x is large enough compared to v to take the asymptotic
 128    // forms above.  From A&S 9.2.28 we require:
 129    //    v < x * eps^1/8
 130    // and from A&S 9.2.29 we require:
 131    //    v^12/10 < 1.5 * x * eps^1/10
 132    // using the former seems to work OK in practice with broadly similar
 133    // error rates either side of the divide for v < 10000.
 134    // At double precision eps^1/8 ~= 0.01.
 135    //
 136    return (std::max)(T(fabs(v)), T(1)) < x * sqrt(boost::math::tools::forth_root_epsilon<T>());
 137 }
 138
 139 }}} // namespaces
 140
 141 #endif // BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP