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
37 /** Multiply a pair of signed 64-bit numbers,
38 * returning a signed 128-bit number.
39 */
40 QofInt128
mult128(gint64 a,gint64 b)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
108 /** Shift right by one bit (i.e. divide by two) */
109 QofInt128
shift128(QofInt128 x)110 shift128 (QofInt128 x)
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
129 /** Shift leftt by one bit (i.e. multiply by two) */
130 QofInt128
shiftleft128(QofInt128 x)131 shiftleft128 (QofInt128 x)
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
151 /** increment a 128-bit number by one */
152 QofInt128
inc128(QofInt128 a)153 inc128 (QofInt128 a)
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
176 /** Divide a signed 128-bit number by a signed 64-bit,
177 * returning a signed 128-bit number.
178 */
179 QofInt128
div128(QofInt128 n,gint64 d)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
214 /** Return the remainder of a signed 128-bit number modulo
215 * a signed 64-bit. That is, return n%d in 128-bit math.
216 * I beleive that ths algo is overflow-free, but should be
217 * audited some more ...
218 */
219 gint64
rem128(QofInt128 n,gint64 d)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
231 /** Return true of two numbers are equal */
232 gboolean
equal128(QofInt128 a,QofInt128 b)233 equal128 (QofInt128 a, QofInt128 b)
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
244 /** Return returns 1 if a>b, -1 if b>a, 0 if a == b */
245 gint
cmp128(QofInt128 a,QofInt128 b)246 cmp128 (QofInt128 a, QofInt128 b)
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
276 /** Return the greatest common factor of two 64-bit numbers */
277 guint64
gcf64(guint64 num,guint64 denom)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
297 /** Return the least common multiple of two 64-bit numbers. */
298 QofInt128
lcm128(guint64 a,guint64 b)299 lcm128 (guint64 a, guint64 b)
300 {
301 guint64 gcf = gcf64 (a, b);
302 b /= gcf;
303 return mult128 (a, b);
304 }
305
306 /** Add a pair of 128-bit numbers, returning a 128-bit number */
307 QofInt128
add128(QofInt128 a,QofInt128 b)308 add128 (QofInt128 a, QofInt128 b)
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
pr(gint64 a,gint64 b)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
prd(gint64 a,gint64 b,gint64 c)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
main()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 =================== */
427