/usr/include/boost/math/special_functions/detail
Edit: /usr/include/boost/math/special_functions/detail/bessel_jy_series.hpp (7241B)
// Copyright (c) 2011 John Maddock
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_BESSEL_JN_SERIES_HPP
#define BOOST_MATH_BESSEL_JN_SERIES_HPP
#ifdef _MSC_VER
#pragma once
#endif
namespace boost { namespace math { namespace detail{
template
struct bessel_j_small_z_series_term
{
typedef T result_type;
bessel_j_small_z_series_term(T v_, T x)
: N(0), v(v_)
{
BOOST_MATH_STD_USING
mult = x / 2;
mult *= -mult;
term = 1;
}
T operator()()
{
T r = term;
++N;
term *= mult / (N * (N + v));
return r;
}
private:
unsigned N;
T v;
T mult;
T term;
};
//
// Series evaluation for BesselJ(v, z) as z -> 0.
// See http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/
// Converges rapidly for all z << v.
//
template
inline T bessel_j_small_z_series(T v, T x, const Policy& pol)
{
BOOST_MATH_STD_USING
T prefix;
if(v < max_factorial::value)
{
prefix = pow(x / 2, v) / boost::math::tgamma(v+1, pol);
}
else
{
prefix = v * log(x / 2) - boost::math::lgamma(v+1, pol);
prefix = exp(prefix);
}
if(0 == prefix)
return prefix;
bessel_j_small_z_series_term s(v, x);
boost::uintmax_t max_iter = policies::get_max_series_iterations();
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T zero = 0;
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter, zero);
#else
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter);
#endif
policies::check_series_iterations("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
return prefix * result;
}
template
struct bessel_y_small_z_series_term_a
{
typedef T result_type;
bessel_y_small_z_series_term_a(T v_, T x)
: N(0), v(v_)
{
BOOST_MATH_STD_USING
mult = x / 2;
mult *= -mult;
term = 1;
}
T operator()()
{
BOOST_MATH_STD_USING
T r = term;
++N;
term *= mult / (N * (N - v));
return r;
}
private:
unsigned N;
T v;
T mult;
T term;
};
template
struct bessel_y_small_z_series_term_b
{
typedef T result_type;
bessel_y_small_z_series_term_b(T v_, T x)
: N(0), v(v_)
{
BOOST_MATH_STD_USING
mult = x / 2;
mult *= -mult;
term = 1;
}
T operator()()
{
T r = term;
++N;
term *= mult / (N * (N + v));
return r;
}
private:
unsigned N;
T v;
T mult;
T term;
};
//
// Series form for BesselY as z -> 0,
// see: http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/01/0003/
// This series is only useful when the second term is small compared to the first
// otherwise we get catastrophic cancellation errors.
//
// Approximating tgamma(v) by v^v, and assuming |tgamma(-z)| < eps we end up requiring:
// eps/2 * v^v(x/2)^-v > (x/2)^v or log(eps/2) > v log((x/2)^2/v)
//
template
inline T bessel_y_small_z_series(T v, T x, T* pscale, const Policy& pol)
{
BOOST_MATH_STD_USING
static const char* function = "bessel_y_small_z_series<%1%>(%1%,%1%)";
T prefix;
T gam;
T p = log(x / 2);
T scale = 1;
bool need_logs = (v >= max_factorial::value) || (tools::log_max_value() / v < fabs(p));
if(!need_logs)
{
gam = boost::math::tgamma(v, pol);
p = pow(x / 2, v);
if(tools::max_value() * p < gam)
{
scale /= gam;
gam = 1;
if(tools::max_value() * p < gam)
{
return -policies::raise_overflow_error(function, 0, pol);
}
}
prefix = -gam / (constants::pi() * p);
}
else
{
gam = boost::math::lgamma(v, pol);
p = v * p;
prefix = gam - log(constants::pi()) - p;
if(tools::log_max_value() < prefix)
{
prefix -= log(tools::max_value() / 4);
scale /= (tools::max_value() / 4);
if(tools::log_max_value() < prefix)
{
return -policies::raise_overflow_error(function, 0, pol);
}
}
prefix = -exp(prefix);
}
bessel_y_small_z_series_term_a s(v, x);
boost::uintmax_t max_iter = policies::get_max_series_iterations();
*pscale = scale;
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T zero = 0;
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter, zero);
#else
T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter);
#endif
policies::check_series_iterations("boost::math::bessel_y_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
result *= prefix;
if(!need_logs)
{
prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v) * p / constants::pi();
}
else
{
int sgn;
prefix = boost::math::lgamma(-v, &sgn, pol) + p;
prefix = exp(prefix) * sgn / constants::pi();
}
bessel_y_small_z_series_term_b s2(v, x);
max_iter = policies::get_max_series_iterations();
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon(), max_iter, zero);
#else
T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon(), max_iter);
#endif
result -= scale * prefix * b;
return result;
}
template
T bessel_yn_small_z(int n, T z, T* scale, const Policy& pol)
{
//
// See http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
//
// Note that when called we assume that x < epsilon and n is a positive integer.
//
BOOST_MATH_STD_USING
BOOST_ASSERT(n >= 0);
BOOST_ASSERT((z < policies::get_epsilon()));
if(n == 0)
{
return (2 / constants::pi()) * (log(z / 2) + constants::euler());
}
else if(n == 1)
{
return (z / constants::pi()) * log(z / 2)
- 2 / (constants::pi() * z)
- (z / (2 * constants::pi())) * (1 - 2 * constants::euler());
}
else if(n == 2)
{
return (z * z) / (4 * constants::pi()) * log(z / 2)
- (4 / (constants::pi() * z * z))
- ((z * z) / (8 * constants::pi())) * (T(3)/2 - 2 * constants::euler());
}
else
{
T p = pow(z / 2, n);
T result = -((boost::math::factorial(n - 1) / constants::pi()));
if(p * tools::max_value() < result)
{
T div = tools::max_value() / 8;
result /= div;
*scale /= div;
if(p * tools::max_value() < result)
{
return -policies::raise_overflow_error("bessel_yn_small_z<%1%>(%1%,%1%)", 0, pol);
}
}
return result / p;
}
}
}}} // namespaces
#endif // BOOST_MATH_BESSEL_JN_SERIES_HPP