/////////////////////////////////////////////////////////////// // Copyright 2012 John Maddock. Distributed under the Boost // Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt #ifndef BOOST_MP_INT_FUNC_HPP #define BOOST_MP_INT_FUNC_HPP #include namespace boost { namespace multiprecision { namespace default_ops { template inline BOOST_MP_CXX14_CONSTEXPR void eval_qr(const Backend& x, const Backend& y, Backend& q, Backend& r) { eval_divide(q, x, y); eval_modulus(r, x, y); } template inline BOOST_MP_CXX14_CONSTEXPR Integer eval_integer_modulus(const Backend& x, Integer val) { BOOST_MP_USING_ABS using default_ops::eval_convert_to; using default_ops::eval_modulus; using int_type = typename boost::multiprecision::detail::canonical::type; Backend t; eval_modulus(t, x, static_cast(val)); Integer result(0); eval_convert_to(&result, t); return abs(result); } template inline BOOST_MP_CXX14_CONSTEXPR void eval_gcd(B& result, const B& a, const B& b) { using default_ops::eval_get_sign; using default_ops::eval_is_zero; using default_ops::eval_lsb; int shift(0); B u(a), v(b); int s = eval_get_sign(u); /* GCD(0,x) := x */ if (s < 0) { u.negate(); } else if (s == 0) { result = v; return; } s = eval_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 = eval_lsb(u); unsigned vs = eval_lsb(v); shift = (std::min)(us, vs); eval_right_shift(u, us); eval_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. */ s = u.compare(v); if (s > 0) u.swap(v); if (s == 0) break; eval_subtract(v, u); vs = eval_lsb(v); eval_right_shift(v, vs); } while (true); result = u; eval_left_shift(result, shift); } template inline BOOST_MP_CXX14_CONSTEXPR void eval_lcm(B& result, const B& a, const B& b) { using ui_type = typename std::tuple_element<0, typename B::unsigned_types>::type; B t; eval_gcd(t, a, b); if (eval_is_zero(t)) { result = static_cast(0); } else { eval_divide(result, a, t); eval_multiply(result, b); } if (eval_get_sign(result) < 0) result.negate(); } } // namespace default_ops template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_integer>::type divide_qr(const number& x, const number& y, number& q, number& r) { using default_ops::eval_qr; eval_qr(x.backend(), y.backend(), q.backend(), r.backend()); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_integer>::type divide_qr(const number& x, const multiprecision::detail::expression& y, number& q, number& r) { divide_qr(x, number(y), q, r); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_integer>::type divide_qr(const multiprecision::detail::expression& x, const number& y, number& q, number& r) { divide_qr(number(x), y, q, r); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_integer>::type divide_qr(const multiprecision::detail::expression& x, const multiprecision::detail::expression& y, number& q, number& r) { divide_qr(number(x), number(y), q, r); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value && (number_category::value == number_kind_integer), Integer>::type integer_modulus(const number& x, Integer val) { using default_ops::eval_integer_modulus; return eval_integer_modulus(x.backend(), val); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value && (number_category::result_type>::value == number_kind_integer), Integer>::type integer_modulus(const multiprecision::detail::expression& x, Integer val) { using result_type = typename multiprecision::detail::expression::result_type; return integer_modulus(result_type(x), val); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_integer, unsigned>::type lsb(const number& x) { using default_ops::eval_lsb; return eval_lsb(x.backend()); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::result_type>::value == number_kind_integer, unsigned>::type lsb(const multiprecision::detail::expression& x) { using number_type = typename multiprecision::detail::expression::result_type; number_type n(x); using default_ops::eval_lsb; return eval_lsb(n.backend()); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_integer, unsigned>::type msb(const number& x) { using default_ops::eval_msb; return eval_msb(x.backend()); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::result_type>::value == number_kind_integer, unsigned>::type msb(const multiprecision::detail::expression& x) { using number_type = typename multiprecision::detail::expression::result_type; number_type n(x); using default_ops::eval_msb; return eval_msb(n.backend()); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_integer, bool>::type bit_test(const number& x, unsigned index) { using default_ops::eval_bit_test; return eval_bit_test(x.backend(), index); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::result_type>::value == number_kind_integer, bool>::type bit_test(const multiprecision::detail::expression& x, unsigned index) { using number_type = typename multiprecision::detail::expression::result_type; number_type n(x); using default_ops::eval_bit_test; return eval_bit_test(n.backend(), index); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_integer, number&>::type bit_set(number& x, unsigned index) { using default_ops::eval_bit_set; eval_bit_set(x.backend(), index); return x; } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_integer, number&>::type bit_unset(number& x, unsigned index) { using default_ops::eval_bit_unset; eval_bit_unset(x.backend(), index); return x; } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_integer, number&>::type bit_flip(number& x, unsigned index) { using default_ops::eval_bit_flip; eval_bit_flip(x.backend(), index); return x; } namespace default_ops { // // Within powm, we need a type with twice as many digits as the argument type, define // a traits class to obtain that type: // template struct double_precision_type { using type = Backend; }; // // If the exponent is a signed integer type, then we need to // check the value is positive: // template inline BOOST_MP_CXX14_CONSTEXPR void check_sign_of_backend(const Backend& v, const std::integral_constant) { if (eval_get_sign(v) < 0) { BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); } } template inline BOOST_MP_CXX14_CONSTEXPR void check_sign_of_backend(const Backend&, const std::integral_constant) {} // // Calculate (a^p)%c: // template BOOST_MP_CXX14_CONSTEXPR void eval_powm(Backend& result, const Backend& a, const Backend& p, const Backend& c) { using default_ops::eval_bit_test; using default_ops::eval_get_sign; using default_ops::eval_modulus; using default_ops::eval_multiply; using default_ops::eval_right_shift; using double_type = typename double_precision_type::type ; using ui_type = typename boost::multiprecision::detail::canonical::type; check_sign_of_backend(p, std::integral_constant >::is_signed>()); double_type x, y(a), b(p), t; x = ui_type(1u); while (eval_get_sign(b) > 0) { if (eval_bit_test(b, 0)) { eval_multiply(t, x, y); eval_modulus(x, t, c); } eval_multiply(t, y, y); eval_modulus(y, t, c); eval_right_shift(b, ui_type(1)); } Backend x2(x); eval_modulus(result, x2, c); } template BOOST_MP_CXX14_CONSTEXPR void eval_powm(Backend& result, const Backend& a, const Backend& p, Integer c) { using double_type = typename double_precision_type::type ; using ui_type = typename boost::multiprecision::detail::canonical::type; using i1_type = typename boost::multiprecision::detail::canonical::type ; using i2_type = typename boost::multiprecision::detail::canonical::type ; using default_ops::eval_bit_test; using default_ops::eval_get_sign; using default_ops::eval_modulus; using default_ops::eval_multiply; using default_ops::eval_right_shift; check_sign_of_backend(p, std::integral_constant >::is_signed>()); if (eval_get_sign(p) < 0) { BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); } double_type x, y(a), b(p), t; x = ui_type(1u); while (eval_get_sign(b) > 0) { if (eval_bit_test(b, 0)) { eval_multiply(t, x, y); eval_modulus(x, t, static_cast(c)); } eval_multiply(t, y, y); eval_modulus(y, t, static_cast(c)); eval_right_shift(b, ui_type(1)); } Backend x2(x); eval_modulus(result, x2, static_cast(c)); } template BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c) { using double_type = typename double_precision_type::type ; using ui_type = typename boost::multiprecision::detail::canonical::type; using default_ops::eval_bit_test; using default_ops::eval_get_sign; using default_ops::eval_modulus; using default_ops::eval_multiply; using default_ops::eval_right_shift; double_type x, y(a), t; x = ui_type(1u); while (b > 0) { if (b & 1) { eval_multiply(t, x, y); eval_modulus(x, t, c); } eval_multiply(t, y, y); eval_modulus(y, t, c); b >>= 1; } Backend x2(x); eval_modulus(result, x2, c); } template BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value && boost::multiprecision::detail::is_integral::value>::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c) { if (b < 0) { BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); } eval_powm(result, a, static_cast::type>(b), c); } template BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c) { using double_type = typename double_precision_type::type ; using ui_type = typename boost::multiprecision::detail::canonical::type; using i1_type = typename boost::multiprecision::detail::canonical::type ; using i2_type = typename boost::multiprecision::detail::canonical::type ; using default_ops::eval_bit_test; using default_ops::eval_get_sign; using default_ops::eval_modulus; using default_ops::eval_multiply; using default_ops::eval_right_shift; double_type x, y(a), t; x = ui_type(1u); while (b > 0) { if (b & 1) { eval_multiply(t, x, y); eval_modulus(x, t, static_cast(c)); } eval_multiply(t, y, y); eval_modulus(y, t, static_cast(c)); b >>= 1; } Backend x2(x); eval_modulus(result, x2, static_cast(c)); } template BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value && boost::multiprecision::detail::is_integral::value>::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c) { if (b < 0) { BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); } eval_powm(result, a, static_cast::type>(b), c); } struct powm_func { template BOOST_MP_CXX14_CONSTEXPR void operator()(T& result, const T& b, const U& p, const V& m) const { eval_powm(result, b, p, m); } template BOOST_MP_CXX14_CONSTEXPR void operator()(R& result, const T& b, const U& p, const V& m) const { T temp; eval_powm(temp, b, p, m); result = std::move(temp); } }; } // namespace default_ops template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if< (number_category::value == number_kind_integer) && (is_number::value || is_number_expression::value) && (is_number::value || is_number_expression::value || boost::multiprecision::detail::is_integral::value) && (is_number::value || is_number_expression::value || boost::multiprecision::detail::is_integral::value), typename std::conditional< is_no_et_number::value, T, typename std::conditional< is_no_et_number::value, U, typename std::conditional< is_no_et_number::value, V, detail::expression >::type>::type>::type>::type powm(const T& b, const U& p, const V& mod) { return detail::expression( default_ops::powm_func(), b, p, mod); } }} // namespace boost::multiprecision #endif