boost/multiprecision/detail/functions/pow.hpp
// Copyright Christopher Kormanyos 2002 - 2013. // Copyright 2011 - 2013 John Maddock. Distributed under the Boost // Distributed under 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) // This work is based on an earlier work: // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations", // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469 // // This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp // #ifdef BOOST_MSVC #pragma warning(push) #pragma warning(disable:6326) // comparison of two constants #endif namespace detail{ template<typename T, typename U> inline void pow_imp(T& result, const T& t, const U& p, const mpl::false_&) { // Compute the pure power of typename T t^p. // Use the S-and-X binary method, as described in // D. E. Knuth, "The Art of Computer Programming", Vol. 2, // Section 4.6.3 . The resulting computational complexity // is order log2[abs(p)]. typedef typename boost::multiprecision::detail::canonical<U, T>::type int_type; if(&result == &t) { T temp; pow_imp(temp, t, p, mpl::false_()); result = temp; return; } // This will store the result. if(U(p % U(2)) != U(0)) { result = t; } else result = int_type(1); U p2(p); // The variable x stores the binary powers of t. T x(t); while(U(p2 /= 2) != U(0)) { // Square x for each binary power. eval_multiply(x, x); const bool has_binary_power = (U(p2 % U(2)) != U(0)); if(has_binary_power) { // Multiply the result with each binary power contained in the exponent. eval_multiply(result, x); } } } template<typename T, typename U> inline void pow_imp(T& result, const T& t, const U& p, const mpl::true_&) { // Signed integer power, just take care of the sign then call the unsigned version: typedef typename boost::multiprecision::detail::canonical<U, T>::type int_type; typedef typename make_unsigned<U>::type ui_type; if(p < 0) { T temp; temp = static_cast<int_type>(1); T denom; pow_imp(denom, t, static_cast<ui_type>(-p), mpl::false_()); eval_divide(result, temp, denom); return; } pow_imp(result, t, static_cast<ui_type>(p), mpl::false_()); } } // namespace detail template<typename T, typename U> inline typename enable_if<is_integral<U> >::type eval_pow(T& result, const T& t, const U& p) { detail::pow_imp(result, t, p, boost::is_signed<U>()); } template <class T> void hyp0F0(T& H0F0, const T& x) { // Compute the series representation of Hypergeometric0F0 taken from // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F0/06/01/ // There are no checks on input range or parameter boundaries. typedef typename mpl::front<typename T::unsigned_types>::type ui_type; BOOST_ASSERT(&H0F0 != &x); long tol = boost::multiprecision::detail::digits2<number<T, et_on> >::value(); T t; T x_pow_n_div_n_fact(x); eval_add(H0F0, x_pow_n_div_n_fact, ui_type(1)); T lim; eval_ldexp(lim, H0F0, 1 - tol); if(eval_get_sign(lim) < 0) lim.negate(); ui_type n; const unsigned series_limit = boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100 ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value(); // Series expansion of hyperg_0f0(; ; x). for(n = 2; n < series_limit; ++n) { eval_multiply(x_pow_n_div_n_fact, x); eval_divide(x_pow_n_div_n_fact, n); eval_add(H0F0, x_pow_n_div_n_fact); bool neg = eval_get_sign(x_pow_n_div_n_fact) < 0; if(neg) x_pow_n_div_n_fact.negate(); if(lim.compare(x_pow_n_div_n_fact) > 0) break; if(neg) x_pow_n_div_n_fact.negate(); } if(n >= series_limit) BOOST_THROW_EXCEPTION(std::runtime_error("H0F0 failed to converge")); } template <class T> void hyp1F0(T& H1F0, const T& a, const T& x) { // Compute the series representation of Hypergeometric1F0 taken from // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F0/06/01/01/ // and also see the corresponding section for the power function (i.e. x^a). // There are no checks on input range or parameter boundaries. typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type; BOOST_ASSERT(&H1F0 != &x); BOOST_ASSERT(&H1F0 != &a); T x_pow_n_div_n_fact(x); T pochham_a (a); T ap (a); eval_multiply(H1F0, pochham_a, x_pow_n_div_n_fact); eval_add(H1F0, si_type(1)); T lim; eval_ldexp(lim, H1F0, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value()); if(eval_get_sign(lim) < 0) lim.negate(); si_type n; T term, part; const si_type series_limit = boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100 ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value(); // Series expansion of hyperg_1f0(a; ; x). for(n = 2; n < series_limit; n++) { eval_multiply(x_pow_n_div_n_fact, x); eval_divide(x_pow_n_div_n_fact, n); eval_increment(ap); eval_multiply(pochham_a, ap); eval_multiply(term, pochham_a, x_pow_n_div_n_fact); eval_add(H1F0, term); if(eval_get_sign(term) < 0) term.negate(); if(lim.compare(term) >= 0) break; } if(n >= series_limit) BOOST_THROW_EXCEPTION(std::runtime_error("H1F0 failed to converge")); } template <class T> void eval_exp(T& result, const T& x) { BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The exp function is only valid for floating point types."); if(&x == &result) { T temp; eval_exp(temp, x); result = temp; return; } typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type; typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type; typedef typename T::exponent_type exp_type; typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type; // Handle special arguments. int type = eval_fpclassify(x); bool isneg = eval_get_sign(x) < 0; if(type == (int)FP_NAN) { result = x; errno = EDOM; return; } else if(type == (int)FP_INFINITE) { if(isneg) result = ui_type(0u); else result = x; return; } else if(type == (int)FP_ZERO) { result = ui_type(1); return; } // Get local copy of argument and force it to be positive. T xx = x; T exp_series; if(isneg) xx.negate(); // Check the range of the argument. if(xx.compare(si_type(1)) <= 0) { // // Use series for exp(x) - 1: // T lim; if(std::numeric_limits<number<T, et_on> >::is_specialized) lim = std::numeric_limits<number<T, et_on> >::epsilon().backend(); else { result = ui_type(1); eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value()); } unsigned k = 2; exp_series = xx; result = si_type(1); if(isneg) eval_subtract(result, exp_series); else eval_add(result, exp_series); eval_multiply(exp_series, xx); eval_divide(exp_series, ui_type(k)); eval_add(result, exp_series); while(exp_series.compare(lim) > 0) { ++k; eval_multiply(exp_series, xx); eval_divide(exp_series, ui_type(k)); if(isneg && (k&1)) eval_subtract(result, exp_series); else eval_add(result, exp_series); } return; } // Check for pure-integer arguments which can be either signed or unsigned. typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type ll; eval_trunc(exp_series, x); eval_convert_to(&ll, exp_series); if(x.compare(ll) == 0) { detail::pow_imp(result, get_constant_e<T>(), ll, mpl::true_()); return; } else if(exp_series.compare(x) == 0) { // We have a value that has no fractional part, but is too large to fit // in a long long, in this situation the code below will fail, so // we're just going to assume that this will overflow: if(isneg) result = ui_type(0); else result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend(); return; } // The algorithm for exp has been taken from MPFUN. // exp(t) = [ (1 + r + r^2/2! + r^3/3! + r^4/4! ...)^p2 ] * 2^n // where p2 is a power of 2 such as 2048, r = t_prime / p2, and // t_prime = t - n*ln2, with n chosen to minimize the absolute // value of t_prime. In the resulting Taylor series, which is // implemented as a hypergeometric function, |r| is bounded by // ln2 / p2. For small arguments, no scaling is done. // Compute the exponential series of the (possibly) scaled argument. eval_divide(result, xx, get_constant_ln2<T>()); exp_type n; eval_convert_to(&n, result); if (n == (std::numeric_limits<exp_type>::max)()) { // Exponent is too large to fit in our exponent type: if (isneg) result = ui_type(0); else result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend(); return; } // The scaling is 2^11 = 2048. const si_type p2 = static_cast<si_type>(si_type(1) << 11); eval_multiply(exp_series, get_constant_ln2<T>(), static_cast<canonical_exp_type>(n)); eval_subtract(exp_series, xx); eval_divide(exp_series, p2); exp_series.negate(); hyp0F0(result, exp_series); detail::pow_imp(exp_series, result, p2, mpl::true_()); result = ui_type(1); eval_ldexp(result, result, n); eval_multiply(exp_series, result); if(isneg) eval_divide(result, ui_type(1), exp_series); else result = exp_series; } template <class T> void eval_log(T& result, const T& arg) { BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types."); // // We use a variation of http://dlmf.nist.gov/4.45#i // using frexp to reduce the argument to x * 2^n, // then let y = x - 1 and compute: // log(x) = log(2) * n + log1p(1 + y) // typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type; typedef typename T::exponent_type exp_type; typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type; typedef typename mpl::front<typename T::float_types>::type fp_type; int s = eval_signbit(arg); switch(eval_fpclassify(arg)) { case FP_NAN: result = arg; errno = EDOM; return; case FP_INFINITE: if(s) break; result = arg; return; case FP_ZERO: result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend(); result.negate(); errno = ERANGE; return; } if(s) { result = std::numeric_limits<number<T> >::quiet_NaN().backend(); errno = EDOM; return; } exp_type e; T t; eval_frexp(t, arg, &e); bool alternate = false; if(t.compare(fp_type(2) / fp_type(3)) <= 0) { alternate = true; eval_ldexp(t, t, 1); --e; } eval_multiply(result, get_constant_ln2<T>(), canonical_exp_type(e)); INSTRUMENT_BACKEND(result); eval_subtract(t, ui_type(1)); /* -0.3 <= t <= 0.3 */ if(!alternate) t.negate(); /* 0 <= t <= 0.33333 */ T pow = t; T lim; T t2; if(alternate) eval_add(result, t); else eval_subtract(result, t); if(std::numeric_limits<number<T, et_on> >::is_specialized) eval_multiply(lim, result, std::numeric_limits<number<T, et_on> >::epsilon().backend()); else eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value()); if(eval_get_sign(lim) < 0) lim.negate(); INSTRUMENT_BACKEND(lim); ui_type k = 1; do { ++k; eval_multiply(pow, t); eval_divide(t2, pow, k); INSTRUMENT_BACKEND(t2); if(alternate && ((k & 1) != 0)) eval_add(result, t2); else eval_subtract(result, t2); INSTRUMENT_BACKEND(result); }while(lim.compare(t2) < 0); } template <class T> const T& get_constant_log10() { static BOOST_MP_THREAD_LOCAL T result; static BOOST_MP_THREAD_LOCAL bool b = false; static BOOST_MP_THREAD_LOCAL long digits = boost::multiprecision::detail::digits2<number<T> >::value(); if(!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value())) { typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type; T ten; ten = ui_type(10u); eval_log(result, ten); b = true; digits = boost::multiprecision::detail::digits2<number<T> >::value(); } constant_initializer<T, &get_constant_log10<T> >::do_nothing(); return result; } template <class T> void eval_log10(T& result, const T& arg) { BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log10 function is only valid for floating point types."); eval_log(result, arg); eval_divide(result, get_constant_log10<T>()); } template <class R, class T> inline void eval_log2(R& result, const T& a) { eval_log(result, a); eval_divide(result, get_constant_ln2<R>()); } template<typename T> inline void eval_pow(T& result, const T& x, const T& a) { BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The pow function is only valid for floating point types."); typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type; typedef typename mpl::front<typename T::float_types>::type fp_type; if((&result == &x) || (&result == &a)) { T t; eval_pow(t, x, a); result = t; return; } if((a.compare(si_type(1)) == 0) || (x.compare(si_type(1)) == 0)) { result = x; return; } if(a.compare(si_type(0)) == 0) { result = si_type(1); return; } int type = eval_fpclassify(x); switch(type) { case FP_ZERO: switch(eval_fpclassify(a)) { case FP_ZERO: result = si_type(1); break; case FP_NAN: result = a; break; case FP_NORMAL: { // Need to check for a an odd integer as a special case: try { typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type i; eval_convert_to(&i, a); if(a.compare(i) == 0) { if(eval_signbit(a)) { if(i & 1) { result = std::numeric_limits<number<T> >::infinity().backend(); if(eval_signbit(x)) result.negate(); errno = ERANGE; } else { result = std::numeric_limits<number<T> >::infinity().backend(); errno = ERANGE; } } else if(i & 1) { result = x; } else result = si_type(0); return; } } catch(const std::exception&) { // fallthrough.. } } default: if(eval_signbit(a)) { result = std::numeric_limits<number<T> >::infinity().backend(); errno = ERANGE; } else result = x; break; } return; case FP_NAN: result = x; errno = ERANGE; return; default: ; } int s = eval_get_sign(a); if(s == 0) { result = si_type(1); return; } if(s < 0) { T t, da; t = a; t.negate(); eval_pow(da, x, t); eval_divide(result, si_type(1), da); return; } typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type an; typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type max_an = std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ? (std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::max)() : static_cast<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>(1) << (sizeof(typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type) * CHAR_BIT - 2); typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type min_an = std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ? (std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::min)() : -min_an; T fa; #ifndef BOOST_NO_EXCEPTIONS try { #endif eval_convert_to(&an, a); if(a.compare(an) == 0) { detail::pow_imp(result, x, an, mpl::true_()); return; } #ifndef BOOST_NO_EXCEPTIONS } catch(const std::exception&) { // conversion failed, just fall through, value is not an integer. an = (std::numeric_limits<boost::intmax_t>::max)(); } #endif if((eval_get_sign(x) < 0)) { typename boost::multiprecision::detail::canonical<boost::uintmax_t, T>::type aun; #ifndef BOOST_NO_EXCEPTIONS try { #endif eval_convert_to(&aun, a); if(a.compare(aun) == 0) { fa = x; fa.negate(); eval_pow(result, fa, a); if(aun & 1u) result.negate(); return; } #ifndef BOOST_NO_EXCEPTIONS } catch(const std::exception&) { // conversion failed, just fall through, value is not an integer. } #endif eval_floor(result, a); // -1^INF is a special case in C99: if((x.compare(si_type(-1)) == 0) && (eval_fpclassify(a) == FP_INFINITE)) { result = si_type(1); } else if(a.compare(result) == 0) { // exponent is so large we have no fractional part: if(x.compare(si_type(-1)) < 0) { result = std::numeric_limits<number<T, et_on> >::infinity().backend(); } else { result = si_type(0); } } else if(type == FP_INFINITE) { result = std::numeric_limits<number<T, et_on> >::infinity().backend(); } else if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN) { result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend(); errno = EDOM; } else { BOOST_THROW_EXCEPTION(std::domain_error("Result of pow is undefined or non-real and there is no NaN for this number type.")); } return; } T t, da; eval_subtract(da, a, an); if((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0) && (an < max_an) && (an > min_an)) { if(a.compare(fp_type(1e-5f)) <= 0) { // Series expansion for small a. eval_log(t, x); eval_multiply(t, a); hyp0F0(result, t); return; } else { // Series expansion for moderately sized x. Note that for large power of a, // the power of the integer part of a is calculated using the pown function. if(an) { da.negate(); t = si_type(1); eval_subtract(t, x); hyp1F0(result, da, t); detail::pow_imp(t, x, an, mpl::true_()); eval_multiply(result, t); } else { da = a; da.negate(); t = si_type(1); eval_subtract(t, x); hyp1F0(result, da, t); } } } else { // Series expansion for pow(x, a). Note that for large power of a, the power // of the integer part of a is calculated using the pown function. if(an) { eval_log(t, x); eval_multiply(t, da); eval_exp(result, t); detail::pow_imp(t, x, an, mpl::true_()); eval_multiply(result, t); } else { eval_log(t, x); eval_multiply(t, a); eval_exp(result, t); } } } template<class T, class A> inline typename enable_if<is_floating_point<A>, void>::type eval_pow(T& result, const T& x, const A& a) { // Note this one is restricted to float arguments since pow.hpp already has a version for // integer powers.... typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type; typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type; cast_type c; c = a; eval_pow(result, x, c); } template<class T, class A> inline typename enable_if<is_arithmetic<A>, void>::type eval_pow(T& result, const A& x, const T& a) { typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type; typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type; cast_type c; c = x; eval_pow(result, c, a); } template <class T> void eval_exp2(T& result, const T& arg) { BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types."); // Check for pure-integer arguments which can be either signed or unsigned. typename boost::multiprecision::detail::canonical<typename T::exponent_type, T>::type i; T temp; try { eval_trunc(temp, arg); eval_convert_to(&i, temp); if(arg.compare(i) == 0) { temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u); eval_ldexp(result, temp, i); return; } } catch(const boost::math::rounding_error&) { /* Fallthrough */ } catch(const std::runtime_error&) { /* Fallthrough */ } temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(2u); eval_pow(result, temp, arg); } namespace detail{ template <class T> void small_sinh_series(T x, T& result) { typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type; bool neg = eval_get_sign(x) < 0; if(neg) x.negate(); T p(x); T mult(x); eval_multiply(mult, x); result = x; ui_type k = 1; T lim(x); eval_ldexp(lim, lim, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value()); do { eval_multiply(p, mult); eval_divide(p, ++k); eval_divide(p, ++k); eval_add(result, p); }while(p.compare(lim) >= 0); if(neg) result.negate(); } template <class T> void sinhcosh(const T& x, T* p_sinh, T* p_cosh) { typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type; typedef typename mpl::front<typename T::float_types>::type fp_type; switch(eval_fpclassify(x)) { case FP_NAN: errno = EDOM; // fallthrough... case FP_INFINITE: if(p_sinh) *p_sinh = x; if(p_cosh) { *p_cosh = x; if(eval_get_sign(x) < 0) p_cosh->negate(); } return; case FP_ZERO: if(p_sinh) *p_sinh = x; if(p_cosh) *p_cosh = ui_type(1); return; default: ; } bool small_sinh = eval_get_sign(x) < 0 ? x.compare(fp_type(-0.5)) > 0 : x.compare(fp_type(0.5)) < 0; if(p_cosh || !small_sinh) { T e_px, e_mx; eval_exp(e_px, x); eval_divide(e_mx, ui_type(1), e_px); if(eval_signbit(e_mx) != eval_signbit(e_px)) e_mx.negate(); // Handles lack of signed zero in some types if(p_sinh) { if(small_sinh) { small_sinh_series(x, *p_sinh); } else { eval_subtract(*p_sinh, e_px, e_mx); eval_ldexp(*p_sinh, *p_sinh, -1); } } if(p_cosh) { eval_add(*p_cosh, e_px, e_mx); eval_ldexp(*p_cosh, *p_cosh, -1); } } else { small_sinh_series(x, *p_sinh); } } } // namespace detail template <class T> inline void eval_sinh(T& result, const T& x) { BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sinh function is only valid for floating point types."); detail::sinhcosh(x, &result, static_cast<T*>(0)); } template <class T> inline void eval_cosh(T& result, const T& x) { BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cosh function is only valid for floating point types."); detail::sinhcosh(x, static_cast<T*>(0), &result); } template <class T> inline void eval_tanh(T& result, const T& x) { BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tanh function is only valid for floating point types."); T c; detail::sinhcosh(x, &result, &c); if((eval_fpclassify(result) == FP_INFINITE) && (eval_fpclassify(c) == FP_INFINITE)) { bool s = eval_signbit(result) != eval_signbit(c); result = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u); if(s) result.negate(); return; } eval_divide(result, c); } #ifdef BOOST_MSVC #pragma warning(pop) #endif