]>
Commit | Line | Data |
---|---|---|
7c673cae | 1 | /////////////////////////////////////////////////////////////// |
1e59de90 TL |
2 | // Copyright 2012-21 John Maddock. |
3 | // Copyright 2021 Iskandarov Lev. Distributed under the Boost | |
7c673cae | 4 | // Software License, Version 1.0. (See accompanying file |
92f5a8d4 | 5 | // LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt |
7c673cae FG |
6 | |
7 | #ifndef BOOST_MP_INTEGER_HPP | |
8 | #define BOOST_MP_INTEGER_HPP | |
9 | ||
1e59de90 | 10 | #include <type_traits> |
7c673cae FG |
11 | #include <boost/multiprecision/cpp_int.hpp> |
12 | #include <boost/multiprecision/detail/bitscan.hpp> | |
1e59de90 TL |
13 | #include <boost/multiprecision/detail/no_exceptions_support.hpp> |
14 | #include <boost/multiprecision/detail/standalone_config.hpp> | |
7c673cae | 15 | |
92f5a8d4 TL |
16 | namespace boost { |
17 | namespace multiprecision { | |
7c673cae FG |
18 | |
19 | template <class Integer, class I2> | |
1e59de90 | 20 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && boost::multiprecision::detail::is_integral<I2>::value, Integer&>::type |
92f5a8d4 | 21 | multiply(Integer& result, const I2& a, const I2& b) |
7c673cae FG |
22 | { |
23 | return result = static_cast<Integer>(a) * static_cast<Integer>(b); | |
24 | } | |
25 | template <class Integer, class I2> | |
1e59de90 | 26 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && boost::multiprecision::detail::is_integral<I2>::value, Integer&>::type |
92f5a8d4 | 27 | add(Integer& result, const I2& a, const I2& b) |
7c673cae FG |
28 | { |
29 | return result = static_cast<Integer>(a) + static_cast<Integer>(b); | |
30 | } | |
31 | template <class Integer, class I2> | |
1e59de90 | 32 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && boost::multiprecision::detail::is_integral<I2>::value, Integer&>::type |
92f5a8d4 | 33 | subtract(Integer& result, const I2& a, const I2& b) |
7c673cae FG |
34 | { |
35 | return result = static_cast<Integer>(a) - static_cast<Integer>(b); | |
36 | } | |
37 | ||
38 | template <class Integer> | |
1e59de90 | 39 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value>::type divide_qr(const Integer& x, const Integer& y, Integer& q, Integer& r) |
7c673cae FG |
40 | { |
41 | q = x / y; | |
42 | r = x % y; | |
43 | } | |
44 | ||
45 | template <class I1, class I2> | |
1e59de90 | 46 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<I1>::value && boost::multiprecision::detail::is_integral<I2>::value, I2>::type integer_modulus(const I1& x, I2 val) |
7c673cae FG |
47 | { |
48 | return static_cast<I2>(x % val); | |
49 | } | |
50 | ||
92f5a8d4 | 51 | namespace detail { |
7c673cae FG |
52 | // |
53 | // Figure out the kind of integer that has twice as many bits as some builtin | |
54 | // integer type I. Use a native type if we can (including types which may not | |
1e59de90 | 55 | // be recognised by boost::int_t because they're larger than long long), |
7c673cae FG |
56 | // otherwise synthesize a cpp_int to do the job. |
57 | // | |
58 | template <class I> | |
59 | struct double_integer | |
60 | { | |
1e59de90 TL |
61 | static constexpr const unsigned int_t_digits = |
62 | 2 * sizeof(I) <= sizeof(long long) ? std::numeric_limits<I>::digits * 2 : 1; | |
7c673cae | 63 | |
1e59de90 TL |
64 | using type = typename std::conditional< |
65 | 2 * sizeof(I) <= sizeof(long long), | |
66 | typename std::conditional< | |
67 | boost::multiprecision::detail::is_signed<I>::value && boost::multiprecision::detail::is_integral<I>::value, | |
68 | typename boost::multiprecision::detail::int_t<int_t_digits>::least, | |
69 | typename boost::multiprecision::detail::uint_t<int_t_digits>::least>::type, | |
70 | typename std::conditional< | |
92f5a8d4 | 71 | 2 * sizeof(I) <= sizeof(double_limb_type), |
1e59de90 TL |
72 | typename std::conditional< |
73 | boost::multiprecision::detail::is_signed<I>::value && boost::multiprecision::detail::is_integral<I>::value, | |
92f5a8d4 TL |
74 | signed_double_limb_type, |
75 | double_limb_type>::type, | |
1e59de90 | 76 | number<cpp_int_backend<sizeof(I) * CHAR_BIT * 2, sizeof(I) * CHAR_BIT * 2, (boost::multiprecision::detail::is_signed<I>::value ? signed_magnitude : unsigned_magnitude), unchecked, void> > >::type>::type; |
7c673cae FG |
77 | }; |
78 | ||
92f5a8d4 | 79 | } // namespace detail |
7c673cae FG |
80 | |
81 | template <class I1, class I2, class I3> | |
1e59de90 | 82 | BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<I1>::value && boost::multiprecision::detail::is_unsigned<I2>::value && boost::multiprecision::detail::is_integral<I3>::value, I1>::type |
92f5a8d4 | 83 | powm(const I1& a, I2 b, I3 c) |
7c673cae | 84 | { |
1e59de90 | 85 | using double_type = typename detail::double_integer<I1>::type; |
7c673cae | 86 | |
92f5a8d4 TL |
87 | I1 x(1), y(a); |
88 | double_type result(0); | |
7c673cae | 89 | |
92f5a8d4 | 90 | while (b > 0) |
7c673cae | 91 | { |
92f5a8d4 | 92 | if (b & 1) |
7c673cae FG |
93 | { |
94 | multiply(result, x, y); | |
95 | x = integer_modulus(result, c); | |
96 | } | |
97 | multiply(result, y, y); | |
98 | y = integer_modulus(result, c); | |
99 | b >>= 1; | |
100 | } | |
101 | return x % c; | |
102 | } | |
103 | ||
104 | template <class I1, class I2, class I3> | |
1e59de90 | 105 | inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<I1>::value && boost::multiprecision::detail::is_signed<I2>::value && boost::multiprecision::detail::is_integral<I2>::value && boost::multiprecision::detail::is_integral<I3>::value, I1>::type |
92f5a8d4 | 106 | powm(const I1& a, I2 b, I3 c) |
7c673cae | 107 | { |
92f5a8d4 | 108 | if (b < 0) |
7c673cae | 109 | { |
1e59de90 | 110 | BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); |
7c673cae | 111 | } |
1e59de90 | 112 | return powm(a, static_cast<typename boost::multiprecision::detail::make_unsigned<I2>::type>(b), c); |
7c673cae FG |
113 | } |
114 | ||
115 | template <class Integer> | |
1e59de90 | 116 | BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, std::size_t>::type lsb(const Integer& val) |
7c673cae | 117 | { |
92f5a8d4 | 118 | if (val <= 0) |
7c673cae | 119 | { |
92f5a8d4 | 120 | if (val == 0) |
7c673cae | 121 | { |
1e59de90 | 122 | BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand.")); |
7c673cae FG |
123 | } |
124 | else | |
125 | { | |
1e59de90 | 126 | BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined.")); |
7c673cae FG |
127 | } |
128 | } | |
129 | return detail::find_lsb(val); | |
130 | } | |
131 | ||
132 | template <class Integer> | |
1e59de90 | 133 | BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, std::size_t>::type msb(Integer val) |
7c673cae | 134 | { |
92f5a8d4 | 135 | if (val <= 0) |
7c673cae | 136 | { |
92f5a8d4 | 137 | if (val == 0) |
7c673cae | 138 | { |
1e59de90 | 139 | BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand.")); |
7c673cae FG |
140 | } |
141 | else | |
142 | { | |
1e59de90 | 143 | BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined.")); |
7c673cae FG |
144 | } |
145 | } | |
146 | return detail::find_msb(val); | |
147 | } | |
148 | ||
149 | template <class Integer> | |
1e59de90 | 150 | BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, bool>::type bit_test(const Integer& val, std::size_t index) |
7c673cae FG |
151 | { |
152 | Integer mask = 1; | |
92f5a8d4 | 153 | if (index >= sizeof(Integer) * CHAR_BIT) |
7c673cae | 154 | return 0; |
92f5a8d4 | 155 | if (index) |
7c673cae FG |
156 | mask <<= index; |
157 | return val & mask ? true : false; | |
158 | } | |
159 | ||
160 | template <class Integer> | |
1e59de90 | 161 | BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer&>::type bit_set(Integer& val, std::size_t index) |
7c673cae FG |
162 | { |
163 | Integer mask = 1; | |
92f5a8d4 | 164 | if (index >= sizeof(Integer) * CHAR_BIT) |
7c673cae | 165 | return val; |
92f5a8d4 | 166 | if (index) |
7c673cae FG |
167 | mask <<= index; |
168 | val |= mask; | |
169 | return val; | |
170 | } | |
171 | ||
172 | template <class Integer> | |
1e59de90 | 173 | BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer&>::type bit_unset(Integer& val, std::size_t index) |
7c673cae FG |
174 | { |
175 | Integer mask = 1; | |
92f5a8d4 | 176 | if (index >= sizeof(Integer) * CHAR_BIT) |
7c673cae | 177 | return val; |
92f5a8d4 | 178 | if (index) |
7c673cae FG |
179 | mask <<= index; |
180 | val &= ~mask; | |
181 | return val; | |
182 | } | |
183 | ||
184 | template <class Integer> | |
1e59de90 | 185 | BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer&>::type bit_flip(Integer& val, std::size_t index) |
7c673cae FG |
186 | { |
187 | Integer mask = 1; | |
92f5a8d4 | 188 | if (index >= sizeof(Integer) * CHAR_BIT) |
7c673cae | 189 | return val; |
92f5a8d4 | 190 | if (index) |
7c673cae FG |
191 | mask <<= index; |
192 | val ^= mask; | |
193 | return val; | |
194 | } | |
195 | ||
1e59de90 TL |
196 | namespace detail { |
197 | ||
198 | template <class Integer> | |
199 | BOOST_MP_CXX14_CONSTEXPR Integer karatsuba_sqrt(const Integer& x, Integer& r, size_t bits) | |
200 | { | |
201 | // | |
202 | // Define the floating point type used for std::sqrt, in our tests, sqrt(double) and sqrt(long double) take | |
203 | // about the same amount of time as long as long double is not an emulated 128-bit type (ie the same type | |
204 | // as __float128 from libquadmath). So only use long double if it's an 80-bit type: | |
205 | // | |
206 | #ifndef __clang__ | |
207 | typedef typename std::conditional<(std::numeric_limits<long double>::digits == 64), long double, double>::type real_cast_type; | |
208 | #else | |
209 | // clang has buggy __int128 -> long double conversion: | |
210 | typedef double real_cast_type; | |
211 | #endif | |
212 | // | |
213 | // As per the Karatsuba sqrt algorithm, the low order bits/4 bits pay no part in the result, only in the remainder, | |
214 | // so define the number of bits our argument must have before passing to std::sqrt is safe, even if doing so | |
215 | // looses a few bits: | |
216 | // | |
217 | constexpr std::size_t cutoff = (std::numeric_limits<real_cast_type>::digits * 4) / 3; | |
218 | // | |
219 | // Type which can hold at least "cutoff" bits: | |
220 | // | |
221 | #ifdef BOOST_HAS_INT128 | |
222 | using cutoff_t = typename std::conditional<(cutoff > 64), uint128_type, std::uint64_t>::type; | |
223 | #else | |
224 | using cutoff_t = std::uint64_t; | |
225 | #endif | |
226 | // | |
227 | // See if we can take the fast path: | |
228 | // | |
229 | if (bits <= cutoff) | |
230 | { | |
231 | constexpr cutoff_t half_bits = (cutoff_t(1u) << ((sizeof(cutoff_t) * CHAR_BIT) / 2)) - 1; | |
232 | cutoff_t val = static_cast<cutoff_t>(x); | |
233 | real_cast_type real_val = static_cast<real_cast_type>(val); | |
234 | cutoff_t s64 = static_cast<cutoff_t>(std::sqrt(real_val)); | |
235 | // converting to long double can loose some precision, and `sqrt` can give eps error, so we'll fix this | |
236 | // this is needed | |
237 | while ((s64 > half_bits) || (s64 * s64 > val)) | |
238 | s64--; | |
239 | // in my tests this never fired, but theoretically this might be needed | |
240 | while ((s64 < half_bits) && ((s64 + 1) * (s64 + 1) <= val)) | |
241 | s64++; | |
242 | r = static_cast<Integer>(val - s64 * s64); | |
243 | return static_cast<Integer>(s64); | |
244 | } | |
245 | // https://hal.inria.fr/file/index/docid/72854/filename/RR-3805.pdf | |
246 | std::size_t b = bits / 4; | |
247 | Integer q = x; | |
248 | q >>= b * 2; | |
249 | Integer s = karatsuba_sqrt(q, r, bits - b * 2); | |
250 | Integer t = 0u; | |
251 | bit_set(t, static_cast<unsigned>(b * 2)); | |
252 | r <<= b; | |
253 | t--; | |
254 | t &= x; | |
255 | t >>= b; | |
256 | t += r; | |
257 | s <<= 1; | |
258 | divide_qr(t, s, q, r); | |
259 | r <<= b; | |
260 | t = 0u; | |
261 | bit_set(t, static_cast<unsigned>(b)); | |
262 | t--; | |
263 | t &= x; | |
264 | r += t; | |
265 | s <<= (b - 1); // we already <<1 it before | |
266 | s += q; | |
267 | q *= q; | |
268 | // we substract after, so it works for unsigned integers too | |
269 | if (r < q) | |
270 | { | |
271 | t = s; | |
272 | t <<= 1; | |
273 | t--; | |
274 | r += t; | |
275 | s--; | |
276 | } | |
277 | r -= q; | |
278 | return s; | |
279 | } | |
280 | ||
7c673cae | 281 | template <class Integer> |
1e59de90 | 282 | BOOST_MP_CXX14_CONSTEXPR Integer bitwise_sqrt(const Integer& x, Integer& r) |
7c673cae FG |
283 | { |
284 | // | |
285 | // This is slow bit-by-bit integer square root, see for example | |
286 | // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29 | |
287 | // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf | |
288 | // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented | |
289 | // at some point. | |
290 | // | |
291 | Integer s = 0; | |
92f5a8d4 | 292 | if (x == 0) |
7c673cae FG |
293 | { |
294 | r = 0; | |
295 | return s; | |
296 | } | |
1e59de90 | 297 | std::ptrdiff_t g = msb(x); |
92f5a8d4 | 298 | if (g == 0) |
7c673cae FG |
299 | { |
300 | r = 1; | |
301 | return s; | |
302 | } | |
92f5a8d4 | 303 | |
7c673cae | 304 | Integer t = 0; |
1e59de90 | 305 | r = x; |
7c673cae FG |
306 | g /= 2; |
307 | bit_set(s, g); | |
308 | bit_set(t, 2 * g); | |
309 | r = x - t; | |
310 | --g; | |
311 | do | |
312 | { | |
313 | t = s; | |
314 | t <<= g + 1; | |
315 | bit_set(t, 2 * g); | |
92f5a8d4 | 316 | if (t <= r) |
7c673cae FG |
317 | { |
318 | bit_set(s, g); | |
319 | r -= t; | |
320 | } | |
321 | --g; | |
92f5a8d4 | 322 | } while (g >= 0); |
7c673cae FG |
323 | return s; |
324 | } | |
325 | ||
1e59de90 TL |
326 | } // namespace detail |
327 | ||
328 | template <class Integer> | |
329 | BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer>::type sqrt(const Integer& x, Integer& r) | |
330 | { | |
331 | #ifndef BOOST_MP_NO_CONSTEXPR_DETECTION | |
332 | // recursive Karatsuba sqrt can cause issues in constexpr context: | |
333 | if (BOOST_MP_IS_CONST_EVALUATED(x)) | |
334 | return detail::bitwise_sqrt(x, r); | |
335 | #endif | |
336 | if (x == 0u) { | |
337 | r = 0u; | |
338 | return 0u; | |
339 | } | |
340 | ||
341 | return detail::karatsuba_sqrt(x, r, msb(x) + 1); | |
342 | } | |
343 | ||
7c673cae | 344 | template <class Integer> |
1e59de90 | 345 | BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer>::type sqrt(const Integer& x) |
7c673cae | 346 | { |
92f5a8d4 | 347 | Integer r(0); |
7c673cae FG |
348 | return sqrt(x, r); |
349 | } | |
350 | ||
92f5a8d4 | 351 | }} // namespace boost::multiprecision |
7c673cae FG |
352 | |
353 | #endif |