+
// Copyright John Maddock 2006-7.
// Copyright Paul A. Bristow 2007.
#endif
#include <boost/config.hpp>
-#ifdef BOOST_MSVC
-# pragma warning(push)
-# pragma warning(disable: 4127 4701)
-// // For lexical_cast, until fixed in 1.35?
-// // conditional expression is constant &
-// // Potentially uninitialized local variable 'name' used
-#endif
-#include <boost/lexical_cast.hpp>
-#ifdef BOOST_MSVC
-# pragma warning(pop)
-#endif
#include <boost/math/tools/series.hpp>
#include <boost/math/tools/fraction.hpp>
#include <boost/math/tools/precision.hpp>
#include <boost/config/no_tr1/cmath.hpp>
#include <algorithm>
-#ifdef BOOST_MATH_INSTRUMENT
-#include <iostream>
-#include <iomanip>
-#include <typeinfo>
-#endif
-
#ifdef BOOST_MSVC
# pragma warning(push)
# pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
//
// tgamma(z), with Lanczos support:
//
-template <class T, class Policy, class L>
-T gamma_imp(T z, const Policy& pol, const L& l)
+template <class T, class Policy, class Lanczos>
+T gamma_imp(T z, const Policy& pol, const Lanczos& l)
{
BOOST_MATH_STD_USING
if(z <= -20)
{
result = gamma_imp(T(-z), pol, l) * sinpx(z);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result = -boost::math::constants::pi<T>() / result;
return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
return result;
}
z += 1;
}
}
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
if((floor(z) == z) && (z < max_factorial<T>::value))
{
result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
{
- result *= L::lanczos_sum(z);
- if(z * log(z) > tools::log_max_value<T>())
+ result *= Lanczos::lanczos_sum(z);
+ T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
+ T lzgh = log(zgh);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
+ BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
+ if(z * lzgh > tools::log_max_value<T>())
{
// we're going to overflow unless this is done with care:
- T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
- if(log(zgh) * z / 2 > tools::log_max_value<T>())
+ BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
+ if(lzgh * z / 2 > tools::log_max_value<T>())
return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
T hp = pow(zgh, (z / 2) - T(0.25));
+ BOOST_MATH_INSTRUMENT_VARIABLE(hp);
result *= hp / exp(zgh);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
if(tools::max_value<T>() / hp < result)
return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result *= hp;
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
{
- T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
+ BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
+ BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>()));
+ BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));
result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
}
return result;
//
// lgamma(z) with Lanczos support:
//
-template <class T, class Policy, class L>
-T lgamma_imp(T z, const Policy& pol, const L& l, int* sign = 0)
+template <class T, class Policy, class Lanczos>
+T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
{
#ifdef BOOST_MATH_INSTRUMENT
static bool b = false;
>,
mpl::int_<113>, mpl::int_<0> >::type
>::type tag_type;
- result = lgamma_small_imp<T>(z, z - 1, z - 2, tag_type(), pol, l);
+ result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);
}
- else if((z >= 3) && (z < 100))
+ else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024))
{
// taking the log of tgamma reduces the error, no danger of overflow here:
result = log(gamma_imp(z, pol, l));
else
{
// regular evaluation:
- T zgh = static_cast<T>(z + L::g() - boost::math::constants::half<T>());
+ T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>());
result = log(zgh) - 1;
result *= z - 0.5f;
- result += log(L::lanczos_sum_expG_scaled(z));
+ result += log(Lanczos::lanczos_sum_expG_scaled(z));
}
if(sign)
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
T factor = policies::get_epsilon<T, Policy>();
T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
- policies::check_series_iterations("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
+ policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
return result;
}
return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
if(z <= -20)
{
- T result = gamma_imp(-z, pol, l) * sinpx(z);
+ T result = gamma_imp(T(-z), pol, l) * sinpx(z);
if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result = -boost::math::constants::pi<T>() / result;
}
else if((z != 1) && (z != 2))
{
- T limit = (std::max)(z+1, T(10));
+ T limit = (std::max)(T(z+1), T(10));
T prefix = z * log(limit) - limit;
T sum = detail::lower_gamma_series(z, limit, pol) / z;
sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policies::get_epsilon<T, Policy>());
// This helper calculates tgamma(dz+1)-1 without cancellation errors,
// used by the upper incomplete gamma with z < 1:
//
-template <class T, class Policy, class L>
-T tgammap1m1_imp(T dz, Policy const& pol, const L& l)
+template <class T, class Policy, class Lanczos>
+T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
{
BOOST_MATH_STD_USING
mpl::greater<precision_type, mpl::int_<113> >
>,
typename mpl::if_<
- is_same<L, lanczos::lanczos24m113>,
+ is_same<Lanczos, lanczos::lanczos24m113>,
mpl::int_<113>,
mpl::int_<0>
>::type,
// algebra isn't easy for the general case....
// Start by subracting 1 from tgamma:
//
- T result = gamma_imp(1 + dz, pol, l) - 1;
+ T result = gamma_imp(T(1 + dz), pol, l) - 1;
BOOST_MATH_INSTRUMENT_CODE(result);
//
// Test the level of cancellation error observed: we loose one bit
// Compute (z^a)(e^-z)/tgamma(a)
// most if the error occurs in this function:
//
-template <class T, class Policy, class L>
-T regularised_gamma_prefix(T a, T z, const Policy& pol, const L& l)
+template <class T, class Policy, class Lanczos>
+T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
{
BOOST_MATH_STD_USING
- T agh = a + static_cast<T>(L::g()) - T(0.5);
+ T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
T prefix;
- T d = ((z - a) - static_cast<T>(L::g()) + T(0.5)) / agh;
+ T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;
if(a < 1)
{
else if((fabs(d*d*a) <= 100) && (a > 150))
{
// special case for large a and a ~ z.
- prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - L::g()) / agh;
+ prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh;
prefix = exp(prefix);
}
else
prefix = pow(z / agh, a) * exp(amz);
}
}
- prefix *= sqrt(agh / boost::math::constants::e<T>()) / L::lanczos_sum_expG_scaled(a);
+ prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a);
return prefix;
}
//
*pderivative = p / (*pgam * exp(x));
T init_value = invert ? *pgam : 0;
result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
- policies::check_series_iterations("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
+ policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
if(invert)
result = -result;
return result;
typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
- T result;
+ T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
BOOST_ASSERT((p_derivative == 0) || (normalised == true));
bool is_int, is_half_int;
- bool is_small_a = (a < 30) && (a <= x + 1);
+ bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());
if(is_small_a)
{
T fa = floor(a);
//
// Ratios of two gamma functions:
//
-template <class T, class Policy, class L>
-T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const L&)
+template <class T, class Policy, class Lanczos>
+T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos&)
{
BOOST_MATH_STD_USING
- T zgh = z + L::g() - constants::half<T>();
+ T zgh = z + Lanczos::g() - constants::half<T>();
T result;
if(fabs(delta) < 10)
{
result = pow(zgh / (zgh + delta), z - constants::half<T>());
}
result *= pow(constants::e<T>() / (zgh + delta), delta);
- result *= L::lanczos_sum(z) / L::lanczos_sum(z + delta);
+ result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
return result;
}
//
return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
}
+template <class T, class Policy>
+T tgamma_ratio_imp(T x, T y, const Policy& pol)
+{
+ BOOST_MATH_STD_USING
+
+ if((x <= tools::min_value<T>()) || (boost::math::isinf)(x))
+ policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", x, pol);
+ if((y <= tools::min_value<T>()) || (boost::math::isinf)(y))
+ policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", y, pol);
+
+ if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
+ {
+ // Rather than subtracting values, lets just call the gamma functions directly:
+ return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+ }
+ T prefix = 1;
+ if(x < 1)
+ {
+ if(y < 2 * max_factorial<T>::value)
+ {
+ // We need to sidestep on x as well, otherwise we'll underflow
+ // before we get to factor in the prefix term:
+ prefix /= x;
+ x += 1;
+ while(y >= max_factorial<T>::value)
+ {
+ y -= 1;
+ prefix /= y;
+ }
+ return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+ }
+ //
+ // result is almost certainly going to underflow to zero, try logs just in case:
+ //
+ return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
+ }
+ if(y < 1)
+ {
+ if(x < 2 * max_factorial<T>::value)
+ {
+ // We need to sidestep on y as well, otherwise we'll overflow
+ // before we get to factor in the prefix term:
+ prefix *= y;
+ y += 1;
+ while(x >= max_factorial<T>::value)
+ {
+ x -= 1;
+ prefix *= x;
+ }
+ return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+ }
+ //
+ // Result will almost certainly overflow, try logs just in case:
+ //
+ return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
+ }
+ //
+ // Regular case, x and y both large and similar in magnitude:
+ //
+ return boost::math::tgamma_delta_ratio(x, y - x, pol);
+}
+
template <class T, class Policy>
T gamma_p_derivative_imp(T a, T x, const Policy& pol)
{
return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma<%1%>(%1%)");
}
+template <class T, class Policy>
+struct igamma_initializer
+{
+ struct init
+ {
+ init()
+ {
+ typedef typename policies::precision<T, Policy>::type precision_type;
+
+ typedef typename mpl::if_<
+ mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
+ mpl::greater<precision_type, mpl::int_<113> > >,
+ mpl::int_<0>,
+ typename mpl::if_<
+ mpl::less_equal<precision_type, mpl::int_<53> >,
+ mpl::int_<53>,
+ typename mpl::if_<
+ mpl::less_equal<precision_type, mpl::int_<64> >,
+ mpl::int_<64>,
+ mpl::int_<113>
+ >::type
+ >::type
+ >::type tag_type;
+
+ do_init(tag_type());
+ }
+ template <int N>
+ static void do_init(const mpl::int_<N>&)
+ {
+ boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
+ }
+ static void do_init(const mpl::int_<53>&){}
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T, class Policy>
+const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;
+
+template <class T, class Policy>
+struct lgamma_initializer
+{
+ struct init
+ {
+ init()
+ {
+ typedef typename policies::precision<T, Policy>::type precision_type;
+ typedef typename mpl::if_<
+ mpl::and_<
+ mpl::less_equal<precision_type, mpl::int_<64> >,
+ mpl::greater<precision_type, mpl::int_<0> >
+ >,
+ mpl::int_<64>,
+ typename mpl::if_<
+ mpl::and_<
+ mpl::less_equal<precision_type, mpl::int_<113> >,
+ mpl::greater<precision_type, mpl::int_<0> >
+ >,
+ mpl::int_<113>, mpl::int_<0> >::type
+ >::type tag_type;
+ do_init(tag_type());
+ }
+ static void do_init(const mpl::int_<64>&)
+ {
+ boost::math::lgamma(static_cast<T>(2.5), Policy());
+ boost::math::lgamma(static_cast<T>(1.25), Policy());
+ boost::math::lgamma(static_cast<T>(1.75), Policy());
+ }
+ static void do_init(const mpl::int_<113>&)
+ {
+ boost::math::lgamma(static_cast<T>(2.5), Policy());
+ boost::math::lgamma(static_cast<T>(1.25), Policy());
+ boost::math::lgamma(static_cast<T>(1.5), Policy());
+ boost::math::lgamma(static_cast<T>(1.75), Policy());
+ }
+ static void do_init(const mpl::int_<0>&)
+ {
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T, class Policy>
+const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;
+
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
- typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+ // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
+
+ igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
return policies::checked_narrowing_cast<result_type, forwarding_policy>(
detail::gamma_incomplete_imp(static_cast<value_type>(a),
static_cast<value_type>(z), false, true,
return tgamma(a, z, policies::policy<>(), tag);
}
+
} // namespace detail
template <class T>
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
+
+ detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::lgamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)");
}
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
- typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+ // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
+ detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
return policies::checked_narrowing_cast<result_type, forwarding_policy>(
detail::gamma_incomplete_imp(static_cast<value_type>(a),
static_cast<value_type>(z), false, false,
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
- typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+ // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
+ detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
return policies::checked_narrowing_cast<result_type, forwarding_policy>(
detail::gamma_incomplete_imp(static_cast<value_type>(a),
static_cast<value_type>(z), true, true,
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
- typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+ // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
+ detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
return policies::checked_narrowing_cast<result_type, forwarding_policy>(
detail::gamma_incomplete_imp(static_cast<value_type>(a),
static_cast<value_type>(z), true, false,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
- return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(static_cast<value_type>(b) - static_cast<value_type>(a)), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
+ return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
projects / rsem.git / blobdiff