From d69424adad4e9b6d28b11708e07e2eb2953a675c Mon Sep 17 00:00:00 2001 From: John Maddock Date: Thu, 26 Jan 2012 10:11:10 +0000 Subject: [PATCH] Fix bug in fixed_int::convert_to with negative numbers. Fix bug in fixed_int shift operator when shifting by 0. Add preliminary gcd/lcm support for integer types. Add static asserts to floating-point only functions. [SVN r76706] --- .../multiprecision/detail/default_ops.hpp | 91 +++++++++++ .../multiprecision/detail/functions/pow.hpp | 8 + .../multiprecision/detail/functions/trig.hpp | 8 + include/boost/multiprecision/fixed_int.hpp | 146 +++++++++++++++++- include/boost/multiprecision/gmp.hpp | 25 +++ include/boost/multiprecision/tommath.hpp | 8 + performance/performance_test.cpp | 15 +- test/test_arithmetic.cpp | 22 +++ test/test_fixed_int.cpp | 3 + 9 files changed, 324 insertions(+), 2 deletions(-) diff --git a/include/boost/multiprecision/detail/default_ops.hpp b/include/boost/multiprecision/detail/default_ops.hpp index 25010500..68a85344 100644 --- a/include/boost/multiprecision/detail/default_ops.hpp +++ b/include/boost/multiprecision/detail/default_ops.hpp @@ -405,6 +405,7 @@ void eval_abs(T& result, const T& arg) template void eval_fabs(T& result, const T& arg) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The fabs function is only valid for floating point types."); typedef typename T::signed_types type_list; typedef typename mpl::front::type front; result = arg; @@ -415,12 +416,14 @@ void eval_fabs(T& result, const T& arg) template inline int eval_fpclassify(const Backend& arg) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The fpclassify function is only valid for floating point types."); return is_zero(arg) ? FP_ZERO : FP_NORMAL; } template inline void eval_fmod(T& result, const T& a, const T& b) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The fmod function is only valid for floating point types."); if((&result == &a) || (&result == &b)) { T temp; @@ -440,6 +443,7 @@ inline void eval_fmod(T& result, const T& a, const T& b) template inline void eval_trunc(T& result, const T& a) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The trunc function is only valid for floating point types."); int c = eval_fpclassify(a); if(c == FP_NAN || c == FP_INFINITE) { @@ -455,6 +459,7 @@ inline void eval_trunc(T& result, const T& a) template inline void eval_round(T& result, const T& a) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The round function is only valid for floating point types."); typedef typename boost::multiprecision::detail::canonical::type fp_type; int c = eval_fpclassify(a); if(c == FP_NAN || c == FP_INFINITE) @@ -474,6 +479,33 @@ inline void eval_round(T& result, const T& a) } } +template +inline typename enable_if >::type eval_gcd(T& result, const T& a, const Arithmetic& b) +{ + typedef typename boost::multiprecision::detail::canonical::type si_type; + T t; + t = static_cast(b); + eval_gcd(result, a, t); +} +template +inline typename enable_if >::type eval_gcd(T& result, const Arithmetic& a, const T& b) +{ + eval_gcd(result, b, a); +} +template +inline typename enable_if >::type eval_lcm(T& result, const T& a, const Arithmetic& b) +{ + typedef typename boost::multiprecision::detail::canonical::type si_type; + T t; + t = static_cast(b); + eval_lcm(result, a, t); +} +template +inline typename enable_if >::type eval_lcm(T& result, const Arithmetic& a, const T& b) +{ + eval_lcm(result, b, a); +} + // // These have to implemented by the backend, declared here so that our macro generated code compiles OK. // @@ -831,6 +863,12 @@ struct BOOST_JOIN(func, _funct)\ using default_ops:: BOOST_JOIN(eval_,func);\ BOOST_JOIN(eval_,func)(result, arg, a);\ }\ + template \ + void operator()(Backend& result, const Arithmetic& arg, const Backend& a)const\ + {\ + using default_ops:: BOOST_JOIN(eval_,func);\ + BOOST_JOIN(eval_,func)(result, arg, a);\ + }\ };\ \ }\ @@ -960,6 +998,52 @@ func(const detail::mp_exp& arg, const Arithmetic& a)\ a\ );\ }\ +template \ +typename enable_if<\ + is_arithmetic,\ + detail::mp_exp<\ + detail::function\ + , detail::BOOST_JOIN(func, _funct) \ + , Arithmetic \ + , mp_number\ + > \ +>::type \ +func(const Arithmetic& arg, const mp_number& a)\ +{\ + return detail::mp_exp<\ + detail::function\ + , detail::BOOST_JOIN(func, _funct) \ + , Arithmetic \ + , mp_number\ + >(\ + detail::BOOST_JOIN(func, _funct)() \ + , arg,\ + a\ + );\ +}\ +template \ +typename enable_if<\ + is_arithmetic,\ + detail::mp_exp<\ + detail::function\ + , detail::BOOST_JOIN(func, _funct) >::type> \ + , Arithmetic \ + , detail::mp_exp\ + > \ +>::type \ +func(const Arithmetic& arg, const detail::mp_exp& a)\ +{\ + return detail::mp_exp<\ + detail::function\ + , detail::BOOST_JOIN(func, _funct) >::type> \ + , Arithmetic \ + , detail::mp_exp\ + >(\ + detail::BOOST_JOIN(func, _funct) >::type>() \ + , arg,\ + a\ + );\ +}\ #define HETERO_BINARY_OP_FUNCTOR(func, Arg2)\ @@ -1042,6 +1126,13 @@ BINARY_OP_FUNCTOR(pow) BINARY_OP_FUNCTOR(fmod) BINARY_OP_FUNCTOR(atan2) +// +// Integer functions: +// +BINARY_OP_FUNCTOR(gcd) +BINARY_OP_FUNCTOR(lcm) + + #undef BINARY_OP_FUNCTOR #undef UNARY_OP_FUNCTOR diff --git a/include/boost/multiprecision/detail/functions/pow.hpp b/include/boost/multiprecision/detail/functions/pow.hpp index 444f7f64..94492bb7 100644 --- a/include/boost/multiprecision/detail/functions/pow.hpp +++ b/include/boost/multiprecision/detail/functions/pow.hpp @@ -86,6 +86,7 @@ inline void pow_imp(T& result, const T& t, const U& p, const mpl::true_&) template inline void eval_pow(T& result, const T& t, const U& p) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The pow function is only valid for floating point types."); typedef typename is_integral::type tag_type; detail::pow_imp(result, t, p, tag_type()); } @@ -184,6 +185,7 @@ void hyp1F0(T& H1F0, const T& a, const T& x) template void eval_exp(T& result, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The ldexp function is only valid for floating point types."); if(&x == &result) { T temp; @@ -323,6 +325,7 @@ void eval_exp(T& result, const T& x) template void eval_log(T& result, const T& arg) { + BOOST_STATIC_ASSERT_MSG(number_category::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, @@ -410,6 +413,7 @@ const T& get_constant_log10() template void eval_log10(T& result, const T& arg) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The fabs function is only valid for floating point types."); eval_log(result, arg); divide(result, get_constant_log10()); } @@ -417,6 +421,7 @@ void eval_log10(T& result, const T& arg) template inline void eval_pow(T& result, const T& x, const T& a) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The pow function is only valid for floating point types."); typedef typename boost::multiprecision::detail::canonical::type si_type; typedef typename boost::multiprecision::detail::canonical::type ui_type; typedef typename T::exponent_type exp_type; @@ -636,18 +641,21 @@ namespace detail{ template inline void eval_sinh(T& result, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The sinh function is only valid for floating point types."); detail::sinhcosh(x, &result, static_cast(0)); } template inline void eval_cosh(T& result, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The cosh function is only valid for floating point types."); detail::sinhcosh(x, static_cast(0), &result); } template inline void eval_tanh(T& result, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The tanh function is only valid for floating point types."); T c; detail::sinhcosh(x, &result, &c); divide(result, c); diff --git a/include/boost/multiprecision/detail/functions/trig.hpp b/include/boost/multiprecision/detail/functions/trig.hpp index 17a7f353..a5c52817 100644 --- a/include/boost/multiprecision/detail/functions/trig.hpp +++ b/include/boost/multiprecision/detail/functions/trig.hpp @@ -68,6 +68,7 @@ void hyp0F1(T& result, const T& b, const T& x) template void eval_sin(T& result, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The sin function is only valid for floating point types."); if(&result == &x) { T temp; @@ -213,6 +214,7 @@ void eval_sin(T& result, const T& x) template void eval_cos(T& result, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The cos function is only valid for floating point types."); if(&result == &x) { T temp; @@ -355,6 +357,7 @@ void eval_cos(T& result, const T& x) template void eval_tan(T& result, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The tan function is only valid for floating point types."); T t; eval_sin(result, x); eval_cos(t, x); @@ -426,6 +429,7 @@ void hyp2F1(T& result, const T& a, const T& b, const T& c, const T& x) template void eval_asin(T& result, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The asin function is only valid for floating point types."); typedef typename boost::multiprecision::detail::canonical::type si_type; typedef typename boost::multiprecision::detail::canonical::type ui_type; typedef typename T::exponent_type exp_type; @@ -535,6 +539,7 @@ void eval_asin(T& result, const T& x) template inline void eval_acos(T& result, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The acos function is only valid for floating point types."); typedef typename boost::multiprecision::detail::canonical::type ui_type; switch(eval_fpclassify(x)) @@ -576,6 +581,7 @@ inline void eval_acos(T& result, const T& x) template void eval_atan(T& result, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The atan function is only valid for floating point types."); typedef typename boost::multiprecision::detail::canonical::type si_type; typedef typename boost::multiprecision::detail::canonical::type ui_type; typedef typename T::exponent_type exp_type; @@ -666,6 +672,7 @@ void eval_atan(T& result, const T& x) template void eval_atan2(T& result, const T& y, const T& x) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The atan2 function is only valid for floating point types."); if(&result == &y) { T temp(y); @@ -763,6 +770,7 @@ void eval_atan2(T& result, const T& y, const T& x) template typename disable_if >::type eval_atan2(T& result, const T& a, const Arithmetic& b) { + BOOST_STATIC_ASSERT_MSG(number_category::value == number_kind_floating_point, "The atan2 function is only valid for floating point types."); T x; x = static_cast::type>(b); eval_atan2(result, a, x); diff --git a/include/boost/multiprecision/fixed_int.hpp b/include/boost/multiprecision/fixed_int.hpp index 449c90d2..a8b99318 100644 --- a/include/boost/multiprecision/fixed_int.hpp +++ b/include/boost/multiprecision/fixed_int.hpp @@ -1254,6 +1254,8 @@ inline void complement(fixed_int& result, const fixed_int inline void left_shift(fixed_int& result, double_limb_type s) { + if(!s) + return; if(s >= Bits) { result = static_cast(0); @@ -1287,6 +1289,8 @@ inline void left_shift(fixed_int& result, double_limb_type s) template inline void right_shift(fixed_int& result, double_limb_type s) { + if(!s) + return; limb_type fill = (Signed && (result.data()[0] & fixed_int::sign_bit_mask)) ? fixed_int::max_limb_value : 0u; if(s >= Bits) { @@ -1327,7 +1331,22 @@ inline void right_shift(fixed_int& result, double_limb_type s) } template -inline typename enable_if, void>::type convert_to(R* result, const fixed_int& backend) +inline typename enable_if, void>::type convert_to(R* result, const fixed_int& backend, const mpl::true_&) +{ + if(backend.data()[0] & fixed_int::sign_bit_mask) + { + fixed_int t(backend); + t.negate(); + convert_to(result, t, mpl::false_()); + *result = -*result; + return; + } + else + convert_to(result, backend, mpl::false_()); +} + +template +inline typename enable_if, void>::type convert_to(R* result, const fixed_int& backend, const mpl::false_&) { unsigned shift = (fixed_int::limb_count - 1) * fixed_int::limb_bits; *result = 0; @@ -1338,9 +1357,24 @@ inline typename enable_if, void>::type convert_to(R* result, cons } } +template +inline typename enable_if, void>::type convert_to(R* result, const fixed_int& backend) +{ + typedef mpl::bool_::is_signed> tag_type; + convert_to(result, backend, tag_type()); +} + template inline typename enable_if, void>::type convert_to(R* result, const fixed_int& backend) { + if(Signed && (backend.data()[0] & fixed_int::sign_bit_mask)) + { + fixed_int t(backend); + t.negate(); + convert_to(result, t); + *result = -*result; + return; + } unsigned shift = (fixed_int::limb_count - 1) * fixed_int::limb_bits; *result = 0; for(unsigned i = 0; i < fixed_int::limb_count; ++i) @@ -1371,6 +1405,116 @@ inline int get_sign(const fixed_int& val) return is_zero(val) ? 0 : val.data()[0] & fixed_int::sign_bit_mask ? -1 : 1; } +namespace detail{ +// +// Get the location of the least-significant-bit: +// +template +inline unsigned get_lsb(const fixed_int& a) +{ + BOOST_ASSERT(get_sign(a) != 0); + + unsigned result = 0; + // + // Find the index of the least significant limb that is non-zero: + // + int index = fixed_int::limb_count - 1; + while(!a.data()[index] && index) + --index; + // + // Find the index of the least significant bit within that limb: + // + limb_type l = a.data()[index]; + while(!(l & 1u)) + { + l >>= 1; + ++result; + } + + return result + (fixed_int::limb_count - 1 - index) * fixed_int::limb_bits; +} + +} + +template +inline void eval_gcd(fixed_int& result, const fixed_int& a, const fixed_int& b) +{ + int shift; + + fixed_int u(a), v(b); + + int s = get_sign(u); + + /* GCD(0,x) := x */ + if(s < 0) + { + u.negate(); + } + else if(s == 0) + { + result = v; + return; + } + s = get_sign(v); + if(s < 0) + { + v.negate(); + } + else if(s == 0) + { + result = u; + return; + } + + /* Let shift := lg K, where K is the greatest power of 2 + dividing both u and v. */ + + unsigned us = detail::get_lsb(u); + unsigned vs = detail::get_lsb(v); + shift = (std::min)(us, vs); + right_shift(u, us); + right_shift(v, vs); + + do + { + /* Now u and v are both odd, so diff(u, v) is even. + Let u = min(u, v), v = diff(u, v)/2. */ + if(u.compare(v) > 0) + u.swap(v); + subtract(v, u); + // Termination condition tries not to do a full compare if possible: + if(!v.data()[fixed_int::limb_count - 1] && is_zero(v)) + break; + vs = detail::get_lsb(v); + right_shift(v, vs); + BOOST_ASSERT((v.data()[fixed_int::limb_count - 1] & 1)); + BOOST_ASSERT((u.data()[fixed_int::limb_count - 1] & 1)); + } + while(true); + + result = u; + left_shift(result, shift); +} + +template +inline void eval_lcm(fixed_int& result, const fixed_int& a, const fixed_int& b) +{ + fixed_int t; + eval_gcd(t, a, b); + + if(is_zero(t)) + { + result = 0; + } + else + { + divide(result, a, t); + multiply(result, b); + } + if(get_sign(result) < 0) + result.negate(); +} + template struct number_category > : public mpl::int_{}; diff --git a/include/boost/multiprecision/gmp.hpp b/include/boost/multiprecision/gmp.hpp index d246b07b..4311e5b0 100644 --- a/include/boost/multiprecision/gmp.hpp +++ b/include/boost/multiprecision/gmp.hpp @@ -1392,6 +1392,31 @@ inline void eval_abs(gmp_int& result, const gmp_int& val) mpz_abs(result.data(), val.data()); } +inline void eval_gcd(gmp_int& result, const gmp_int& a, const gmp_int& b) +{ + mpz_gcd(result.data(), a.data(), b.data()); +} +inline void eval_lcm(gmp_int& result, const gmp_int& a, const gmp_int& b) +{ + mpz_lcm(result.data(), a.data(), b.data()); +} +inline void eval_gcd(gmp_int& result, const gmp_int& a, const unsigned long b) +{ + mpz_gcd_ui(result.data(), a.data(), b); +} +inline void eval_lcm(gmp_int& result, const gmp_int& a, const unsigned long b) +{ + mpz_lcm_ui(result.data(), a.data(), b); +} +inline void eval_gcd(gmp_int& result, const gmp_int& a, const long b) +{ + mpz_gcd_ui(result.data(), a.data(), std::abs(b)); +} +inline void eval_lcm(gmp_int& result, const gmp_int& a, const long b) +{ + mpz_lcm_ui(result.data(), a.data(), std::abs(b)); +} + struct gmp_rational; void add(gmp_rational& t, const gmp_rational& o); diff --git a/include/boost/multiprecision/tommath.hpp b/include/boost/multiprecision/tommath.hpp index 41757e32..41eeaf76 100644 --- a/include/boost/multiprecision/tommath.hpp +++ b/include/boost/multiprecision/tommath.hpp @@ -432,6 +432,14 @@ inline void eval_abs(tommath_int& result, const tommath_int& val) { detail::check_tommath_result(mp_abs(const_cast< ::mp_int*>(&val.data()), &result.data())); } +inline void eval_gcd(tommath_int& result, const tommath_int& a, const tommath_int& b) +{ + detail::check_tommath_result(mp_gcd(const_cast< ::mp_int*>(&a.data()), const_cast< ::mp_int*>(&b.data()), const_cast< ::mp_int*>(&result.data()))); +} +inline void eval_lcm(tommath_int& result, const tommath_int& a, const tommath_int& b) +{ + detail::check_tommath_result(mp_lcm(const_cast< ::mp_int*>(&a.data()), const_cast< ::mp_int*>(&b.data()), const_cast< ::mp_int*>(&result.data()))); +} template<> diff --git a/performance/performance_test.cpp b/performance/performance_test.cpp index a365ee82..bb7ddc92 100644 --- a/performance/performance_test.cpp +++ b/performance/performance_test.cpp @@ -15,6 +15,7 @@ && !defined(TEST_FIXED_INT) # define TEST_MPF # define TEST_MPZ +# define TEST_MPQ # define TEST_MPFR # define TEST_CPP_FLOAT # define TEST_MPQ @@ -299,6 +300,16 @@ struct tester } return boost::chrono::duration_cast >(w.elapsed()).count(); } + double test_gcd() + { + stopwatch w; + for(unsigned i = 0; i < 1000; ++i) + { + for(unsigned i = 0; i < b.size(); ++i) + a[i] = gcd(b[i], c[i]); + } + return boost::chrono::duration_cast >(w.elapsed()).count(); + } private: T generate_random() { @@ -445,6 +456,7 @@ void test_int_ops(tester& t, const char* type, unsigned precision, co report_result(cat, type, "|(int)", precision, t.test_or_int()); report_result(cat, type, "&(int)", precision, t.test_and_int()); report_result(cat, type, "^(int)", precision, t.test_xor_int()); + report_result(cat, type, "gcd", precision, t.test_gcd()); } template void test_int_ops(tester& t, const char* type, unsigned precision, const U&) @@ -564,7 +576,8 @@ int main() test("gmp_int", 256); test("gmp_int", 512); test("gmp_int", 1024); - +#endif +#ifdef TEST_MPQ test("mpq_rational", 64); test("mpq_rational", 128); test("mpq_rational", 256); diff --git a/test/test_arithmetic.cpp b/test/test_arithmetic.cpp index 282baac6..03da9c63 100644 --- a/test/test_arithmetic.cpp +++ b/test/test_arithmetic.cpp @@ -9,6 +9,7 @@ #include #include +#include #if !defined(TEST_MPF_50) && !defined(TEST_MPF) && !defined(TEST_BACKEND) && !defined(TEST_MPZ) && \ !defined(TEST_CPP_FLOAT) && !defined(TEST_MPFR) && !defined(TEST_MPFR_50) && !defined(TEST_MPQ) \ @@ -345,6 +346,14 @@ void test_integer_ops(const boost::mpl::int_::is_signed) { a = -20; @@ -355,6 +364,15 @@ void test_integer_ops(const boost::mpl::int_() == n1); + BOOST_TEST(Real(n2).template convert_to() == n2); + BOOST_TEST(Real(n3).template convert_to() == n3); + BOOST_TEST(Real(n4).template convert_to() == n4); #if defined(TEST_MPFR) || defined(TEST_MPFR_50) Num tol = 10 * std::numeric_limits::epsilon(); #else diff --git a/test/test_fixed_int.cpp b/test/test_fixed_int.cpp index cbf2458e..e69e2e29 100644 --- a/test/test_fixed_int.cpp +++ b/test/test_fixed_int.cpp @@ -119,6 +119,8 @@ int main() BOOST_CHECK_EQUAL(mpz_int(si|a).str(), packed_type(si|a1).str()); BOOST_CHECK_EQUAL(mpz_int(si&a).str(), packed_type(si&a1).str()); BOOST_CHECK_EQUAL(mpz_int(si^a).str(), packed_type(si^a1).str()); + BOOST_CHECK_EQUAL(mpz_int(gcd(a, b)).str(), packed_type(gcd(a1, b1)).str()); + BOOST_CHECK_EQUAL(mpz_int(lcm(c, d)).str(), packed_type(lcm(c1, d1)).str()); if(last_error_count != boost::detail::test_errors()) { @@ -135,6 +137,7 @@ int main() std::cout << "d1 = " << d1 << std::endl; std::cout << "a + b = " << a+b << std::endl; std::cout << "a1 + b1 = " << a1+b1 << std::endl; + std::cout << std::dec; std::cout << "a - b = " << a-b << std::endl; std::cout << "a1 - b1 = " << a1-b1 << std::endl; std::cout << "-a + b = " << mpz_int(-a)+b << std::endl;