boost/math/special_functions/detail/bessel_jy_series.hpp

   1 //  Copyright (c) 2011 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 #ifndef BOOST_MATH_BESSEL_JN_SERIES_HPP
   7 #define BOOST_MATH_BESSEL_JN_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_small_z_series_term
  17 {
  18    typedef T result_type;
  19
  20    bessel_j_small_z_series_term(T v_, T x)
  21       : N(0), v(v_)
  22    {
  23       BOOST_MATH_STD_USING
  24       mult = x / 2;
  25       mult *= -mult;
  26       term = 1;
  27    }
  28    T operator()()
  29    {
  30       T r = term;
  31       ++N;
  32       term *= mult / (N * (N + v));
  33       return r;
  34    }
  35 private:
  36    unsigned N;
  37    T v;
  38    T mult;
  39    T term;
  40 };
  41 //
  42 // Series evaluation for BesselJ(v, z) as z -> 0.
  43 // See http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/
  44 // Converges rapidly for all z << v.
  45 //
  46 template <class T, class Policy>
  47 inline T bessel_j_small_z_series(T v, T x, const Policy& pol)
  48 {
  49    BOOST_MATH_STD_USING
  50    T prefix;
  51    if(v < max_factorial<T>::value)
  52    {
  53       prefix = pow(x / 2, v) / boost::math::tgamma(v+1, pol);
  54    }
  55    else
  56    {
  57       prefix = v * log(x / 2) - boost::math::lgamma(v+1, pol);
  58       prefix = exp(prefix);
  59    }
  60    if(0 == prefix)
  61       return prefix;
  62
  63    bessel_j_small_z_series_term<T, Policy> s(v, x);
  64    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
  65 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
  66    T zero = 0;
  67    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
  68 #else
  69    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
  70 #endif
  71    policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
  72    return prefix * result;
  73 }
  74
  75 template <class T, class Policy>
  76 struct bessel_y_small_z_series_term_a
  77 {
  78    typedef T result_type;
  79
  80    bessel_y_small_z_series_term_a(T v_, T x)
  81       : N(0), v(v_)
  82    {
  83       BOOST_MATH_STD_USING
  84       mult = x / 2;
  85       mult *= -mult;
  86       term = 1;
  87    }
  88    T operator()()
  89    {
  90       BOOST_MATH_STD_USING
  91       T r = term;
  92       ++N;
  93       term *= mult / (N * (N - v));
  94       return r;
  95    }
  96 private:
  97    unsigned N;
  98    T v;
  99    T mult;
 100    T term;
 101 };
 102
 103 template <class T, class Policy>
 104 struct bessel_y_small_z_series_term_b
 105 {
 106    typedef T result_type;
 107
 108    bessel_y_small_z_series_term_b(T v_, T x)
 109       : N(0), v(v_)
 110    {
 111       BOOST_MATH_STD_USING
 112       mult = x / 2;
 113       mult *= -mult;
 114       term = 1;
 115    }
 116    T operator()()
 117    {
 118       T r = term;
 119       ++N;
 120       term *= mult / (N * (N + v));
 121       return r;
 122    }
 123 private:
 124    unsigned N;
 125    T v;
 126    T mult;
 127    T term;
 128 };
 129 //
 130 // Series form for BesselY as z -> 0,
 131 // see: http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/01/0003/
 132 // This series is only useful when the second term is small compared to the first
 133 // otherwise we get catestrophic cancellation errors.
 134 //
 135 // Approximating tgamma(v) by v^v, and assuming |tgamma(-z)| < eps we end up requiring:
 136 // eps/2 * v^v(x/2)^-v > (x/2)^v or log(eps/2) > v log((x/2)^2/v)
 137 //
 138 template <class T, class Policy>
 139 inline T bessel_y_small_z_series(T v, T x, T* pscale, const Policy& pol)
 140 {
 141    BOOST_MATH_STD_USING
 142    static const char* function = "bessel_y_small_z_series<%1%>(%1%,%1%)";
 143    T prefix;
 144    T gam;
 145    T p = log(x / 2);
 146    T scale = 1;
 147    bool need_logs = (v >= max_factorial<T>::value) || (tools::log_max_value<T>() / v < fabs(p));
 148    if(!need_logs)
 149    {
 150       gam = boost::math::tgamma(v, pol);
 151       p = pow(x / 2, v);
 152       if(tools::max_value<T>() * p < gam)
 153       {
 154          scale /= gam;
 155          gam = 1;
 156          if(tools::max_value<T>() * p < gam)
 157          {
 158             return -policies::raise_overflow_error<T>(function, 0, pol);
 159          }
 160       }
 161       prefix = -gam / (constants::pi<T>() * p);
 162    }
 163    else
 164    {
 165       gam = boost::math::lgamma(v, pol);
 166       p = v * p;
 167       prefix = gam - log(constants::pi<T>()) - p;
 168       if(tools::log_max_value<T>() < prefix)
 169       {
 170          prefix -= log(tools::max_value<T>() / 4);
 171          scale /= (tools::max_value<T>() / 4);
 172          if(tools::log_max_value<T>() < prefix)
 173          {
 174             return -policies::raise_overflow_error<T>(function, 0, pol);
 175          }
 176       }
 177       prefix = -exp(prefix);
 178    }
 179    bessel_y_small_z_series_term_a<T, Policy> s(v, x);
 180    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 181    *pscale = scale;
 182 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 183    T zero = 0;
 184    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
 185 #else
 186    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 187 #endif
 188    policies::check_series_iterations<T>("boost::math::bessel_y_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
 189    result *= prefix;
 190
 191    if(!need_logs)
 192    {
 193       prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v) * p / constants::pi<T>();
 194    }
 195    else
 196    {
 197       int sgn;
 198       prefix = boost::math::lgamma(-v, &sgn, pol) + p;
 199       prefix = exp(prefix) * sgn / constants::pi<T>();
 200    }
 201    bessel_y_small_z_series_term_b<T, Policy> s2(v, x);
 202    max_iter = policies::get_max_series_iterations<Policy>();
 203 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 204    T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
 205 #else
 206    T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 207 #endif
 208    result -= scale * prefix * b;
 209    return result;
 210 }
 211
 212 template <class T, class Policy>
 213 T bessel_yn_small_z(int n, T z, T* scale, const Policy& pol)
 214 {
 215    //
 216    // See http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
 217    //
 218    // Note that when called we assume that x < epsilon and n is a positive integer.
 219    //
 220    BOOST_MATH_STD_USING
 221    BOOST_ASSERT(n >= 0);
 222    BOOST_ASSERT((z < policies::get_epsilon<T, Policy>()));
 223
 224    if(n == 0)
 225    {
 226       return (2 / constants::pi<T>()) * (log(z / 2) +  constants::euler<T>());
 227    }
 228    else if(n == 1)
 229    {
 230       return (z / constants::pi<T>()) * log(z / 2)
 231          - 2 / (constants::pi<T>() * z)
 232          - (z / (2 * constants::pi<T>())) * (1 - 2 * constants::euler<T>());
 233    }
 234    else if(n == 2)
 235    {
 236       return (z * z) / (4 * constants::pi<T>()) * log(z / 2)
 237          - (4 / (constants::pi<T>() * z * z))
 238          - ((z * z) / (8 * constants::pi<T>())) * (T(3)/2 - 2 * constants::euler<T>());
 239    }
 240    else
 241    {
 242       T p = pow(z / 2, n);
 243       T result = -((boost::math::factorial<T>(n - 1) / constants::pi<T>()));
 244       if(p * tools::max_value<T>() < result)
 245       {
 246          T div = tools::max_value<T>() / 8;
 247          result /= div;
 248          *scale /= div;
 249          if(p * tools::max_value<T>() < result)
 250          {
 251             return -policies::raise_overflow_error<T>("bessel_yn_small_z<%1%>(%1%,%1%)", 0, pol);
 252          }
 253       }
 254       return result / p;
 255    }
 256 }
 257
 258 }}} // namespaces
 259
 260 #endif // BOOST_MATH_BESSEL_JN_SERIES_HPP
 261