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