ceph/src/boost/boost/math/special_functions/detail/bessel_jy_asym.hpp

   1 //  Copyright (c) 2007 John Maddock
   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 Bessel J(v,x) and Y(v,x)
  10 // functions, as x -> INF.
  11 //
  12 #ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP
  13 #define BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP
  14
  15 #ifdef _MSC_VER
  16 #pragma once
  17 #endif
  18
  19 #include <boost/math/special_functions/factorials.hpp>
  20
  21 namespace boost{ namespace math{ namespace detail{
  22
  23 template <class T>
  24 inline T asymptotic_bessel_amplitude(T v, T x)
  25 {
  26    // Calculate the amplitude of J(v, x) and Y(v, x) for large
  27    // x: see A&S 9.2.28.
  28    BOOST_MATH_STD_USING
  29    T s = 1;
  30    T mu = 4 * v * v;
  31    T txq = 2 * x;
  32    txq *= txq;
  33
  34    s += (mu - 1) / (2 * txq);
  35    s += 3 * (mu - 1) * (mu - 9) / (txq * txq * 8);
  36    s += 15 * (mu - 1) * (mu - 9) * (mu - 25) / (txq * txq * txq * 8 * 6);
  37
  38    return sqrt(s * 2 / (constants::pi<T>() * x));
  39 }
  40
  41 template <class T>
  42 T asymptotic_bessel_phase_mx(T v, T x)
  43 {
  44    //
  45    // Calculate the phase of J(v, x) and Y(v, x) for large x.
  46    // See A&S 9.2.29.
  47    // Note that the result returned is the phase less (x - PI(v/2 + 1/4))
  48    // which we'll factor in later when we calculate the sines/cosines of the result:
  49    //
  50    T mu = 4 * v * v;
  51    T denom = 4 * x;
  52    T denom_mult = denom * denom;
  53
  54    T s = 0;
  55    s += (mu - 1) / (2 * denom);
  56    denom *= denom_mult;
  57    s += (mu - 1) * (mu - 25) / (6 * denom);
  58    denom *= denom_mult;
  59    s += (mu - 1) * (mu * mu - 114 * mu + 1073) / (5 * denom);
  60    denom *= denom_mult;
  61    s += (mu - 1) * (5 * mu * mu * mu - 1535 * mu * mu + 54703 * mu - 375733) / (14 * denom);
  62    return s;
  63 }
  64
  65 template <class T>
  66 inline T asymptotic_bessel_y_large_x_2(T v, T x)
  67 {
  68    // See A&S 9.2.19.
  69    BOOST_MATH_STD_USING
  70    // Get the phase and amplitude:
  71    T ampl = asymptotic_bessel_amplitude(v, x);
  72    T phase = asymptotic_bessel_phase_mx(v, x);
  73    BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
  74    BOOST_MATH_INSTRUMENT_VARIABLE(phase);
  75    //
  76    // Calculate the sine of the phase, using
  77    // sine/cosine addition rules to factor in
  78    // the x - PI(v/2 + 1/4) term not added to the
  79    // phase when we calculated it.
  80    //
  81    T cx = cos(x);
  82    T sx = sin(x);
  83    T ci = cos_pi(v / 2 + 0.25f);
  84    T si = sin_pi(v / 2 + 0.25f);
  85    T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si);
  86    BOOST_MATH_INSTRUMENT_CODE(sin(phase));
  87    BOOST_MATH_INSTRUMENT_CODE(cos(x));
  88    BOOST_MATH_INSTRUMENT_CODE(cos(phase));
  89    BOOST_MATH_INSTRUMENT_CODE(sin(x));
  90    return sin_phase * ampl;
  91 }
  92
  93 template <class T>
  94 inline T asymptotic_bessel_j_large_x_2(T v, T x)
  95 {
  96    // See A&S 9.2.19.
  97    BOOST_MATH_STD_USING
  98    // Get the phase and amplitude:
  99    T ampl = asymptotic_bessel_amplitude(v, x);
 100    T phase = asymptotic_bessel_phase_mx(v, x);
 101    BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
 102    BOOST_MATH_INSTRUMENT_VARIABLE(phase);
 103    //
 104    // Calculate the sine of the phase, using
 105    // sine/cosine addition rules to factor in
 106    // the x - PI(v/2 + 1/4) term not added to the
 107    // phase when we calculated it.
 108    //
 109    BOOST_MATH_INSTRUMENT_CODE(cos(phase));
 110    BOOST_MATH_INSTRUMENT_CODE(cos(x));
 111    BOOST_MATH_INSTRUMENT_CODE(sin(phase));
 112    BOOST_MATH_INSTRUMENT_CODE(sin(x));
 113    T cx = cos(x);
 114    T sx = sin(x);
 115    T ci = cos_pi(v / 2 + 0.25f);
 116    T si = sin_pi(v / 2 + 0.25f);
 117    T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si);
 118    BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase);
 119    return sin_phase * ampl;
 120 }
 121
 122 template <class T>
 123 inline bool asymptotic_bessel_large_x_limit(int v, const T& x)
 124 {
 125    BOOST_MATH_STD_USING
 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       BOOST_ASSERT(v >= 0);
 137       return (v ? v : 1) < x * 0.004f;
 138 }
 139
 140 template <class T>
 141 inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
 142 {
 143    BOOST_MATH_STD_USING
 144    //
 145    // Determines if x is large enough compared to v to take the asymptotic
 146    // forms above.  From A&S 9.2.28 we require:
 147    //    v < x * eps^1/8
 148    // and from A&S 9.2.29 we require:
 149    //    v^12/10 < 1.5 * x * eps^1/10
 150    // using the former seems to work OK in practice with broadly similar
 151    // error rates either side of the divide for v < 10000.
 152    // At double precision eps^1/8 ~= 0.01.
 153    //
 154    return (std::max)(T(fabs(v)), T(1)) < x * sqrt(tools::forth_root_epsilon<T>());
 155 }
 156
 157 template <class T, class Policy>
 158 void temme_asyptotic_y_small_x(T v, T x, T* Y, T* Y1, const Policy& pol)
 159 {
 160    T c = 1;
 161    T p = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, -v) / boost::math::tgamma(1 - v, pol);
 162    T q = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, v) / boost::math::tgamma(1 + v, pol);
 163    T f = (p - q) / v;
 164    T g_prefix = boost::math::sin_pi(v / 2, pol);
 165    g_prefix *= g_prefix * 2 / v;
 166    T g = f + g_prefix * q;
 167    T h = p;
 168    T c_mult = -x * x / 4;
 169
 170    T y(c * g), y1(c * h);
 171
 172    for(int k = 1; k < policies::get_max_series_iterations<Policy>(); ++k)
 173    {
 174       f = (k * f + p + q) / (k*k - v*v);
 175       p /= k - v;
 176       q /= k + v;
 177       c *= c_mult / k;
 178       T c1 = pow(-x * x / 4, k) / factorial<T>(k, pol);
 179       g = f + g_prefix * q;
 180       h = -k * g + p;
 181       y += c * g;
 182       y1 += c * h;
 183       if(c * g / tools::epsilon<T>() < y)
 184          break;
 185    }
 186
 187    *Y = -y;
 188    *Y1 = (-2 / x) * y1;
 189 }
 190
 191 template <class T, class Policy>
 192 T asymptotic_bessel_i_large_x(T v, T x, const Policy& pol)
 193 {
 194    BOOST_MATH_STD_USING  // ADL of std names
 195    T s = 1;
 196    T mu = 4 * v * v;
 197    T ex = 8 * x;
 198    T num = mu - 1;
 199    T denom = ex;
 200
 201    s -= num / denom;
 202
 203    num *= mu - 9;
 204    denom *= ex * 2;
 205    s += num / denom;
 206
 207    num *= mu - 25;
 208    denom *= ex * 3;
 209    s -= num / denom;
 210
 211    // Try and avoid overflow to the last minute:
 212    T e = exp(x/2);
 213
 214    s = e * (e * s / sqrt(2 * x * constants::pi<T>()));
 215
 216    return (boost::math::isfinite)(s) ?
 217       s : policies::raise_overflow_error<T>("boost::math::asymptotic_bessel_i_large_x<%1%>(%1%,%1%)", 0, pol);
 218 }
 219
 220 }}} // namespaces
 221
 222 #endif
 223