Branch data Line data Source code
1 : : /********************************************************************
2 : : * qofmath128.c -- an 128-bit integer library *
3 : : * Copyright (C) 2004 Linas Vepstas <linas@linas.org> *
4 : : * Copyright (C) 2014 John Ralls <jralls@ceridwen.us> *
5 : : * *
6 : : * This program is free software; you can redistribute it and/or *
7 : : * modify it under the terms of the GNU General Public License as *
8 : : * published by the Free Software Foundation; either version 2 of *
9 : : * the License, or (at your option) any later version. *
10 : : * *
11 : : * This program is distributed in the hope that it will be useful, *
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 : : * GNU General Public License for more details. *
15 : : * *
16 : : * You should have received a copy of the GNU General Public License*
17 : : * along with this program; if not, contact: *
18 : : * *
19 : : * Free Software Foundation Voice: +1-617-542-5942 *
20 : : * 51 Franklin Street, Fifth Floor Fax: +1-617-542-2652 *
21 : : * Boston, MA 02110-1301, USA gnu@gnu.org *
22 : : * *
23 : : *******************************************************************/
24 : : #include <config.h>
25 : :
26 : : #include "gnc-int128.hpp"
27 : :
28 : : #include <iomanip>
29 : : #include <utility>
30 : : #include <cassert>
31 : : #include <cstdio>
32 : : #include <sstream>
33 : :
34 : : /* All algorithms from Donald E. Knuth, "The Art of Computer
35 : : * Programming, Volume 2: Seminumerical Algorithms", 3rd Ed.,
36 : : * Addison-Wesley, 1998.
37 : : */
38 : :
39 : : namespace {
40 : : static const unsigned int upper_num_bits = 61;
41 : : static const unsigned int sublegs = GncInt128::numlegs * 2;
42 : : static const unsigned int sublegbits = GncInt128::legbits / 2;
43 : : static const uint64_t sublegmask = (UINT64_C(1) << sublegbits) - 1;
44 : : static const uint64_t flagmask = UINT64_C(0xe000000000000000);
45 : : static const uint64_t nummask = UINT64_C(0x1fffffffffffffff);
46 : : /* Assumes that the high bits represent old flags to be replaced.
47 : : * Any overflow must be detected first!
48 : : */
49 : 24419525 : static inline uint64_t set_flags(uint64_t leg, uint8_t flags)
50 : : {
51 : 24419525 : auto flag_part = static_cast<uint64_t>(flags) << upper_num_bits;
52 : 24419525 : return flag_part + (leg & nummask);
53 : : }
54 : 76319081 : static inline uint8_t get_flags(uint64_t leg)
55 : : {
56 : 76319081 : return (leg & flagmask) >> upper_num_bits;
57 : : }
58 : 33720461 : static inline uint64_t get_num(uint64_t leg)
59 : : {
60 : 33720461 : return leg & nummask;
61 : : }
62 : : }
63 : :
64 : 2484462 : GncInt128::GncInt128 () : m_hi {0}, m_lo {0}{}
65 : :
66 : 3991610 : GncInt128::GncInt128 (int64_t upper, int64_t lower, uint8_t flags) :
67 : 3991610 : m_hi {static_cast<uint64_t>(upper < 0 ? -upper : upper)},
68 : 3991610 : m_lo {static_cast<uint64_t>(lower < 0 ? -lower : lower)}
69 : : {
70 : 3991610 : if ((upper < 0 && lower > 0) || (upper > 0 && lower < 0))
71 : 0 : m_lo = (m_hi << 63) - m_lo;
72 : : else
73 : 3991610 : m_lo += (m_hi << 63);
74 : :
75 : 3991610 : m_hi >>= 1;
76 : 3991610 : if (m_hi & flagmask)
77 : : {
78 : 0 : std::ostringstream ss;
79 : 0 : ss << "Constructing GncInt128 with int64_t " << upper
80 : 0 : << " which is too big.";
81 : 0 : throw std::overflow_error(ss.str());
82 : 0 : }
83 : 3991610 : flags ^= (upper < 0 ? neg : upper == 0 && lower < 0 ? neg : pos);
84 : 3991610 : m_hi = set_flags(m_hi, flags);
85 : 3991610 : }
86 : :
87 : 0 : GncInt128::GncInt128 (int64_t upper, uint64_t lower, uint8_t flags) :
88 : 0 : m_hi {static_cast<uint64_t>(upper < 0 ? -upper : upper)}, m_lo {lower}
89 : : {
90 : 0 : if (m_hi & flagmask)
91 : : {
92 : 0 : std::ostringstream ss;
93 : 0 : ss << "Constructing GncInt128 with int64_t " << upper
94 : 0 : << " which is too big when lower is unsigned.";
95 : 0 : throw std::overflow_error(ss.str());
96 : 0 : }
97 : 0 : flags ^= (upper < 0 ? neg : pos);
98 : 0 : m_hi = set_flags(m_hi, flags);
99 : 0 : }
100 : :
101 : 3969017 : GncInt128::GncInt128 (uint64_t upper, uint64_t lower, uint8_t flags) :
102 : 3969017 : m_hi {upper}, m_lo {lower}
103 : : {
104 : : /* somewhere in the bowels of gnc_module.scm compilation this gets called
105 : : * with upper=INT64_MAX, which would otherwise throw. Make it our max value
106 : : * instead.
107 : : */
108 : 3969017 : if (m_hi == UINT64_MAX)
109 : 1338 : m_hi = nummask;
110 : 3969017 : if (m_hi & flagmask)
111 : : {
112 : 0 : std::ostringstream ss;
113 : 0 : ss << "Constructing GncInt128 with uint64_t " << upper
114 : 0 : << " which is too big.";
115 : 0 : throw std::overflow_error(ss.str());
116 : 0 : }
117 : :
118 : 3969017 : m_hi = set_flags(m_hi, flags);
119 : 3969017 : }
120 : :
121 : : GncInt128&
122 : 2484462 : GncInt128::zero () noexcept
123 : : {
124 : 2484462 : m_lo = m_hi = UINT64_C(0);
125 : 2484462 : return *this;
126 : : }
127 : :
128 : 814878 : GncInt128::operator int64_t() const
129 : : {
130 : 814878 : auto flags = get_flags(m_hi);
131 : 814878 : if ((flags & neg) && isBig())
132 : 0 : throw std::underflow_error ("Negative value too large to represent as int64_t");
133 : 814878 : if ((flags & (overflow | NaN)) || isBig())
134 : 15 : throw std::overflow_error ("Value too large to represent as int64_t");
135 : 814863 : int64_t retval = static_cast<int64_t>(m_lo);
136 : 814863 : return flags & neg ? -retval : retval;
137 : : }
138 : :
139 : 0 : GncInt128::operator uint64_t() const
140 : : {
141 : 0 : auto flags = get_flags(m_hi);
142 : 0 : if ((flags & neg) && !isZero()) // exclude negative zero
143 : 0 : throw std::underflow_error ("Can't represent negative value as uint64_t");
144 : 0 : if ((flags & (overflow | NaN)) || (m_hi || m_lo > UINT64_MAX))
145 : 0 : throw std::overflow_error ("Value to large to represent as uint64_t");
146 : 0 : return m_lo;
147 : : }
148 : :
149 : :
150 : : int
151 : 2578048 : GncInt128::cmp (const GncInt128& b) const noexcept
152 : : {
153 : 2578048 : auto flags = get_flags(m_hi);
154 : 2578048 : if (flags & (overflow | NaN))
155 : 0 : return -1;
156 : 2578048 : if (b.isOverflow () || b.isNan ())
157 : 0 : return 1;
158 : 2578048 : auto hi = get_num(m_hi);
159 : 2578048 : auto bhi = get_num(b.m_hi);
160 : 2578048 : if (isZero() && b.isZero()) return 0;
161 : 2556197 : if (flags & neg)
162 : : {
163 : 864 : if (!b.isNeg()) return -1;
164 : 212 : if (hi > bhi) return -1;
165 : 212 : if (hi < bhi) return 1;
166 : 212 : if (m_lo > b.m_lo) return -1;
167 : 80 : if (m_lo < b.m_lo) return 1;
168 : 10 : return 0;
169 : : }
170 : 2555333 : if (b.isNeg()) return 1;
171 : 2555320 : if (hi < bhi) return -1;
172 : 2554008 : if (hi > bhi) return 1;
173 : 2543708 : if (m_lo < b.m_lo) return -1;
174 : 1856410 : if (m_lo > b.m_lo) return 1;
175 : 866424 : return 0;
176 : : }
177 : :
178 : : /* Knuth 4.5.3 Algo B, recommended by GMP as much faster than Algo A (Euclidean
179 : : * method).
180 : : */
181 : : GncInt128
182 : 362122 : GncInt128::gcd(GncInt128 b) const noexcept
183 : : {
184 : 362122 : if (b.isZero())
185 : 252 : return *this;
186 : 361870 : if (isZero())
187 : 0 : return b;
188 : :
189 : 361870 : if (b.isOverflow() || b.isNan())
190 : 0 : return b;
191 : 361870 : if (isOverflow() || isNan())
192 : 0 : return *this;
193 : :
194 : 361870 : GncInt128 a (isNeg() ? -(*this) : *this);
195 : 361870 : if (b.isNeg()) b = -b;
196 : :
197 : 361870 : unsigned int k {};
198 : 361870 : const uint64_t one {1};
199 : 753773 : while (!((a & one) || (b & one))) //B1
200 : : {
201 : 391903 : a >>= 1;
202 : 391903 : b >>= 1;
203 : 391903 : ++k;
204 : : }
205 : 361870 : GncInt128 t {a & one ? -b : a}; //B2
206 : 761645 : while (a != b)
207 : : {
208 : 1213552 : while (t && ((t & one) ^ one)) t >>= 1; //B3 & B4
209 : 399775 : if (t.isNeg()) //B5
210 : 320932 : b = -t;
211 : : else
212 : 78843 : a = t;
213 : 399775 : t = a - b; //B6
214 : : }
215 : 361870 : return a << k;
216 : : }
217 : :
218 : : /* Since u * v = gcd(u, v) * lcm(u, v), we find lcm by u / gcd * v. */
219 : :
220 : : GncInt128
221 : 257841 : GncInt128::lcm(const GncInt128& b) const noexcept
222 : : {
223 : 257841 : auto common = gcd(b);
224 : 257841 : return *this / common * b.abs(); //Preserve our sign, discard the other's.
225 : : }
226 : :
227 : : /* Knuth section 4.6.3 */
228 : : GncInt128
229 : 0 : GncInt128::pow(unsigned int b) const noexcept
230 : : {
231 : 0 : if (isZero() || (m_lo == 1 && m_hi == 0) || isNan() || isOverflow())
232 : 0 : return *this;
233 : 0 : if (b == 0)
234 : 0 : return GncInt128 (1);
235 : 0 : GncInt128 retval (1), squares = *this;
236 : 0 : while (b && !retval.isOverflow())
237 : : {
238 : 0 : if (b & 1)
239 : 0 : retval *= squares;
240 : 0 : squares *= squares;
241 : 0 : b >>= 1;
242 : : }
243 : 0 : return retval;
244 : : }
245 : :
246 : : bool
247 : 13804756 : GncInt128::isNeg () const noexcept
248 : : {
249 : 13804756 : return (get_flags(m_hi) & neg);
250 : : }
251 : :
252 : : bool
253 : 1599249 : GncInt128::isBig () const noexcept
254 : : {
255 : 1599249 : return (get_num(m_hi) || m_lo > INT64_MAX);
256 : : }
257 : :
258 : : bool
259 : 17456522 : GncInt128::isOverflow () const noexcept
260 : : {
261 : 17456522 : return (get_flags(m_hi) & overflow);
262 : : }
263 : :
264 : : bool
265 : 17583454 : GncInt128::isNan () const noexcept
266 : : {
267 : 17583454 : return (get_flags(m_hi) & NaN);
268 : : }
269 : :
270 : : bool
271 : 1615030 : GncInt128::valid() const noexcept
272 : : {
273 : 1615030 : return !(get_flags(m_hi) & (overflow | NaN));
274 : : }
275 : :
276 : : bool
277 : 10387796 : GncInt128::isZero() const noexcept
278 : : {
279 : 20775592 : return ((get_flags(m_hi) & (overflow | NaN)) == 0 &&
280 : 20775592 : get_num(m_hi) == 0 && m_lo == 0);
281 : : }
282 : :
283 : : GncInt128
284 : 3718565 : GncInt128::abs() const noexcept
285 : : {
286 : 3718565 : if (isNeg())
287 : 411536 : return operator-();
288 : :
289 : 3307029 : return *this;
290 : : }
291 : :
292 : : unsigned int
293 : 1656116 : GncInt128::bits() const noexcept
294 : : {
295 : 1656116 : auto hi = get_num(m_hi);
296 : 1656116 : unsigned int bits {static_cast<unsigned int>(hi == 0 ? 0 : 64)};
297 : 1656116 : uint64_t temp {(hi == 0 ? m_lo : hi)};
298 : 14306078 : for (;temp > 0; temp >>= 1)
299 : 12649962 : ++bits;
300 : 1656116 : return bits;
301 : : }
302 : :
303 : :
304 : : GncInt128
305 : 1165593 : GncInt128::operator-() const noexcept
306 : : {
307 : 1165593 : auto retval = *this;
308 : 1165593 : auto flags = get_flags(retval.m_hi);
309 : 1165593 : if (isNeg())
310 : 784468 : flags ^= neg;
311 : : else
312 : 381125 : flags |= neg;
313 : 1165593 : retval.m_hi = set_flags(retval.m_hi, flags);
314 : 1165593 : return retval;
315 : : }
316 : :
317 : 3970139 : GncInt128::operator bool() const noexcept
318 : : {
319 : 3970139 : return ! isZero ();
320 : : }
321 : :
322 : : GncInt128&
323 : 0 : GncInt128::operator++ () noexcept
324 : : {
325 : 0 : return operator+=(UINT64_C(1));
326 : : }
327 : :
328 : : GncInt128&
329 : 0 : GncInt128::operator++ (int) noexcept
330 : : {
331 : 0 : return operator+=(UINT64_C(1));
332 : : }
333 : :
334 : : GncInt128&
335 : 0 : GncInt128::operator-- () noexcept
336 : : {
337 : 0 : return operator-=(UINT64_C(1));
338 : : }
339 : :
340 : : GncInt128&
341 : 0 : GncInt128::operator-- (int) noexcept
342 : : {
343 : 0 : return operator-=(UINT64_C(1));
344 : : }
345 : :
346 : : GncInt128&
347 : 257581 : GncInt128::operator+= (const GncInt128& b) noexcept
348 : : {
349 : 257581 : auto flags = get_flags(m_hi);
350 : 257581 : if (b.isOverflow())
351 : 0 : flags |= overflow;
352 : 257581 : if (b.isNan())
353 : 0 : flags |= NaN;
354 : 257581 : m_hi = set_flags(m_hi, flags);
355 : 257581 : if (isOverflow() || isNan())
356 : 0 : return *this;
357 : 257581 : if ((isNeg () && !b.isNeg ()) || (!isNeg () && b.isNeg ()))
358 : 87673 : return this->operator-= (-b);
359 : 169908 : uint64_t result = m_lo + b.m_lo;
360 : 169908 : uint64_t carry = static_cast<int64_t>(result < m_lo); //Wrapping
361 : 169908 : m_lo = result;
362 : 169908 : auto hi = get_num(m_hi);
363 : 169908 : auto bhi = get_num(b.m_hi);
364 : 169908 : result = hi + bhi + carry;
365 : 169908 : if (result < hi || result & flagmask)
366 : 0 : flags |= overflow;
367 : 169908 : m_hi = set_flags(result, flags);
368 : 169908 : return *this;
369 : : }
370 : :
371 : : GncInt128&
372 : 361870 : GncInt128::operator<<= (unsigned int i) noexcept
373 : : {
374 : 361870 : auto flags = get_flags(m_hi);
375 : 361870 : if (i == 0)
376 : 170452 : return *this;
377 : 191418 : if (i > maxbits)
378 : : {
379 : 0 : flags &= 0xfe;
380 : 0 : m_hi = set_flags(0, flags);
381 : 0 : m_lo = 0;
382 : 0 : return *this;
383 : : }
384 : 191418 : auto hi = get_num(m_hi);
385 : 191418 : if (i < legbits)
386 : : {
387 : 191418 : uint64_t carry {(m_lo & (((UINT64_C(1) << i) - 1) << (legbits - i))) >> (legbits - i)};
388 : 191418 : m_lo <<= i;
389 : 191418 : hi <<= i;
390 : 191418 : hi += carry;
391 : 191418 : m_hi = set_flags(hi, flags);
392 : 191418 : return *this;
393 : : }
394 : 0 : hi = m_lo << (i - legbits);
395 : 0 : m_hi = set_flags(hi, flags);
396 : 0 : m_lo = 0;
397 : 0 : return *this;
398 : : }
399 : :
400 : : GncInt128&
401 : 1597650 : GncInt128::operator>>= (unsigned int i) noexcept
402 : : {
403 : 1597650 : auto flags = get_flags(m_hi);
404 : 1597650 : if (i > maxbits)
405 : : {
406 : 0 : flags &= 0xfe;
407 : 0 : m_hi = set_flags(0, flags);
408 : 0 : m_lo = 0;
409 : 0 : return *this;
410 : : }
411 : 1597650 : auto hi = get_num(m_hi);
412 : 1597650 : if (i < legbits)
413 : : {
414 : 1597622 : uint64_t carry {(hi & ((UINT64_C(1) << i) - 1))};
415 : 1597622 : m_lo >>= i;
416 : 1597622 : hi >>= i;
417 : 1597622 : m_lo += (carry << (legbits - i));
418 : 1597622 : m_hi = set_flags(hi, flags);
419 : 1597622 : return *this;
420 : : }
421 : 28 : m_lo = hi >> (i - legbits);
422 : 28 : m_hi = set_flags(0, flags);
423 : 28 : return *this;
424 : : }
425 : :
426 : : GncInt128&
427 : 487448 : GncInt128::operator-= (const GncInt128& b) noexcept
428 : : {
429 : 487448 : auto flags = get_flags(m_hi);
430 : 487448 : if (b.isOverflow())
431 : 0 : flags |= overflow;
432 : 487448 : if (b.isNan())
433 : 0 : flags |= NaN;
434 : 487448 : m_hi = set_flags(m_hi, flags);
435 : :
436 : 487448 : if (isOverflow() || isNan())
437 : 0 : return *this;
438 : :
439 : 487448 : if ((!isNeg() && b.isNeg()) || (isNeg() && !b.isNeg()))
440 : 0 : return this->operator+= (-b);
441 : 487448 : bool operand_bigger {abs().cmp (b.abs()) < 0};
442 : 487448 : auto hi = get_num(m_hi);
443 : 487448 : auto far_hi = get_num(b.m_hi);
444 : 487448 : if (operand_bigger)
445 : : {
446 : 248814 : flags ^= neg; // ^= flips the bit
447 : 248814 : if (m_lo > b.m_lo)
448 : : {
449 : : /* The + 1 on the end is because we really want to use 2^64, or
450 : : * UINT64_MAX + 1, but that can't be represented in a uint64_t.
451 : : */
452 : 0 : m_lo = UINT64_MAX - m_lo + b.m_lo + 1;
453 : 0 : --far_hi; //borrow
454 : : }
455 : : else
456 : 248814 : m_lo = b.m_lo - m_lo;
457 : 248814 : hi = far_hi - hi;
458 : 248814 : m_hi = set_flags(hi, flags);
459 : 248814 : return *this;
460 : : }
461 : 238634 : if (m_lo < b.m_lo)
462 : : {
463 : 0 : m_lo = UINT64_MAX - b.m_lo + m_lo + 1; //See UINT64_MAX comment above
464 : 0 : --hi; //borrow
465 : : }
466 : : else
467 : 238634 : m_lo -= b.m_lo;
468 : :
469 : 238634 : hi -= far_hi;
470 : 238634 : m_hi = set_flags(hi, flags);
471 : 238634 : return *this;
472 : : }
473 : :
474 : : GncInt128&
475 : 886218 : GncInt128::operator*= (const GncInt128& b) noexcept
476 : : {
477 : : /* Test for 0 first */
478 : 886218 : auto flags = get_flags(m_hi);
479 : : /* Handle the sign; ^ flips if b is negative */
480 : 886218 : flags ^= (get_flags(b.m_hi) & neg);
481 : 886218 : if (isZero() || b.isZero())
482 : : {
483 : 58160 : m_lo = 0;
484 : 58160 : m_hi = set_flags(0, flags);
485 : 58160 : return *this;
486 : : }
487 : 828058 : if (b.isOverflow())
488 : 0 : flags |= overflow;
489 : 828058 : if (b.isNan())
490 : 0 : flags |= NaN;
491 : 828058 : m_hi = set_flags(m_hi, flags);
492 : 828058 : if (isOverflow() || isNan())
493 : 0 : return *this;
494 : :
495 : : /* Test for overflow before spending time on the calculation */
496 : 828058 : auto hi = get_num(m_hi);
497 : 828058 : auto bhi = get_num(b.m_hi);
498 : 828058 : if (hi && bhi)
499 : : {
500 : 0 : flags |= overflow;
501 : 0 : m_hi = set_flags(hi, flags);
502 : 0 : return *this;
503 : : }
504 : :
505 : 828058 : unsigned int abits {bits()}, bbits {b.bits()};
506 : : /* If the product of the high bytes < 7fff then the result will have abits +
507 : : * bbits -1 bits and won't actually overflow. It's not worth the effort to
508 : : * work that out for this corner-case, so we'll take the risk and calculate
509 : : * the product and catch the overflow later.
510 : : */
511 : 828058 : if (abits + bbits - 1 > maxbits)
512 : : {
513 : 0 : flags |= overflow;
514 : 0 : m_hi = set_flags(m_hi, flags);
515 : 0 : return *this;
516 : : }
517 : :
518 : : /* The trivial case */
519 : 828058 : if (abits + bbits <= legbits)
520 : : {
521 : 825399 : m_lo *= b.m_lo;
522 : 825399 : m_hi = set_flags(m_hi, flags);
523 : 825399 : return *this;
524 : : }
525 : :
526 : : /* This is Knuth's "classical" multi-precision multiplication algorithm
527 : : * truncated to a GncInt128 result with the loop unrolled for clarity and with
528 : : * overflow and zero checks beforehand to save time. See Donald Knuth, "The Art
529 : : * of Computer Programming Volume 2: Seminumerical Algorithms", Addison-Wesley,
530 : : * 1998, p 268.
531 : : *
532 : : * There are potentially faster algorithms (O(n^1.6) instead of O(n^2) for the
533 : : * full precision), but this is already close to that fast because of truncation
534 : : * and it's not clear that the truncation is applicable to the faster
535 : : * algorithms.
536 : : */
537 : :
538 : 2659 : uint64_t av[sublegs] {(m_lo & sublegmask), (m_lo >> sublegbits),
539 : 2659 : (hi & sublegmask), (hi >> sublegbits)};
540 : 2659 : uint64_t bv[sublegs] {(b.m_lo & sublegmask), (b.m_lo >> sublegbits),
541 : 2659 : (bhi & sublegmask), (bhi >> sublegbits)};
542 : 2659 : uint64_t rv[sublegs] {};
543 : 2659 : uint64_t carry {}, scratch {};
544 : :
545 : 2659 : rv[0] = av[0] * bv[0];
546 : :
547 : 2659 : rv[1] = av[1] * bv [0];
548 : 2659 : scratch = rv[1] + av[0] * bv[1];
549 : 2659 : carry = rv[1] > scratch ? 1 : 0;
550 : 2659 : rv[1] = scratch;
551 : :
552 : 2659 : rv[2] = av[2] * bv[0] + carry; //0xffff^2 + 1 < 0xffffffff, can't overflow
553 : 2659 : scratch = rv[2] + av[1] * bv[1];
554 : 2659 : carry = rv[2] > scratch ? 1 : 0;
555 : 2659 : rv[2] = scratch + av[0] * bv[2];
556 : 2659 : carry += scratch > rv[2] ? 1 : 0;
557 : :
558 : 2659 : rv[3] = av[3] * bv[0] + carry;
559 : 2659 : scratch = rv[3] + av[2] * bv[1];
560 : 2659 : carry = rv[3] > scratch ? 1 : 0;
561 : 2659 : rv[3] = scratch + av[1] * bv[2];
562 : 2659 : carry += scratch > rv[3] ? 1 : 0;
563 : 2659 : scratch = rv[3] + av[0] * bv[3];
564 : 2659 : carry += rv[3] > scratch ? 1 : 0;
565 : 2659 : rv[3] = scratch;
566 : :
567 : 2659 : if (carry) //Shouldn't happen because of the checks above
568 : : {
569 : 0 : flags |= overflow;
570 : 0 : m_hi = set_flags(m_hi, flags);
571 : 0 : return *this;
572 : : }
573 : :
574 : 2659 : m_lo = rv[0] + (rv[1] << sublegbits);
575 : 2659 : carry = rv[1] >> sublegbits;
576 : 2659 : carry += (rv[1] << sublegbits) > m_lo || rv[0] > m_lo ? 1 : 0;
577 : 2659 : hi = rv[2] + (rv[3] << sublegbits) + carry;
578 : 2659 : if ((rv[3] << sublegbits) > hi || rv[2] > hi || (rv[3] >> sublegbits) ||
579 : : hi & flagmask)
580 : : {
581 : 0 : flags |= overflow;
582 : 0 : m_hi = set_flags(hi, flags);
583 : 0 : return *this;
584 : : }
585 : 2659 : m_hi = set_flags(hi, flags);
586 : 2659 : return *this;
587 : : }
588 : :
589 : : namespace {
590 : : /* Algorithm from Knuth (full citation at operator*=) p272ff. Again, there
591 : : * are faster algorithms out there, but they require much larger numbers to
592 : : * be of benefit.
593 : : */
594 : : /* We're using arrays here instead of vectors to avoid an alloc. */
595 : : void
596 : 82 : div_multi_leg (uint64_t* u, size_t m, uint64_t* v, size_t n,
597 : : GncInt128& q, GncInt128& r) noexcept
598 : : {
599 : : /* D1, Normalization */
600 : 82 : uint64_t qv[sublegs] {};
601 : 82 : uint64_t d {(UINT64_C(1) << sublegbits)/(v[n - 1] + UINT64_C(1))};
602 : 82 : uint64_t carry {UINT64_C(0)};
603 : 82 : bool negative {q.isNeg()};
604 : 82 : bool rnegative {r.isNeg()};
605 : 359 : for (size_t i = 0; i < m; ++i)
606 : : {
607 : 277 : u[i] = u[i] * d + carry;
608 : 277 : if (u[i] > sublegmask)
609 : : {
610 : 226 : carry = u[i] >> sublegbits;
611 : 226 : u[i] &= sublegmask;
612 : : }
613 : : else
614 : 51 : carry = UINT64_C(0);
615 : 277 : assert (u[i] <= sublegmask);
616 : : }
617 : 82 : if (carry)
618 : : {
619 : 39 : u[m++] = carry;
620 : 39 : carry = UINT64_C(0);
621 : : }
622 : 250 : for (size_t i = 0; i < n; ++i)
623 : : {
624 : 168 : v[i] = v[i] * d + carry;
625 : 168 : if (v[i] > sublegmask)
626 : : {
627 : 82 : carry = v[i] >> sublegbits;
628 : 82 : v[i] &= sublegmask;
629 : : }
630 : : else
631 : 86 : carry = UINT64_C(0);
632 : 168 : assert (v[i] < sublegmask);
633 : : }
634 : 82 : assert (carry == UINT64_C(0));
635 : 312 : for (int j = m - n; j >= 0; j--) //D3
636 : : {
637 : : uint64_t qhat, rhat;
638 : 230 : qhat = ((u[j + n] << sublegbits) + u[j + n - 1]) / v[n - 1];
639 : 230 : rhat = ((u[j + n] << sublegbits) + u[j + n - 1]) % v[n - 1];
640 : :
641 : 257 : while (qhat > sublegmask ||
642 : 236 : (rhat <= sublegmask &&
643 : 236 : ((qhat * v[n - 2]) > ((rhat << sublegbits) + u[j + n - 2]))))
644 : : {
645 : 27 : --qhat;
646 : 27 : rhat += v[n - 1];
647 : : }
648 : 230 : carry = UINT64_C(0);
649 : 230 : uint64_t borrow {};
650 : 697 : for (size_t k = 0; k < n; ++k) //D4
651 : : {
652 : 467 : auto subend = qhat * v[k] + carry;
653 : 467 : carry = subend >> sublegbits;
654 : 467 : subend &= sublegmask;
655 : 467 : if (u[j + k] >= subend)
656 : : {
657 : 345 : u[j + k] = u[j + k] - subend;
658 : 345 : borrow = UINT64_C(0);
659 : : }
660 : : else
661 : : {
662 : 122 : if (u[j + k + 1] > 0)
663 : 122 : --u[j + k + 1];
664 : : else
665 : 0 : ++borrow;
666 : 122 : u[j + k] = u[j + k] + sublegmask + 1 - subend;
667 : 122 : u[j + k] &= sublegmask;
668 : : }
669 : : }
670 : 230 : u[j + n] -= carry;
671 : 230 : qv[j] = qhat;
672 : 230 : if (borrow) //D5
673 : : { //D6
674 : 0 : --qv[j];
675 : 0 : carry = UINT64_C(0);
676 : 0 : for (size_t k = 0; k < n; ++k)
677 : : {
678 : 0 : u[j + k] += v[k] + carry;
679 : 0 : if (u[j + k] > sublegmask)
680 : : {
681 : 0 : carry = u[j + k] >> sublegbits;
682 : 0 : u[j + k] &= sublegmask;
683 : : }
684 : : }
685 : 0 : u[j + n] += carry;
686 : : }
687 : : }//D7
688 : 82 : q = GncInt128 ((qv[3] << sublegbits) + qv[2], (qv[1] << sublegbits) + qv[0]);
689 : 82 : r = GncInt128 ((u[3] << sublegbits) + u[2], (u[1] << sublegbits) + u[0]);
690 : 82 : r /= d;
691 : 82 : if (negative) q = -q;
692 : 82 : if (rnegative) r = -r;
693 : 82 : }
694 : :
695 : : void
696 : 7938 : div_single_leg (uint64_t* u, size_t m, uint64_t v,
697 : : GncInt128& q, GncInt128& r) noexcept
698 : : {
699 : 7938 : uint64_t qv[sublegs] {};
700 : 7938 : bool negative {q.isNeg()};
701 : 7938 : bool rnegative {r.isNeg()};
702 : 31766 : for (int i = m - 1; i >= 0; --i)
703 : : {
704 : 23828 : qv[i] = u[i] / v;
705 : 23828 : if (i > 0)
706 : : {
707 : 15890 : u[i - 1] += ((u[i] % v) << sublegbits);
708 : 15890 : u[i] = UINT64_C(0);
709 : : }
710 : : else
711 : 7938 : u[i] %= v;
712 : : }
713 : :
714 : 7938 : q = GncInt128 ((qv[3] << sublegbits) + qv[2], (qv[1] << sublegbits) + qv[0]);
715 : 7938 : r = GncInt128 ((u[3] << sublegbits) + u[2], (u[1] << sublegbits) + u[0]);
716 : 7938 : if (negative) q = -q;
717 : 7938 : if (rnegative) r = -r;
718 : 7938 : }
719 : :
720 : : }// namespace
721 : :
722 : : void
723 : 1242231 : GncInt128::div (const GncInt128& b, GncInt128& q, GncInt128& r) const noexcept
724 : : {
725 : : // clear remainder and quotient
726 : 1242231 : r = GncInt128(0);
727 : 1242231 : q = GncInt128(0);
728 : :
729 : 1242231 : auto qflags = get_flags(q.m_hi);
730 : 1242231 : auto rflags = get_flags(r.m_hi);
731 : 1242231 : if (isOverflow() || b.isOverflow())
732 : : {
733 : 0 : qflags |= overflow;
734 : 0 : rflags |= overflow;
735 : 0 : q.m_hi = set_flags(q.m_hi, qflags);
736 : 0 : r.m_hi = set_flags(r.m_hi, rflags);
737 : 0 : return;
738 : : }
739 : :
740 : 1242231 : if (isNan() || b.isNan())
741 : : {
742 : 0 : qflags |= NaN;
743 : 0 : rflags |= NaN;
744 : 0 : q.m_hi = set_flags(q.m_hi, qflags);
745 : 0 : r.m_hi = set_flags(r.m_hi, rflags);
746 : 0 : return;
747 : : }
748 : 1242231 : assert (&q != this);
749 : 1242231 : assert (&r != this);
750 : 1242231 : assert (&q != &b);
751 : 1242231 : assert (&r != &b);
752 : :
753 : 1242231 : q.zero(), r.zero();
754 : 1242231 : qflags = rflags = 0;
755 : 1242231 : if (b.isZero())
756 : : {
757 : 0 : qflags |= NaN;
758 : 0 : rflags |= NaN;
759 : 0 : q.m_hi = set_flags(q.m_hi, qflags);
760 : 0 : r.m_hi = set_flags(r.m_hi, rflags);
761 : 0 : return;
762 : : }
763 : :
764 : 1242231 : if (isNeg())
765 : : {
766 : 339235 : qflags |= neg;
767 : 339235 : rflags |= neg; // the remainder inherits the dividend's sign
768 : : }
769 : :
770 : 1242231 : if (b.isNeg())
771 : 0 : qflags ^= neg;
772 : :
773 : 1242231 : q.m_hi = set_flags(q.m_hi, qflags);
774 : 1242231 : r.m_hi = set_flags(r.m_hi, rflags);
775 : :
776 : 1242231 : if (abs() < b.abs())
777 : : {
778 : 113134 : r = *this;
779 : 113134 : return;
780 : : }
781 : 1129097 : auto hi = get_num(m_hi);
782 : 1129097 : auto bhi = get_num(b.m_hi);
783 : :
784 : 1129097 : if (hi == 0 && bhi == 0) //let the hardware do it
785 : : {
786 : 1121077 : assert (b.m_lo != 0); // b.m_hi is 0 but b isn't or we didn't get here.
787 : 1121077 : q.m_lo = m_lo / b.m_lo;
788 : 1121077 : r.m_lo = m_lo % b.m_lo;
789 : 1121077 : return;
790 : : }
791 : :
792 : 8020 : uint64_t u[sublegs + 2] {(m_lo & sublegmask), (m_lo >> sublegbits),
793 : 8020 : (hi & sublegmask), (hi >> sublegbits), 0, 0};
794 : 8020 : uint64_t v[sublegs] {(b.m_lo & sublegmask), (b.m_lo >> sublegbits),
795 : 8020 : (bhi & sublegmask), (bhi >> sublegbits)};
796 : 8020 : auto m = u[3] ? 4 : u[2] ? 3 : u[1] ? 2 : u[0] ? 1 : 0;
797 : 8020 : auto n = v[3] ? 4 : v[2] ? 3 : v[1] ? 2 : v[0] ? 1 : 0;
798 : 8020 : if (m == 0 || n == 0) //Shouldn't happen
799 : 0 : return;
800 : 8020 : if (n == 1)
801 : 7938 : return div_single_leg (u, m, v[0], q, r);
802 : :
803 : 82 : return div_multi_leg (u, m, v, n, q, r);
804 : : }
805 : :
806 : : GncInt128&
807 : 1109667 : GncInt128::operator/= (const GncInt128& b) noexcept
808 : : {
809 : 1109667 : GncInt128 q {}, r {};
810 : 1109667 : div(b, q, r);
811 : 1109667 : std::swap (*this, q);
812 : 1109667 : return *this;
813 : : }
814 : :
815 : : GncInt128&
816 : 132555 : GncInt128::operator%= (const GncInt128& b) noexcept
817 : : {
818 : 132555 : GncInt128 q {}, r {};
819 : 132555 : div(b, q, r);
820 : 132555 : std::swap (*this, r);
821 : 132555 : if (q.isNan())
822 : 0 : m_hi = set_flags(m_hi, (get_flags(m_hi) | NaN));
823 : 132555 : return *this;
824 : : }
825 : :
826 : : GncInt128&
827 : 2738005 : GncInt128::operator&= (const GncInt128& b) noexcept
828 : : {
829 : 2738005 : auto flags = get_flags(m_hi);
830 : 2738005 : if (b.isOverflow())
831 : 0 : flags |= overflow;
832 : 2738005 : if (b.isNan())
833 : 0 : flags |= NaN;
834 : 2738005 : m_hi = set_flags(m_hi, flags);
835 : 2738005 : if (isOverflow() || isNan())
836 : 0 : return *this;
837 : 2738005 : auto hi = get_num(m_hi);
838 : 2738005 : hi &= get_num(b.m_hi);
839 : 2738005 : m_lo &= b.m_lo;
840 : 2738005 : m_hi = set_flags(hi, flags);
841 : 2738005 : return *this;
842 : : }
843 : :
844 : : GncInt128&
845 : 0 : GncInt128::operator|= (const GncInt128& b) noexcept
846 : : {
847 : 0 : auto flags = get_flags(m_hi);
848 : 0 : auto hi = get_num(m_hi);
849 : 0 : hi ^= get_num(b.m_hi);
850 : 0 : m_hi = set_flags(hi, flags);
851 : 0 : m_lo ^= b.m_lo;
852 : 0 : return *this;
853 : : }
854 : :
855 : : GncInt128&
856 : 1213552 : GncInt128::operator^= (const GncInt128& b) noexcept
857 : : {
858 : 1213552 : auto flags = get_flags(m_hi);
859 : 1213552 : if (b.isOverflow())
860 : 0 : flags |= overflow;
861 : 1213552 : if (b.isNan())
862 : 0 : flags |= NaN;
863 : 1213552 : m_hi = set_flags(m_hi, flags);
864 : 1213552 : if (isOverflow() || isNan())
865 : 0 : return *this;
866 : 1213552 : auto hi = get_num(m_hi);
867 : 1213552 : hi ^= get_num(b.m_hi);
868 : 1213552 : m_hi = set_flags(hi, flags);
869 : 1213552 : m_lo ^= b.m_lo;
870 : 1213552 : return *this;
871 : : }
872 : :
873 : : static const uint8_t dec_array_size {5};
874 : : /* Convert a uint128 represented as a binary number into 4 10-digit decimal
875 : : * equivalents. Adapted from Douglas W. Jones, "Binary to Decimal Conversion in
876 : : * Limited Precision", http://homepage.cs.uiowa.edu/~jones/bcd/decimal.html,
877 : : * accessed 28 Oct 2014. */
878 : : static void
879 : 0 : decimal_from_binary (uint64_t d[dec_array_size], uint64_t hi, uint64_t lo)
880 : : {
881 : : /* Coefficients are the values of 2^96, 2^64, and 2^32 divided into 8-digit
882 : : * segments:
883 : : * 2^96 = 79228,16251426,43375935,43950336
884 : : * 2^64 = 1844,67440737,09551616
885 : : * 2^32 = 42,94967296
886 : : */
887 : 0 : const uint8_t coeff_array_size = dec_array_size - 1;
888 : 0 : const uint32_t coeff_3 [coeff_array_size] {79228, 16251426, 43375935, 43950336};
889 : 0 : const uint32_t coeff_2 [coeff_array_size] {0, 1844, 67440737, 9551616};
890 : 0 : const uint32_t coeff_1 [coeff_array_size] {0, 0, 42, 94967296};
891 : 0 : const uint32_t bin_mask {0xffffffff};
892 : 0 : const uint32_t dec_div {UINT32_C(100000000)};
893 : 0 : const uint8_t last {dec_array_size - 1};
894 : :
895 : 0 : d[0] = lo & bin_mask;
896 : 0 : d[1] = (lo >> 32) & bin_mask;
897 : 0 : d[2] = hi & bin_mask;
898 : 0 : d[3] = (hi >> 32) & bin_mask;
899 : :
900 : 0 : d[0] += coeff_3[3] * d[3] + coeff_2[3] * d[2] + coeff_1[3] * d[1];
901 : 0 : uint64_t q {d[0] / dec_div};
902 : 0 : d[0] %= dec_div;
903 : :
904 : 0 : for (int i {1}; i < coeff_array_size; ++i)
905 : : {
906 : 0 : int j = coeff_array_size - i - 1;
907 : 0 : d[i] = q + coeff_3[j] * d[3] + coeff_2[j] * d[2] + coeff_1[j] * d[1];
908 : 0 : q = d[i] / dec_div;
909 : 0 : d[i] %= dec_div;
910 : : }
911 : 0 : d[last] = q;
912 : 0 : return;
913 : : }
914 : :
915 : : static const uint8_t char_buf_size {41}; //39 digits plus sign and trailing null
916 : :
917 : : char*
918 : 0 : GncInt128::asCharBufR(char* buf, uint32_t size) const noexcept
919 : : {
920 : 0 : if (isOverflow())
921 : : {
922 : 0 : snprintf (buf, size, "%s", "Overflow");
923 : 0 : return buf;
924 : : }
925 : 0 : if (isNan())
926 : : {
927 : 0 : snprintf (buf, size, "%s", "NaN");
928 : 0 : return buf;
929 : : }
930 : 0 : if (isZero())
931 : : {
932 : 0 : snprintf (buf, size, "%d", 0);
933 : 0 : return buf;
934 : : }
935 : 0 : uint64_t d[dec_array_size] {};
936 : 0 : decimal_from_binary(d, get_num(m_hi), m_lo);
937 : 0 : char* next = buf;
938 : 0 : char neg {'-'};
939 : :
940 : 0 : if (isNeg()) *(next++) = neg;
941 : 0 : bool trailing {false};
942 : 0 : for (unsigned int i {dec_array_size}; i; --i)
943 : 0 : if (d[i - 1] || trailing)
944 : : {
945 : 0 : uint32_t new_size = size - (next - buf);
946 : 0 : if (trailing)
947 : 0 : next += snprintf (next, new_size, "%8.8" PRIu64, d[i - 1]);
948 : : else
949 : 0 : next += snprintf (next, new_size, "%" PRIu64, d[i - 1]);
950 : :
951 : 0 : trailing = true;
952 : : }
953 : :
954 : 0 : return buf;
955 : : }
956 : :
957 : : std::ostream&
958 : 0 : operator<< (std::ostream& stream, const GncInt128& a) noexcept
959 : : {
960 : 0 : char buf[char_buf_size] {};
961 : 0 : stream << a.asCharBufR (buf, char_buf_size - 1);
962 : 0 : return stream;
963 : : }
964 : :
965 : : bool
966 : 48223 : operator== (const GncInt128& a, const GncInt128& b) noexcept
967 : : {
968 : 48223 : return a.cmp(b) == 0;
969 : : }
970 : :
971 : : bool
972 : 761645 : operator!= (const GncInt128& a, const GncInt128& b) noexcept
973 : : {
974 : 761645 : return a.cmp(b) != 0;
975 : : }
976 : :
977 : : bool
978 : 1246457 : operator< (const GncInt128& a, const GncInt128& b) noexcept
979 : : {
980 : 1246457 : return a.cmp(b) < 0;
981 : : }
982 : :
983 : : bool
984 : 34213 : operator> (const GncInt128& a, const GncInt128& b) noexcept
985 : : {
986 : 34213 : return a.cmp(b) > 0;
987 : : }
988 : :
989 : : bool
990 : 0 : operator<= (const GncInt128& a, const GncInt128& b) noexcept
991 : : {
992 : 0 : return a.cmp(b) <= 0;
993 : : }
994 : :
995 : : bool
996 : 62 : operator>= (const GncInt128& a, const GncInt128& b) noexcept
997 : : {
998 : 62 : return a.cmp(b) >= 0;
999 : : }
1000 : :
1001 : : GncInt128
1002 : 257572 : operator+ (GncInt128 a, const GncInt128& b) noexcept
1003 : : {
1004 : 257572 : a += b;
1005 : 257572 : return a;
1006 : : }
1007 : :
1008 : : GncInt128
1009 : 399775 : operator- (GncInt128 a, const GncInt128& b) noexcept
1010 : : {
1011 : 399775 : a -= b;
1012 : 399775 : return a;
1013 : : }
1014 : : GncInt128
1015 : 886218 : operator* (GncInt128 a, const GncInt128& b) noexcept
1016 : : {
1017 : 886218 : a *= b;
1018 : 886218 : return a;
1019 : : }
1020 : :
1021 : : GncInt128
1022 : 980740 : operator/ (GncInt128 a, const GncInt128& b) noexcept
1023 : : {
1024 : 980740 : a /= b;
1025 : 980740 : return a;
1026 : : }
1027 : :
1028 : : GncInt128
1029 : 132555 : operator% (GncInt128 a, const GncInt128& b) noexcept
1030 : : {
1031 : 132555 : a %= b;
1032 : 132555 : return a;
1033 : : }
1034 : :
1035 : : GncInt128
1036 : 2738005 : operator& (GncInt128 a, const GncInt128& b) noexcept
1037 : : {
1038 : 2738005 : a &= b;
1039 : 2738005 : return a;
1040 : : }
1041 : :
1042 : : GncInt128
1043 : 0 : operator| (GncInt128 a, const GncInt128& b) noexcept
1044 : : {
1045 : 0 : a |= b;
1046 : 0 : return a;
1047 : : }
1048 : :
1049 : : GncInt128
1050 : 1213552 : operator^ (GncInt128 a, const GncInt128& b) noexcept
1051 : : {
1052 : 1213552 : a ^= b;
1053 : 1213552 : return a;
1054 : : }
1055 : :
1056 : : GncInt128
1057 : 361870 : operator<< (GncInt128 a, unsigned int b) noexcept
1058 : : {
1059 : 361870 : a <<= b;
1060 : 361870 : return a;
1061 : : }
1062 : :
1063 : : GncInt128
1064 : 67 : operator>> (GncInt128 a, unsigned int b) noexcept
1065 : : {
1066 : 67 : a >>= b;
1067 : 67 : return a;
1068 : : }
|
