From 0d1d5016b57c59817fa0619e5cf4c6df046b28c3 Mon Sep 17 00:00:00 2001 From: Mark Benvenuto Date: Fri, 15 Aug 2014 17:06:16 -0400 Subject: SERVER-8994: Boost 1.56.0 Initial import of source Libraries: algorithm array asio bind chrono config container date_time filesystem function integer intrusive noncopyable optional program_options random smart_ptr static_assert thread unordered utility --- .../boost/math/special_functions/factorials.hpp | 268 +++++++++++++++++++++ 1 file changed, 268 insertions(+) create mode 100644 src/third_party/boost-1.56.0/boost/math/special_functions/factorials.hpp (limited to 'src/third_party/boost-1.56.0/boost/math/special_functions/factorials.hpp') diff --git a/src/third_party/boost-1.56.0/boost/math/special_functions/factorials.hpp b/src/third_party/boost-1.56.0/boost/math/special_functions/factorials.hpp new file mode 100644 index 00000000000..de24642ac41 --- /dev/null +++ b/src/third_party/boost-1.56.0/boost/math/special_functions/factorials.hpp @@ -0,0 +1,268 @@ +// Copyright John Maddock 2006, 2010. +// 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_SP_FACTORIALS_HPP +#define BOOST_MATH_SP_FACTORIALS_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include +#include +#include +#include +#ifdef BOOST_MSVC +#pragma warning(push) // Temporary until lexical cast fixed. +#pragma warning(disable: 4127 4701) +#endif +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif +#include + +namespace boost { namespace math +{ + +template +inline T factorial(unsigned i, const Policy& pol) +{ + BOOST_STATIC_ASSERT(!boost::is_integral::value); + // factorial(n) is not implemented + // because it would overflow integral type T for too small n + // to be useful. Use instead a floating-point type, + // and convert to an unsigned type if essential, for example: + // unsigned int nfac = static_cast(factorial(n)); + // See factorial documentation for more detail. + + BOOST_MATH_STD_USING // Aid ADL for floor. + + if(i <= max_factorial::value) + return unchecked_factorial(i); + T result = boost::math::tgamma(static_cast(i+1), pol); + if(result > tools::max_value()) + return result; // Overflowed value! (But tgamma will have signalled the error already). + return floor(result + 0.5f); +} + +template +inline T factorial(unsigned i) +{ + return factorial(i, policies::policy<>()); +} +/* +// Can't have these in a policy enabled world? +template<> +inline float factorial(unsigned i) +{ + if(i <= max_factorial::value) + return unchecked_factorial(i); + return tools::overflow_error(BOOST_CURRENT_FUNCTION); +} + +template<> +inline double factorial(unsigned i) +{ + if(i <= max_factorial::value) + return unchecked_factorial(i); + return tools::overflow_error(BOOST_CURRENT_FUNCTION); +} +*/ +template +T double_factorial(unsigned i, const Policy& pol) +{ + BOOST_STATIC_ASSERT(!boost::is_integral::value); + BOOST_MATH_STD_USING // ADL lookup of std names + if(i & 1) + { + // odd i: + if(i < max_factorial::value) + { + unsigned n = (i - 1) / 2; + return ceil(unchecked_factorial(i) / (ldexp(T(1), (int)n) * unchecked_factorial(n)) - 0.5f); + } + // + // Fallthrough: i is too large to use table lookup, try the + // gamma function instead. + // + T result = boost::math::tgamma(static_cast(i) / 2 + 1, pol) / sqrt(constants::pi()); + if(ldexp(tools::max_value(), -static_cast(i+1) / 2) > result) + return ceil(result * ldexp(T(1), static_cast(i+1) / 2) - 0.5f); + } + else + { + // even i: + unsigned n = i / 2; + T result = factorial(n, pol); + if(ldexp(tools::max_value(), -(int)n) > result) + return result * ldexp(T(1), (int)n); + } + // + // If we fall through to here then the result is infinite: + // + return policies::raise_overflow_error("boost::math::double_factorial<%1%>(unsigned)", 0, pol); +} + +template +inline T double_factorial(unsigned i) +{ + return double_factorial(i, policies::policy<>()); +} + +namespace detail{ + +template +T rising_factorial_imp(T x, int n, const Policy& pol) +{ + BOOST_STATIC_ASSERT(!boost::is_integral::value); + if(x < 0) + { + // + // For x less than zero, we really have a falling + // factorial, modulo a possible change of sign. + // + // Note that the falling factorial isn't defined + // for negative n, so we'll get rid of that case + // first: + // + bool inv = false; + if(n < 0) + { + x += n; + n = -n; + inv = true; + } + T result = ((n&1) ? -1 : 1) * falling_factorial(-x, n, pol); + if(inv) + result = 1 / result; + return result; + } + if(n == 0) + return 1; + if(x == 0) + { + if(n < 0) + return -boost::math::tgamma_delta_ratio(x + 1, static_cast(-n), pol); + else + return 0; + } + if((x < 1) && (x + n < 0)) + { + T val = boost::math::tgamma_delta_ratio(1 - x, static_cast(-n), pol); + return (n & 1) ? -val : val; + } + // + // We don't optimise this for small n, because + // tgamma_delta_ratio is alreay optimised for that + // use case: + // + return 1 / boost::math::tgamma_delta_ratio(x, static_cast(n), pol); +} + +template +inline T falling_factorial_imp(T x, unsigned n, const Policy& pol) +{ + BOOST_STATIC_ASSERT(!boost::is_integral::value); + BOOST_MATH_STD_USING // ADL of std names + if((x == 0) && (n >= 0)) + return 0; + if(x < 0) + { + // + // For x < 0 we really have a rising factorial + // modulo a possible change of sign: + // + return (n&1 ? -1 : 1) * rising_factorial(-x, n, pol); + } + if(n == 0) + return 1; + if(x < 0.5f) + { + // + // 1 + x below will throw away digits, so split up calculation: + // + if(n > max_factorial::value - 2) + { + // If the two end of the range are far apart we have a ratio of two very large + // numbers, split the calculation up into two blocks: + T t1 = x * boost::math::falling_factorial(x - 1, max_factorial::value - 2); + T t2 = boost::math::falling_factorial(x - max_factorial::value + 1, n - max_factorial::value + 1); + if(tools::max_value() / fabs(t1) < fabs(t2)) + return boost::math::sign(t1) * boost::math::sign(t2) * policies::raise_overflow_error("boost::math::falling_factorial<%1%>", 0, pol); + return t1 * t2; + } + return x * boost::math::falling_factorial(x - 1, n - 1); + } + if(x <= n - 1) + { + // + // x+1-n will be negative and tgamma_delta_ratio won't + // handle it, split the product up into three parts: + // + T xp1 = x + 1; + unsigned n2 = itrunc((T)floor(xp1), pol); + if(n2 == xp1) + return 0; + T result = boost::math::tgamma_delta_ratio(xp1, -static_cast(n2), pol); + x -= n2; + result *= x; + ++n2; + if(n2 < n) + result *= falling_factorial(x - 1, n - n2, pol); + return result; + } + // + // Simple case: just the ratio of two + // (positive argument) gamma functions. + // Note that we don't optimise this for small n, + // because tgamma_delta_ratio is alreay optimised + // for that use case: + // + return boost::math::tgamma_delta_ratio(x + 1, -static_cast(n), pol); +} + +} // namespace detail + +template +inline typename tools::promote_args::type + falling_factorial(RT x, unsigned n) +{ + typedef typename tools::promote_args::type result_type; + return detail::falling_factorial_imp( + static_cast(x), n, policies::policy<>()); +} + +template +inline typename tools::promote_args::type + falling_factorial(RT x, unsigned n, const Policy& pol) +{ + typedef typename tools::promote_args::type result_type; + return detail::falling_factorial_imp( + static_cast(x), n, pol); +} + +template +inline typename tools::promote_args::type + rising_factorial(RT x, int n) +{ + typedef typename tools::promote_args::type result_type; + return detail::rising_factorial_imp( + static_cast(x), n, policies::policy<>()); +} + +template +inline typename tools::promote_args::type + rising_factorial(RT x, int n, const Policy& pol) +{ + typedef typename tools::promote_args::type result_type; + return detail::rising_factorial_imp( + static_cast(x), n, pol); +} + +} // namespace math +} // namespace boost + +#endif // BOOST_MATH_SP_FACTORIALS_HPP + -- cgit v1.2.1