--- /dev/null
+// 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 <class T, class Policy>
+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 <class T, class Policy>
+inline T bessel_j_small_z_series(T v, T x, const Policy& pol)
+{
+ BOOST_MATH_STD_USING
+ T prefix;
+ if(v < max_factorial<T>::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<T, Policy> s(v, x);
+ boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+ T zero = 0;
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
+#else
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
+#endif
+ policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
+ return prefix * result;
+}
+
+template <class T, class Policy>
+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 <class T, class Policy>
+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 catestrophic 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 <class T, class Policy>
+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<T>::value) || (tools::log_max_value<T>() / v < fabs(p));
+ if(!need_logs)
+ {
+ gam = boost::math::tgamma(v, pol);
+ p = pow(x / 2, v);
+ if(tools::max_value<T>() * p < gam)
+ {
+ scale /= gam;
+ gam = 1;
+ if(tools::max_value<T>() * p < gam)
+ {
+ return -policies::raise_overflow_error<T>(function, 0, pol);
+ }
+ }
+ prefix = -gam / (constants::pi<T>() * p);
+ }
+ else
+ {
+ gam = boost::math::lgamma(v, pol);
+ p = v * p;
+ prefix = gam - log(constants::pi<T>()) - p;
+ if(tools::log_max_value<T>() < prefix)
+ {
+ prefix -= log(tools::max_value<T>() / 4);
+ scale /= (tools::max_value<T>() / 4);
+ if(tools::log_max_value<T>() < prefix)
+ {
+ return -policies::raise_overflow_error<T>(function, 0, pol);
+ }
+ }
+ prefix = -exp(prefix);
+ }
+ bessel_y_small_z_series_term_a<T, Policy> s(v, x);
+ boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+ *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<T, Policy>(), max_iter, zero);
+#else
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
+#endif
+ policies::check_series_iterations<T>("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<T>();
+ }
+ else
+ {
+ int sgn;
+ prefix = boost::math::lgamma(-v, &sgn, pol) + p;
+ prefix = exp(prefix) * sgn / constants::pi<T>();
+ }
+ bessel_y_small_z_series_term_b<T, Policy> s2(v, x);
+ max_iter = policies::get_max_series_iterations<Policy>();
+#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+ T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
+#else
+ T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
+#endif
+ result -= scale * prefix * b;
+ return result;
+}
+
+template <class T, class Policy>
+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<T, Policy>()));
+
+ if(n == 0)
+ {
+ return (2 / constants::pi<T>()) * (log(z / 2) + constants::euler<T>());
+ }
+ else if(n == 1)
+ {
+ return (z / constants::pi<T>()) * log(z / 2)
+ - 2 / (constants::pi<T>() * z)
+ - (z / (2 * constants::pi<T>())) * (1 - 2 * constants::euler<T>());
+ }
+ else if(n == 2)
+ {
+ return (z * z) / (4 * constants::pi<T>()) * log(z / 2)
+ - (4 / (constants::pi<T>() * z * z))
+ - ((z * z) / (8 * constants::pi<T>())) * (T(3)/2 - 2 * constants::euler<T>());
+ }
+ else
+ {
+ T p = pow(z / 2, n);
+ T result = -((boost::math::factorial<T>(n - 1) / constants::pi<T>()));
+ if(p * tools::max_value<T>() < result)
+ {
+ T div = tools::max_value<T>() / 8;
+ result /= div;
+ *scale /= div;
+ if(p * tools::max_value<T>() < result)
+ {
+ return -policies::raise_overflow_error<T>("bessel_yn_small_z<%1%>(%1%,%1%)", 0, pol);
+ }
+ }
+ return result / p;
+ }
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_BESSEL_JN_SERIES_HPP
+
projects / rsem.git / blobdiff