ceph/src/boost/libs/math/include/boost/math/special_functions/detail/bessel_jy_derivatives_series.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_JY_DERIVATIVES_SERIES_HPP
   7 #define BOOST_MATH_BESSEL_JY_DERIVATIVES_SERIES_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 namespace boost{ namespace math{ namespace detail{
  14
  15 template <class T, class Policy>
  16 struct bessel_j_derivative_small_z_series_term
  17 {
  18    typedef T result_type;
  19
  20    bessel_j_derivative_small_z_series_term(T v_, T x)
  21       : N(0), v(v_), term(1), mult(x / 2)
  22    {
  23       mult *= -mult;
  24       // iterate if v == 0; otherwise result of
  25       // first term is 0 and tools::sum_series stops
  26       if (v == 0)
  27          iterate();
  28    }
  29    T operator()()
  30    {
  31       T r = term * (v + 2 * N);
  32       iterate();
  33       return r;
  34    }
  35 private:
  36    void iterate()
  37    {
  38       ++N;
  39       term *= mult / (N * (N + v));
  40    }
  41    unsigned N;
  42    T v;
  43    T term;
  44    T mult;
  45 };
  46 //
  47 // Series evaluation for BesselJ'(v, z) as z -> 0.
  48 // It's derivative of http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/
  49 // Converges rapidly for all z << v.
  50 //
  51 template <class T, class Policy>
  52 inline T bessel_j_derivative_small_z_series(T v, T x, const Policy& pol)
  53 {
  54    BOOST_MATH_STD_USING
  55    T prefix;
  56    if (v < boost::math::max_factorial<T>::value)
  57    {
  58       prefix = pow(x / 2, v - 1) / 2 / boost::math::tgamma(v + 1, pol);
  59    }
  60    else
  61    {
  62       prefix = (v - 1) * log(x / 2) - constants::ln_two<T>() - boost::math::lgamma(v + 1, pol);
  63       prefix = exp(prefix);
  64    }
  65    if (0 == prefix)
  66       return prefix;
  67
  68    bessel_j_derivative_small_z_series_term<T, Policy> s(v, x);
  69    boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
  70 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
  71    T zero = 0;
  72    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
  73 #else
  74    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
  75 #endif
  76    boost::math::policies::check_series_iterations<T>("boost::math::bessel_j_derivative_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
  77    return prefix * result;
  78 }
  79
  80 template <class T, class Policy>
  81 struct bessel_y_derivative_small_z_series_term_a
  82 {
  83    typedef T result_type;
  84
  85    bessel_y_derivative_small_z_series_term_a(T v_, T x)
  86       : N(0), v(v_)
  87    {
  88       mult = x / 2;
  89       mult *= -mult;
  90       term = 1;
  91    }
  92    T operator()()
  93    {
  94       T r = term * (-v + 2 * N);
  95       ++N;
  96       term *= mult / (N * (N - v));
  97       return r;
  98    }
  99 private:
 100    unsigned N;
 101    T v;
 102    T mult;
 103    T term;
 104 };
 105
 106 template <class T, class Policy>
 107 struct bessel_y_derivative_small_z_series_term_b
 108 {
 109    typedef T result_type;
 110
 111    bessel_y_derivative_small_z_series_term_b(T v_, T x)
 112       : N(0), v(v_)
 113    {
 114       mult = x / 2;
 115       mult *= -mult;
 116       term = 1;
 117    }
 118    T operator()()
 119    {
 120       T r = term * (v + 2 * N);
 121       ++N;
 122       term *= mult / (N * (N + v));
 123       return r;
 124    }
 125 private:
 126    unsigned N;
 127    T v;
 128    T mult;
 129    T term;
 130 };
 131 //
 132 // Series form for BesselY' as z -> 0,
 133 // It's derivative of http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/01/0003/
 134 // This series is only useful when the second term is small compared to the first
 135 // otherwise we get catestrophic cancellation errors.
 136 //
 137 // Approximating tgamma(v) by v^v, and assuming |tgamma(-z)| < eps we end up requiring:
 138 // eps/2 * v^v(x/2)^-v > (x/2)^v or log(eps/2) > v log((x/2)^2/v)
 139 //
 140 template <class T, class Policy>
 141 inline T bessel_y_derivative_small_z_series(T v, T x, const Policy& pol)
 142 {
 143    BOOST_MATH_STD_USING
 144    static const char* function = "bessel_y_derivative_small_z_series<%1%>(%1%,%1%)";
 145    T prefix;
 146    T gam;
 147    T p = log(x / 2);
 148    T scale = 1;
 149    bool need_logs = (v >= boost::math::max_factorial<T>::value) || (boost::math::tools::log_max_value<T>() / v < fabs(p));
 150    if (!need_logs)
 151    {
 152       gam = boost::math::tgamma(v, pol);
 153       p = pow(x / 2, v + 1) * 2;
 154       if (boost::math::tools::max_value<T>() * p < gam)
 155       {
 156          scale /= gam;
 157          gam = 1;
 158          if (boost::math::tools::max_value<T>() * p < gam)
 159          {
 160             // This term will overflow to -INF, when combined with the series below it becomes +INF:
 161             return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
 162          }
 163       }
 164       prefix = -gam / (boost::math::constants::pi<T>() * p);
 165    }
 166    else
 167    {
 168       gam = boost::math::lgamma(v, pol);
 169       p = (v + 1) * p + constants::ln_two<T>();
 170       prefix = gam - log(boost::math::constants::pi<T>()) - p;
 171       if (boost::math::tools::log_max_value<T>() < prefix)
 172       {
 173          prefix -= log(boost::math::tools::max_value<T>() / 4);
 174          scale /= (boost::math::tools::max_value<T>() / 4);
 175          if (boost::math::tools::log_max_value<T>() < prefix)
 176          {
 177             return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
 178          }
 179       }
 180       prefix = -exp(prefix);
 181    }
 182    bessel_y_derivative_small_z_series_term_a<T, Policy> s(v, x);
 183    boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 184 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 185    T zero = 0;
 186    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
 187 #else
 188    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 189 #endif
 190    boost::math::policies::check_series_iterations<T>("boost::math::bessel_y_derivative_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
 191    result *= prefix;
 192
 193    p = pow(x / 2, v - 1) / 2;
 194    if (!need_logs)
 195    {
 196       prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v) * p / boost::math::constants::pi<T>();
 197    }
 198    else
 199    {
 200       int sgn;
 201       prefix = boost::math::lgamma(-v, &sgn, pol) + (v - 1) * log(x / 2) - constants::ln_two<T>();
 202       prefix = exp(prefix) * sgn / boost::math::constants::pi<T>();
 203    }
 204    bessel_y_derivative_small_z_series_term_b<T, Policy> s2(v, x);
 205    max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 206 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 207    T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
 208 #else
 209    T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 210 #endif
 211    result += scale * prefix * b;
 212    return result;
 213 }
 214
 215 // Calculating of BesselY'(v,x) with small x (x < epsilon) and integer x using derivatives
 216 // of formulas in http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
 217 // seems to lose precision. Instead using linear combination of regular Bessel is preferred.
 218
 219 }}} // namespaces
 220
 221 #endif // BOOST_MATH_BESSEL_JY_DERIVATVIES_SERIES_HPP