[Gnucash-changes] 128-bit division algorithm was badly broken,
replace it with standard
Linas Vepstas
linas at cvs.gnucash.org
Mon Jun 28 23:45:09 EDT 2004
Log Message:
-----------
128-bit division algorithm was badly broken,
replace it with standard grade-school long division
Modified Files:
--------------
gnucash/src/engine:
qofmath128.c
qofmath128.h
Revision Data
-------------
Index: qofmath128.h
===================================================================
RCS file: /home/cvs/cvsroot/gnucash/src/engine/qofmath128.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -Lsrc/engine/qofmath128.h -Lsrc/engine/qofmath128.h -u -r1.2 -r1.3
--- src/engine/qofmath128.h
+++ src/engine/qofmath128.h
@@ -39,6 +39,21 @@
short isbig; /**< sizeflag -- T if number won't fit in signed 64-bit */
} qofint128;
+/** Return true of two numbers are equal */
+inline gboolean equal128 (qofint128 a, qofint128 b);
+
+/** Shift right by one bit (i.e. divide by two) */
+inline qofint128 shift128 (qofint128 x);
+
+/** Shift left by one bit (i.e. multiply by two) */
+inline qofint128 shiftleft128 (qofint128 x);
+
+/** Increment by one */
+inline qofint128 inc128 (qofint128 a);
+
+/** Add a pair of 128-bit numbers, returning a 128-bit number */
+inline qofint128 add128 (qofint128 a, qofint128 b);
+
/** Multiply a pair of signed 64-bit numbers,
* returning a signed 128-bit number.
*/
@@ -56,24 +71,12 @@
*/
inline gint64 rem128 (qofint128 n, gint64 d);
-/** Return the ratio n/d reduced so that there are no common factors. */
-inline gnc_numeric reduce128(qofint128 n, gint64 d);
-
-/** Return true of two numbers are equal */
-inline gboolean equal128 (qofint128 a, qofint128 b);
-
/** Return the greatest common factor of two 64-bit numbers */
inline guint64 gcf64(guint64 num, guint64 denom);
/** Return the least common multiple of two 64-bit numbers. */
inline qofint128 lcm128 (guint64 a, guint64 b);
-/** Add a pair of 128-bit numbers, returning a 128-bit number */
-inline qofint128 add128 (qofint128 a, qofint128 b);
-
-/** Shift right by one bit (i.e. divide by two) */
-inline qofint128 shift128 (qofint128 x);
-
#endif
/** @} */
Index: qofmath128.c
===================================================================
RCS file: /home/cvs/cvsroot/gnucash/src/engine/qofmath128.c,v
retrieving revision 1.4
retrieving revision 1.5
diff -Lsrc/engine/qofmath128.c -Lsrc/engine/qofmath128.c -u -r1.4 -r1.5
--- src/engine/qofmath128.c
+++ src/engine/qofmath128.c
@@ -34,6 +34,8 @@
* work, and have been tested, but not comprehensively tested.
*/
+#define HIBIT (0x8000000000000000ULL)
+
/** Multiply a pair of signed 64-bit numbers,
* returning a signed 128-bit number.
*/
@@ -98,6 +100,73 @@
return prod;
}
+/** Shift right by one bit (i.e. divide by two) */
+inline qofint128
+shift128 (qofint128 x)
+{
+ guint64 sbit = x.hi & 0x1;
+ x.hi >>= 1;
+ x.lo >>= 1;
+ x.isbig = 0;
+ if (sbit)
+ {
+ x.lo |= HIBIT;
+ x.isbig = 1;
+ return x;
+ }
+ if (x.hi)
+ {
+ x.isbig = 1;
+ }
+ return x;
+}
+
+/** Shift leftt by one bit (i.e. multiply by two) */
+inline qofint128
+shiftleft128 (qofint128 x)
+{
+ guint64 sbit = x.lo & HIBIT;
+ x.hi <<= 1;
+ x.lo <<= 1;
+ x.isbig = 0;
+ if (sbit)
+ {
+ x.hi |= 1;
+ x.isbig = 1;
+ return x;
+ }
+ if (x.hi)
+ {
+ x.isbig = 1;
+ }
+ return x;
+}
+
+/** increment a 128-bit number by one */
+inline qofint128
+inc128 (qofint128 a)
+{
+ if (0 == a.isneg)
+ {
+ a.lo ++;
+ if (0 == a.lo)
+ {
+ a.hi ++;
+ }
+ }
+ else
+ {
+ if (0 == a.lo)
+ {
+ a.hi --;
+ }
+ a.lo --;
+ }
+
+ a.isbig = (a.hi != 0) || (a.lo >> 63);
+ return a;
+}
+
/** Divide a signed 128-bit number by a signed 64-bit,
* returning a signed 128-bit number.
*/
@@ -105,43 +174,28 @@
div128 (qofint128 n, gint64 d)
{
qofint128 quotient;
- guint64 hirem; /* hi remainder */
- guint64 qlo;
+ guint64 remainder = 0;
- quotient.isneg = n.isneg;
+ quotient = n;
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;
+ /* Use grade-school long division algorithm */
+ int i;
+ for (i=0; i<128; i++)
+ {
+ guint64 sbit = HIBIT & quotient.hi;
+ remainder <<= 1;
+ if (sbit) remainder |= 1;
+ quotient = shiftleft128 (quotient);
+ if (remainder >= d)
+ {
+ remainder -= d;
+ quotient.lo |= 1;
+ }
}
/* compute the carry situation */
@@ -167,39 +221,6 @@
return nn - rr;
}
-/** Return the ratio n/d reduced so that there are no common factors. */
-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 */
inline gboolean
equal128 (qofint128 a, qofint128 b)
@@ -278,28 +299,6 @@
return sum;
}
-/** Shift right by one bit (i.e. divide by two) */
-inline qofint128
-shift128 (qofint128 x)
-{
- guint64 sbit = x.hi & 0x1;
- x.hi >>= 1;
- x.lo >>= 1;
- x.isbig = 0;
- if (sbit)
- {
- sbit = 1<<30; /* in two step to avoid 1ULL<<63 */
- sbit <<= 33;
- x.lo |= sbit;
- x.isbig = 1;
- return x;
- }
- if (x.hi)
- {
- x.isbig = 1;
- }
- return x;
-}
#ifdef TEST_128_BIT_MULT
static void pr (gint64 a, gint64 b)
@@ -356,6 +355,21 @@
prd (777,x,7);
prd (1111,x,11);
+ /* Really test division */
+ qofint128 n;
+ n.hi = 0xdd91;
+ n.lo = 0x6c5abefbb9e13480ULL;
+
+ gint64 d = 0x2ae79964d3ae1d04ULL;
+
+ int i;
+ for (i=0; i<20; i++) {
+
+ qofint128 quot = div128 (n, d);
+ printf ("%d result = %llx %llx\n", i, quot.hi, quot.lo);
+ d >>=1;
+ n = shift128 (n);
+ }
return 0;
}
More information about the gnucash-changes
mailing list