QOF 0.8.2
|
00001 /******************************************************************** 00002 * qofmath128.c -- an 128-bit integer library * 00003 * Copyright (C) 2004 Linas Vepstas <linas@linas.org> * 00004 * * 00005 * This program is free software; you can redistribute it and/or * 00006 * modify it under the terms of the GNU General Public License as * 00007 * published by the Free Software Foundation; either version 2 of * 00008 * the License, or (at your option) any later version. * 00009 * * 00010 * This program is distributed in the hope that it will be useful, * 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 00013 * GNU General Public License for more details. * 00014 * * 00015 * You should have received a copy of the GNU General Public License* 00016 * along with this program; if not, contact: * 00017 * * 00018 * Free Software Foundation Voice: +1-617-542-5942 * 00019 * 51 Franklin Street, Fifth Floor Fax: +1-617-542-2652 * 00020 * Boston, MA 02110-1301, USA gnu@gnu.org * 00021 * * 00022 *******************************************************************/ 00023 00024 #include "config.h" 00025 #include "qofmath128.h" 00026 00027 #include <glib.h> 00028 00029 /* =============================================================== */ 00030 /* 00031 * Quick-n-dirty 128-bit integer math lib. Things seem to mostly 00032 * work, and have been tested, but not comprehensively tested. 00033 */ 00034 00035 #define HIBIT (0x8000000000000000ULL) 00036 00040 QofInt128 00041 mult128 (gint64 a, gint64 b) 00042 { 00043 QofInt128 prod; 00044 guint64 a0, a1; 00045 guint64 b0, b1; 00046 guint64 d, d0, d1; 00047 guint64 e, e0, e1; 00048 guint64 f, f0, f1; 00049 guint64 g, g0, g1; 00050 guint64 sum, carry, roll, pmax; 00051 00052 prod.isneg = 0; 00053 if (0 > a) 00054 { 00055 prod.isneg = !prod.isneg; 00056 a = -a; 00057 } 00058 00059 if (0 > b) 00060 { 00061 prod.isneg = !prod.isneg; 00062 b = -b; 00063 } 00064 00065 a1 = a >> 32; 00066 a0 = a - (a1 << 32); 00067 00068 b1 = b >> 32; 00069 b0 = b - (b1 << 32); 00070 00071 d = a0 * b0; 00072 d1 = d >> 32; 00073 d0 = d - (d1 << 32); 00074 00075 e = a0 * b1; 00076 e1 = e >> 32; 00077 e0 = e - (e1 << 32); 00078 00079 f = a1 * b0; 00080 f1 = f >> 32; 00081 f0 = f - (f1 << 32); 00082 00083 g = a1 * b1; 00084 g1 = g >> 32; 00085 g0 = g - (g1 << 32); 00086 00087 sum = d1 + e0 + f0; 00088 carry = 0; 00089 /* Can't say 1<<32 cause cpp will goof it up; 1ULL<<32 might work */ 00090 roll = 1 << 30; 00091 roll <<= 2; 00092 00093 pmax = roll - 1; 00094 while (pmax < sum) 00095 { 00096 sum -= roll; 00097 carry++; 00098 } 00099 00100 prod.lo = d0 + (sum << 32); 00101 prod.hi = carry + e1 + f1 + g0 + (g1 << 32); 00102 // prod.isbig = (prod.hi || (sum >> 31)); 00103 prod.isbig = prod.hi || (prod.lo >> 63); 00104 00105 return prod; 00106 } 00107 00109 QofInt128 00110 shift128 (QofInt128 x) 00111 { 00112 guint64 sbit = x.hi & 0x1; 00113 x.hi >>= 1; 00114 x.lo >>= 1; 00115 x.isbig = 0; 00116 if (sbit) 00117 { 00118 x.lo |= HIBIT; 00119 x.isbig = 1; 00120 return x; 00121 } 00122 if (x.hi) 00123 { 00124 x.isbig = 1; 00125 } 00126 return x; 00127 } 00128 00130 QofInt128 00131 shiftleft128 (QofInt128 x) 00132 { 00133 guint64 sbit; 00134 sbit = x.lo & HIBIT; 00135 x.hi <<= 1; 00136 x.lo <<= 1; 00137 x.isbig = 0; 00138 if (sbit) 00139 { 00140 x.hi |= 1; 00141 x.isbig = 1; 00142 return x; 00143 } 00144 if (x.hi) 00145 { 00146 x.isbig = 1; 00147 } 00148 return x; 00149 } 00150 00152 QofInt128 00153 inc128 (QofInt128 a) 00154 { 00155 if (0 == a.isneg) 00156 { 00157 a.lo++; 00158 if (0 == a.lo) 00159 { 00160 a.hi++; 00161 } 00162 } 00163 else 00164 { 00165 if (0 == a.lo) 00166 { 00167 a.hi--; 00168 } 00169 a.lo--; 00170 } 00171 00172 a.isbig = (a.hi != 0) || (a.lo >> 63); 00173 return a; 00174 } 00175 00179 QofInt128 00180 div128 (QofInt128 n, gint64 d) 00181 { 00182 QofInt128 quotient; 00183 int i; 00184 gint64 remainder = 0; 00185 00186 quotient = n; 00187 if (0 > d) 00188 { 00189 d = -d; 00190 quotient.isneg = !quotient.isneg; 00191 } 00192 00193 /* Use grade-school long division algorithm */ 00194 for (i = 0; i < 128; i++) 00195 { 00196 guint64 sbit = HIBIT & quotient.hi; 00197 remainder <<= 1; 00198 if (sbit) 00199 remainder |= 1; 00200 quotient = shiftleft128 (quotient); 00201 if (remainder >= d) 00202 { 00203 remainder -= d; 00204 quotient.lo |= 1; 00205 } 00206 } 00207 00208 /* compute the carry situation */ 00209 quotient.isbig = (quotient.hi || (quotient.lo >> 63)); 00210 00211 return quotient; 00212 } 00213 00219 gint64 00220 rem128 (QofInt128 n, gint64 d) 00221 { 00222 QofInt128 quotient = div128 (n, d); 00223 00224 QofInt128 mu = mult128 (quotient.lo, d); 00225 00226 gint64 nn = 0x7fffffffffffffffULL & n.lo; 00227 gint64 rr = 0x7fffffffffffffffULL & mu.lo; 00228 return nn - rr; 00229 } 00230 00232 gboolean 00233 equal128 (QofInt128 a, QofInt128 b) 00234 { 00235 if (a.lo != b.lo) 00236 return 0; 00237 if (a.hi != b.hi) 00238 return 0; 00239 if (a.isneg != b.isneg) 00240 return 0; 00241 return 1; 00242 } 00243 00245 gint 00246 cmp128 (QofInt128 a, QofInt128 b) 00247 { 00248 if ((0 == a.isneg) && b.isneg) 00249 return 1; 00250 if (a.isneg && (0 == b.isneg)) 00251 return -1; 00252 if (0 == a.isneg) 00253 { 00254 if (a.hi > b.hi) 00255 return 1; 00256 if (a.hi < b.hi) 00257 return -1; 00258 if (a.lo > b.lo) 00259 return 1; 00260 if (a.lo < b.lo) 00261 return -1; 00262 return 0; 00263 } 00264 00265 if (a.hi > b.hi) 00266 return -1; 00267 if (a.hi < b.hi) 00268 return 1; 00269 if (a.lo > b.lo) 00270 return -1; 00271 if (a.lo < b.lo) 00272 return 1; 00273 return 0; 00274 } 00275 00277 guint64 00278 gcf64 (guint64 num, guint64 denom) 00279 { 00280 guint64 t; 00281 00282 t = num % denom; 00283 num = denom; 00284 denom = t; 00285 00286 /* The strategy is to use Euclid's algorithm */ 00287 while (0 != denom) 00288 { 00289 t = num % denom; 00290 num = denom; 00291 denom = t; 00292 } 00293 /* num now holds the GCD (Greatest Common Divisor) */ 00294 return num; 00295 } 00296 00298 QofInt128 00299 lcm128 (guint64 a, guint64 b) 00300 { 00301 guint64 gcf = gcf64 (a, b); 00302 b /= gcf; 00303 return mult128 (a, b); 00304 } 00305 00307 QofInt128 00308 add128 (QofInt128 a, QofInt128 b) 00309 { 00310 QofInt128 sum; 00311 if (a.isneg == b.isneg) 00312 { 00313 sum.isneg = a.isneg; 00314 sum.hi = a.hi + b.hi; 00315 sum.lo = a.lo + b.lo; 00316 if ((sum.lo < a.lo) || (sum.lo < b.lo)) 00317 { 00318 sum.hi++; 00319 } 00320 sum.isbig = sum.hi || (sum.lo >> 63); 00321 return sum; 00322 } 00323 if ((b.hi > a.hi) || ((b.hi == a.hi) && (b.lo > a.lo))) 00324 { 00325 QofInt128 tmp = a; 00326 a = b; 00327 b = tmp; 00328 } 00329 00330 sum.isneg = a.isneg; 00331 sum.hi = a.hi - b.hi; 00332 sum.lo = a.lo - b.lo; 00333 00334 if (sum.lo > a.lo) 00335 { 00336 sum.hi--; 00337 } 00338 00339 sum.isbig = sum.hi || (sum.lo >> 63); 00340 return sum; 00341 } 00342 00343 00344 #ifdef TEST_128_BIT_MULT 00345 static void 00346 pr (gint64 a, gint64 b) 00347 { 00348 QofInt128 prod = mult128 (a, b); 00349 printf ("%" G_GINT64_FORMAT " * %" G_GINT64_FORMAT " = %" 00350 G_GUINT64_FORMAT " %" G_GUINT64_FORMAT " (0x%llx %llx) %hd\n", a, 00351 b, prod.hi, prod.lo, prod.hi, prod.lo, prod.isbig); 00352 } 00353 00354 static void 00355 prd (gint64 a, gint64 b, gint64 c) 00356 { 00357 QofInt128 prod = mult128 (a, b); 00358 QofInt128 quot = div128 (prod, c); 00359 gint64 rem = rem128 (prod, c); 00360 printf ("%" G_GINT64_FORMAT " * %" G_GINT64_FORMAT " / %" 00361 G_GINT64_FORMAT " = %" G_GUINT64_FORMAT " %" G_GUINT64_FORMAT 00362 " + %" G_GINT64_FORMAT " (0x%llx %llx) %hd\n", a, b, c, quot.hi, 00363 quot.lo, rem, quot.hi, quot.lo, quot.isbig); 00364 } 00365 00366 int 00367 main () 00368 { 00369 pr (2, 2); 00370 00371 gint64 x = 1 << 30; 00372 x <<= 2; 00373 00374 pr (x, x); 00375 pr (x + 1, x); 00376 pr (x + 1, x + 1); 00377 00378 pr (x, -x); 00379 pr (-x, -x); 00380 pr (x - 1, x); 00381 pr (x - 1, x - 1); 00382 pr (x - 2, x - 2); 00383 00384 x <<= 1; 00385 pr (x, x); 00386 pr (x, -x); 00387 00388 pr (1000000, 10000000000000); 00389 00390 prd (x, x, 2); 00391 prd (x, x, 3); 00392 prd (x, x, 4); 00393 prd (x, x, 5); 00394 prd (x, x, 6); 00395 00396 x <<= 29; 00397 prd (3, x, 3); 00398 prd (6, x, 3); 00399 prd (99, x, 3); 00400 prd (100, x, 5); 00401 prd (540, x, 5); 00402 prd (777, x, 7); 00403 prd (1111, x, 11); 00404 00405 /* Really test division */ 00406 QofInt128 n; 00407 n.hi = 0xdd91; 00408 n.lo = 0x6c5abefbb9e13480ULL; 00409 00410 gint64 d = 0x2ae79964d3ae1d04ULL; 00411 00412 int i; 00413 for (i = 0; i < 20; i++) 00414 { 00415 00416 QofInt128 quot = div128 (n, d); 00417 printf ("%d result = %llx %llx\n", i, quot.hi, quot.lo); 00418 d >>= 1; 00419 n = shift128 (n); 00420 } 00421 return 0; 00422 } 00423 00424 #endif /* TEST_128_BIT_MULT */ 00425 00426 /* ======================== END OF FILE =================== */