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