QOF  0.8.7
qofmath128.c
1 /********************************************************************
2  * qofmath128.c -- an 128-bit integer library *
3  * Copyright (C) 2004 Linas Vepstas <linas@linas.org> *
4  * *
5  * This program is free software; you can redistribute it and/or *
6  * modify it under the terms of the GNU General Public License as *
7  * published by the Free Software Foundation; either version 2 of *
8  * the License, or (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public License*
16  * along with this program; if not, contact: *
17  * *
18  * Free Software Foundation Voice: +1-617-542-5942 *
19  * 51 Franklin Street, Fifth Floor Fax: +1-617-542-2652 *
20  * Boston, MA 02110-1301, USA gnu@gnu.org *
21  * *
22  *******************************************************************/
23 
24 #include "config.h"
25 #include "qofmath128.h"
26 
27 #include <glib.h>
28 
29 /* =============================================================== */
30 /*
31  * Quick-n-dirty 128-bit integer math lib. Things seem to mostly
32  * work, and have been tested, but not comprehensively tested.
33  */
34 
35 #define HIBIT (0x8000000000000000ULL)
36 
41 mult128 (gint64 a, gint64 b)
42 {
43  QofInt128 prod;
44  guint64 a0, a1;
45  guint64 b0, b1;
46  guint64 d, d0, d1;
47  guint64 e, e0, e1;
48  guint64 f, f0, f1;
49  guint64 g, g0, g1;
50  guint64 sum, carry, roll, pmax;
51 
52  prod.isneg = 0;
53  if (0 > a)
54  {
55  prod.isneg = !prod.isneg;
56  a = -a;
57  }
58 
59  if (0 > b)
60  {
61  prod.isneg = !prod.isneg;
62  b = -b;
63  }
64 
65  a1 = a >> 32;
66  a0 = a - (a1 << 32);
67 
68  b1 = b >> 32;
69  b0 = b - (b1 << 32);
70 
71  d = a0 * b0;
72  d1 = d >> 32;
73  d0 = d - (d1 << 32);
74 
75  e = a0 * b1;
76  e1 = e >> 32;
77  e0 = e - (e1 << 32);
78 
79  f = a1 * b0;
80  f1 = f >> 32;
81  f0 = f - (f1 << 32);
82 
83  g = a1 * b1;
84  g1 = g >> 32;
85  g0 = g - (g1 << 32);
86 
87  sum = d1 + e0 + f0;
88  carry = 0;
89  /* Can't say 1<<32 cause cpp will goof it up; 1ULL<<32 might work */
90  roll = 1 << 30;
91  roll <<= 2;
92 
93  pmax = roll - 1;
94  while (pmax < sum)
95  {
96  sum -= roll;
97  carry++;
98  }
99 
100  prod.lo = d0 + (sum << 32);
101  prod.hi = carry + e1 + f1 + g0 + (g1 << 32);
102  // prod.isbig = (prod.hi || (sum >> 31));
103  prod.isbig = prod.hi || (prod.lo >> 63);
104 
105  return prod;
106 }
107 
109 QofInt128
111 {
112  guint64 sbit = x.hi & 0x1;
113  x.hi >>= 1;
114  x.lo >>= 1;
115  x.isbig = 0;
116  if (sbit)
117  {
118  x.lo |= HIBIT;
119  x.isbig = 1;
120  return x;
121  }
122  if (x.hi)
123  {
124  x.isbig = 1;
125  }
126  return x;
127 }
128 
130 QofInt128
132 {
133  guint64 sbit;
134  sbit = x.lo & HIBIT;
135  x.hi <<= 1;
136  x.lo <<= 1;
137  x.isbig = 0;
138  if (sbit)
139  {
140  x.hi |= 1;
141  x.isbig = 1;
142  return x;
143  }
144  if (x.hi)
145  {
146  x.isbig = 1;
147  }
148  return x;
149 }
150 
152 QofInt128
154 {
155  if (0 == a.isneg)
156  {
157  a.lo++;
158  if (0 == a.lo)
159  {
160  a.hi++;
161  }
162  }
163  else
164  {
165  if (0 == a.lo)
166  {
167  a.hi--;
168  }
169  a.lo--;
170  }
171 
172  a.isbig = (a.hi != 0) || (a.lo >> 63);
173  return a;
174 }
175 
179 QofInt128
180 div128 (QofInt128 n, gint64 d)
181 {
182  QofInt128 quotient;
183  int i;
184  gint64 remainder = 0;
185 
186  quotient = n;
187  if (0 > d)
188  {
189  d = -d;
190  quotient.isneg = !quotient.isneg;
191  }
192 
193  /* Use grade-school long division algorithm */
194  for (i = 0; i < 128; i++)
195  {
196  guint64 sbit = HIBIT & quotient.hi;
197  remainder <<= 1;
198  if (sbit)
199  remainder |= 1;
200  quotient = shiftleft128 (quotient);
201  if (remainder >= d)
202  {
203  remainder -= d;
204  quotient.lo |= 1;
205  }
206  }
207 
208  /* compute the carry situation */
209  quotient.isbig = (quotient.hi || (quotient.lo >> 63));
210 
211  return quotient;
212 }
213 
219 gint64
220 rem128 (QofInt128 n, gint64 d)
221 {
222  QofInt128 quotient = div128 (n, d);
223 
224  QofInt128 mu = mult128 (quotient.lo, d);
225 
226  gint64 nn = 0x7fffffffffffffffULL & n.lo;
227  gint64 rr = 0x7fffffffffffffffULL & mu.lo;
228  return nn - rr;
229 }
230 
232 gboolean
234 {
235  if (a.lo != b.lo)
236  return 0;
237  if (a.hi != b.hi)
238  return 0;
239  if (a.isneg != b.isneg)
240  return 0;
241  return 1;
242 }
243 
245 gint
247 {
248  if ((0 == a.isneg) && b.isneg)
249  return 1;
250  if (a.isneg && (0 == b.isneg))
251  return -1;
252  if (0 == a.isneg)
253  {
254  if (a.hi > b.hi)
255  return 1;
256  if (a.hi < b.hi)
257  return -1;
258  if (a.lo > b.lo)
259  return 1;
260  if (a.lo < b.lo)
261  return -1;
262  return 0;
263  }
264 
265  if (a.hi > b.hi)
266  return -1;
267  if (a.hi < b.hi)
268  return 1;
269  if (a.lo > b.lo)
270  return -1;
271  if (a.lo < b.lo)
272  return 1;
273  return 0;
274 }
275 
277 guint64
278 gcf64 (guint64 num, guint64 denom)
279 {
280  guint64 t;
281 
282  t = num % denom;
283  num = denom;
284  denom = t;
285 
286  /* The strategy is to use Euclid's algorithm */
287  while (0 != denom)
288  {
289  t = num % denom;
290  num = denom;
291  denom = t;
292  }
293  /* num now holds the GCD (Greatest Common Divisor) */
294  return num;
295 }
296 
298 QofInt128
299 lcm128 (guint64 a, guint64 b)
300 {
301  guint64 gcf = gcf64 (a, b);
302  b /= gcf;
303  return mult128 (a, b);
304 }
305 
307 QofInt128
309 {
310  QofInt128 sum;
311  if (a.isneg == b.isneg)
312  {
313  sum.isneg = a.isneg;
314  sum.hi = a.hi + b.hi;
315  sum.lo = a.lo + b.lo;
316  if ((sum.lo < a.lo) || (sum.lo < b.lo))
317  {
318  sum.hi++;
319  }
320  sum.isbig = sum.hi || (sum.lo >> 63);
321  return sum;
322  }
323  if ((b.hi > a.hi) || ((b.hi == a.hi) && (b.lo > a.lo)))
324  {
325  QofInt128 tmp = a;
326  a = b;
327  b = tmp;
328  }
329 
330  sum.isneg = a.isneg;
331  sum.hi = a.hi - b.hi;
332  sum.lo = a.lo - b.lo;
333 
334  if (sum.lo > a.lo)
335  {
336  sum.hi--;
337  }
338 
339  sum.isbig = sum.hi || (sum.lo >> 63);
340  return sum;
341 }
342 
343 
344 #ifdef TEST_128_BIT_MULT
345 static void
346 pr (gint64 a, gint64 b)
347 {
348  QofInt128 prod = mult128 (a, b);
349  printf ("%" G_GINT64_FORMAT " * %" G_GINT64_FORMAT " = %"
350  G_GUINT64_FORMAT " %" G_GUINT64_FORMAT " (0x%llx %llx) %hd\n", a,
351  b, prod.hi, prod.lo, prod.hi, prod.lo, prod.isbig);
352 }
353 
354 static void
355 prd (gint64 a, gint64 b, gint64 c)
356 {
357  QofInt128 prod = mult128 (a, b);
358  QofInt128 quot = div128 (prod, c);
359  gint64 rem = rem128 (prod, c);
360  printf ("%" G_GINT64_FORMAT " * %" G_GINT64_FORMAT " / %"
361  G_GINT64_FORMAT " = %" G_GUINT64_FORMAT " %" G_GUINT64_FORMAT
362  " + %" G_GINT64_FORMAT " (0x%llx %llx) %hd\n", a, b, c, quot.hi,
363  quot.lo, rem, quot.hi, quot.lo, quot.isbig);
364 }
365 
366 int
367 main ()
368 {
369  pr (2, 2);
370 
371  gint64 x = 1 << 30;
372  x <<= 2;
373 
374  pr (x, x);
375  pr (x + 1, x);
376  pr (x + 1, x + 1);
377 
378  pr (x, -x);
379  pr (-x, -x);
380  pr (x - 1, x);
381  pr (x - 1, x - 1);
382  pr (x - 2, x - 2);
383 
384  x <<= 1;
385  pr (x, x);
386  pr (x, -x);
387 
388  pr (1000000, 10000000000000);
389 
390  prd (x, x, 2);
391  prd (x, x, 3);
392  prd (x, x, 4);
393  prd (x, x, 5);
394  prd (x, x, 6);
395 
396  x <<= 29;
397  prd (3, x, 3);
398  prd (6, x, 3);
399  prd (99, x, 3);
400  prd (100, x, 5);
401  prd (540, x, 5);
402  prd (777, x, 7);
403  prd (1111, x, 11);
404 
405  /* Really test division */
406  QofInt128 n;
407  n.hi = 0xdd91;
408  n.lo = 0x6c5abefbb9e13480ULL;
409 
410  gint64 d = 0x2ae79964d3ae1d04ULL;
411 
412  int i;
413  for (i = 0; i < 20; i++)
414  {
415 
416  QofInt128 quot = div128 (n, d);
417  printf ("%d result = %llx %llx\n", i, quot.hi, quot.lo);
418  d >>= 1;
419  n = shift128 (n);
420  }
421  return 0;
422 }
423 
424 #endif /* TEST_128_BIT_MULT */
425 
426 /* ======================== END OF FILE =================== */
guint64 gcf64(guint64 num, guint64 denom)
Definition: qofmath128.c:278
gshort isneg
Definition: qofmath128.h:42
QofInt128 inc128(QofInt128 a)
Definition: qofmath128.c:153
gint cmp128(QofInt128 a, QofInt128 b)
Definition: qofmath128.c:246
gint64 rem128(QofInt128 n, gint64 d)
Definition: qofmath128.c:220
QofInt128 shift128(QofInt128 x)
Definition: qofmath128.c:110
QofInt128 div128(QofInt128 n, gint64 d)
Definition: qofmath128.c:180
QofInt128 add128(QofInt128 a, QofInt128 b)
Definition: qofmath128.c:308
gboolean equal128(QofInt128 a, QofInt128 b)
Definition: qofmath128.c:233
QofInt128 shiftleft128(QofInt128 x)
Definition: qofmath128.c:131
gshort isbig
Definition: qofmath128.h:43
QofInt128 mult128(gint64 a, gint64 b)
Definition: qofmath128.c:41
QofInt128 lcm128(guint64 a, guint64 b)
Definition: qofmath128.c:299