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