[Gnucash-changes] initial checkin of 18-bit integer math lib

Linas Vepstas linas at cvs.gnucash.org
Sun Jun 27 00:05:08 EDT 2004


Log Message:
-----------
initial checkin of 18-bit integer math lib

Added Files:
-----------
    gnucash/src/engine:
        qofmath128.c

Revision Data
-------------
--- /dev/null
+++ src/engine/qofmath128.c
@@ -0,0 +1,346 @@
+/********************************************************************
+ * qofmath128.c -- an 128-bit integer library                       *
+ * Copyright (C) 2004 Linas Vepstas <linas at linas.org>               *
+ *                                                                  *
+ * This program is free software; you can redistribute it and/or    *
+ * modify it under the terms of the GNU General Public License as   *
+ * published by the Free Software Foundation; either version 2 of   *
+ * the License, or (at your option) any later version.              *
+ *                                                                  *
+ * This program is distributed in the hope that it will be useful,  *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of   *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the    *
+ * GNU General Public License for more details.                     *
+ *                                                                  *
+ * You should have received a copy of the GNU General Public License*
+ * along with this program; if not, contact:                        *
+ *                                                                  *
+ * Free Software Foundation           Voice:  +1-617-542-5942       *
+ * 59 Temple Place - Suite 330        Fax:    +1-617-542-2652       *
+ * Boston, MA  02111-1307,  USA       gnu at gnu.org                   *
+ *                                                                  *
+ *******************************************************************/
+
+#define _GNU_SOURCE
+
+#include "config.h"
+
+#include <glib.h>
+
+/* =============================================================== */
+/* Quick-n-dirty 128-bit math lib. The mult128 routine should work 
+ * great; I think that div128 works, but its not really tested.
+ */
+
+typedef struct {
+  guint64 hi;
+  guint64 lo;
+  short isneg;    /* sign-bit -- T if number is negative */
+  short isbig;    /* sizeflag -- T if number won't fit in signed 64-bit */
+} qofint128;
+
+/** Multiply a pair of signed 64-bit numbers, 
+ *  returning a signed 128-bit number.
+ */
+static inline qofint128
+mult128 (gint64 a, gint64 b)
+{
+  qofint128 prod;
+
+  prod.isneg = 0;
+  if (0>a)
+  {
+    prod.isneg = !prod.isneg;
+    a = -a;
+  }
+
+  if (0>b)
+  {
+    prod.isneg = !prod.isneg;
+    b = -b;
+  }
+
+  guint64 a1 = a >> 32;
+  guint64 a0 = a - (a1<<32);
+
+  guint64 b1 = b >> 32;
+  guint64 b0 = b - (b1<<32);
+
+  guint64 d = a0*b0;
+  guint64 d1 = d >> 32;
+  guint64 d0 = d - (d1<<32);
+
+  guint64 e = a0*b1;
+  guint64 e1 = e >> 32;
+  guint64 e0 = e - (e1<<32);
+
+  guint64 f = a1*b0;
+  guint64 f1 = f >> 32;
+  guint64 f0 = f - (f1<<32);
+
+  guint64 g = a1*b1;
+  guint64 g1 = g >> 32;
+  guint64 g0 = g - (g1<<32);
+
+  guint64 sum = d1+e0+f0;
+  guint64 carry = 0;
+  /* Can't say 1<<32 cause cpp will goof it up; 1ULL<<32 might work */
+  guint64 roll = 1<<30;
+  roll <<= 2;
+
+  guint64 pmax = roll-1;
+  while (pmax < sum)
+  {
+    sum -= roll;
+    carry ++;
+  }
+
+  prod.lo = d0 + (sum<<32);
+  prod.hi = carry + e1 + f1 + g0 + (g1<<32);
+  // prod.isbig = (prod.hi || (sum >> 31));
+  prod.isbig = prod.hi || (prod.lo >> 63);
+
+  return prod;
+}
+
+/** Divide a signed 128-bit number by a signed 64-bit,
+ *  returning a signed 128-bit number.
+ */
+static inline qofint128
+div128 (qofint128 n, gint64 d)
+{
+  qofint128 quotient;
+  guint64 hirem;   /* hi remainder */
+  guint64 qlo;
+
+  quotient.isneg = n.isneg;
+  if (0 > d)
+  {
+    d = -d;
+    quotient.isneg = !quotient.isneg;
+  }
+
+  quotient.hi = n.hi / d;
+  hirem = n.hi - quotient.hi * d;
+  
+  guint64 lo = 1<<30;
+  lo <<= 33;
+  lo /= d;
+  lo <<= 1;
+
+  lo *= hirem; 
+  quotient.lo = lo + n.lo/d;
+
+  /* Deal with low remainder bits.
+   * Is there a more efficient way of doing this?
+   */
+  qofint128 mu = mult128 (quotient.lo, d);
+
+  gint64 nn = 0x7fffffffffffffffULL & n.lo;
+  gint64 rr = 0x7fffffffffffffffULL & mu.lo;
+  gint64 rnd = nn - rr;
+  rnd /= d;   
+
+  /* ?? will this ever overflow ? */
+  qlo = quotient.lo;
+  quotient.lo += rnd;
+  if (qlo > quotient.lo)
+  {
+    quotient.hi += 1;
+  }
+
+  /* compute the carry situation */
+  quotient.isbig = (quotient.hi || (quotient.lo >> 63));
+
+  return quotient;
+}
+
+/** Return the remainder of a signed 128-bit number modulo 
+ *  a signed 64-bit.  That is, return n%d in 128-bit math.
+ *  I beleive that ths algo is overflow-free, but should be 
+ *  audited some more ... 
+ */
+static inline gint64
+rem128 (qofint128 n, gint64 d)
+{
+  qofint128 quotient = div128 (n,d);
+
+  qofint128 mu = mult128 (quotient.lo, d);
+
+  gint64 nn = 0x7fffffffffffffffULL & n.lo;
+  gint64 rr = 0x7fffffffffffffffULL & mu.lo;
+  return nn - rr;
+}
+
+/** Return the ratio n/d reduced so that there are no common factors. */
+static inline gnc_numeric
+reduce128(qofint128 n, gint64 d)
+{
+  gint64   t;
+  gint64   num;
+  gint64   denom;
+  gnc_numeric out;
+
+  t =  rem128 (n, d);
+  num = d;
+  denom = t;
+
+  /* The strategy is to use Euclid's algorithm */
+  while (denom > 0) 
+  {
+    t = num % denom;
+    num = denom;
+    denom = t;
+  }
+  /* num now holds the GCD (Greatest Common Divisor) */
+
+  qofint128 red = div128 (n, num);
+  if (red.isbig)
+  {
+    return gnc_numeric_error (GNC_ERROR_OVERFLOW);
+  }
+  out.num   = red.lo;
+  if (red.isneg) out.num = -out.num;
+  out.denom = d / num;
+  return out;
+}
+
+/** Return true of two numbers are equal */
+static inline gboolean
+equal128 (qofint128 a, qofint128 b)
+{
+	if (a.lo != b.lo) return 0;
+	if (a.hi != b.hi) return 0;
+	if (a.isneg != b.isneg) return 0;
+	return 1;
+}
+
+/** Return the greatest common factor of two 64-bit numbers */
+static inline guint64
+gcf64(guint64 num, guint64 denom)
+{
+  guint64   t;
+
+  t =  num % denom;
+  num = denom;
+  denom = t;
+
+  /* The strategy is to use Euclid's algorithm */
+  while (0 != denom) 
+  {
+    t = num % denom;
+    num = denom;
+    denom = t;
+  }
+  /* num now holds the GCD (Greatest Common Divisor) */
+  return num;
+}
+
+/** Return the least common multiple of two 64-bit numbers. */
+static inline qofint128
+lcm128 (guint64 a, guint64 b)
+{
+  guint64 gcf = gcf64 (a,b);
+  b /= gcf;
+  return mult128 (a,b);
+}
+
+/** Add a pair of 128-bit numbers, returning a 128-bit number */
+static inline qofint128
+add128 (qofint128 a, qofint128 b)
+{
+  qofint128 sum;
+  if (a.isneg == b.isneg)
+  {
+    sum.isneg = a.isneg;
+    sum.hi = a.hi + b.hi;
+    sum.lo = a.lo + b.lo;
+    if ((sum.lo < a.lo) || (sum.lo < b.lo))
+    {
+     sum.hi ++;
+    }
+    sum.isbig = sum.hi || (sum.lo >> 63);
+    return sum;
+  }
+  if ((b.hi > a.hi) ||
+     ((b.hi == a.hi) && (b.lo > a.lo)))
+  {
+    qofint128 tmp = a;
+    a = b;
+    b = tmp;
+  }
+
+  sum.isneg = a.isneg;
+  sum.hi = a.hi - b.hi;
+  sum.lo = a.lo - b.lo;
+
+  if (sum.lo > a.lo)
+  {
+    sum.hi --;
+  }
+
+  sum.isbig = sum.hi || (sum.lo >> 63);
+  return sum;
+}
+
+#ifdef TEST_128_BIT_MULT
+static void pr (gint64 a, gint64 b)
+{
+   qofint128 prod = mult128 (a,b);
+   printf ("%lld * %lld = %lld %llu (0x%llx %llx) %hd\n",
+	   a, b, prod.hi, prod.lo, prod.hi, prod.lo, prod.isbig);
+}
+
+static void prd (gint64 a, gint64 b, gint64 c)
+{
+   qofint128 prod = mult128 (a,b);
+   qofint128 quot = div128 (prod, c);
+   gint64 rem = rem128 (prod, c);
+   printf ("%lld * %lld / %lld = %lld %llu + %lld (0x%llx %llx) %hd\n",
+	   a, b, c, quot.hi, quot.lo, rem, quot.hi, quot.lo, quot.isbig);
+}
+
+int main ()
+{
+  pr (2,2);
+
+  gint64 x = 1<<30;
+  x <<= 2;
+
+  pr (x,x);
+  pr (x+1,x);
+  pr (x+1,x+1);
+
+  pr (x,-x);
+  pr (-x,-x);
+  pr (x-1,x);
+  pr (x-1,x-1);
+  pr (x-2,x-2);
+
+  x <<= 1;
+  pr (x,x);
+  pr (x,-x);
+
+  pr (1000000, 10000000000000);
+
+  prd (x,x,2);
+  prd (x,x,3);
+  prd (x,x,4);
+  prd (x,x,5);
+  prd (x,x,6);
+
+  x <<= 29;
+  prd (3,x,3);
+  prd (6,x,3);
+  prd (99,x,3);
+  prd (100,x,5);
+  prd (540,x,5);
+  prd (777,x,7);
+  prd (1111,x,11);
+
+  return 0;
+}
+
+#endif /* TEST_128_BIT_MULT */
+
+/* ======================== END OF FILE =================== */


More information about the gnucash-changes mailing list