diff options
Diffstat (limited to 'boost/multiprecision')
-rw-r--r-- | boost/multiprecision/cpp_bin_float.hpp | 253 | ||||
-rw-r--r-- | boost/multiprecision/cpp_bin_float/transcendental.hpp | 1 | ||||
-rw-r--r-- | boost/multiprecision/cpp_dec_float.hpp | 41 | ||||
-rw-r--r-- | boost/multiprecision/cpp_int/add.hpp | 4 | ||||
-rw-r--r-- | boost/multiprecision/cpp_int/bitwise.hpp | 12 | ||||
-rw-r--r-- | boost/multiprecision/cpp_int/divide.hpp | 2 | ||||
-rw-r--r-- | boost/multiprecision/cpp_int/import_export.hpp | 8 | ||||
-rw-r--r-- | boost/multiprecision/cpp_int/limits.hpp | 4 | ||||
-rw-r--r-- | boost/multiprecision/cpp_int/multiply.hpp | 4 | ||||
-rw-r--r-- | boost/multiprecision/debug_adaptor.hpp | 9 | ||||
-rw-r--r-- | boost/multiprecision/detail/bitscan.hpp | 32 | ||||
-rw-r--r-- | boost/multiprecision/detail/default_ops.hpp | 96 | ||||
-rw-r--r-- | boost/multiprecision/detail/functions/pow.hpp | 164 | ||||
-rw-r--r-- | boost/multiprecision/detail/functions/trig.hpp | 65 | ||||
-rw-r--r-- | boost/multiprecision/float128.hpp | 20 | ||||
-rw-r--r-- | boost/multiprecision/logged_adaptor.hpp | 7 | ||||
-rw-r--r-- | boost/multiprecision/mpfi.hpp | 8 | ||||
-rw-r--r-- | boost/multiprecision/mpfr.hpp | 32 | ||||
-rw-r--r-- | boost/multiprecision/number.hpp | 76 |
19 files changed, 604 insertions, 234 deletions
diff --git a/boost/multiprecision/cpp_bin_float.hpp b/boost/multiprecision/cpp_bin_float.hpp index 38f5387528..f3e593a142 100644 --- a/boost/multiprecision/cpp_bin_float.hpp +++ b/boost/multiprecision/cpp_bin_float.hpp @@ -21,6 +21,10 @@ #include <boost/math/special_functions/expm1.hpp> #include <boost/math/special_functions/gamma.hpp> +#ifdef BOOST_HAS_FLOAT128 +#include <quadmath.h> +#endif + namespace boost{ namespace multiprecision{ namespace backends{ enum digit_base_type @@ -163,6 +167,64 @@ public: return assign_float(f); } +#ifdef BOOST_HAS_FLOAT128 + cpp_bin_float& operator=(__float128 f) + { + using default_ops::eval_add; + typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type; + + if(f == 0) + { + m_data = limb_type(0); + m_sign = (signbitq(f) > 0); + m_exponent = exponent_zero; + return *this; + } + else if(isnanq(f)) + { + m_data = limb_type(0); + m_sign = false; + m_exponent = exponent_nan; + return *this; + } + else if(isinfq(f)) + { + m_data = limb_type(0); + m_sign = (f < 0); + m_exponent = exponent_infinity; + return *this; + } + if(f < 0) + { + *this = -f; + this->negate(); + return *this; + } + + typedef typename mpl::front<unsigned_types>::type ui_type; + m_data = static_cast<ui_type>(0u); + m_sign = false; + m_exponent = 0; + + static const int bits = sizeof(int) * CHAR_BIT - 1; + int e; + f = frexpq(f, &e); + while(f) + { + f = ldexpq(f, bits); + e -= bits; + int ipart = (int)truncq(f); + f -= ipart; + m_exponent += bits; + cpp_bin_float t; + t = static_cast<bf_int_type>(ipart); + eval_add(*this, t); + } + m_exponent += static_cast<Exponent>(e); + return *this; + } +#endif + template <class Float> typename boost::enable_if_c<is_floating_point<Float>::value, cpp_bin_float&>::type assign_float(Float f) { @@ -414,10 +476,10 @@ public: #endif template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class Int> -inline void copy_and_round(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, Int &arg) +inline void copy_and_round(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, Int &arg, int bits_to_keep = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) { // Precondition: exponent of res must have been set before this function is called - // as we may need to adjust it based on how many cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in arg are set. + // as we may need to adjust it based on how many bits_to_keep in arg are set. using default_ops::eval_msb; using default_ops::eval_lsb; using default_ops::eval_left_shift; @@ -435,45 +497,70 @@ inline void copy_and_round(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, return; } int msb = eval_msb(arg); - if(static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) > msb + 1) + if(static_cast<int>(bits_to_keep) > msb + 1) { - // Must have had cancellation in subtraction, shift left and copy: + // Must have had cancellation in subtraction, + // or be converting from a narrower type, so shift left: res.bits() = arg; - eval_left_shift(res.bits(), cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb - 1); - res.exponent() -= static_cast<Exponent>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb - 1); + eval_left_shift(res.bits(), bits_to_keep - msb - 1); + res.exponent() -= static_cast<Exponent>(bits_to_keep - msb - 1); } - else if(static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) < msb + 1) + else if(static_cast<int>(bits_to_keep) < msb + 1) { - // We have more cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count than we need, so round as required, + // We have more bits_to_keep than we need, so round as required, // first get the rounding bit: - bool roundup = eval_bit_test(arg, msb - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count); + bool roundup = eval_bit_test(arg, msb - bits_to_keep); // Then check for a tie: - if(roundup && (msb - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count == eval_lsb(arg))) + if(roundup && (msb - bits_to_keep == (int)eval_lsb(arg))) { // Ties round towards even: - if(!eval_bit_test(arg, msb - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1)) + if(!eval_bit_test(arg, msb - bits_to_keep + 1)) roundup = false; } - // Shift off the cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count we don't need: - eval_right_shift(arg, msb - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1); - res.exponent() += static_cast<Exponent>(msb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1); + // Shift off the bits_to_keep we don't need: + eval_right_shift(arg, msb - bits_to_keep + 1); + res.exponent() += static_cast<Exponent>(msb - bits_to_keep + 1); if(roundup) { eval_increment(arg); - if(eval_bit_test(arg, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)) + if(bits_to_keep) + { + if(eval_bit_test(arg, bits_to_keep)) + { + // This happens very very rairly, all the bits left after + // truncation must be 1's and we're rounding up an order of magnitude: + eval_right_shift(arg, 1u); + ++res.exponent(); + } + } + else { - // This happens very very rairly: - eval_right_shift(arg, 1u); - ++res.exponent(); + // We get here when bits_to_keep is zero but we're rounding up, + // as a result we end up with a single digit that is a 1: + ++bits_to_keep; } } + if(bits_to_keep != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) + { + // Normalize result when we're rounding to fewer bits than we can hold, only happens in conversions + // to narrower types: + eval_left_shift(arg, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - bits_to_keep); + res.exponent() -= static_cast<Exponent>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - bits_to_keep); + } res.bits() = arg; } else { res.bits() = arg; } - BOOST_ASSERT((eval_msb(res.bits()) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)); + if(!bits_to_keep && !res.bits().limbs()[0]) + { + // We're keeping zero bits and did not round up, so result is zero: + res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero; + return; + } + // Result must be normalized: + BOOST_ASSERT(((int)eval_msb(res.bits()) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)); if(res.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) { @@ -571,6 +658,7 @@ inline void do_eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponen { using default_ops::eval_subtract; using default_ops::eval_bit_test; + using default_ops::eval_decrement; typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt; @@ -626,19 +714,41 @@ inline void do_eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponen res.exponent() = a.exponent() - e_diff; eval_subtract(dt, b.bits()); } + else if(a.exponent() == (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent() + 1) + { + if(eval_lsb(b.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1) + { + eval_left_shift(dt, 1); + eval_decrement(dt); + res.exponent() = a.exponent() - 1; + } + else + res.exponent() = a.exponent(); + } else res.exponent() = a.exponent(); } else { dt = b.bits(); - if(b.exponent() <= a.exponent() + (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) + if(b.exponent() <= (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + a.exponent()) { typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type e_diff = a.exponent() - b.exponent(); eval_left_shift(dt, -e_diff); res.exponent() = b.exponent() + e_diff; eval_subtract(dt, a.bits()); } + else if(b.exponent() == (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + a.exponent() + 1) + { + if(eval_lsb(a.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1) + { + eval_left_shift(dt, 1); + eval_decrement(dt); + res.exponent() = b.exponent() - 1; + } + else + res.exponent() = b.exponent(); + } else res.exponent() = b.exponent(); s = !s; @@ -1232,50 +1342,77 @@ inline void eval_convert_to(boost::ulong_long_type *res, const cpp_bin_float<Dig template <class Float, unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> inline typename boost::enable_if_c<boost::is_float<Float>::value>::type eval_convert_to(Float *res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &original_arg) { - // - // Perform rounding first, then afterwards extract the digits: - // typedef cpp_bin_float<std::numeric_limits<Float>::digits, digit_base_2, void, Exponent, MinE, MaxE> conv_type; typedef typename common_type<typename conv_type::exponent_type, int>::type common_exp_type; - conv_type arg(original_arg); - switch(arg.exponent()) + // + // Special cases first: + // + switch(original_arg.exponent()) { - case conv_type::exponent_zero: + case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: *res = 0; - if(arg.sign()) + if(original_arg.sign()) *res = -*res; return; - case conv_type::exponent_nan: + case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: *res = std::numeric_limits<Float>::quiet_NaN(); return; - case conv_type::exponent_infinity: + case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: *res = (std::numeric_limits<Float>::infinity)(); - if(arg.sign()) + if(original_arg.sign()) *res = -*res; return; } - common_exp_type e = arg.exponent(); - static const common_exp_type min_exp_limit = std::numeric_limits<Float>::min_exponent - - (common_exp_type)cpp_bin_float<std::numeric_limits<Float>::digits, digit_base_2, void, Exponent, MinE, MaxE>::bit_count - std::numeric_limits<Float>::digits - 2; - e -= cpp_bin_float<std::numeric_limits<Float>::digits, digit_base_2, void, Exponent, MinE, MaxE>::bit_count - 1; - if(e < min_exp_limit) + // + // Check for super large exponent that must be converted to infinity: + // + if(original_arg.exponent() > std::numeric_limits<Float>::max_exponent) { - *res = 0; + *res = std::numeric_limits<Float>::has_infinity ? std::numeric_limits<Float>::infinity() : (std::numeric_limits<Float>::max)(); + if(original_arg.sign()) + *res = -*res; return; } - if(e > std::numeric_limits<Float>::max_exponent) + // + // Figure out how many digits we will have in our result, + // allowing for a possibly denormalized result: + // + common_exp_type digits_to_round_to = std::numeric_limits<Float>::digits; + if(original_arg.exponent() < std::numeric_limits<Float>::min_exponent - 1) { - *res = std::numeric_limits<Float>::has_infinity ? std::numeric_limits<Float>::infinity() : (std::numeric_limits<Float>::max)(); - return; + common_exp_type diff = original_arg.exponent(); + diff -= std::numeric_limits<Float>::min_exponent - 1; + digits_to_round_to += diff; } - - *res = std::ldexp(static_cast<Float>(*arg.bits().limbs()), static_cast<int>(e)); - for(unsigned i = 1; i < arg.bits().size(); ++i) + if(digits_to_round_to < 0) { + // Result must be zero: + *res = 0; + if(original_arg.sign()) + *res = -*res; + return; + } + // + // Perform rounding first, then afterwards extract the digits: + // + cpp_bin_float<std::numeric_limits<Float>::digits, digit_base_2, Allocator, Exponent, MinE, MaxE> arg; + typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type bits(original_arg.bits()); + arg.exponent() = original_arg.exponent(); + copy_and_round(arg, bits, (int)digits_to_round_to); + common_exp_type e = arg.exponent(); + e -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1; + static const unsigned limbs_needed = std::numeric_limits<Float>::digits / (sizeof(*arg.bits().limbs()) * CHAR_BIT) + + (std::numeric_limits<Float>::digits % (sizeof(*arg.bits().limbs()) * CHAR_BIT) ? 1 : 0); + unsigned first_limb_needed = arg.bits().size() - limbs_needed; + *res = 0; + e += first_limb_needed * sizeof(*arg.bits().limbs()) * CHAR_BIT; + while(first_limb_needed < arg.bits().size()) + { + *res += std::ldexp(static_cast<Float>(arg.bits().limbs()[first_limb_needed]), static_cast<int>(e)); + ++first_limb_needed; e += sizeof(*arg.bits().limbs()) * CHAR_BIT; - *res += std::ldexp(static_cast<Float>(arg.bits().limbs()[i]), static_cast<int>(e)); } - if(arg.sign()) + if(original_arg.sign()) *res = -*res; } @@ -1401,13 +1538,18 @@ inline void eval_sqrt(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE using default_ops::eval_increment; switch(arg.exponent()) { - case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: + errno = EDOM; + // fallthrough... + case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: res = arg; return; case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: if(arg.sign()) + { res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); + errno = EDOM; + } else res = arg; return; @@ -1415,6 +1557,7 @@ inline void eval_sqrt(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE if(arg.sign()) { res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); + errno = EDOM; return; } @@ -1443,8 +1586,10 @@ inline void eval_floor(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, Min using default_ops::eval_increment; switch(arg.exponent()) { - case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: + errno = EDOM; + // fallthrough... + case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: res = arg; return; @@ -1484,9 +1629,11 @@ inline void eval_ceil(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE using default_ops::eval_increment; switch(arg.exponent()) { + case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: + errno = EDOM; + // fallthrough... case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: - case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: res = arg; return; } @@ -1521,6 +1668,12 @@ inline void eval_ceil(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE } template<unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2> +int eval_signbit(const cpp_bin_float<D1, B1, A1, E1, M1, M2>& val) +{ + return val.sign(); +} + +template<unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2> inline std::size_t hash_value(const cpp_bin_float<D1, B1, A1, E1, M1, M2>& val) { std::size_t result = hash_value(val.bits()); @@ -1545,12 +1698,6 @@ struct is_explicitly_convertible<FloatT, backends::cpp_bin_float<D2, B2, A2, E2, #endif template<unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Exponent, Exponent MinE, Exponent MaxE, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates> -inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>& arg) -{ - return arg.backend().sign(); -} - -template<unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Exponent, Exponent MinE, Exponent MaxE, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates> inline boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION( const boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>& a, diff --git a/boost/multiprecision/cpp_bin_float/transcendental.hpp b/boost/multiprecision/cpp_bin_float/transcendental.hpp index 9037dd3982..066bc458e4 100644 --- a/boost/multiprecision/cpp_bin_float/transcendental.hpp +++ b/boost/multiprecision/cpp_bin_float/transcendental.hpp @@ -70,6 +70,7 @@ void eval_exp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> if(type == (int)FP_NAN) { res = arg; + errno = EDOM; return; } else if(type == (int)FP_INFINITE) diff --git a/boost/multiprecision/cpp_dec_float.hpp b/boost/multiprecision/cpp_dec_float.hpp index 9e0c2900ae..f9e6a22189 100644 --- a/boost/multiprecision/cpp_dec_float.hpp +++ b/boost/multiprecision/cpp_dec_float.hpp @@ -1254,9 +1254,15 @@ cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, Expone { // Compute the square root of *this. + if((isinf)() && !isneg()) + { + return *this; + } + if(isneg() || (!(isfinite)())) { *this = nan(); + errno = EDOM; return *this; } @@ -2312,6 +2318,11 @@ void cpp_dec_float<Digits10, ExponentType, Allocator>::from_unsigned_long_long(c fpclass = cpp_dec_float_finite; prec_elem = cpp_dec_float_elem_number; + if(u == 0) + { + return; + } + std::size_t i =static_cast<std::size_t>(0u); boost::ulong_long_type uu = u; @@ -2822,6 +2833,8 @@ inline void eval_floor(cpp_dec_float<Digits10, ExponentType, Allocator>& result, result = x; if(!(x.isfinite)() || x.isint()) { + if((x.isnan)()) + errno = EDOM; return; } @@ -2836,6 +2849,8 @@ inline void eval_ceil(cpp_dec_float<Digits10, ExponentType, Allocator>& result, result = x; if(!(x.isfinite)() || x.isint()) { + if((x.isnan)()) + errno = EDOM; return; } @@ -2847,14 +2862,11 @@ inline void eval_ceil(cpp_dec_float<Digits10, ExponentType, Allocator>& result, template <unsigned Digits10, class ExponentType, class Allocator> inline void eval_trunc(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x) { - if(!(x.isfinite)()) - { - result = boost::math::policies::raise_rounding_error("boost::multiprecision::trunc<%1%>(%1%)", 0, number<cpp_dec_float<Digits10, ExponentType, Allocator> >(x), number<cpp_dec_float<Digits10, ExponentType, Allocator> >(x), boost::math::policies::policy<>()).backend(); - return; - } - else if(x.isint()) + if(x.isint() || !(x.isfinite)()) { result = x; + if((x.isnan)()) + errno = EDOM; return; } result = x.extract_integer_part(); @@ -2863,6 +2875,16 @@ inline void eval_trunc(cpp_dec_float<Digits10, ExponentType, Allocator>& result, template <unsigned Digits10, class ExponentType, class Allocator> inline ExponentType eval_ilogb(const cpp_dec_float<Digits10, ExponentType, Allocator>& val) { + if(val.iszero()) + return (std::numeric_limits<ExponentType>::min)(); + if((val.isinf)()) + return INT_MAX; + if((val.isnan)()) +#ifdef FP_ILOGBNAN + return FP_ILOGBNAN; +#else + return INT_MAX; +#endif // Set result, to the exponent of val: return val.order(); } @@ -2897,15 +2919,16 @@ template <unsigned Digits10, class ExponentType, class Allocator> inline void eval_frexp(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x, ExponentType* e) { result = x; - if(result.isneg()) - result.negate(); - if(result.iszero()) + if(result.iszero() || (result.isinf)() || (result.isnan)()) { *e = 0; return; } + if(result.isneg()) + result.negate(); + ExponentType t = result.order(); BOOST_MP_USING_ABS if(abs(t) < ((std::numeric_limits<ExponentType>::max)() / 1000)) diff --git a/boost/multiprecision/cpp_int/add.hpp b/boost/multiprecision/cpp_int/add.hpp index 22b8c186dc..2fca9b9597 100644 --- a/boost/multiprecision/cpp_int/add.hpp +++ b/boost/multiprecision/cpp_int/add.hpp @@ -84,7 +84,7 @@ inline void add_unsigned(CppInt1& result, const CppInt2& a, const CppInt3& b) BO { // We overflowed, need to add one more limb: result.resize(x + 1, x + 1); - if(CppInt1::variable || (result.size() > x)) + if(result.size() > x) result.limbs()[x] = static_cast<limb_type>(carry); } result.normalize(); @@ -126,7 +126,7 @@ inline void add_unsigned(CppInt1& result, const CppInt2& a, const limb_type& o) // We overflowed, need to add one more limb: unsigned x = result.size(); result.resize(x + 1, x + 1); - if(CppInt1::variable || (result.size() > x)) + if(result.size() > x) result.limbs()[x] = static_cast<limb_type>(carry); } result.normalize(); diff --git a/boost/multiprecision/cpp_int/bitwise.hpp b/boost/multiprecision/cpp_int/bitwise.hpp index 748d088056..168544c9bb 100644 --- a/boost/multiprecision/cpp_int/bitwise.hpp +++ b/boost/multiprecision/cpp_int/bitwise.hpp @@ -326,7 +326,7 @@ inline void left_shift_byte(Int& result, double_limb_type s) if(rs != ors) pr[rs - 1] = 0u; std::size_t bytes = static_cast<std::size_t>(s / CHAR_BIT); - std::size_t len = std::min(ors * sizeof(limb_type), rs * sizeof(limb_type) - bytes); + std::size_t len = (std::min)(ors * sizeof(limb_type), rs * sizeof(limb_type) - bytes); if(bytes >= rs * sizeof(limb_type)) result = static_cast<limb_type>(0u); else @@ -411,14 +411,14 @@ inline void left_shift_generic(Int& result, double_limb_type s) ++i; } } - for(; ors > 1 + i; ++i) + for(; rs - i >= 2 + offset; ++i) { - pr[rs - 1 - i] = pr[ors - 1 - i] << shift; - pr[rs - 1 - i] |= pr[ors - 2 - i] >> (Int::limb_bits - shift); + pr[rs - 1 - i] = pr[rs - 1 - i - offset] << shift; + pr[rs - 1 - i] |= pr[rs - 2 - i - offset] >> (Int::limb_bits - shift); } - if(ors >= 1 + i) + if(rs - i >= 1 + offset) { - pr[rs - 1 - i] = pr[ors - 1 - i] << shift; + pr[rs - 1 - i] = pr[rs - 1 - i - offset] << shift; ++i; } for(; i < rs; ++i) diff --git a/boost/multiprecision/cpp_int/divide.hpp b/boost/multiprecision/cpp_int/divide.hpp index 2f83d1ba71..81583ec384 100644 --- a/boost/multiprecision/cpp_int/divide.hpp +++ b/boost/multiprecision/cpp_int/divide.hpp @@ -220,7 +220,7 @@ void divide_unsigned_helper( // double_limb_type carry = 0; t.resize(y.size() + shift + 1, y.size() + shift); - bool truncated_t = !CppInt1::variable && (t.size() != y.size() + shift + 1); + bool truncated_t = (t.size() != y.size() + shift + 1); typename CppInt1::limb_pointer pt = t.limbs(); for(unsigned i = 0; i < shift; ++i) pt[i] = 0; diff --git a/boost/multiprecision/cpp_int/import_export.hpp b/boost/multiprecision/cpp_int/import_export.hpp index 0fe3249ddc..7bd843fb32 100644 --- a/boost/multiprecision/cpp_int/import_export.hpp +++ b/boost/multiprecision/cpp_int/import_export.hpp @@ -134,7 +134,7 @@ namespace boost { cpp_int_backend<MinBits, MaxBits, SignType, Checked, Allocator>& result = val.backend(); result.resize(static_cast<unsigned>(limb_len), static_cast<unsigned>(limb_len)); // checked types may throw here if they're not large enough to hold the data! result.limbs()[result.size() - 1] = 0u; - std::memcpy(result.limbs(), i, std::min(byte_len, result.size() * sizeof(limb_type))); + std::memcpy(result.limbs(), i, (std::min)(byte_len, result.size() * sizeof(limb_type))); result.normalize(); // In case data has leading zeros. return val; } @@ -150,7 +150,7 @@ namespace boost { ++limb_len; result.limbs()[0] = 0u; result.resize(static_cast<unsigned>(limb_len), static_cast<unsigned>(limb_len)); // checked types may throw here if they're not large enough to hold the data! - std::memcpy(result.limbs(), i, std::min(byte_len, result.size() * sizeof(result.limbs()[0]))); + std::memcpy(result.limbs(), i, (std::min)(byte_len, result.size() * sizeof(result.limbs()[0]))); result.normalize(); // In case data has leading zeros. return val; } @@ -199,8 +199,8 @@ namespace boost { template <class Backend> inline boost::uintmax_t extract_bits(const Backend& val, unsigned location, unsigned count, const mpl::true_&) { - boost::uintmax_t result = *val.limbs(); - boost::uintmax_t mask = count == std::numeric_limits<boost::uintmax_t>::digits ? ~static_cast<boost::uintmax_t>(0) : (static_cast<boost::uintmax_t>(1u) << count) - 1; + typename Backend::local_limb_type result = *val.limbs(); + typename Backend::local_limb_type mask = count >= std::numeric_limits<typename Backend::local_limb_type>::digits ? ~static_cast<typename Backend::local_limb_type>(0) : (static_cast<typename Backend::local_limb_type>(1u) << count) - 1; return (result >> location) & mask; } diff --git a/boost/multiprecision/cpp_int/limits.hpp b/boost/multiprecision/cpp_int/limits.hpp index 518bb9f4e0..b19e1ebbe6 100644 --- a/boost/multiprecision/cpp_int/limits.hpp +++ b/boost/multiprecision/cpp_int/limits.hpp @@ -23,7 +23,7 @@ inline boost::multiprecision::number<boost::multiprecision::cpp_int_backend<MinB { // Bounded and signed. typedef boost::multiprecision::number<boost::multiprecision::cpp_int_backend<MinBits, MaxBits, SignType, Checked, Allocator>, ExpressionTemplates> result_type; - typedef boost::multiprecision::number<boost::multiprecision::cpp_int_backend<MinBits, MaxBits, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked, Allocator>, ExpressionTemplates> ui_type; + typedef boost::multiprecision::number<boost::multiprecision::cpp_int_backend<MaxBits, MaxBits, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked>, ExpressionTemplates> ui_type; static const result_type val = -result_type(~ui_type(0)); return val; } @@ -62,7 +62,7 @@ inline boost::multiprecision::number<boost::multiprecision::cpp_int_backend<MinB { // Bounded and signed. typedef boost::multiprecision::number<boost::multiprecision::cpp_int_backend<MinBits, MaxBits, SignType, Checked, Allocator>, ExpressionTemplates> result_type; - typedef boost::multiprecision::number<boost::multiprecision::cpp_int_backend<MinBits, MaxBits, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked, Allocator>, ExpressionTemplates> ui_type; + typedef boost::multiprecision::number<boost::multiprecision::cpp_int_backend<MaxBits, MaxBits, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked>, ExpressionTemplates> ui_type; static const result_type val = ~ui_type(0); return val; } diff --git a/boost/multiprecision/cpp_int/multiply.hpp b/boost/multiprecision/cpp_int/multiply.hpp index 88264c47be..afe1db7a9f 100644 --- a/boost/multiprecision/cpp_int/multiply.hpp +++ b/boost/multiprecision/cpp_int/multiply.hpp @@ -48,7 +48,7 @@ inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBit { unsigned i = result.size(); result.resize(i + 1, i + 1); - if(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::variable || (result.size() > i)) + if(result.size() > i) result.limbs()[i] = static_cast<limb_type>(carry); } result.sign(a.sign()); @@ -157,7 +157,7 @@ inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBit BOOST_ASSERT(carry <= (cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::max_limb_value)); } resize_for_carry(result, as + bs); // May throw if checking is enabled - if(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::variable || (i + bs < result.size())) + if(i + bs < result.size()) pr[i + bs] = static_cast<limb_type>(carry); carry = 0; } diff --git a/boost/multiprecision/debug_adaptor.hpp b/boost/multiprecision/debug_adaptor.hpp index 1b69ba4aa6..1e395035f7 100644 --- a/boost/multiprecision/debug_adaptor.hpp +++ b/boost/multiprecision/debug_adaptor.hpp @@ -282,6 +282,7 @@ inline void eval_ldexp(debug_adaptor<Backend>& result, const debug_adaptor<Backe template <class Backend, class Exp> inline void eval_scalbn(debug_adaptor<Backend>& result, const debug_adaptor<Backend>& arg, Exp exp) { + using default_ops::eval_scalbn; eval_scalbn(result.value(), arg.value(), exp); result.update_view(); } @@ -289,6 +290,7 @@ inline void eval_scalbn(debug_adaptor<Backend>& result, const debug_adaptor<Back template <class Backend> inline typename Backend::exponent_type eval_ilogb(const debug_adaptor<Backend>& arg) { + using default_ops::eval_ilogb; return eval_ilogb(arg.value()); } @@ -456,12 +458,17 @@ NON_MEMBER_OP3(pow, "pow"); NON_MEMBER_OP3(atan2, "atan2"); template <class Backend> +int eval_signbit(const debug_adaptor<Backend>& val) +{ + return eval_signbit(val.value()); +} + +template <class Backend> std::size_t hash_value(const debug_adaptor<Backend>& val) { return hash_value(val.value()); } - } // namespace backends using backends::debug_adaptor; diff --git a/boost/multiprecision/detail/bitscan.hpp b/boost/multiprecision/detail/bitscan.hpp index a21fd58651..21e8b8a3a9 100644 --- a/boost/multiprecision/detail/bitscan.hpp +++ b/boost/multiprecision/detail/bitscan.hpp @@ -8,6 +8,8 @@ #ifndef BOOST_MP_DETAIL_BITSCAN_HPP #define BOOST_MP_DETAIL_BITSCAN_HPP +#include <boost/detail/endian.hpp> + #if (defined(BOOST_MSVC) || (defined(__clang__) && defined(__c2__)) || (defined(BOOST_INTEL) && defined(_MSC_VER))) && (defined(_M_IX86) || defined(_M_X64)) #include <intrin.h> #endif @@ -140,6 +142,36 @@ BOOST_FORCEINLINE unsigned find_msb(boost::ulong_long_type mask, mpl::int_<3> co { return sizeof(boost::ulong_long_type) * CHAR_BIT - 1 - __builtin_clzll(mask); } +#ifdef BOOST_HAS_INT128 +BOOST_FORCEINLINE unsigned find_msb(unsigned __int128 mask, mpl::int_<0> const&) +{ + union { unsigned __int128 v; boost::uint64_t sv[2]; } val; + val.v = mask; +#ifdef BOOST_LITTLE_ENDIAN + if(val.sv[1]) + return find_msb(val.sv[1], mpl::int_<3>()) + 64; + return find_msb(val.sv[0], mpl::int_<3>()); +#else + if(val.sv[0]) + return find_msb(val.sv[0], mpl::int_<3>()) + 64; + return find_msb(val.sv[1], mpl::int_<3>()); +#endif +} +BOOST_FORCEINLINE unsigned find_lsb(unsigned __int128 mask, mpl::int_<0> const&) +{ + union { unsigned __int128 v; boost::uint64_t sv[2]; } val; + val.v = mask; +#ifdef BOOST_LITTLE_ENDIAN + if(val.sv[0] == 0) + return find_lsb(val.sv[1], mpl::int_<3>()) + 64; + return find_lsb(val.sv[0], mpl::int_<3>()); +#else + if(val.sv[1] == 0) + return find_lsb(val.sv[0], mpl::int_<3>()) + 64; + return find_lsb(val.sv[1], mpl::int_<3>()); +#endif +} +#endif template <class Unsigned> BOOST_FORCEINLINE unsigned find_lsb(Unsigned mask) diff --git a/boost/multiprecision/detail/default_ops.hpp b/boost/multiprecision/detail/default_ops.hpp index fe8d5acd16..334db06a64 100644 --- a/boost/multiprecision/detail/default_ops.hpp +++ b/boost/multiprecision/detail/default_ops.hpp @@ -983,6 +983,25 @@ inline void eval_fmod(T& result, const T& a, const T& b) result = temp; return; } + switch(eval_fpclassify(a)) + { + case FP_ZERO: + result = a; + return; + case FP_INFINITE: + case FP_NAN: + result = std::numeric_limits<number<T> >::quiet_NaN().backend(); + errno = EDOM; + return; + } + switch(eval_fpclassify(b)) + { + case FP_ZERO: + case FP_NAN: + result = std::numeric_limits<number<T> >::quiet_NaN().backend(); + errno = EDOM; + return; + } T n; eval_divide(result, a, b); if(eval_get_sign(result) < 0) @@ -1162,10 +1181,14 @@ template <class T> inline void eval_trunc(T& result, const T& a) { BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The trunc function is only valid for floating point types."); - int c = eval_fpclassify(a); - if(c == (int)FP_NAN || c == (int)FP_INFINITE) + switch(eval_fpclassify(a)) { - result = boost::math::policies::raise_rounding_error("boost::multiprecision::trunc<%1%>(%1%)", 0, number<T>(a), number<T>(a), boost::math::policies::policy<>()).backend(); + case FP_NAN: + errno = EDOM; + // fallthrough... + case FP_ZERO: + case FP_INFINITE: + result = a; return; } if(eval_get_sign(a) < 0) @@ -1212,12 +1235,17 @@ inline void eval_round(T& result, const T& a) BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The round function is only valid for floating point types."); typedef typename boost::multiprecision::detail::canonical<float, T>::type fp_type; int c = eval_fpclassify(a); - if((c == (int)FP_NAN) || (c == (int)FP_INFINITE)) + if(c == (int)FP_NAN) { - result = boost::math::policies::raise_rounding_error("boost::multiprecision::round<%1%>(%1%)", 0, number<T>(a), number<T>(a), boost::math::policies::policy<>()).backend(); + result = a; + errno = EDOM; return; } - if(eval_get_sign(a) < 0) + if((c == FP_ZERO) || (c == (int)FP_INFINITE)) + { + result = a; + } + else if(eval_get_sign(a) < 0) { eval_subtract(result, a, fp_type(0.5f)); eval_ceil(result, result); @@ -1382,9 +1410,10 @@ void eval_integer_sqrt(B& s, B& r, const B& x) return; } int g = eval_msb(x); - if(g == 0) + if(g <= 1) { - r = ui_type(1); + s = ui_type(1); + eval_subtract(r, x, s); return; } @@ -1408,6 +1437,7 @@ void eval_integer_sqrt(B& s, B& r, const B& x) eval_bit_set(t, 2 * g); if(t.compare(r) <= 0) { + BOOST_ASSERT(g >= 0); eval_bit_set(s, g); eval_subtract(r, t); if(eval_get_sign(r) == 0) @@ -1448,7 +1478,11 @@ inline typename B::exponent_type eval_ilogb(const B& val) switch(eval_fpclassify(val)) { case FP_NAN: - return (std::numeric_limits<typename B::exponent_type>::min)(); +#ifdef FP_ILOGBNAN + return FP_ILOGBNAN; +#else + return (std::numeric_limits<typename B::exponent_type>::max)(); +#endif case FP_INFINITE: return (std::numeric_limits<typename B::exponent_type>::max)(); case FP_ZERO: @@ -1458,9 +1492,30 @@ inline typename B::exponent_type eval_ilogb(const B& val) eval_frexp(result, val, &e); return e - 1; } + +template <class T> +int eval_signbit(const T& val); + template <class B> inline void eval_logb(B& result, const B& val) { + switch(eval_fpclassify(val)) + { + case FP_NAN: + result = val; + errno = EDOM; + return; + case FP_ZERO: + result = std::numeric_limits<number<B> >::infinity().backend(); + result.negate(); + errno = ERANGE; + return; + case FP_INFINITE: + result = val; + if(eval_signbit(val)) + result.negate(); + return; + } typedef typename boost::mpl::if_c<boost::is_same<boost::intmax_t, long>::value, boost::long_long_type, boost::intmax_t>::type max_t; result = static_cast<max_t>(eval_ilogb(val)); } @@ -1598,6 +1653,12 @@ inline void eval_rint(R& result, const T& a) eval_nearbyint(result, a); } +template <class T> +inline int eval_signbit(const T& val) +{ + return eval_get_sign(val) < 0 ? 1 : 0; +} + // // These functions are implemented in separate files, but expanded inline here, // DO NOT CHANGE THE ORDER OF THESE INCLUDES: @@ -1687,7 +1748,8 @@ inline int sign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::e template <class Backend, multiprecision::expression_template_option ExpressionTemplates> inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg) { - return arg.sign() < 0; + using default_ops::eval_signbit; + return eval_signbit(arg.backend()); } template <class tag, class A1, class A2, class A3, class A4> inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg) @@ -1837,7 +1899,14 @@ namespace multiprecision{ template <class Backend, multiprecision::expression_template_option ExpressionTemplates> inline multiprecision::number<Backend, ExpressionTemplates> lgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg) { - return boost::math::lgamma(arg, c99_error_policy()); + multiprecision::number<Backend, ExpressionTemplates> result; + result = boost::math::lgamma(arg, c99_error_policy()); + if((boost::multiprecision::isnan)(result) && !(boost::multiprecision::isnan)(arg)) + { + result = std::numeric_limits<multiprecision::number<Backend, ExpressionTemplates> >::infinity(); + errno = ERANGE; + } + return result; } template <class tag, class A1, class A2, class A3, class A4> inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type lgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg) @@ -1848,6 +1917,11 @@ namespace multiprecision{ template <class Backend, multiprecision::expression_template_option ExpressionTemplates> inline multiprecision::number<Backend, ExpressionTemplates> tgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg) { + if((arg == 0) && std::numeric_limits<multiprecision::number<Backend, ExpressionTemplates> >::has_infinity) + { + errno = ERANGE; + return 1 / arg; + } return boost::math::tgamma(arg, c99_error_policy()); } template <class tag, class A1, class A2, class A3, class A4> diff --git a/boost/multiprecision/detail/functions/pow.hpp b/boost/multiprecision/detail/functions/pow.hpp index b244a18c8b..460658e034 100644 --- a/boost/multiprecision/detail/functions/pow.hpp +++ b/boost/multiprecision/detail/functions/pow.hpp @@ -208,11 +208,11 @@ void eval_exp(T& result, const T& x) if(type == (int)FP_NAN) { result = x; + errno = EDOM; return; } else if(type == (int)FP_INFINITE) { - result = x; if(isneg) result = ui_type(0u); else @@ -277,6 +277,17 @@ void eval_exp(T& result, const T& x) 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 @@ -326,6 +337,29 @@ void eval_log(T& result, const T& arg) 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; @@ -427,19 +461,21 @@ inline void eval_pow(T& result, const T& x, const T& a) return; } - if(a.compare(si_type(1)) == 0) + 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_INFINITE: - result = x; - return; case FP_ZERO: switch(eval_fpclassify(a)) { @@ -449,13 +485,58 @@ inline void eval_pow(T& result, const T& x, const T& a) 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: - result = x; + 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: ; } @@ -478,6 +559,16 @@ inline void eval_pow(T& result, const T& x, const T& a) } 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 @@ -521,8 +612,33 @@ inline void eval_pow(T& result, const T& x, const T& a) // conversion failed, just fall through, value is not an integer. } #endif - if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN) + 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.")); @@ -534,7 +650,7 @@ inline void eval_pow(T& result, const T& x, const T& a) eval_subtract(da, a, an); - if((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0)) + 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) { @@ -618,14 +734,20 @@ void eval_exp2(T& result, const T& arg) // 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; - 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; + 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); @@ -669,6 +791,8 @@ namespace detail{ switch(eval_fpclassify(x)) { case FP_NAN: + errno = EDOM; + // fallthrough... case FP_INFINITE: if(p_sinh) *p_sinh = x; @@ -695,6 +819,8 @@ namespace detail{ 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) { @@ -742,6 +868,14 @@ 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); } diff --git a/boost/multiprecision/detail/functions/trig.hpp b/boost/multiprecision/detail/functions/trig.hpp index c84a639418..319b708b0c 100644 --- a/boost/multiprecision/detail/functions/trig.hpp +++ b/boost/multiprecision/detail/functions/trig.hpp @@ -91,12 +91,15 @@ void eval_sin(T& result, const T& x) case FP_INFINITE: case FP_NAN: 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 is undefined or complex and there is no NaN for this number type.")); return; case FP_ZERO: - result = ui_type(0); + result = x; return; default: ; } @@ -238,7 +241,10 @@ void eval_cos(T& result, const T& x) case FP_INFINITE: case FP_NAN: 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 is undefined or complex and there is no NaN for this number type.")); return; @@ -422,12 +428,15 @@ void eval_asin(T& result, const T& x) case FP_NAN: case FP_INFINITE: 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 is undefined or complex and there is no NaN for this number type.")); return; case FP_ZERO: - result = ui_type(0); + result = x; return; default: ; } @@ -442,7 +451,10 @@ void eval_asin(T& result, const T& x) if(c > 0) { 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 is undefined or complex and there is no NaN for this number type.")); return; @@ -515,10 +527,8 @@ void eval_asin(T& result, const T& x) eval_divide(sine, cosine); eval_subtract(result, sine); current_precision = eval_ilogb(sine); -#ifdef FP_ILOGB0 - if(current_precision == FP_ILOGB0) + if(current_precision <= (std::numeric_limits<typename T::exponent_type>::min)() + 1) break; -#endif } if(b_neg) result.negate(); @@ -535,7 +545,10 @@ inline void eval_acos(T& result, const T& x) case FP_NAN: case FP_INFINITE: 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 is undefined or complex and there is no NaN for this number type.")); return; @@ -551,7 +564,10 @@ inline void eval_acos(T& result, const T& x) if(c > 0) { 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 is undefined or complex and there is no NaN for this number type.")); return; @@ -584,9 +600,10 @@ void eval_atan(T& result, const T& x) { case FP_NAN: result = x; + errno = EDOM; return; case FP_ZERO: - result = ui_type(0); + result = x; return; case FP_INFINITE: if(eval_get_sign(x) < 0) @@ -662,10 +679,8 @@ void eval_atan(T& result, const T& x) eval_multiply(s, t, c); eval_add(result, s); current_precision = eval_ilogb(s); -#ifdef FP_ILOGB0 - if(current_precision == FP_ILOGB0) + if(current_precision <= (std::numeric_limits<typename T::exponent_type>::min)() + 1) break; -#endif } if(b_neg) result.negate(); @@ -694,24 +709,41 @@ void eval_atan2(T& result, const T& y, const T& x) { case FP_NAN: result = y; + errno = EDOM; return; case FP_ZERO: { - int c = eval_get_sign(x); - if(c < 0) + if(eval_signbit(x)) + { result = get_constant_pi<T>(); - else if(c >= 0) - result = ui_type(0); // Note we allow atan2(0,0) to be zero, even though it's mathematically undefined + if(eval_signbit(y)) + result.negate(); + } + else + { + result = y; // Note we allow atan2(0,0) to be +-zero, even though it's mathematically undefined + } return; } case FP_INFINITE: { if(eval_fpclassify(x) == FP_INFINITE) { - if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN) - result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend(); + if(eval_signbit(x)) + { + // 3Pi/4 + eval_ldexp(result, get_constant_pi<T>(), -2); + eval_subtract(result, get_constant_pi<T>()); + if(eval_get_sign(y) >= 0) + result.negate(); + } else - BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type.")); + { + // Pi/4 + eval_ldexp(result, get_constant_pi<T>(), -2); + if(eval_get_sign(y) < 0) + result.negate(); + } } else { @@ -727,6 +759,7 @@ void eval_atan2(T& result, const T& y, const T& x) { case FP_NAN: result = x; + errno = EDOM; return; case FP_ZERO: { diff --git a/boost/multiprecision/float128.hpp b/boost/multiprecision/float128.hpp index 702d488cd6..8de3af6519 100644 --- a/boost/multiprecision/float128.hpp +++ b/boost/multiprecision/float128.hpp @@ -399,15 +399,6 @@ inline void eval_fabs(float128_backend& result, const float128_backend& arg) inline void eval_trunc(float128_backend& result, const float128_backend& arg) { - if(isnanq(arg.value()) || isinfq(arg.value())) - { - result = boost::math::policies::raise_rounding_error( - "boost::multiprecision::trunc<%1%>(%1%)", 0, - number<float128_backend, et_off>(arg), - number<float128_backend, et_off>(arg), - boost::math::policies::policy<>()).backend(); - return; - } result.value() = truncq(arg.value()); } /* @@ -494,6 +485,11 @@ inline void eval_multiply_add(float128_backend& result, const float128_backend& result.value() = fmaq(a.value(), b.value(), c.value()); } +inline int eval_signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const float128_backend& arg) +{ + return ::signbitq(arg.value()); +} + inline std::size_t hash_value(const float128_backend& val) { return boost::hash_value(static_cast<double>(val.value())); @@ -553,12 +549,6 @@ inline std::size_t hash_value(const float128_backend& val) } template <multiprecision::expression_template_option ExpressionTemplates> - inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates>& arg) - { - return ::signbitq(arg.backend().value()); - } - - template <multiprecision::expression_template_option ExpressionTemplates> inline boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates>& a, const boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates>& b) { return ::copysignq(a.backend().value(), b.backend().value()); diff --git a/boost/multiprecision/logged_adaptor.hpp b/boost/multiprecision/logged_adaptor.hpp index 2f916dd101..7cb33ee57a 100644 --- a/boost/multiprecision/logged_adaptor.hpp +++ b/boost/multiprecision/logged_adaptor.hpp @@ -507,6 +507,13 @@ NON_MEMBER_OP3(pow, "pow"); NON_MEMBER_OP3(atan2, "atan2"); template <class Backend> +int eval_signbit(const logged_adaptor<Backend>& val) +{ + using default_ops::eval_signbit; + return eval_signbit(val.value()); +} + +template <class Backend> std::size_t hash_value(const logged_adaptor<Backend>& val) { return hash_value(val.value()); diff --git a/boost/multiprecision/mpfi.hpp b/boost/multiprecision/mpfi.hpp index b863e31352..a06e4a39d4 100644 --- a/boost/multiprecision/mpfi.hpp +++ b/boost/multiprecision/mpfi.hpp @@ -1042,7 +1042,7 @@ inline std::size_t hash_value(const mpfi_float_backend<Digits10>& val) std::size_t len = val.left_data()[0]._mpfr_prec / mp_bits_per_limb; if(val.left_data()[0]._mpfr_prec % mp_bits_per_limb) ++len; - for(int i = 0; i < len; ++i) + for(std::size_t i = 0; i < len; ++i) boost::hash_combine(result, val.left_data()[0]._mpfr_d[i]); boost::hash_combine(result, val.left_data()[0]._mpfr_exp); boost::hash_combine(result, val.left_data()[0]._mpfr_sign); @@ -1050,7 +1050,7 @@ inline std::size_t hash_value(const mpfi_float_backend<Digits10>& val) len = val.right_data()[0]._mpfr_prec / mp_bits_per_limb; if(val.right_data()[0]._mpfr_prec % mp_bits_per_limb) ++len; - for(int i = 0; i < len; ++i) + for(std::size_t i = 0; i < len; ++i) boost::hash_combine(result, val.right_data()[0]._mpfr_d[i]); boost::hash_combine(result, val.right_data()[0]._mpfr_exp); boost::hash_combine(result, val.right_data()[0]._mpfr_sign); @@ -1420,7 +1420,7 @@ struct constant_pi<boost::multiprecision::number<boost::multiprecision::mpfi_flo { typedef boost::multiprecision::number<boost::multiprecision::mpfi_float_backend<Digits10>, ExpressionTemplates> result_type; template<int N> - static inline result_type const& get(const mpl::int_<N>&) + static inline const result_type& get(const mpl::int_<N>&) { mpfi_initializer<result_type>::force_instantiate(); static result_type result; @@ -1444,7 +1444,7 @@ struct constant_ln_two<boost::multiprecision::number<boost::multiprecision::mpfi { typedef boost::multiprecision::number<boost::multiprecision::mpfi_float_backend<Digits10>, ExpressionTemplates> result_type; template<int N> - static inline result_type get(const mpl::int_<N>&) + static inline const result_type& get(const mpl::int_<N>&) { mpfi_initializer<result_type>::force_instantiate(); static result_type result; diff --git a/boost/multiprecision/mpfr.hpp b/boost/multiprecision/mpfr.hpp index 440517ff86..2759a8a8de 100644 --- a/boost/multiprecision/mpfr.hpp +++ b/boost/multiprecision/mpfr.hpp @@ -1286,11 +1286,6 @@ inline void eval_floor(mpfr_float_backend<Digits10, AllocateType>& result, const template <unsigned Digits10, mpfr_allocation_type AllocateType> inline void eval_trunc(mpfr_float_backend<Digits10, AllocateType>& result, const mpfr_float_backend<Digits10, AllocateType>& val) { - if(0 == mpfr_number_p(val.data())) - { - result = boost::math::policies::raise_rounding_error("boost::multiprecision::trunc<%1%>(%1%)", 0, number<mpfr_float_backend<Digits10, AllocateType> >(val), number<mpfr_float_backend<Digits10, AllocateType> >(val), boost::math::policies::policy<>()).backend(); - return; - } mpfr_trunc(result.data(), val.data()); } template <unsigned Digits10, mpfr_allocation_type AllocateType> @@ -1327,7 +1322,12 @@ inline int eval_fpclassify(const mpfr_float_backend<Digits10, AllocateType>& val template <unsigned Digits10, mpfr_allocation_type AllocateType> inline void eval_pow(mpfr_float_backend<Digits10, AllocateType>& result, const mpfr_float_backend<Digits10, AllocateType>& b, const mpfr_float_backend<Digits10, AllocateType>& e) { - mpfr_pow(result.data(), b.data(), e.data(), GMP_RNDN); + if(mpfr_zero_p(b.data()) && mpfr_integer_p(e.data()) && (mpfr_signbit(e.data()) == 0) && mpfr_fits_ulong_p(e.data(), GMP_RNDN) && (mpfr_get_ui(e.data(), GMP_RNDN) & 1)) + { + mpfr_set(result.data(), b.data(), GMP_RNDN); + } + else + mpfr_pow(result.data(), b.data(), e.data(), GMP_RNDN); } #ifdef BOOST_MSVC @@ -1505,13 +1505,19 @@ inline void eval_multiply_subtract(mpfr_float_backend<Digits10, AllocateType>& r } template <unsigned Digits10, mpfr_allocation_type AllocateType> +inline int eval_signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const mpfr_float_backend<Digits10, AllocateType>& arg) +{ + return (arg.data()[0]._mpfr_sign < 0) ? 1 : 0; +} + +template <unsigned Digits10, mpfr_allocation_type AllocateType> inline std::size_t hash_value(const mpfr_float_backend<Digits10, AllocateType>& val) { std::size_t result = 0; std::size_t len = val.data()[0]._mpfr_prec / mp_bits_per_limb; if(val.data()[0]._mpfr_prec % mp_bits_per_limb) ++len; - for(int i = 0; i < len; ++i) + for(std::size_t i = 0; i < len; ++i) boost::hash_combine(result, val.data()[0]._mpfr_d[i]); boost::hash_combine(result, val.data()[0]._mpfr_exp); boost::hash_combine(result, val.data()[0]._mpfr_sign); @@ -1546,24 +1552,12 @@ typedef number<mpfr_float_backend<50, allocate_stack> > static_mpfr_float_50; typedef number<mpfr_float_backend<100, allocate_stack> > static_mpfr_float_100; template<unsigned Digits10, boost::multiprecision::mpfr_allocation_type AllocateType, boost::multiprecision::expression_template_option ExpressionTemplates> -inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType>, ExpressionTemplates>& arg) -{ - return (arg.backend().data()[0]._mpfr_sign < 0) ? 1 : 0; -} - -template<unsigned Digits10, boost::multiprecision::mpfr_allocation_type AllocateType, boost::multiprecision::expression_template_option ExpressionTemplates> inline boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType>, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType>, ExpressionTemplates>& a, const boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType>, ExpressionTemplates>& b) { return (boost::multiprecision::signbit)(a) != (boost::multiprecision::signbit)(b) ? boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType>, ExpressionTemplates>(-a) : a; } template<unsigned Digits10, boost::multiprecision::mpfr_allocation_type AllocateType, boost::multiprecision::expression_template_option ExpressionTemplates> -inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType> >, ExpressionTemplates>& arg) -{ - return (arg.backend().value().data()[0]._mpfr_sign < 0) ? 1 : 0; -} - -template<unsigned Digits10, boost::multiprecision::mpfr_allocation_type AllocateType, boost::multiprecision::expression_template_option ExpressionTemplates> inline boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType> >, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType> >, ExpressionTemplates>& a, const boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType> >, ExpressionTemplates>& b) { return (boost::multiprecision::signbit)(a) != (boost::multiprecision::signbit)(b) ? boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType> >, ExpressionTemplates>(-a) : a; diff --git a/boost/multiprecision/number.hpp b/boost/multiprecision/number.hpp index aaf6532e13..c3819dc561 100644 --- a/boost/multiprecision/number.hpp +++ b/boost/multiprecision/number.hpp @@ -1732,7 +1732,7 @@ inline std::string read_string_while(std::istream& is, std::string const& permit for(;; c = is.rdbuf()->snextc()) if(std::istream::traits_type::eq_int_type(std::istream::traits_type::eof(), c)) - { // end of file: + { // end of file: state |= std::ios_base::eofbit; break; } @@ -1743,7 +1743,7 @@ inline std::string read_string_while(std::istream& is, std::string const& permit break; } else - { + { result.append(1, std::istream::traits_type::to_char_type(c)); } } @@ -1856,78 +1856,6 @@ inline std::istream& operator >> (std::istream& is, rational<multiprecision::num return is; } -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator == (const rational<multiprecision::number<T, ExpressionTemplates> >& a, const Arithmetic& b) -{ - return a == multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator == (const Arithmetic& b, const rational<multiprecision::number<T, ExpressionTemplates> >& a) -{ - return a == multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator != (const rational<multiprecision::number<T, ExpressionTemplates> >& a, const Arithmetic& b) -{ - return a != multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator != (const Arithmetic& b, const rational<multiprecision::number<T, ExpressionTemplates> >& a) -{ - return a != multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator < (const rational<multiprecision::number<T, ExpressionTemplates> >& a, const Arithmetic& b) -{ - return a < multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator < (const Arithmetic& b, const rational<multiprecision::number<T, ExpressionTemplates> >& a) -{ - return a > multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator <= (const rational<multiprecision::number<T, ExpressionTemplates> >& a, const Arithmetic& b) -{ - return a <= multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator <= (const Arithmetic& b, const rational<multiprecision::number<T, ExpressionTemplates> >& a) -{ - return a >= multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator > (const rational<multiprecision::number<T, ExpressionTemplates> >& a, const Arithmetic& b) -{ - return a > multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator > (const Arithmetic& b, const rational<multiprecision::number<T, ExpressionTemplates> >& a) -{ - return a < multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator >= (const rational<multiprecision::number<T, ExpressionTemplates> >& a, const Arithmetic& b) -{ - return a >= multiprecision::number<T, ExpressionTemplates>(b); -} - -template <class T, multiprecision::expression_template_option ExpressionTemplates, class Arithmetic> -typename boost::enable_if<boost::is_arithmetic<Arithmetic>, bool>::type operator >= (const Arithmetic& b, const rational<multiprecision::number<T, ExpressionTemplates> >& a) -{ - return a <= multiprecision::number<T, ExpressionTemplates>(b); -} - template <class T, multiprecision::expression_template_option ExpressionTemplates> inline multiprecision::number<T, ExpressionTemplates> numerator(const rational<multiprecision::number<T, ExpressionTemplates> >& a) { |