qofmath128.c

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 #define _GNU_SOURCE
00025 
00026 #include "config.h"
00027 #include "qofmath128.h"
00028 
00029 #include <glib.h>
00030 
00031 /* =============================================================== */
00032 /*
00033  *  Quick-n-dirty 128-bit integer math lib.   Things seem to mostly
00034  *  work, and have been tested, but not comprehensively tested.
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   /* Can't say 1<<32 cause cpp will goof it up; 1ULL<<32 might work */
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   // prod.isbig = (prod.hi || (sum >> 31));
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   /* Use grade-school long division algorithm */
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   /* compute the carry situation */
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   /* The strategy is to use Euclid's algorithm */
00275   while (0 != denom) 
00276   {
00277     t = num % denom;
00278     num = denom;
00279     denom = t;
00280   }
00281   /* num now holds the GCD (Greatest Common Divisor) */
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   /* Really test division */
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 /* TEST_128_BIT_MULT */
00407 
00408 /* ======================== END OF FILE =================== */

Generated on Fri May 12 17:57:20 2006 for QOF by  doxygen 1.4.4