]>
Commit | Line | Data |
---|---|---|
7c673cae | 1 | /////////////////////////////////////////////////////////////// |
20effc67 TL |
2 | // Copyright 2012-2020 John Maddock. |
3 | // Copyright 2020 Madhur Chauhan. | |
1e59de90 | 4 | // Copyright 2021 Matt Borland. |
20effc67 TL |
5 | // Distributed under the Boost Software License, Version 1.0. |
6 | // (See accompanying file LICENSE_1_0.txt or copy at | |
7 | // https://www.boost.org/LICENSE_1_0.txt) | |
7c673cae FG |
8 | // |
9 | // Comparison operators for cpp_int_backend: | |
10 | // | |
11 | #ifndef BOOST_MP_CPP_INT_MISC_HPP | |
12 | #define BOOST_MP_CPP_INT_MISC_HPP | |
13 | ||
1e59de90 TL |
14 | #include <boost/multiprecision/detail/standalone_config.hpp> |
15 | #include <boost/multiprecision/detail/number_base.hpp> | |
16 | #include <boost/multiprecision/cpp_int/cpp_int_config.hpp> | |
17 | #include <boost/multiprecision/detail/float128_functions.hpp> | |
18 | #include <boost/multiprecision/detail/assert.hpp> | |
92f5a8d4 | 19 | #include <boost/multiprecision/detail/constexpr.hpp> |
7c673cae | 20 | #include <boost/multiprecision/detail/bitscan.hpp> // lsb etc |
1e59de90 TL |
21 | #include <boost/multiprecision/detail/hash.hpp> |
22 | #include <boost/multiprecision/detail/no_exceptions_support.hpp> | |
20effc67 | 23 | #include <numeric> // std::gcd |
1e59de90 TL |
24 | #include <type_traits> |
25 | #include <stdexcept> | |
26 | #include <cmath> | |
27 | ||
28 | #ifndef BOOST_MP_STANDALONE | |
29 | #include <boost/integer/common_factor_rt.hpp> | |
30 | #endif | |
31 | ||
32 | #ifdef BOOST_MP_MATH_AVAILABLE | |
33 | #include <boost/math/special_functions/next.hpp> | |
34 | #endif | |
7c673cae FG |
35 | |
36 | #ifdef BOOST_MSVC | |
37 | #pragma warning(push) | |
92f5a8d4 TL |
38 | #pragma warning(disable : 4702) |
39 | #pragma warning(disable : 4127) // conditional expression is constant | |
40 | #pragma warning(disable : 4146) // unary minus operator applied to unsigned type, result still unsigned | |
7c673cae FG |
41 | #endif |
42 | ||
1e59de90 TL |
43 | // Forward decleration of gcd and lcm functions |
44 | namespace boost { namespace multiprecision { namespace detail { | |
45 | ||
46 | template <typename T> | |
47 | inline BOOST_CXX14_CONSTEXPR T constexpr_gcd(T a, T b) noexcept; | |
48 | ||
49 | template <typename T> | |
50 | inline BOOST_CXX14_CONSTEXPR T constexpr_lcm(T a, T b) noexcept; | |
51 | ||
52 | }}} // namespace boost::multiprecision::detail | |
53 | ||
92f5a8d4 | 54 | namespace boost { namespace multiprecision { namespace backends { |
7c673cae | 55 | |
20effc67 TL |
56 | template <class T, bool has_limits = std::numeric_limits<T>::is_specialized> |
57 | struct numeric_limits_workaround : public std::numeric_limits<T> | |
58 | { | |
59 | }; | |
60 | template <class R> | |
61 | struct numeric_limits_workaround<R, false> | |
62 | { | |
1e59de90 TL |
63 | static constexpr unsigned digits = ~static_cast<R>(0) < 0 ? sizeof(R) * CHAR_BIT - 1 : sizeof(R) * CHAR_BIT; |
64 | static constexpr R (min)(){ return (static_cast<R>(-1) < 0) ? static_cast<R>(1) << digits : 0; } | |
65 | static constexpr R (max)() { return (static_cast<R>(-1) < 0) ? ~(static_cast<R>(1) << digits) : ~static_cast<R>(0); } | |
20effc67 TL |
66 | }; |
67 | ||
7c673cae | 68 | template <class R, class CppInt> |
1e59de90 | 69 | BOOST_MP_CXX14_CONSTEXPR void check_in_range(const CppInt& val, const std::integral_constant<int, checked>&) |
7c673cae | 70 | { |
1e59de90 | 71 | using cast_type = typename boost::multiprecision::detail::canonical<R, CppInt>::type; |
20effc67 | 72 | |
92f5a8d4 | 73 | if (val.sign()) |
7c673cae | 74 | { |
1e59de90 TL |
75 | BOOST_IF_CONSTEXPR (boost::multiprecision::detail::is_signed<R>::value == false) |
76 | BOOST_MP_THROW_EXCEPTION(std::range_error("Attempt to assign a negative value to an unsigned type.")); | |
20effc67 | 77 | if (val.compare(static_cast<cast_type>((numeric_limits_workaround<R>::min)())) < 0) |
1e59de90 | 78 | BOOST_MP_THROW_EXCEPTION(std::overflow_error("Could not convert to the target type - -value is out of range.")); |
7c673cae FG |
79 | } |
80 | else | |
81 | { | |
20effc67 | 82 | if (val.compare(static_cast<cast_type>((numeric_limits_workaround<R>::max)())) > 0) |
1e59de90 | 83 | BOOST_MP_THROW_EXCEPTION(std::overflow_error("Could not convert to the target type - -value is out of range.")); |
7c673cae FG |
84 | } |
85 | } | |
86 | template <class R, class CppInt> | |
1e59de90 | 87 | inline BOOST_MP_CXX14_CONSTEXPR void check_in_range(const CppInt& /*val*/, const std::integral_constant<int, unchecked>&) noexcept {} |
7c673cae | 88 | |
1e59de90 TL |
89 | inline BOOST_MP_CXX14_CONSTEXPR void check_is_negative(const std::integral_constant<bool, true>&) noexcept {} |
90 | inline void check_is_negative(const std::integral_constant<bool, false>&) | |
7c673cae | 91 | { |
1e59de90 | 92 | BOOST_MP_THROW_EXCEPTION(std::range_error("Attempt to assign a negative value to an unsigned type.")); |
7c673cae FG |
93 | } |
94 | ||
95 | template <class Integer> | |
1e59de90 | 96 | inline BOOST_MP_CXX14_CONSTEXPR Integer negate_integer(Integer i, const std::integral_constant<bool, true>&) noexcept |
7c673cae FG |
97 | { |
98 | return -i; | |
99 | } | |
100 | template <class Integer> | |
1e59de90 | 101 | inline BOOST_MP_CXX14_CONSTEXPR Integer negate_integer(Integer i, const std::integral_constant<bool, false>&) noexcept |
7c673cae | 102 | { |
92f5a8d4 | 103 | return ~(i - 1); |
7c673cae FG |
104 | } |
105 | ||
1e59de90 TL |
106 | template <class R, std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
107 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<R>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, void>::type | |
92f5a8d4 | 108 | eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& backend) |
7c673cae | 109 | { |
1e59de90 | 110 | using checked_type = std::integral_constant<int, Checked1>; |
7c673cae FG |
111 | check_in_range<R>(backend, checked_type()); |
112 | ||
20effc67 | 113 | BOOST_IF_CONSTEXPR(numeric_limits_workaround<R>::digits < cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits) |
11fdf7f2 | 114 | { |
1e59de90 | 115 | if ((backend.sign() && boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value) && (1 + static_cast<boost::multiprecision::limb_type>((std::numeric_limits<R>::max)()) <= backend.limbs()[0])) |
11fdf7f2 | 116 | { |
20effc67 | 117 | *result = (numeric_limits_workaround<R>::min)(); |
11fdf7f2 TL |
118 | return; |
119 | } | |
1e59de90 | 120 | else if (boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value && !backend.sign() && static_cast<boost::multiprecision::limb_type>((std::numeric_limits<R>::max)()) <= backend.limbs()[0]) |
11fdf7f2 | 121 | { |
20effc67 | 122 | *result = (numeric_limits_workaround<R>::max)(); |
11fdf7f2 TL |
123 | return; |
124 | } | |
125 | else | |
126 | *result = static_cast<R>(backend.limbs()[0]); | |
127 | } | |
128 | else | |
129 | *result = static_cast<R>(backend.limbs()[0]); | |
1e59de90 TL |
130 | std::size_t shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; |
131 | std::size_t i = 1; | |
20effc67 | 132 | BOOST_IF_CONSTEXPR(numeric_limits_workaround<R>::digits > cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits) |
7c673cae | 133 | { |
20effc67 | 134 | while ((i < backend.size()) && (shift < static_cast<unsigned>(numeric_limits_workaround<R>::digits - cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits))) |
11fdf7f2 TL |
135 | { |
136 | *result += static_cast<R>(backend.limbs()[i]) << shift; | |
137 | shift += cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
138 | ++i; | |
139 | } | |
140 | // | |
141 | // We have one more limb to extract, but may not need all the bits, so treat this as a special case: | |
142 | // | |
143 | if (i < backend.size()) | |
144 | { | |
20effc67 | 145 | const limb_type mask = numeric_limits_workaround<R>::digits - shift == cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits ? ~static_cast<limb_type>(0) : (static_cast<limb_type>(1u) << (numeric_limits_workaround<R>::digits - shift)) - 1; |
11fdf7f2 TL |
146 | *result += (static_cast<R>(backend.limbs()[i]) & mask) << shift; |
147 | if ((static_cast<R>(backend.limbs()[i]) & static_cast<limb_type>(~mask)) || (i + 1 < backend.size())) | |
148 | { | |
149 | // Overflow: | |
150 | if (backend.sign()) | |
151 | { | |
1e59de90 | 152 | check_is_negative(boost::multiprecision::detail::is_signed<R>()); |
20effc67 | 153 | *result = (numeric_limits_workaround<R>::min)(); |
11fdf7f2 | 154 | } |
1e59de90 | 155 | else if (boost::multiprecision::detail::is_signed<R>::value) |
20effc67 | 156 | *result = (numeric_limits_workaround<R>::max)(); |
11fdf7f2 TL |
157 | return; |
158 | } | |
159 | } | |
160 | } | |
161 | else if (backend.size() > 1) | |
162 | { | |
163 | // Overflow: | |
164 | if (backend.sign()) | |
165 | { | |
1e59de90 | 166 | check_is_negative(boost::multiprecision::detail::is_signed<R>()); |
20effc67 | 167 | *result = (numeric_limits_workaround<R>::min)(); |
11fdf7f2 | 168 | } |
1e59de90 | 169 | else if (boost::multiprecision::detail::is_signed<R>::value) |
20effc67 | 170 | *result = (numeric_limits_workaround<R>::max)(); |
11fdf7f2 | 171 | return; |
7c673cae | 172 | } |
92f5a8d4 | 173 | if (backend.sign()) |
7c673cae | 174 | { |
1e59de90 TL |
175 | check_is_negative(std::integral_constant<bool, boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value>()); |
176 | *result = negate_integer(*result, std::integral_constant<bool, boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value>()); | |
7c673cae FG |
177 | } |
178 | } | |
179 | ||
1e59de90 TL |
180 | template <class R, std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
181 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<std::is_floating_point<R>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, void>::type | |
182 | eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& backend) noexcept(boost::multiprecision::detail::is_arithmetic<R>::value) | |
7c673cae | 183 | { |
1e59de90 TL |
184 | BOOST_MP_FLOAT128_USING using std::ldexp; |
185 | if (eval_is_zero(backend)) | |
186 | { | |
187 | *result = 0.0f; | |
188 | return; | |
189 | } | |
190 | ||
191 | #ifdef BOOST_HAS_FLOAT128 | |
192 | std::ptrdiff_t bits_to_keep = std::is_same<R, float128_type>::value ? 113 : std::numeric_limits<R>::digits; | |
193 | #else | |
194 | std::ptrdiff_t bits_to_keep = std::numeric_limits<R>::digits; | |
195 | #endif | |
196 | std::ptrdiff_t bits = eval_msb_imp(backend) + 1; | |
197 | ||
198 | if (bits > bits_to_keep) | |
199 | { | |
200 | // Extract the bits we need, and then manually round the result: | |
201 | *result = 0.0f; | |
202 | typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::const_limb_pointer p = backend.limbs(); | |
203 | limb_type mask = ~static_cast<limb_type>(0u); | |
204 | std::size_t index = backend.size() - 1; | |
205 | std::size_t shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits * index; | |
206 | while (bits_to_keep > 0) | |
207 | { | |
208 | if (bits_to_keep < (std::ptrdiff_t)cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits) | |
209 | { | |
210 | if(index != backend.size() - 1) | |
211 | mask <<= cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits - bits_to_keep; | |
212 | else | |
213 | { | |
214 | std::ptrdiff_t bits_in_first_limb = bits % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
215 | if (bits_in_first_limb == 0) | |
216 | bits_in_first_limb = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
217 | if (bits_in_first_limb > bits_to_keep) | |
218 | mask <<= bits_in_first_limb - bits_to_keep; | |
219 | } | |
220 | } | |
221 | *result += ldexp(static_cast<R>(p[index] & mask), static_cast<int>(shift)); | |
222 | shift -= cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
223 | bits_to_keep -= (index == backend.size() - 1) && (bits % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits) | |
224 | ? bits % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits | |
225 | : cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
226 | --index; | |
227 | } | |
228 | // Perform rounding: | |
229 | bits -= 1 + std::numeric_limits<R>::digits; | |
230 | if (eval_bit_test(backend, static_cast<unsigned>(bits))) | |
231 | { | |
232 | if ((eval_lsb_imp(backend) < static_cast<std::size_t>(bits)) || eval_bit_test(backend, static_cast<std::size_t>(bits + 1))) | |
233 | { | |
234 | #ifdef BOOST_MP_MATH_AVAILABLE | |
235 | *result = boost::math::float_next(*result); | |
236 | #else | |
237 | *result = std::nextafter(*result, *result * 2); | |
238 | #endif | |
239 | } | |
240 | } | |
241 | } | |
242 | else | |
7c673cae | 243 | { |
1e59de90 TL |
244 | typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::const_limb_pointer p = backend.limbs(); |
245 | std::size_t shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
246 | *result = static_cast<R>(*p); | |
247 | for (std::size_t i = 1; i < backend.size(); ++i) | |
248 | { | |
249 | *result += static_cast<R>(ldexp(static_cast<long double>(p[i]), static_cast<int>(shift))); | |
250 | shift += cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
251 | } | |
7c673cae | 252 | } |
92f5a8d4 | 253 | if (backend.sign()) |
7c673cae FG |
254 | *result = -*result; |
255 | } | |
256 | ||
1e59de90 TL |
257 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
258 | BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, bool>::type | |
259 | eval_is_zero(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept | |
7c673cae FG |
260 | { |
261 | return (val.size() == 1) && (val.limbs()[0] == 0); | |
262 | } | |
1e59de90 TL |
263 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
264 | BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, int>::type | |
265 | eval_get_sign(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept | |
7c673cae FG |
266 | { |
267 | return eval_is_zero(val) ? 0 : val.sign() ? -1 : 1; | |
268 | } | |
1e59de90 TL |
269 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
270 | BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
271 | eval_abs(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value)) | |
7c673cae FG |
272 | { |
273 | result = val; | |
274 | result.sign(false); | |
275 | } | |
276 | ||
277 | // | |
278 | // Get the location of the least-significant-bit: | |
279 | // | |
1e59de90 TL |
280 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
281 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type | |
282 | eval_lsb_imp(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a) | |
7c673cae | 283 | { |
7c673cae FG |
284 | // |
285 | // Find the index of the least significant limb that is non-zero: | |
286 | // | |
1e59de90 | 287 | std::size_t index = 0; |
92f5a8d4 | 288 | while (!a.limbs()[index] && (index < a.size())) |
7c673cae FG |
289 | ++index; |
290 | // | |
291 | // Find the index of the least significant bit within that limb: | |
292 | // | |
1e59de90 | 293 | std::size_t result = boost::multiprecision::detail::find_lsb(a.limbs()[index]); |
7c673cae FG |
294 | |
295 | return result + index * cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
296 | } | |
297 | ||
1e59de90 TL |
298 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
299 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type | |
300 | eval_lsb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a) | |
301 | { | |
302 | using default_ops::eval_get_sign; | |
303 | if (eval_get_sign(a) == 0) | |
304 | { | |
305 | BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand.")); | |
306 | } | |
307 | if (a.sign()) | |
308 | { | |
309 | BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined.")); | |
310 | } | |
311 | return eval_lsb_imp(a); | |
312 | } | |
313 | ||
7c673cae FG |
314 | // |
315 | // Get the location of the most-significant-bit: | |
316 | // | |
1e59de90 TL |
317 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
318 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type | |
7c673cae FG |
319 | eval_msb_imp(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a) |
320 | { | |
321 | // | |
322 | // Find the index of the most significant bit that is non-zero: | |
323 | // | |
324 | return (a.size() - 1) * cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits + boost::multiprecision::detail::find_msb(a.limbs()[a.size() - 1]); | |
325 | } | |
326 | ||
1e59de90 TL |
327 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
328 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type | |
92f5a8d4 | 329 | eval_msb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a) |
7c673cae FG |
330 | { |
331 | using default_ops::eval_get_sign; | |
92f5a8d4 | 332 | if (eval_get_sign(a) == 0) |
7c673cae | 333 | { |
1e59de90 | 334 | BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand.")); |
7c673cae | 335 | } |
92f5a8d4 | 336 | if (a.sign()) |
7c673cae | 337 | { |
1e59de90 | 338 | BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined.")); |
7c673cae FG |
339 | } |
340 | return eval_msb_imp(a); | |
341 | } | |
342 | ||
1e59de90 TL |
343 | #ifdef BOOST_GCC |
344 | // | |
345 | // We really shouldn't need to be disabling this warning, but it really does appear to be | |
346 | // spurious. The warning appears only when in release mode, and asserts are on. | |
347 | // | |
348 | #pragma GCC diagnostic push | |
349 | #pragma GCC diagnostic ignored "-Warray-bounds" | |
350 | #endif | |
351 | ||
352 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> | |
353 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, bool>::type | |
354 | eval_bit_test(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, std::size_t index) noexcept | |
7c673cae | 355 | { |
1e59de90 TL |
356 | std::size_t offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; |
357 | std::size_t shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
92f5a8d4 TL |
358 | limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u); |
359 | if (offset >= val.size()) | |
7c673cae FG |
360 | return false; |
361 | return val.limbs()[offset] & mask ? true : false; | |
362 | } | |
363 | ||
1e59de90 TL |
364 | #ifdef BOOST_GCC |
365 | #pragma GCC diagnostic pop | |
366 | #endif | |
367 | ||
368 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> | |
369 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
370 | eval_bit_set(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, std::size_t index) | |
7c673cae | 371 | { |
1e59de90 TL |
372 | std::size_t offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; |
373 | std::size_t shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
92f5a8d4 TL |
374 | limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u); |
375 | if (offset >= val.size()) | |
7c673cae | 376 | { |
1e59de90 | 377 | std::size_t os = val.size(); |
7c673cae | 378 | val.resize(offset + 1, offset + 1); |
92f5a8d4 TL |
379 | if (offset >= val.size()) |
380 | return; // fixed precision overflow | |
1e59de90 | 381 | for (std::size_t i = os; i <= offset; ++i) |
7c673cae FG |
382 | val.limbs()[i] = 0; |
383 | } | |
384 | val.limbs()[offset] |= mask; | |
385 | } | |
386 | ||
1e59de90 TL |
387 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
388 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
389 | eval_bit_unset(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, std::size_t index) noexcept | |
7c673cae | 390 | { |
1e59de90 TL |
391 | std::size_t offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; |
392 | std::size_t shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
92f5a8d4 TL |
393 | limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u); |
394 | if (offset >= val.size()) | |
7c673cae FG |
395 | return; |
396 | val.limbs()[offset] &= ~mask; | |
397 | val.normalize(); | |
398 | } | |
399 | ||
1e59de90 TL |
400 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
401 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
402 | eval_bit_flip(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, std::size_t index) | |
7c673cae | 403 | { |
1e59de90 TL |
404 | std::size_t offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; |
405 | std::size_t shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits; | |
92f5a8d4 TL |
406 | limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u); |
407 | if (offset >= val.size()) | |
7c673cae | 408 | { |
1e59de90 | 409 | std::size_t os = val.size(); |
7c673cae | 410 | val.resize(offset + 1, offset + 1); |
92f5a8d4 TL |
411 | if (offset >= val.size()) |
412 | return; // fixed precision overflow | |
1e59de90 | 413 | for (std::size_t i = os; i <= offset; ++i) |
7c673cae FG |
414 | val.limbs()[i] = 0; |
415 | } | |
416 | val.limbs()[offset] ^= mask; | |
417 | val.normalize(); | |
418 | } | |
419 | ||
1e59de90 TL |
420 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
421 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
92f5a8d4 TL |
422 | eval_qr( |
423 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x, | |
424 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& y, | |
425 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q, | |
1e59de90 | 426 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value)) |
7c673cae FG |
427 | { |
428 | divide_unsigned_helper(&q, x, y, r); | |
429 | q.sign(x.sign() != y.sign()); | |
430 | r.sign(x.sign()); | |
431 | } | |
432 | ||
1e59de90 TL |
433 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
434 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
92f5a8d4 TL |
435 | eval_qr( |
436 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x, | |
437 | limb_type y, | |
438 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q, | |
1e59de90 | 439 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value)) |
7c673cae FG |
440 | { |
441 | divide_unsigned_helper(&q, x, y, r); | |
442 | q.sign(x.sign()); | |
443 | r.sign(x.sign()); | |
444 | } | |
445 | ||
1e59de90 TL |
446 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class U> |
447 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<U>::value>::type eval_qr( | |
92f5a8d4 TL |
448 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x, |
449 | U y, | |
450 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q, | |
1e59de90 | 451 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value)) |
7c673cae FG |
452 | { |
453 | using default_ops::eval_qr; | |
454 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t; | |
455 | t = y; | |
456 | eval_qr(x, t, q, r); | |
457 | } | |
458 | ||
1e59de90 TL |
459 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer> |
460 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_unsigned<Integer>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, Integer>::type | |
20effc67 | 461 | eval_integer_modulus(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, Integer mod) |
7c673cae | 462 | { |
1e59de90 | 463 | BOOST_IF_CONSTEXPR (sizeof(Integer) <= sizeof(limb_type)) |
7c673cae | 464 | { |
1e59de90 TL |
465 | if (mod <= (std::numeric_limits<limb_type>::max)()) |
466 | { | |
467 | const std::ptrdiff_t n = a.size(); | |
468 | const double_limb_type two_n_mod = static_cast<limb_type>(1u) + (~static_cast<limb_type>(0u) - mod) % mod; | |
469 | limb_type res = a.limbs()[n - 1] % mod; | |
20effc67 | 470 | |
1e59de90 TL |
471 | for (std::ptrdiff_t i = n - 2; i >= 0; --i) |
472 | res = static_cast<limb_type>((res * two_n_mod + a.limbs()[i]) % mod); | |
473 | return res; | |
474 | } | |
475 | else | |
476 | return default_ops::eval_integer_modulus(a, mod); | |
7c673cae FG |
477 | } |
478 | else | |
479 | { | |
20effc67 | 480 | return default_ops::eval_integer_modulus(a, mod); |
7c673cae FG |
481 | } |
482 | } | |
483 | ||
1e59de90 TL |
484 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer> |
485 | BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_signed<Integer>::value && boost::multiprecision::detail::is_integral<Integer>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, Integer>::type | |
92f5a8d4 | 486 | eval_integer_modulus(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x, Integer val) |
7c673cae FG |
487 | { |
488 | return eval_integer_modulus(x, boost::multiprecision::detail::unsigned_abs(val)); | |
489 | } | |
490 | ||
20effc67 | 491 | BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR limb_type eval_gcd(limb_type u, limb_type v) |
7c673cae | 492 | { |
20effc67 TL |
493 | // boundary cases |
494 | if (!u || !v) | |
495 | return u | v; | |
496 | #if __cpp_lib_gcd_lcm >= 201606L | |
497 | return std::gcd(u, v); | |
498 | #else | |
1e59de90 | 499 | std::size_t shift = boost::multiprecision::detail::find_lsb(u | v); |
20effc67 | 500 | u >>= boost::multiprecision::detail::find_lsb(u); |
7c673cae FG |
501 | do |
502 | { | |
20effc67 | 503 | v >>= boost::multiprecision::detail::find_lsb(v); |
92f5a8d4 TL |
504 | if (u > v) |
505 | std_constexpr::swap(u, v); | |
7c673cae | 506 | v -= u; |
20effc67 TL |
507 | } while (v); |
508 | return u << shift; | |
509 | #endif | |
7c673cae FG |
510 | } |
511 | ||
20effc67 | 512 | inline BOOST_MP_CXX14_CONSTEXPR double_limb_type eval_gcd(double_limb_type u, double_limb_type v) |
7c673cae | 513 | { |
20effc67 TL |
514 | #if (__cpp_lib_gcd_lcm >= 201606L) && (!defined(BOOST_HAS_INT128) || !defined(__STRICT_ANSI__)) |
515 | return std::gcd(u, v); | |
516 | #else | |
1e59de90 TL |
517 | if (u == 0) |
518 | return v; | |
519 | ||
520 | std::size_t shift = boost::multiprecision::detail::find_lsb(u | v); | |
20effc67 | 521 | u >>= boost::multiprecision::detail::find_lsb(u); |
7c673cae FG |
522 | do |
523 | { | |
20effc67 | 524 | v >>= boost::multiprecision::detail::find_lsb(v); |
92f5a8d4 TL |
525 | if (u > v) |
526 | std_constexpr::swap(u, v); | |
7c673cae | 527 | v -= u; |
20effc67 TL |
528 | } while (v); |
529 | return u << shift; | |
7c673cae | 530 | #endif |
7c673cae FG |
531 | } |
532 | ||
1e59de90 TL |
533 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
534 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
92f5a8d4 TL |
535 | eval_gcd( |
536 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, | |
537 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, | |
20effc67 | 538 | limb_type b) |
7c673cae | 539 | { |
20effc67 TL |
540 | int s = eval_get_sign(a); |
541 | if (!b || !s) | |
7c673cae | 542 | { |
20effc67 TL |
543 | result = a; |
544 | *result.limbs() |= b; | |
7c673cae | 545 | } |
20effc67 | 546 | else |
7c673cae | 547 | { |
20effc67 TL |
548 | eval_modulus(result, a, b); |
549 | limb_type& res = *result.limbs(); | |
550 | res = eval_gcd(res, b); | |
7c673cae | 551 | } |
20effc67 TL |
552 | result.sign(false); |
553 | } | |
554 | ||
1e59de90 TL |
555 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
556 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
20effc67 TL |
557 | eval_gcd( |
558 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, | |
559 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, | |
560 | double_limb_type b) | |
561 | { | |
562 | int s = eval_get_sign(a); | |
563 | if (!b || !s) | |
7c673cae | 564 | { |
20effc67 TL |
565 | if (!s) |
566 | result = b; | |
567 | else | |
568 | result = a; | |
7c673cae FG |
569 | return; |
570 | } | |
20effc67 TL |
571 | double_limb_type res = 0; |
572 | if(a.sign() == 0) | |
573 | res = eval_integer_modulus(a, b); | |
574 | else | |
7c673cae | 575 | { |
20effc67 TL |
576 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(a); |
577 | t.negate(); | |
578 | res = eval_integer_modulus(t, b); | |
579 | } | |
580 | res = eval_gcd(res, b); | |
581 | result = res; | |
582 | result.sign(false); | |
7c673cae | 583 | } |
1e59de90 TL |
584 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
585 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
20effc67 TL |
586 | eval_gcd( |
587 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, | |
588 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, | |
589 | signed_double_limb_type v) | |
590 | { | |
591 | eval_gcd(result, a, static_cast<double_limb_type>(v < 0 ? -v : v)); | |
592 | } | |
593 | // | |
594 | // These 2 overloads take care of gcd against an (unsigned) short etc: | |
595 | // | |
1e59de90 TL |
596 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer> |
597 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_unsigned<Integer>::value && (sizeof(Integer) <= sizeof(limb_type)) && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
92f5a8d4 TL |
598 | eval_gcd( |
599 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, | |
600 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, | |
601 | const Integer& v) | |
7c673cae FG |
602 | { |
603 | eval_gcd(result, a, static_cast<limb_type>(v)); | |
604 | } | |
1e59de90 TL |
605 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer> |
606 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_signed<Integer>::value && boost::multiprecision::detail::is_integral<Integer>::value && (sizeof(Integer) <= sizeof(limb_type)) && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
92f5a8d4 TL |
607 | eval_gcd( |
608 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, | |
609 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, | |
610 | const Integer& v) | |
7c673cae FG |
611 | { |
612 | eval_gcd(result, a, static_cast<limb_type>(v < 0 ? -v : v)); | |
613 | } | |
20effc67 TL |
614 | // |
615 | // What follows is Lehmer's GCD algorithm: | |
616 | // Essentially this uses the leading digit(s) of U and V | |
617 | // only to run a "simulated" Euclid algorithm. It stops | |
618 | // when the calculated quotient differs from what would have been | |
619 | // the true quotient. At that point the cosequences are used to | |
620 | // calculate the new U and V. A nice lucid description appears | |
621 | // in "An Analysis of Lehmer's Euclidean GCD Algorithm", | |
622 | // by Jonathan Sorenson. https://www.researchgate.net/publication/2424634_An_Analysis_of_Lehmer%27s_Euclidean_GCD_Algorithm | |
623 | // DOI: 10.1145/220346.220378. | |
624 | // | |
625 | // There are two versions of this algorithm here, and both are "double digit" | |
626 | // variations: which is to say if there are k bits per limb, then they extract | |
627 | // 2k bits into a double_limb_type and then run the algorithm on that. The first | |
628 | // version is a straightforward version of the algorithm, and is designed for | |
629 | // situations where double_limb_type is a native integer (for example where | |
630 | // limb_type is a 32-bit integer on a 64-bit machine). For 32-bit limbs it | |
631 | // reduces the size of U by about 30 bits per call. The second is a more complex | |
632 | // version for situations where double_limb_type is a synthetic type: for example | |
633 | // __int128. For 64 bit limbs it reduces the size of U by about 62 bits per call. | |
634 | // | |
635 | // The complexity of the algorithm given by Sorenson is roughly O(ln^2(N)) for | |
636 | // two N bit numbers. | |
637 | // | |
638 | // The original double-digit version of the algorithm is described in: | |
639 | // | |
640 | // "A Double Digit Lehmer-Euclid Algorithm for Finding the GCD of Long Integers", | |
641 | // Tudor Jebelean, J Symbolic Computation, 1995 (19), 145. | |
642 | // | |
643 | #ifndef BOOST_HAS_INT128 | |
644 | // | |
645 | // When double_limb_type is a native integer type then we should just use it and not worry about the consequences. | |
646 | // This can eliminate approximately a full limb with each call. | |
647 | // | |
1e59de90 TL |
648 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Storage> |
649 | void eval_gcd_lehmer(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& U, cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& V, std::size_t lu, Storage& storage) | |
20effc67 TL |
650 | { |
651 | // | |
652 | // Extract the leading 2 * bits_per_limb bits from U and V: | |
653 | // | |
1e59de90 | 654 | std::size_t h = lu % bits_per_limb; |
20effc67 TL |
655 | double_limb_type u = (static_cast<double_limb_type>((U.limbs()[U.size() - 1])) << bits_per_limb) | U.limbs()[U.size() - 2]; |
656 | double_limb_type v = (static_cast<double_limb_type>((V.size() < U.size() ? 0 : V.limbs()[V.size() - 1])) << bits_per_limb) | V.limbs()[U.size() - 2]; | |
657 | if (h) | |
658 | { | |
659 | u <<= bits_per_limb - h; | |
660 | u |= U.limbs()[U.size() - 3] >> h; | |
661 | v <<= bits_per_limb - h; | |
662 | v |= V.limbs()[U.size() - 3] >> h; | |
663 | } | |
664 | // | |
665 | // Co-sequences x an y: we need only the last 3 values of these, | |
666 | // the first 2 values are known correct, the third gets checked | |
667 | // in each loop operation, and we terminate when they go wrong. | |
668 | // | |
669 | // x[i+0] is positive for even i. | |
670 | // y[i+0] is positive for odd i. | |
671 | // | |
672 | // However we track only absolute values here: | |
673 | // | |
674 | double_limb_type x[3] = {1, 0}; | |
675 | double_limb_type y[3] = {0, 1}; | |
1e59de90 | 676 | std::size_t i = 0; |
20effc67 TL |
677 | |
678 | #ifdef BOOST_MP_GCD_DEBUG | |
679 | cpp_int UU, VV; | |
680 | UU = U; | |
681 | VV = V; | |
682 | #endif | |
683 | ||
684 | while (true) | |
685 | { | |
686 | double_limb_type q = u / v; | |
687 | x[2] = x[0] + q * x[1]; | |
688 | y[2] = y[0] + q * y[1]; | |
689 | double_limb_type tu = u; | |
690 | u = v; | |
691 | v = tu - q * v; | |
692 | ++i; | |
693 | // | |
694 | // We must make sure that y[2] occupies a single limb otherwise | |
695 | // the multiprecision multiplications below would be much more expensive. | |
696 | // This can sometimes lose us one iteration, but is worth it for improved | |
697 | // calculation efficiency. | |
698 | // | |
699 | if (y[2] >> bits_per_limb) | |
700 | break; | |
701 | // | |
702 | // These are Jebelean's exact termination conditions: | |
703 | // | |
704 | if ((i & 1u) == 0) | |
705 | { | |
1e59de90 | 706 | BOOST_MP_ASSERT(u > v); |
20effc67 TL |
707 | if ((v < x[2]) || ((u - v) < (y[2] + y[1]))) |
708 | break; | |
709 | } | |
710 | else | |
711 | { | |
1e59de90 | 712 | BOOST_MP_ASSERT(u > v); |
20effc67 TL |
713 | if ((v < y[2]) || ((u - v) < (x[2] + x[1]))) |
714 | break; | |
715 | } | |
716 | #ifdef BOOST_MP_GCD_DEBUG | |
1e59de90 | 717 | BOOST_MP_ASSERT(q == UU / VV); |
20effc67 TL |
718 | UU %= VV; |
719 | UU.swap(VV); | |
720 | #endif | |
721 | x[0] = x[1]; | |
722 | x[1] = x[2]; | |
723 | y[0] = y[1]; | |
724 | y[1] = y[2]; | |
725 | } | |
726 | if (i == 1) | |
727 | { | |
728 | // No change to U and V we've stalled! | |
729 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t; | |
730 | eval_modulus(t, U, V); | |
731 | U.swap(V); | |
732 | V.swap(t); | |
733 | return; | |
734 | } | |
735 | // | |
736 | // Update U and V. | |
737 | // We have: | |
738 | // | |
739 | // U = x[0]U + y[0]V and | |
740 | // V = x[1]U + y[1]V. | |
741 | // | |
742 | // But since we track only absolute values of x and y | |
743 | // we have to take account of the implied signs and perform | |
744 | // the appropriate subtraction depending on the whether i is | |
745 | // even or odd: | |
746 | // | |
1e59de90 | 747 | std::size_t ts = U.size() + 1; |
20effc67 TL |
748 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t1(storage, ts), t2(storage, ts), t3(storage, ts); |
749 | eval_multiply(t1, U, static_cast<limb_type>(x[0])); | |
750 | eval_multiply(t2, V, static_cast<limb_type>(y[0])); | |
751 | eval_multiply(t3, U, static_cast<limb_type>(x[1])); | |
752 | if ((i & 1u) == 0) | |
753 | { | |
754 | if (x[0] == 0) | |
755 | U = t2; | |
756 | else | |
757 | { | |
1e59de90 | 758 | BOOST_MP_ASSERT(t2.compare(t1) >= 0); |
20effc67 | 759 | eval_subtract(U, t2, t1); |
1e59de90 | 760 | BOOST_MP_ASSERT(U.sign() == false); |
20effc67 TL |
761 | } |
762 | } | |
763 | else | |
764 | { | |
1e59de90 | 765 | BOOST_MP_ASSERT(t1.compare(t2) >= 0); |
20effc67 | 766 | eval_subtract(U, t1, t2); |
1e59de90 | 767 | BOOST_MP_ASSERT(U.sign() == false); |
20effc67 TL |
768 | } |
769 | eval_multiply(t2, V, static_cast<limb_type>(y[1])); | |
770 | if (i & 1u) | |
771 | { | |
772 | if (x[1] == 0) | |
773 | V = t2; | |
774 | else | |
775 | { | |
1e59de90 | 776 | BOOST_MP_ASSERT(t2.compare(t3) >= 0); |
20effc67 | 777 | eval_subtract(V, t2, t3); |
1e59de90 | 778 | BOOST_MP_ASSERT(V.sign() == false); |
20effc67 TL |
779 | } |
780 | } | |
781 | else | |
782 | { | |
1e59de90 | 783 | BOOST_MP_ASSERT(t3.compare(t2) >= 0); |
20effc67 | 784 | eval_subtract(V, t3, t2); |
1e59de90 | 785 | BOOST_MP_ASSERT(V.sign() == false); |
20effc67 | 786 | } |
1e59de90 TL |
787 | BOOST_MP_ASSERT(U.compare(V) >= 0); |
788 | BOOST_MP_ASSERT(lu > eval_msb(U)); | |
20effc67 TL |
789 | #ifdef BOOST_MP_GCD_DEBUG |
790 | ||
1e59de90 TL |
791 | BOOST_MP_ASSERT(UU == U); |
792 | BOOST_MP_ASSERT(VV == V); | |
20effc67 | 793 | |
1e59de90 TL |
794 | extern std::size_t total_lehmer_gcd_calls; |
795 | extern std::size_t total_lehmer_gcd_bits_saved; | |
796 | extern std::size_t total_lehmer_gcd_cycles; | |
20effc67 TL |
797 | |
798 | ++total_lehmer_gcd_calls; | |
799 | total_lehmer_gcd_bits_saved += lu - eval_msb(U); | |
800 | total_lehmer_gcd_cycles += i; | |
801 | #endif | |
802 | if (lu < 2048) | |
803 | { | |
804 | // | |
805 | // Since we have stripped all common powers of 2 from U and V at the start | |
806 | // if either are even at this point, we can remove stray powers of 2 now. | |
807 | // Note that it is not possible for *both* U and V to be even at this point. | |
808 | // | |
809 | // This has an adverse effect on performance for high bit counts, but has | |
810 | // a significant positive effect for smaller counts. | |
811 | // | |
812 | if ((U.limbs()[0] & 1u) == 0) | |
813 | { | |
814 | eval_right_shift(U, eval_lsb(U)); | |
815 | if (U.compare(V) < 0) | |
816 | U.swap(V); | |
817 | } | |
818 | else if ((V.limbs()[0] & 1u) == 0) | |
819 | { | |
820 | eval_right_shift(V, eval_lsb(V)); | |
821 | } | |
822 | } | |
823 | storage.deallocate(ts * 3); | |
824 | } | |
825 | ||
826 | #else | |
827 | // | |
828 | // This branch is taken when double_limb_type is a synthetic type with no native hardware support. | |
829 | // For example __int128. The assumption is that add/subtract/multiply of double_limb_type are efficient, | |
830 | // but that division is very slow. | |
831 | // | |
832 | // We begin with a specialized routine for division. | |
1e59de90 | 833 | // We know that most of the time this is called the result will be 1. |
20effc67 TL |
834 | // For small limb counts, this almost doubles the performance of Lehmer's routine! |
835 | // | |
1e59de90 | 836 | BOOST_FORCEINLINE void divide_subtract(double_limb_type& q, double_limb_type& u, const double_limb_type& v) |
20effc67 | 837 | { |
1e59de90 | 838 | BOOST_MP_ASSERT(q == 1); // precondition on entry. |
20effc67 TL |
839 | u -= v; |
840 | while (u >= v) | |
841 | { | |
842 | u -= v; | |
843 | if (++q > 30) | |
844 | { | |
1e59de90 | 845 | double_limb_type t = u / v; |
20effc67 TL |
846 | u -= t * v; |
847 | q += t; | |
848 | } | |
849 | } | |
850 | } | |
851 | ||
1e59de90 TL |
852 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Storage> |
853 | void eval_gcd_lehmer(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& U, cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& V, std::size_t lu, Storage& storage) | |
20effc67 TL |
854 | { |
855 | // | |
856 | // Extract the leading 2*bits_per_limb bits from U and V: | |
857 | // | |
1e59de90 | 858 | std::size_t h = lu % bits_per_limb; |
20effc67 TL |
859 | double_limb_type u, v; |
860 | if (h) | |
861 | { | |
862 | u = (static_cast<double_limb_type>((U.limbs()[U.size() - 1])) << bits_per_limb) | U.limbs()[U.size() - 2]; | |
863 | v = (static_cast<double_limb_type>((V.size() < U.size() ? 0 : V.limbs()[V.size() - 1])) << bits_per_limb) | V.limbs()[U.size() - 2]; | |
864 | u <<= bits_per_limb - h; | |
865 | u |= U.limbs()[U.size() - 3] >> h; | |
866 | v <<= bits_per_limb - h; | |
867 | v |= V.limbs()[U.size() - 3] >> h; | |
868 | } | |
869 | else | |
870 | { | |
871 | u = (static_cast<double_limb_type>(U.limbs()[U.size() - 1]) << bits_per_limb) | U.limbs()[U.size() - 2]; | |
872 | v = (static_cast<double_limb_type>(V.limbs()[U.size() - 1]) << bits_per_limb) | V.limbs()[U.size() - 2]; | |
873 | } | |
874 | // | |
875 | // Cosequences are stored as limb_types, we take care not to overflow these: | |
876 | // | |
877 | // x[i+0] is positive for even i. | |
878 | // y[i+0] is positive for odd i. | |
879 | // | |
880 | // However we track only absolute values here: | |
881 | // | |
882 | limb_type x[3] = { 1, 0 }; | |
883 | limb_type y[3] = { 0, 1 }; | |
1e59de90 | 884 | std::size_t i = 0; |
20effc67 TL |
885 | |
886 | #ifdef BOOST_MP_GCD_DEBUG | |
887 | cpp_int UU, VV; | |
888 | UU = U; | |
889 | VV = V; | |
890 | #endif | |
891 | // | |
892 | // We begine by running a single digit version of Lehmer's algorithm, we still have | |
893 | // to track u and v at double precision, but this adds only a tiny performance penalty. | |
894 | // What we gain is fast division, and fast termination testing. | |
895 | // When you see static_cast<limb_type>(u >> bits_per_limb) here, this is really just | |
896 | // a direct access to the upper bits_per_limb of the double limb type. For __int128 | |
897 | // this is simple a load of the upper 64 bits and the "shift" is optimised away. | |
898 | // | |
899 | double_limb_type old_u, old_v; | |
900 | while (true) | |
901 | { | |
902 | limb_type q = static_cast<limb_type>(u >> bits_per_limb) / static_cast<limb_type>(v >> bits_per_limb); | |
903 | x[2] = x[0] + q * x[1]; | |
904 | y[2] = y[0] + q * y[1]; | |
905 | double_limb_type tu = u; | |
906 | old_u = u; | |
907 | old_v = v; | |
908 | u = v; | |
909 | double_limb_type t = q * v; | |
910 | if (tu < t) | |
911 | { | |
912 | ++i; | |
913 | break; | |
914 | } | |
915 | v = tu - t; | |
916 | ++i; | |
1e59de90 TL |
917 | BOOST_MP_ASSERT((u <= v) || (t / q == old_v)); |
918 | if (u <= v) | |
919 | { | |
920 | // We've gone terribly wrong, probably numeric overflow: | |
921 | break; | |
922 | } | |
20effc67 TL |
923 | if ((i & 1u) == 0) |
924 | { | |
20effc67 TL |
925 | if ((static_cast<limb_type>(v >> bits_per_limb) < x[2]) || ((static_cast<limb_type>(u >> bits_per_limb) - static_cast<limb_type>(v >> bits_per_limb)) < (y[2] + y[1]))) |
926 | break; | |
927 | } | |
928 | else | |
929 | { | |
20effc67 TL |
930 | if ((static_cast<limb_type>(v >> bits_per_limb) < y[2]) || ((static_cast<limb_type>(u >> bits_per_limb) - static_cast<limb_type>(v >> bits_per_limb)) < (x[2] + x[1]))) |
931 | break; | |
932 | } | |
933 | #ifdef BOOST_MP_GCD_DEBUG | |
1e59de90 | 934 | BOOST_MP_ASSERT(q == UU / VV); |
20effc67 TL |
935 | UU %= VV; |
936 | UU.swap(VV); | |
937 | #endif | |
938 | x[0] = x[1]; | |
939 | x[1] = x[2]; | |
940 | y[0] = y[1]; | |
941 | y[1] = y[2]; | |
942 | } | |
943 | // | |
944 | // We get here when the single digit algorithm has gone wrong, back up i, u and v: | |
945 | // | |
946 | --i; | |
947 | u = old_u; | |
948 | v = old_v; | |
949 | // | |
950 | // Now run the full double-digit algorithm: | |
951 | // | |
952 | while (true) | |
953 | { | |
1e59de90 | 954 | double_limb_type q = 1; |
20effc67 TL |
955 | double_limb_type tt = u; |
956 | divide_subtract(q, u, v); | |
957 | std::swap(u, v); | |
958 | tt = y[0] + q * static_cast<double_limb_type>(y[1]); | |
959 | // | |
960 | // If calculation of y[2] would overflow a single limb, then we *must* terminate. | |
961 | // Note that x[2] < y[2] so there is no need to check that as well: | |
962 | // | |
963 | if (tt >> bits_per_limb) | |
964 | { | |
965 | ++i; | |
966 | break; | |
967 | } | |
968 | x[2] = x[0] + q * x[1]; | |
969 | y[2] = tt; | |
970 | ++i; | |
971 | if ((i & 1u) == 0) | |
972 | { | |
1e59de90 | 973 | BOOST_MP_ASSERT(u > v); |
20effc67 TL |
974 | if ((v < x[2]) || ((u - v) < (static_cast<double_limb_type>(y[2]) + y[1]))) |
975 | break; | |
976 | } | |
977 | else | |
978 | { | |
1e59de90 | 979 | BOOST_MP_ASSERT(u > v); |
20effc67 TL |
980 | if ((v < y[2]) || ((u - v) < (static_cast<double_limb_type>(x[2]) + x[1]))) |
981 | break; | |
982 | } | |
983 | #ifdef BOOST_MP_GCD_DEBUG | |
1e59de90 | 984 | BOOST_MP_ASSERT(q == UU / VV); |
20effc67 TL |
985 | UU %= VV; |
986 | UU.swap(VV); | |
987 | #endif | |
988 | x[0] = x[1]; | |
989 | x[1] = x[2]; | |
990 | y[0] = y[1]; | |
991 | y[1] = y[2]; | |
992 | } | |
993 | if (i == 1) | |
994 | { | |
995 | // No change to U and V we've stalled! | |
996 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t; | |
997 | eval_modulus(t, U, V); | |
998 | U.swap(V); | |
999 | V.swap(t); | |
1000 | return; | |
1001 | } | |
1002 | // | |
1003 | // Update U and V. | |
1004 | // We have: | |
1005 | // | |
1006 | // U = x[0]U + y[0]V and | |
1007 | // V = x[1]U + y[1]V. | |
1008 | // | |
1009 | // But since we track only absolute values of x and y | |
1010 | // we have to take account of the implied signs and perform | |
1011 | // the appropriate subtraction depending on the whether i is | |
1012 | // even or odd: | |
1013 | // | |
1e59de90 | 1014 | std::size_t ts = U.size() + 1; |
20effc67 TL |
1015 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t1(storage, ts), t2(storage, ts), t3(storage, ts); |
1016 | eval_multiply(t1, U, x[0]); | |
1017 | eval_multiply(t2, V, y[0]); | |
1018 | eval_multiply(t3, U, x[1]); | |
1019 | if ((i & 1u) == 0) | |
1020 | { | |
1021 | if (x[0] == 0) | |
1022 | U = t2; | |
1023 | else | |
1024 | { | |
1e59de90 | 1025 | BOOST_MP_ASSERT(t2.compare(t1) >= 0); |
20effc67 | 1026 | eval_subtract(U, t2, t1); |
1e59de90 | 1027 | BOOST_MP_ASSERT(U.sign() == false); |
20effc67 TL |
1028 | } |
1029 | } | |
1030 | else | |
1031 | { | |
1e59de90 | 1032 | BOOST_MP_ASSERT(t1.compare(t2) >= 0); |
20effc67 | 1033 | eval_subtract(U, t1, t2); |
1e59de90 | 1034 | BOOST_MP_ASSERT(U.sign() == false); |
20effc67 TL |
1035 | } |
1036 | eval_multiply(t2, V, y[1]); | |
1037 | if (i & 1u) | |
1038 | { | |
1039 | if (x[1] == 0) | |
1040 | V = t2; | |
1041 | else | |
1042 | { | |
1e59de90 | 1043 | BOOST_MP_ASSERT(t2.compare(t3) >= 0); |
20effc67 | 1044 | eval_subtract(V, t2, t3); |
1e59de90 | 1045 | BOOST_MP_ASSERT(V.sign() == false); |
20effc67 TL |
1046 | } |
1047 | } | |
1048 | else | |
1049 | { | |
1e59de90 | 1050 | BOOST_MP_ASSERT(t3.compare(t2) >= 0); |
20effc67 | 1051 | eval_subtract(V, t3, t2); |
1e59de90 | 1052 | BOOST_MP_ASSERT(V.sign() == false); |
20effc67 | 1053 | } |
1e59de90 TL |
1054 | BOOST_MP_ASSERT(U.compare(V) >= 0); |
1055 | BOOST_MP_ASSERT(lu > eval_msb(U)); | |
20effc67 TL |
1056 | #ifdef BOOST_MP_GCD_DEBUG |
1057 | ||
1e59de90 TL |
1058 | BOOST_MP_ASSERT(UU == U); |
1059 | BOOST_MP_ASSERT(VV == V); | |
20effc67 | 1060 | |
1e59de90 TL |
1061 | extern std::size_t total_lehmer_gcd_calls; |
1062 | extern std::size_t total_lehmer_gcd_bits_saved; | |
1063 | extern std::size_t total_lehmer_gcd_cycles; | |
20effc67 TL |
1064 | |
1065 | ++total_lehmer_gcd_calls; | |
1066 | total_lehmer_gcd_bits_saved += lu - eval_msb(U); | |
1067 | total_lehmer_gcd_cycles += i; | |
1068 | #endif | |
1069 | if (lu < 2048) | |
1070 | { | |
1071 | // | |
1072 | // Since we have stripped all common powers of 2 from U and V at the start | |
1073 | // if either are even at this point, we can remove stray powers of 2 now. | |
1074 | // Note that it is not possible for *both* U and V to be even at this point. | |
1075 | // | |
1076 | // This has an adverse effect on performance for high bit counts, but has | |
1077 | // a significant positive effect for smaller counts. | |
1078 | // | |
1079 | if ((U.limbs()[0] & 1u) == 0) | |
1080 | { | |
1081 | eval_right_shift(U, eval_lsb(U)); | |
1082 | if (U.compare(V) < 0) | |
1083 | U.swap(V); | |
1084 | } | |
1085 | else if ((V.limbs()[0] & 1u) == 0) | |
1086 | { | |
1087 | eval_right_shift(V, eval_lsb(V)); | |
1088 | } | |
1089 | } | |
1090 | storage.deallocate(ts * 3); | |
1091 | } | |
1092 | ||
1093 | #endif | |
7c673cae | 1094 | |
1e59de90 TL |
1095 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
1096 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
92f5a8d4 | 1097 | eval_gcd( |
20effc67 TL |
1098 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, |
1099 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, | |
1100 | const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b) | |
7c673cae | 1101 | { |
7c673cae | 1102 | using default_ops::eval_get_sign; |
92f5a8d4 TL |
1103 | using default_ops::eval_is_zero; |
1104 | using default_ops::eval_lsb; | |
7c673cae | 1105 | |
92f5a8d4 | 1106 | if (a.size() == 1) |
7c673cae FG |
1107 | { |
1108 | eval_gcd(result, b, *a.limbs()); | |
1109 | return; | |
1110 | } | |
92f5a8d4 | 1111 | if (b.size() == 1) |
7c673cae FG |
1112 | { |
1113 | eval_gcd(result, a, *b.limbs()); | |
1114 | return; | |
1115 | } | |
1e59de90 | 1116 | std::size_t temp_size = (std::max)(a.size(), b.size()) + 1; |
20effc67 | 1117 | typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::scoped_shared_storage storage(a, temp_size * 6); |
7c673cae | 1118 | |
20effc67 TL |
1119 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> U(storage, temp_size); |
1120 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> V(storage, temp_size); | |
1121 | cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(storage, temp_size); | |
1122 | U = a; | |
1123 | V = b; | |
7c673cae | 1124 | |
20effc67 | 1125 | int s = eval_get_sign(U); |
7c673cae FG |
1126 | |
1127 | /* GCD(0,x) := x */ | |
92f5a8d4 | 1128 | if (s < 0) |
7c673cae | 1129 | { |
20effc67 | 1130 | U.negate(); |
7c673cae | 1131 | } |
92f5a8d4 | 1132 | else if (s == 0) |
7c673cae | 1133 | { |
20effc67 | 1134 | result = V; |
7c673cae FG |
1135 | return; |
1136 | } | |
20effc67 | 1137 | s = eval_get_sign(V); |
92f5a8d4 | 1138 | if (s < 0) |
7c673cae | 1139 | { |
20effc67 | 1140 | V.negate(); |
7c673cae | 1141 | } |
92f5a8d4 | 1142 | else if (s == 0) |
7c673cae | 1143 | { |
20effc67 | 1144 | result = U; |
7c673cae FG |
1145 | return; |
1146 | } | |
20effc67 TL |
1147 | // |
1148 | // Remove common factors of 2: | |
1149 | // | |
1e59de90 TL |
1150 | std::size_t us = eval_lsb(U); |
1151 | std::size_t vs = eval_lsb(V); | |
1152 | std::size_t shift = (std::min)(us, vs); | |
20effc67 TL |
1153 | if (us) |
1154 | eval_right_shift(U, us); | |
1155 | if (vs) | |
1156 | eval_right_shift(V, vs); | |
7c673cae | 1157 | |
20effc67 TL |
1158 | if (U.compare(V) < 0) |
1159 | U.swap(V); | |
7c673cae | 1160 | |
20effc67 | 1161 | while (!eval_is_zero(V)) |
7c673cae | 1162 | { |
20effc67 | 1163 | if (U.size() <= 2) |
7c673cae | 1164 | { |
20effc67 TL |
1165 | // |
1166 | // Special case: if V has no more than 2 limbs | |
1167 | // then we can reduce U and V to a pair of integers and perform | |
1168 | // direct integer gcd: | |
1169 | // | |
1170 | if (U.size() == 1) | |
1171 | U = eval_gcd(*V.limbs(), *U.limbs()); | |
7c673cae FG |
1172 | else |
1173 | { | |
20effc67 TL |
1174 | double_limb_type i = U.limbs()[0] | (static_cast<double_limb_type>(U.limbs()[1]) << sizeof(limb_type) * CHAR_BIT); |
1175 | double_limb_type j = (V.size() == 1) ? *V.limbs() : V.limbs()[0] | (static_cast<double_limb_type>(V.limbs()[1]) << sizeof(limb_type) * CHAR_BIT); | |
1176 | U = eval_gcd(i, j); | |
7c673cae FG |
1177 | } |
1178 | break; | |
1179 | } | |
1e59de90 TL |
1180 | std::size_t lu = eval_msb(U) + 1; |
1181 | std::size_t lv = eval_msb(V) + 1; | |
20effc67 TL |
1182 | #ifndef BOOST_MP_NO_CONSTEXPR_DETECTION |
1183 | if (!BOOST_MP_IS_CONST_EVALUATED(lu) && (lu - lv <= bits_per_limb / 2)) | |
1184 | #else | |
1185 | if (lu - lv <= bits_per_limb / 2) | |
1186 | #endif | |
1187 | { | |
1188 | eval_gcd_lehmer(U, V, lu, storage); | |
1189 | } | |
1190 | else | |
1191 | { | |
1192 | eval_modulus(t, U, V); | |
1193 | U.swap(V); | |
1194 | V.swap(t); | |
1195 | } | |
1196 | } | |
1197 | result = U; | |
1198 | if (shift) | |
1199 | eval_left_shift(result, shift); | |
7c673cae FG |
1200 | } |
1201 | // | |
1202 | // Now again for trivial backends: | |
1203 | // | |
1e59de90 TL |
1204 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
1205 | BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type | |
1206 | eval_gcd(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b) noexcept | |
7c673cae | 1207 | { |
1e59de90 TL |
1208 | *result.limbs() = boost::multiprecision::detail::constexpr_gcd(*a.limbs(), *b.limbs()); |
1209 | result.sign(false); | |
7c673cae FG |
1210 | } |
1211 | // This one is only enabled for unchecked cpp_int's, for checked int's we need the checking in the default version: | |
1e59de90 TL |
1212 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
1213 | BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && (Checked1 == unchecked)>::type | |
1214 | eval_lcm(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value)) | |
7c673cae | 1215 | { |
1e59de90 | 1216 | *result.limbs() = boost::multiprecision::detail::constexpr_lcm(*a.limbs(), *b.limbs()); |
7c673cae | 1217 | result.normalize(); // result may overflow the specified number of bits |
1e59de90 | 1218 | result.sign(false); |
7c673cae FG |
1219 | } |
1220 | ||
1e59de90 TL |
1221 | inline void conversion_overflow(const std::integral_constant<int, checked>&) |
1222 | { | |
1223 | BOOST_MP_THROW_EXCEPTION(std::overflow_error("Overflow in conversion to narrower type")); | |
1224 | } | |
1225 | inline BOOST_MP_CXX14_CONSTEXPR void conversion_overflow(const std::integral_constant<int, unchecked>&) {} | |
1226 | ||
1227 | #if defined(__clang__) && defined(__MINGW32__) | |
1228 | // | |
1229 | // clang-11 on Mingw segfaults on conversion of __int128 -> float. | |
1230 | // See: https://bugs.llvm.org/show_bug.cgi?id=48941 | |
1231 | // These workarounds pass everything through an intermediate uint64_t. | |
1232 | // | |
1233 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> | |
1234 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if< | |
1235 | is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && std::is_same<typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type, double_limb_type>::value>::type | |
1236 | eval_convert_to(float* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) | |
1237 | { | |
1238 | float f = static_cast<std::uint64_t>((*val.limbs()) >> 64); | |
1239 | *result = std::ldexp(f, 64); | |
1240 | *result += static_cast<std::uint64_t>((*val.limbs())); | |
1241 | if(val.sign()) | |
1242 | *result = -*result; | |
1243 | } | |
1244 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> | |
1245 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if< | |
1246 | is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && std::is_same<typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type, double_limb_type>::value>::type | |
1247 | eval_convert_to(double* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) | |
1248 | { | |
1249 | float f = static_cast<std::uint64_t>((*val.limbs()) >> 64); | |
1250 | *result = std::ldexp(f, 64); | |
1251 | *result += static_cast<std::uint64_t>((*val.limbs())); | |
1252 | if(val.sign()) | |
1253 | *result = -*result; | |
1254 | } | |
1255 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> | |
1256 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if< | |
1257 | is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && std::is_same<typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type, double_limb_type>::value>::type | |
1258 | eval_convert_to(long double* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) | |
7c673cae | 1259 | { |
1e59de90 TL |
1260 | float f = static_cast<std::uint64_t>((*val.limbs()) >> 64); |
1261 | *result = std::ldexp(f, 64); | |
1262 | *result += static_cast<std::uint64_t>((*val.limbs())); | |
1263 | if(val.sign()) | |
1264 | *result = -*result; | |
7c673cae | 1265 | } |
1e59de90 | 1266 | #endif |
7c673cae | 1267 | |
1e59de90 TL |
1268 | template <class R, std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
1269 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if< | |
1270 | is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && std::is_convertible<typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type, R>::value>::type | |
92f5a8d4 | 1271 | eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) |
7c673cae | 1272 | { |
1e59de90 TL |
1273 | using common_type = typename std::common_type<R, typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type>::type; |
1274 | BOOST_IF_CONSTEXPR(std::numeric_limits<R>::is_specialized) | |
7c673cae | 1275 | { |
1e59de90 | 1276 | if (static_cast<common_type>(*val.limbs()) > static_cast<common_type>((std::numeric_limits<R>::max)())) |
7c673cae | 1277 | { |
1e59de90 TL |
1278 | if (val.isneg()) |
1279 | { | |
1280 | check_is_negative(std::integral_constant < bool, (boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value) || (number_category<R>::value == number_kind_floating_point) > ()); | |
1281 | if (static_cast<common_type>(*val.limbs()) > -static_cast<common_type>((std::numeric_limits<R>::min)())) | |
1282 | conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type()); | |
1283 | *result = (std::numeric_limits<R>::min)(); | |
1284 | } | |
1285 | else | |
1286 | { | |
7c673cae | 1287 | conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type()); |
1e59de90 TL |
1288 | *result = boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value ? (std::numeric_limits<R>::max)() : static_cast<R>(*val.limbs()); |
1289 | } | |
7c673cae FG |
1290 | } |
1291 | else | |
1292 | { | |
1e59de90 TL |
1293 | *result = static_cast<R>(*val.limbs()); |
1294 | if (val.isneg()) | |
1295 | { | |
1296 | check_is_negative(std::integral_constant < bool, (boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value) || (number_category<R>::value == number_kind_floating_point) > ()); | |
1297 | *result = negate_integer(*result, std::integral_constant < bool, is_signed_number<R>::value || (number_category<R>::value == number_kind_floating_point) > ()); | |
1298 | } | |
7c673cae FG |
1299 | } |
1300 | } | |
1301 | else | |
1302 | { | |
1303 | *result = static_cast<R>(*val.limbs()); | |
92f5a8d4 | 1304 | if (val.isneg()) |
7c673cae | 1305 | { |
1e59de90 TL |
1306 | check_is_negative(std::integral_constant<bool, (boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value) || (number_category<R>::value == number_kind_floating_point) > ()); |
1307 | *result = negate_integer(*result, std::integral_constant<bool, is_signed_number<R>::value || (number_category<R>::value == number_kind_floating_point) > ()); | |
7c673cae FG |
1308 | } |
1309 | } | |
1310 | } | |
1311 | ||
1e59de90 TL |
1312 | template <class R, std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
1313 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if< | |
1314 | is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_unsigned_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && std::is_convertible<typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type, R>::value>::type | |
92f5a8d4 | 1315 | eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) |
7c673cae | 1316 | { |
1e59de90 TL |
1317 | using common_type = typename std::common_type<R, typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type>::type; |
1318 | BOOST_IF_CONSTEXPR(std::numeric_limits<R>::is_specialized) | |
7c673cae | 1319 | { |
1e59de90 TL |
1320 | if(static_cast<common_type>(*val.limbs()) > static_cast<common_type>((std::numeric_limits<R>::max)())) |
1321 | { | |
1322 | conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type()); | |
1323 | *result = boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value ? (std::numeric_limits<R>::max)() : static_cast<R>(*val.limbs()); | |
1324 | } | |
1325 | else | |
1326 | *result = static_cast<R>(*val.limbs()); | |
7c673cae FG |
1327 | } |
1328 | else | |
1329 | *result = static_cast<R>(*val.limbs()); | |
1330 | } | |
1331 | ||
1e59de90 TL |
1332 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
1333 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type | |
92f5a8d4 | 1334 | eval_lsb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a) |
7c673cae FG |
1335 | { |
1336 | using default_ops::eval_get_sign; | |
92f5a8d4 | 1337 | if (eval_get_sign(a) == 0) |
7c673cae | 1338 | { |
1e59de90 | 1339 | BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand.")); |
7c673cae | 1340 | } |
92f5a8d4 | 1341 | if (a.sign()) |
7c673cae | 1342 | { |
1e59de90 | 1343 | BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined.")); |
7c673cae FG |
1344 | } |
1345 | // | |
1346 | // Find the index of the least significant bit within that limb: | |
1347 | // | |
1348 | return boost::multiprecision::detail::find_lsb(*a.limbs()); | |
1349 | } | |
1350 | ||
1e59de90 TL |
1351 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
1352 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type | |
7c673cae FG |
1353 | eval_msb_imp(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a) |
1354 | { | |
1355 | // | |
1356 | // Find the index of the least significant bit within that limb: | |
1357 | // | |
1358 | return boost::multiprecision::detail::find_msb(*a.limbs()); | |
1359 | } | |
1360 | ||
1e59de90 TL |
1361 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
1362 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type | |
92f5a8d4 | 1363 | eval_msb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a) |
7c673cae FG |
1364 | { |
1365 | using default_ops::eval_get_sign; | |
92f5a8d4 | 1366 | if (eval_get_sign(a) == 0) |
7c673cae | 1367 | { |
1e59de90 | 1368 | BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand.")); |
7c673cae | 1369 | } |
92f5a8d4 | 1370 | if (a.sign()) |
7c673cae | 1371 | { |
1e59de90 | 1372 | BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined.")); |
7c673cae FG |
1373 | } |
1374 | return eval_msb_imp(a); | |
1375 | } | |
1376 | ||
1e59de90 TL |
1377 | template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1> |
1378 | inline BOOST_MP_CXX14_CONSTEXPR std::size_t hash_value(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept | |
7c673cae FG |
1379 | { |
1380 | std::size_t result = 0; | |
1e59de90 | 1381 | for (std::size_t i = 0; i < val.size(); ++i) |
7c673cae | 1382 | { |
1e59de90 | 1383 | boost::multiprecision::detail::hash_combine(result, val.limbs()[i]); |
7c673cae | 1384 | } |
1e59de90 | 1385 | boost::multiprecision::detail::hash_combine(result, val.sign()); |
7c673cae FG |
1386 | return result; |
1387 | } | |
1388 | ||
1389 | #ifdef BOOST_MSVC | |
1390 | #pragma warning(pop) | |
1391 | #endif | |
1392 | ||
1e59de90 TL |
1393 | } // Namespace backends |
1394 | ||
1395 | namespace detail { | |
1396 | ||
1397 | #ifndef BOOST_MP_STANDALONE | |
1398 | template <typename T> | |
1399 | inline BOOST_CXX14_CONSTEXPR T constexpr_gcd(T a, T b) noexcept | |
1400 | { | |
1401 | return boost::integer::gcd(a, b); | |
1402 | } | |
1403 | ||
1404 | template <typename T> | |
1405 | inline BOOST_CXX14_CONSTEXPR T constexpr_lcm(T a, T b) noexcept | |
1406 | { | |
1407 | return boost::integer::lcm(a, b); | |
1408 | } | |
1409 | ||
1410 | #else | |
1411 | ||
1412 | template <typename T> | |
1413 | inline BOOST_CXX14_CONSTEXPR T constexpr_gcd(T a, T b) noexcept | |
1414 | { | |
1415 | return boost::multiprecision::backends::eval_gcd(a, b); | |
1416 | } | |
1417 | ||
1418 | template <typename T> | |
1419 | inline BOOST_CXX14_CONSTEXPR T constexpr_lcm(T a, T b) noexcept | |
1420 | { | |
1421 | const T ab_gcd = boost::multiprecision::detail::constexpr_gcd(a, b); | |
1422 | return (a * b) / ab_gcd; | |
1423 | } | |
1424 | ||
1425 | #endif // BOOST_MP_STANDALONE | |
1426 | ||
1427 | } | |
1428 | ||
1429 | }} // Namespace boost::multiprecision | |
7c673cae FG |
1430 | |
1431 | #endif |