mirror of
https://github.com/Gnucash/gnucash.git
synced 2025-02-25 18:55:30 -06:00
initial checkin of 18-bit integer math lib
git-svn-id: svn+ssh://svn.gnucash.org/repo/gnucash/trunk@10107 57a11ea4-9604-0410-9ed3-97b8803252fd
This commit is contained in:
Linas Vepstas
parent
fab4674a0e
commit
c24d87eeb3
346
src/engine/qofmath128.c
Normal file
346
src/engine/qofmath128.c
Normal file
@ -0,0 +1,346 @@
|
||||
/********************************************************************
|
||||
* qofmath128.c -- an 128-bit integer library *
|
||||
* Copyright (C) 2004 Linas Vepstas <linas@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@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 =================== */
|
||||
Loading…
Reference in New Issue
Block a user