00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #define _GNU_SOURCE
00025
00026 #include "config.h"
00027 #include "qofmath128.h"
00028
00029 #include <glib.h>
00030
00031
00032
00033
00034
00035
00036
00037 #define HIBIT (0x8000000000000000ULL)
00038
00042 inline qofint128
00043 mult128 (gint64 a, gint64 b)
00044 {
00045 qofint128 prod;
00046 guint64 a0, a1;
00047 guint64 b0, b1;
00048 guint64 d, d0, d1;
00049 guint64 e, e0, e1;
00050 guint64 f, f0, f1;
00051 guint64 g, g0, g1;
00052 guint64 sum, carry, roll, pmax;
00053
00054 prod.isneg = 0;
00055 if (0>a)
00056 {
00057 prod.isneg = !prod.isneg;
00058 a = -a;
00059 }
00060
00061 if (0>b)
00062 {
00063 prod.isneg = !prod.isneg;
00064 b = -b;
00065 }
00066
00067 a1 = a >> 32;
00068 a0 = a - (a1<<32);
00069
00070 b1 = b >> 32;
00071 b0 = b - (b1<<32);
00072
00073 d = a0*b0;
00074 d1 = d >> 32;
00075 d0 = d - (d1<<32);
00076
00077 e = a0*b1;
00078 e1 = e >> 32;
00079 e0 = e - (e1<<32);
00080
00081 f = a1*b0;
00082 f1 = f >> 32;
00083 f0 = f - (f1<<32);
00084
00085 g = a1*b1;
00086 g1 = g >> 32;
00087 g0 = g - (g1<<32);
00088
00089 sum = d1+e0+f0;
00090 carry = 0;
00091
00092 roll = 1<<30;
00093 roll <<= 2;
00094
00095 pmax = roll-1;
00096 while (pmax < sum)
00097 {
00098 sum -= roll;
00099 carry ++;
00100 }
00101
00102 prod.lo = d0 + (sum<<32);
00103 prod.hi = carry + e1 + f1 + g0 + (g1<<32);
00104
00105 prod.isbig = prod.hi || (prod.lo >> 63);
00106
00107 return prod;
00108 }
00109
00111 inline qofint128
00112 shift128 (qofint128 x)
00113 {
00114 guint64 sbit = x.hi & 0x1;
00115 x.hi >>= 1;
00116 x.lo >>= 1;
00117 x.isbig = 0;
00118 if (sbit)
00119 {
00120 x.lo |= HIBIT;
00121 x.isbig = 1;
00122 return x;
00123 }
00124 if (x.hi)
00125 {
00126 x.isbig = 1;
00127 }
00128 return x;
00129 }
00130
00132 inline qofint128
00133 shiftleft128 (qofint128 x)
00134 {
00135 guint64 sbit;
00136 sbit = x.lo & HIBIT;
00137 x.hi <<= 1;
00138 x.lo <<= 1;
00139 x.isbig = 0;
00140 if (sbit)
00141 {
00142 x.hi |= 1;
00143 x.isbig = 1;
00144 return x;
00145 }
00146 if (x.hi)
00147 {
00148 x.isbig = 1;
00149 }
00150 return x;
00151 }
00152
00154 inline qofint128
00155 inc128 (qofint128 a)
00156 {
00157 if (0 == a.isneg)
00158 {
00159 a.lo ++;
00160 if (0 == a.lo)
00161 {
00162 a.hi ++;
00163 }
00164 }
00165 else
00166 {
00167 if (0 == a.lo)
00168 {
00169 a.hi --;
00170 }
00171 a.lo --;
00172 }
00173
00174 a.isbig = (a.hi != 0) || (a.lo >> 63);
00175 return a;
00176 }
00177
00181 inline qofint128
00182 div128 (qofint128 n, gint64 d)
00183 {
00184 qofint128 quotient;
00185 int i;
00186 guint64 remainder = 0;
00187
00188 quotient = n;
00189 if (0 > d)
00190 {
00191 d = -d;
00192 quotient.isneg = !quotient.isneg;
00193 }
00194
00195
00196 for (i=0; i<128; i++)
00197 {
00198 guint64 sbit = HIBIT & quotient.hi;
00199 remainder <<= 1;
00200 if (sbit) remainder |= 1;
00201 quotient = shiftleft128 (quotient);
00202 if (remainder >= d)
00203 {
00204 remainder -= d;
00205 quotient.lo |= 1;
00206 }
00207 }
00208
00209
00210 quotient.isbig = (quotient.hi || (quotient.lo >> 63));
00211
00212 return quotient;
00213 }
00214
00220 inline gint64
00221 rem128 (qofint128 n, gint64 d)
00222 {
00223 qofint128 quotient = div128 (n,d);
00224
00225 qofint128 mu = mult128 (quotient.lo, d);
00226
00227 gint64 nn = 0x7fffffffffffffffULL & n.lo;
00228 gint64 rr = 0x7fffffffffffffffULL & mu.lo;
00229 return nn - rr;
00230 }
00231
00233 inline gboolean
00234 equal128 (qofint128 a, qofint128 b)
00235 {
00236 if (a.lo != b.lo) return 0;
00237 if (a.hi != b.hi) return 0;
00238 if (a.isneg != b.isneg) return 0;
00239 return 1;
00240 }
00241
00243 inline int
00244 cmp128 (qofint128 a, qofint128 b)
00245 {
00246 if ((0 == a.isneg) && b.isneg) return 1;
00247 if (a.isneg && (0 == b.isneg)) return -1;
00248 if (0 == a.isneg)
00249 {
00250 if (a.hi > b.hi) return 1;
00251 if (a.hi < b.hi) return -1;
00252 if (a.lo > b.lo) return 1;
00253 if (a.lo < b.lo) return -1;
00254 return 0;
00255 }
00256
00257 if (a.hi > b.hi) return -1;
00258 if (a.hi < b.hi) return 1;
00259 if (a.lo > b.lo) return -1;
00260 if (a.lo < b.lo) return 1;
00261 return 0;
00262 }
00263
00265 inline guint64
00266 gcf64(guint64 num, guint64 denom)
00267 {
00268 guint64 t;
00269
00270 t = num % denom;
00271 num = denom;
00272 denom = t;
00273
00274
00275 while (0 != denom)
00276 {
00277 t = num % denom;
00278 num = denom;
00279 denom = t;
00280 }
00281
00282 return num;
00283 }
00284
00286 inline qofint128
00287 lcm128 (guint64 a, guint64 b)
00288 {
00289 guint64 gcf = gcf64 (a,b);
00290 b /= gcf;
00291 return mult128 (a,b);
00292 }
00293
00295 inline qofint128
00296 add128 (qofint128 a, qofint128 b)
00297 {
00298 qofint128 sum;
00299 if (a.isneg == b.isneg)
00300 {
00301 sum.isneg = a.isneg;
00302 sum.hi = a.hi + b.hi;
00303 sum.lo = a.lo + b.lo;
00304 if ((sum.lo < a.lo) || (sum.lo < b.lo))
00305 {
00306 sum.hi ++;
00307 }
00308 sum.isbig = sum.hi || (sum.lo >> 63);
00309 return sum;
00310 }
00311 if ((b.hi > a.hi) ||
00312 ((b.hi == a.hi) && (b.lo > a.lo)))
00313 {
00314 qofint128 tmp = a;
00315 a = b;
00316 b = tmp;
00317 }
00318
00319 sum.isneg = a.isneg;
00320 sum.hi = a.hi - b.hi;
00321 sum.lo = a.lo - b.lo;
00322
00323 if (sum.lo > a.lo)
00324 {
00325 sum.hi --;
00326 }
00327
00328 sum.isbig = sum.hi || (sum.lo >> 63);
00329 return sum;
00330 }
00331
00332
00333 #ifdef TEST_128_BIT_MULT
00334 static void pr (gint64 a, gint64 b)
00335 {
00336 qofint128 prod = mult128 (a,b);
00337 printf ("%" G_GINT64_FORMAT " * %" G_GINT64_FORMAT " = %" G_GUINT64_FORMAT " %" G_GUINT64_FORMAT " (0x%llx %llx) %hd\n",
00338 a, b, prod.hi, prod.lo, prod.hi, prod.lo, prod.isbig);
00339 }
00340
00341 static void prd (gint64 a, gint64 b, gint64 c)
00342 {
00343 qofint128 prod = mult128 (a,b);
00344 qofint128 quot = div128 (prod, c);
00345 gint64 rem = rem128 (prod, c);
00346 printf ("%" G_GINT64_FORMAT " * %" G_GINT64_FORMAT " / %" G_GINT64_FORMAT " = %" G_GUINT64_FORMAT " %" G_GUINT64_FORMAT " + %" G_GINT64_FORMAT " (0x%llx %llx) %hd\n",
00347 a, b, c, quot.hi, quot.lo, rem, quot.hi, quot.lo, quot.isbig);
00348 }
00349
00350 int main ()
00351 {
00352 pr (2,2);
00353
00354 gint64 x = 1<<30;
00355 x <<= 2;
00356
00357 pr (x,x);
00358 pr (x+1,x);
00359 pr (x+1,x+1);
00360
00361 pr (x,-x);
00362 pr (-x,-x);
00363 pr (x-1,x);
00364 pr (x-1,x-1);
00365 pr (x-2,x-2);
00366
00367 x <<= 1;
00368 pr (x,x);
00369 pr (x,-x);
00370
00371 pr (1000000, 10000000000000);
00372
00373 prd (x,x,2);
00374 prd (x,x,3);
00375 prd (x,x,4);
00376 prd (x,x,5);
00377 prd (x,x,6);
00378
00379 x <<= 29;
00380 prd (3,x,3);
00381 prd (6,x,3);
00382 prd (99,x,3);
00383 prd (100,x,5);
00384 prd (540,x,5);
00385 prd (777,x,7);
00386 prd (1111,x,11);
00387
00388
00389 qofint128 n;
00390 n.hi = 0xdd91;
00391 n.lo = 0x6c5abefbb9e13480ULL;
00392
00393 gint64 d = 0x2ae79964d3ae1d04ULL;
00394
00395 int i;
00396 for (i=0; i<20; i++) {
00397
00398 qofint128 quot = div128 (n, d);
00399 printf ("%d result = %llx %llx\n", i, quot.hi, quot.lo);
00400 d >>=1;
00401 n = shift128 (n);
00402 }
00403 return 0;
00404 }
00405
00406 #endif
00407
00408