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 14841620 : static inline uint64_t set_flags(uint64_t leg, uint8_t flags)
50 : {
51 14841620 : auto flag_part = static_cast<uint64_t>(flags) << upper_num_bits;
52 14841620 : return flag_part + (leg & nummask);
53 : }
54 46580379 : static inline uint8_t get_flags(uint64_t leg)
55 : {
56 46580379 : return (leg & flagmask) >> upper_num_bits;
57 : }
58 20873466 : static inline uint64_t get_num(uint64_t leg)
59 : {
60 20873466 : return leg & nummask;
61 : }
62 : }
63 :
64 1544160 : GncInt128::GncInt128 () : m_hi {0}, m_lo {0}{}
65 :
66 2555229 : GncInt128::GncInt128 (int64_t upper, int64_t lower, uint8_t flags) :
67 2555229 : m_hi {static_cast<uint64_t>(upper < 0 ? -upper : upper)},
68 2555229 : m_lo {static_cast<uint64_t>(lower < 0 ? -lower : lower)}
69 : {
70 2555229 : if ((upper < 0 && lower > 0) || (upper > 0 && lower < 0))
71 0 : m_lo = (m_hi << 63) - m_lo;
72 : else
73 2555229 : m_lo += (m_hi << 63);
74 :
75 2555229 : m_hi >>= 1;
76 2555229 : 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 2555229 : flags ^= (upper < 0 ? neg : upper == 0 && lower < 0 ? neg : pos);
84 2555229 : m_hi = set_flags(m_hi, flags);
85 2555229 : }
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 2293609 : GncInt128::GncInt128 (uint64_t upper, uint64_t lower, uint8_t flags) :
102 2293609 : 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 2293609 : if (m_hi == UINT64_MAX)
109 1624 : m_hi = nummask;
110 2293609 : 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 2293609 : m_hi = set_flags(m_hi, flags);
119 2293609 : }
120 :
121 : GncInt128&
122 1544160 : GncInt128::zero () noexcept
123 : {
124 1544160 : m_lo = m_hi = UINT64_C(0);
125 1544160 : return *this;
126 : }
127 :
128 500750 : GncInt128::operator int64_t() const
129 : {
130 500750 : auto flags = get_flags(m_hi);
131 500750 : if ((flags & neg) && isBig())
132 0 : throw std::underflow_error ("Negative value too large to represent as int64_t");
133 500750 : if ((flags & (overflow | NaN)) || isBig())
134 15 : throw std::overflow_error ("Value too large to represent as int64_t");
135 500735 : int64_t retval = static_cast<int64_t>(m_lo);
136 500735 : 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 1495570 : GncInt128::cmp (const GncInt128& b) const noexcept
152 : {
153 1495570 : auto flags = get_flags(m_hi);
154 1495570 : if (flags & (overflow | NaN))
155 0 : return -1;
156 1495570 : if (b.isOverflow () || b.isNan ())
157 0 : return 1;
158 1495570 : auto hi = get_num(m_hi);
159 1495570 : auto bhi = get_num(b.m_hi);
160 1495570 : if (isZero() && b.isZero()) return 0;
161 1482740 : if (flags & neg)
162 : {
163 757 : if (!b.isNeg()) return -1;
164 214 : if (hi > bhi) return -1;
165 214 : if (hi < bhi) return 1;
166 214 : if (m_lo > b.m_lo) return -1;
167 82 : if (m_lo < b.m_lo) return 1;
168 12 : return 0;
169 : }
170 1481983 : if (b.isNeg()) return 1;
171 1481981 : if (hi < bhi) return -1;
172 1480669 : if (hi > bhi) return 1;
173 1473801 : if (m_lo < b.m_lo) return -1;
174 1231804 : if (m_lo > b.m_lo) return 1;
175 523530 : 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 235652 : GncInt128::gcd(GncInt128 b) const noexcept
183 : {
184 235652 : if (b.isZero())
185 268 : return *this;
186 235384 : if (isZero())
187 0 : return b;
188 :
189 235384 : if (b.isOverflow() || b.isNan())
190 0 : return b;
191 235384 : if (isOverflow() || isNan())
192 0 : return *this;
193 :
194 235384 : GncInt128 a (isNeg() ? -(*this) : *this);
195 235384 : if (b.isNeg()) b = -b;
196 :
197 235384 : unsigned int k {};
198 235384 : const uint64_t one {1};
199 535963 : while (!((a & one) || (b & one))) //B1
200 : {
201 300579 : a >>= 1;
202 300579 : b >>= 1;
203 300579 : ++k;
204 : }
205 235384 : GncInt128 t {a & one ? -b : a}; //B2
206 429432 : while (a != b)
207 : {
208 597399 : while (t && ((t & one) ^ one)) t >>= 1; //B3 & B4
209 194048 : if (t.isNeg()) //B5
210 119930 : b = -t;
211 : else
212 74118 : a = t;
213 194048 : t = a - b; //B6
214 : }
215 235384 : 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 198070 : GncInt128::lcm(const GncInt128& b) const noexcept
222 : {
223 198070 : auto common = gcd(b);
224 198070 : 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 8284020 : GncInt128::isNeg () const noexcept
248 : {
249 8284020 : return (get_flags(m_hi) & neg);
250 : }
251 :
252 : bool
253 1049204 : GncInt128::isBig () const noexcept
254 : {
255 1049204 : return (get_num(m_hi) || m_lo > INT64_MAX);
256 : }
257 :
258 : bool
259 10621227 : GncInt128::isOverflow () const noexcept
260 : {
261 10621227 : return (get_flags(m_hi) & overflow);
262 : }
263 :
264 : bool
265 10671482 : GncInt128::isNan () const noexcept
266 : {
267 10671482 : return (get_flags(m_hi) & NaN);
268 : }
269 :
270 : bool
271 1236106 : GncInt128::valid() const noexcept
272 : {
273 1236106 : return !(get_flags(m_hi) & (overflow | NaN));
274 : }
275 :
276 : bool
277 6347411 : GncInt128::isZero() const noexcept
278 : {
279 12694822 : return ((get_flags(m_hi) & (overflow | NaN)) == 0 &&
280 12694822 : get_num(m_hi) == 0 && m_lo == 0);
281 : }
282 :
283 : GncInt128
284 2220810 : GncInt128::abs() const noexcept
285 : {
286 2220810 : if (isNeg())
287 260633 : return operator-();
288 :
289 1960177 : return *this;
290 : }
291 :
292 : unsigned int
293 1247692 : GncInt128::bits() const noexcept
294 : {
295 1247692 : auto hi = get_num(m_hi);
296 1247692 : unsigned int bits {static_cast<unsigned int>(hi == 0 ? 0 : 64)};
297 1247692 : uint64_t temp {(hi == 0 ? m_lo : hi)};
298 10962787 : for (;temp > 0; temp >>= 1)
299 9715095 : ++bits;
300 1247692 : return bits;
301 : }
302 :
303 :
304 : GncInt128
305 646748 : GncInt128::operator-() const noexcept
306 : {
307 646748 : auto retval = *this;
308 646748 : auto flags = get_flags(retval.m_hi);
309 646748 : if (isNeg())
310 414291 : flags ^= neg;
311 : else
312 232457 : flags |= neg;
313 646748 : retval.m_hi = set_flags(retval.m_hi, flags);
314 646748 : return retval;
315 : }
316 :
317 2291801 : GncInt128::operator bool() const noexcept
318 : {
319 2291801 : 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 197725 : GncInt128::operator+= (const GncInt128& b) noexcept
348 : {
349 197725 : auto flags = get_flags(m_hi);
350 197725 : if (b.isOverflow())
351 0 : flags |= overflow;
352 197725 : if (b.isNan())
353 0 : flags |= NaN;
354 197725 : m_hi = set_flags(m_hi, flags);
355 197725 : if (isOverflow() || isNan())
356 0 : return *this;
357 197725 : if ((isNeg () && !b.isNeg ()) || (!isNeg () && b.isNeg ()))
358 44549 : return this->operator-= (-b);
359 153176 : uint64_t result = m_lo + b.m_lo;
360 153176 : uint64_t carry = static_cast<int64_t>(result < m_lo); //Wrapping
361 153176 : m_lo = result;
362 153176 : auto hi = get_num(m_hi);
363 153176 : auto bhi = get_num(b.m_hi);
364 153176 : result = hi + bhi + carry;
365 153176 : if (result < hi || result & flagmask)
366 0 : flags |= overflow;
367 153176 : m_hi = set_flags(result, flags);
368 153176 : return *this;
369 : }
370 :
371 : GncInt128&
372 235384 : GncInt128::operator<<= (unsigned int i) noexcept
373 : {
374 235384 : auto flags = get_flags(m_hi);
375 235384 : if (i == 0)
376 89567 : return *this;
377 145817 : 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 145817 : auto hi = get_num(m_hi);
385 145817 : if (i < legbits)
386 : {
387 145817 : uint64_t carry {(m_lo & (((UINT64_C(1) << i) - 1) << (legbits - i))) >> (legbits - i)};
388 145817 : m_lo <<= i;
389 145817 : hi <<= i;
390 145817 : hi += carry;
391 145817 : m_hi = set_flags(hi, flags);
392 145817 : 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 1004576 : GncInt128::operator>>= (unsigned int i) noexcept
402 : {
403 1004576 : auto flags = get_flags(m_hi);
404 1004576 : 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 1004576 : auto hi = get_num(m_hi);
412 1004576 : if (i < legbits)
413 : {
414 1004548 : uint64_t carry {(hi & ((UINT64_C(1) << i) - 1))};
415 1004548 : m_lo >>= i;
416 1004548 : hi >>= i;
417 1004548 : m_lo += (carry << (legbits - i));
418 1004548 : m_hi = set_flags(hi, flags);
419 1004548 : 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 238597 : GncInt128::operator-= (const GncInt128& b) noexcept
428 : {
429 238597 : auto flags = get_flags(m_hi);
430 238597 : if (b.isOverflow())
431 0 : flags |= overflow;
432 238597 : if (b.isNan())
433 0 : flags |= NaN;
434 238597 : m_hi = set_flags(m_hi, flags);
435 :
436 238597 : if (isOverflow() || isNan())
437 0 : return *this;
438 :
439 238597 : if ((!isNeg() && b.isNeg()) || (isNeg() && !b.isNeg()))
440 0 : return this->operator+= (-b);
441 238597 : bool operand_bigger {abs().cmp (b.abs()) < 0};
442 238597 : auto hi = get_num(m_hi);
443 238597 : auto far_hi = get_num(b.m_hi);
444 238597 : if (operand_bigger)
445 : {
446 94673 : flags ^= neg; // ^= flips the bit
447 94673 : 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 94673 : m_lo = b.m_lo - m_lo;
457 94673 : hi = far_hi - hi;
458 94673 : m_hi = set_flags(hi, flags);
459 94673 : return *this;
460 : }
461 143924 : 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 143924 : m_lo -= b.m_lo;
468 :
469 143924 : hi -= far_hi;
470 143924 : m_hi = set_flags(hi, flags);
471 143924 : return *this;
472 : }
473 :
474 : GncInt128&
475 637902 : GncInt128::operator*= (const GncInt128& b) noexcept
476 : {
477 : /* Test for 0 first */
478 637902 : auto flags = get_flags(m_hi);
479 : /* Handle the sign; ^ flips if b is negative */
480 637902 : flags ^= (get_flags(b.m_hi) & neg);
481 637902 : if (isZero() || b.isZero())
482 : {
483 14056 : m_lo = 0;
484 14056 : m_hi = set_flags(0, flags);
485 14056 : return *this;
486 : }
487 623846 : if (b.isOverflow())
488 0 : flags |= overflow;
489 623846 : if (b.isNan())
490 0 : flags |= NaN;
491 623846 : m_hi = set_flags(m_hi, flags);
492 623846 : if (isOverflow() || isNan())
493 0 : return *this;
494 :
495 : /* Test for overflow before spending time on the calculation */
496 623846 : auto hi = get_num(m_hi);
497 623846 : auto bhi = get_num(b.m_hi);
498 623846 : if (hi && bhi)
499 : {
500 0 : flags |= overflow;
501 0 : m_hi = set_flags(hi, flags);
502 0 : return *this;
503 : }
504 :
505 623846 : 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 623846 : 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 623846 : if (abits + bbits <= legbits)
520 : {
521 621472 : m_lo *= b.m_lo;
522 621472 : m_hi = set_flags(m_hi, flags);
523 621472 : 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 2374 : uint64_t av[sublegs] {(m_lo & sublegmask), (m_lo >> sublegbits),
539 2374 : (hi & sublegmask), (hi >> sublegbits)};
540 2374 : uint64_t bv[sublegs] {(b.m_lo & sublegmask), (b.m_lo >> sublegbits),
541 2374 : (bhi & sublegmask), (bhi >> sublegbits)};
542 2374 : uint64_t rv[sublegs] {};
543 2374 : uint64_t carry {}, scratch {};
544 :
545 2374 : rv[0] = av[0] * bv[0];
546 :
547 2374 : rv[1] = av[1] * bv [0];
548 2374 : scratch = rv[1] + av[0] * bv[1];
549 2374 : carry = rv[1] > scratch ? 1 : 0;
550 2374 : rv[1] = scratch;
551 :
552 2374 : rv[2] = av[2] * bv[0] + carry; //0xffff^2 + 1 < 0xffffffff, can't overflow
553 2374 : scratch = rv[2] + av[1] * bv[1];
554 2374 : carry = rv[2] > scratch ? 1 : 0;
555 2374 : rv[2] = scratch + av[0] * bv[2];
556 2374 : carry += scratch > rv[2] ? 1 : 0;
557 :
558 2374 : rv[3] = av[3] * bv[0] + carry;
559 2374 : scratch = rv[3] + av[2] * bv[1];
560 2374 : carry = rv[3] > scratch ? 1 : 0;
561 2374 : rv[3] = scratch + av[1] * bv[2];
562 2374 : carry += scratch > rv[3] ? 1 : 0;
563 2374 : scratch = rv[3] + av[0] * bv[3];
564 2374 : carry += rv[3] > scratch ? 1 : 0;
565 2374 : rv[3] = scratch;
566 :
567 2374 : 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 2374 : m_lo = rv[0] + (rv[1] << sublegbits);
575 2374 : carry = rv[1] >> sublegbits;
576 2374 : carry += (rv[1] << sublegbits) > m_lo || rv[0] > m_lo ? 1 : 0;
577 2374 : hi = rv[2] + (rv[3] << sublegbits) + carry;
578 2374 : 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 2374 : m_hi = set_flags(hi, flags);
586 2374 : 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 5460 : div_single_leg (uint64_t* u, size_t m, uint64_t v,
697 : GncInt128& q, GncInt128& r) noexcept
698 : {
699 5460 : uint64_t qv[sublegs] {};
700 5460 : bool negative {q.isNeg()};
701 5460 : bool rnegative {r.isNeg()};
702 21854 : for (int i = m - 1; i >= 0; --i)
703 : {
704 16394 : qv[i] = u[i] / v;
705 16394 : if (i > 0)
706 : {
707 10934 : u[i - 1] += ((u[i] % v) << sublegbits);
708 10934 : u[i] = UINT64_C(0);
709 : }
710 : else
711 5460 : u[i] %= v;
712 : }
713 :
714 5460 : q = GncInt128 ((qv[3] << sublegbits) + qv[2], (qv[1] << sublegbits) + qv[0]);
715 5460 : r = GncInt128 ((u[3] << sublegbits) + u[2], (u[1] << sublegbits) + u[0]);
716 5460 : if (negative) q = -q;
717 5460 : if (rnegative) r = -r;
718 5460 : }
719 :
720 : }// namespace
721 :
722 : void
723 772080 : GncInt128::div (const GncInt128& b, GncInt128& q, GncInt128& r) const noexcept
724 : {
725 : // clear remainder and quotient
726 772080 : r = GncInt128(0);
727 772080 : q = GncInt128(0);
728 :
729 772080 : auto qflags = get_flags(q.m_hi);
730 772080 : auto rflags = get_flags(r.m_hi);
731 772080 : 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 772080 : 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 772080 : assert (&q != this);
749 772080 : assert (&r != this);
750 772080 : assert (&q != &b);
751 772080 : assert (&r != &b);
752 :
753 772080 : q.zero(), r.zero();
754 772080 : qflags = rflags = 0;
755 772080 : 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 772080 : if (isNeg())
765 : {
766 238330 : qflags |= neg;
767 238330 : rflags |= neg; // the remainder inherits the dividend's sign
768 : }
769 :
770 772080 : if (b.isNeg())
771 0 : qflags ^= neg;
772 :
773 772080 : q.m_hi = set_flags(q.m_hi, qflags);
774 772080 : r.m_hi = set_flags(r.m_hi, rflags);
775 :
776 772080 : if (abs() < b.abs())
777 : {
778 24705 : r = *this;
779 24705 : return;
780 : }
781 747375 : auto hi = get_num(m_hi);
782 747375 : auto bhi = get_num(b.m_hi);
783 :
784 747375 : if (hi == 0 && bhi == 0) //let the hardware do it
785 : {
786 741833 : assert (b.m_lo != 0); // b.m_hi is 0 but b isn't or we didn't get here.
787 741833 : q.m_lo = m_lo / b.m_lo;
788 741833 : r.m_lo = m_lo % b.m_lo;
789 741833 : return;
790 : }
791 :
792 5542 : uint64_t u[sublegs + 2] {(m_lo & sublegmask), (m_lo >> sublegbits),
793 5542 : (hi & sublegmask), (hi >> sublegbits), 0, 0};
794 5542 : uint64_t v[sublegs] {(b.m_lo & sublegmask), (b.m_lo >> sublegbits),
795 5542 : (bhi & sublegmask), (bhi >> sublegbits)};
796 5542 : auto m = u[3] ? 4 : u[2] ? 3 : u[1] ? 2 : u[0] ? 1 : 0;
797 5542 : auto n = v[3] ? 4 : v[2] ? 3 : v[1] ? 2 : v[0] ? 1 : 0;
798 5542 : if (m == 0 || n == 0) //Shouldn't happen
799 0 : return;
800 5542 : if (n == 1)
801 5460 : 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 718913 : GncInt128::operator/= (const GncInt128& b) noexcept
808 : {
809 718913 : GncInt128 q {}, r {};
810 718913 : div(b, q, r);
811 718913 : std::swap (*this, q);
812 718913 : return *this;
813 : }
814 :
815 : GncInt128&
816 53158 : GncInt128::operator%= (const GncInt128& b) noexcept
817 : {
818 53158 : GncInt128 q {}, r {};
819 53158 : div(b, q, r);
820 53158 : std::swap (*this, r);
821 53158 : if (q.isNan())
822 0 : m_hi = set_flags(m_hi, (get_flags(m_hi) | NaN));
823 53158 : return *this;
824 : }
825 :
826 : GncInt128&
827 1683420 : GncInt128::operator&= (const GncInt128& b) noexcept
828 : {
829 1683420 : auto flags = get_flags(m_hi);
830 1683420 : if (b.isOverflow())
831 0 : flags |= overflow;
832 1683420 : if (b.isNan())
833 0 : flags |= NaN;
834 1683420 : m_hi = set_flags(m_hi, flags);
835 1683420 : if (isOverflow() || isNan())
836 0 : return *this;
837 1683420 : auto hi = get_num(m_hi);
838 1683420 : hi &= get_num(b.m_hi);
839 1683420 : m_lo &= b.m_lo;
840 1683420 : m_hi = set_flags(hi, flags);
841 1683420 : 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 597399 : GncInt128::operator^= (const GncInt128& b) noexcept
857 : {
858 597399 : auto flags = get_flags(m_hi);
859 597399 : if (b.isOverflow())
860 0 : flags |= overflow;
861 597399 : if (b.isNan())
862 0 : flags |= NaN;
863 597399 : m_hi = set_flags(m_hi, flags);
864 597399 : if (isOverflow() || isNan())
865 0 : return *this;
866 597399 : auto hi = get_num(m_hi);
867 597399 : hi ^= get_num(b.m_hi);
868 597399 : m_hi = set_flags(hi, flags);
869 597399 : m_lo ^= b.m_lo;
870 597399 : 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 30746 : operator== (const GncInt128& a, const GncInt128& b) noexcept
967 : {
968 30746 : return a.cmp(b) == 0;
969 : }
970 :
971 : bool
972 429432 : operator!= (const GncInt128& a, const GncInt128& b) noexcept
973 : {
974 429432 : return a.cmp(b) != 0;
975 : }
976 :
977 : bool
978 776246 : operator< (const GncInt128& a, const GncInt128& b) noexcept
979 : {
980 776246 : return a.cmp(b) < 0;
981 : }
982 :
983 : bool
984 20477 : operator> (const GncInt128& a, const GncInt128& b) noexcept
985 : {
986 20477 : 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 72 : operator>= (const GncInt128& a, const GncInt128& b) noexcept
997 : {
998 72 : return a.cmp(b) >= 0;
999 : }
1000 :
1001 : GncInt128
1002 197716 : operator+ (GncInt128 a, const GncInt128& b) noexcept
1003 : {
1004 197716 : a += b;
1005 197716 : return a;
1006 : }
1007 :
1008 : GncInt128
1009 194048 : operator- (GncInt128 a, const GncInt128& b) noexcept
1010 : {
1011 194048 : a -= b;
1012 194048 : return a;
1013 : }
1014 : GncInt128
1015 637902 : operator* (GncInt128 a, const GncInt128& b) noexcept
1016 : {
1017 637902 : a *= b;
1018 637902 : return a;
1019 : }
1020 :
1021 : GncInt128
1022 667851 : operator/ (GncInt128 a, const GncInt128& b) noexcept
1023 : {
1024 667851 : a /= b;
1025 667851 : return a;
1026 : }
1027 :
1028 : GncInt128
1029 53158 : operator% (GncInt128 a, const GncInt128& b) noexcept
1030 : {
1031 53158 : a %= b;
1032 53158 : return a;
1033 : }
1034 :
1035 : GncInt128
1036 1683420 : operator& (GncInt128 a, const GncInt128& b) noexcept
1037 : {
1038 1683420 : a &= b;
1039 1683420 : 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 597399 : operator^ (GncInt128 a, const GncInt128& b) noexcept
1051 : {
1052 597399 : a ^= b;
1053 597399 : return a;
1054 : }
1055 :
1056 : GncInt128
1057 235384 : operator<< (GncInt128 a, unsigned int b) noexcept
1058 : {
1059 235384 : a <<= b;
1060 235384 : 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 : }
|
