Make mantissa type of mp_no configurable

The mantissa of mp_no is intended to take only integral values.  This
is a relatively good choice for powerpc due to its 4 fpus, but not for
other architectures, which suffer due to this choice.  This change
makes the default mantissa a long integer and allows powerpc to
override it.  Additionally, some operations have been optimized for
integer manipulation, resulting in a significant improvement in
performance.

1 files changed, 38 insertions, 58 deletions

diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
index 0766476544..c238ccf2af 100644
--- a/sysdeps/ieee754/dbl-64/mpa.c
+++ b/sysdeps/ieee754/dbl-64/mpa.c
@@ -125,7 +125,8 @@ norm (const mp_no *x, double *y, int p)
 {
 #define R RADIXI
   long i;
-  double a, c, u, v, z[5];
+  double c;
+  mantissa_t a, u, v, z[5];
   if (p < 5)
     {
       if (p == 1)
@@ -147,17 +148,14 @@ norm (const mp_no *x, double *y, int p)
 
       for (i = 2; i < 5; i++)
 	{
-	  z[i] = X[i] * a;
-	  u = (z[i] + CUTTER) - CUTTER;
-	  if (u > z[i])
-	    u -= RADIX;
-	  z[i] -= u;
-	  z[i - 1] += u * RADIXI;
+	  mantissa_store_t d, r;
+	  d = X[i] * (mantissa_store_t) a;
+	  DIV_RADIX (d, r);
+	  z[i] = r;
+	  z[i - 1] += d;
 	}
 
-      u = (z[3] + TWO71) - TWO71;
-      if (u > z[3])
-	u -= TWO19;
+      u = ALIGN_DOWN_TO (z[3], TWO19);
       v = z[3] - u;
 
       if (v == TWO18)
@@ -200,7 +198,8 @@ denorm (const mp_no *x, double *y, int p)
 {
   long i, k;
   long p2 = p;
-  double c, u, z[5];
+  double c;
+  mantissa_t u, z[5];
 
 #define R RADIXI
   if (EX < -44 || (EX == -44 && X[1] < TWO5))
@@ -280,9 +279,7 @@ denorm (const mp_no *x, double *y, int p)
       z[3] = X[k];
     }
 
-  u = (z[3] + TWO57) - TWO57;
-  if (u > z[3])
-    u -= TWO5;
+  u = ALIGN_DOWN_TO (z[3], TWO5);
 
   if (u == z[3])
     {
@@ -330,7 +327,6 @@ __dbl_mp (double x, mp_no *y, int p)
 {
   long i, n;
   long p2 = p;
-  double u;
 
   /* Sign.  */
   if (x == ZERO)
@@ -356,11 +352,7 @@ __dbl_mp (double x, mp_no *y, int p)
   n = MIN (p2, 4);
   for (i = 1; i <= n; i++)
     {
-      u = (x + TWO52) - TWO52;
-      if (u > x)
-	u -= ONE;
-      Y[i] = u;
-      x -= u;
+      INTEGER_OF (x, Y[i]);
       x *= RADIX;
     }
   for (; i <= p2; i++)
@@ -377,7 +369,7 @@ add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
   long i, j, k;
   long p2 = p;
-  double zk;
+  mantissa_t zk;
 
   EZ = EX;
 
@@ -445,7 +437,7 @@ sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
   long i, j, k;
   long p2 = p;
-  double zk;
+  mantissa_t zk;
 
   EZ = EX;
   i = p2;
@@ -621,9 +613,9 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
   long i, j, k, ip, ip2;
   long p2 = p;
-  double u, zk;
+  mantissa_store_t zk;
   const mp_no *a;
-  double *diag;
+  mantissa_store_t *diag;
 
   /* Is z=0?  */
   if (__glibc_unlikely (X[0] * Y[0] == ZERO))
@@ -680,11 +672,11 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
 
   /* Precompute sums of diagonal elements so that we can directly use them
      later.  See the next comment to know we why need them.  */
-  diag = alloca (k * sizeof (double));
-  double d = ZERO;
+  diag = alloca (k * sizeof (mantissa_store_t));
+  mantissa_store_t d = ZERO;
   for (i = 1; i <= ip; i++)
     {
-      d += X[i] * Y[i];
+      d += X[i] * (mantissa_store_t) Y[i];
       diag[i] = d;
     }
   while (i < k)
@@ -697,18 +689,15 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
       if (k % 2 == 0)
 	/* We want to add this only once, but since we subtract it in the sum
 	   of products above, we add twice.  */
-	zk += 2 * X[lim] * Y[lim];
+	zk += 2 * X[lim] * (mantissa_store_t) Y[lim];
 
       for (i = k - p2, j = p2; i < j; i++, j--)
-	zk += (X[i] + X[j]) * (Y[i] + Y[j]);
+	zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]);
 
       zk -= diag[k - 1];
 
-      u = (zk + CUTTER) - CUTTER;
-      if (u > zk)
-	u -= RADIX;
-      Z[k--] = zk - u;
-      zk = u * RADIXI;
+      DIV_RADIX (zk, Z[k]);
+      k--;
     }
 
   /* The real deal.  Mantissa digit Z[k] is the sum of all X[i] * Y[j] where i
@@ -731,18 +720,15 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
       if (k % 2 == 0)
 	/* We want to add this only once, but since we subtract it in the sum
 	   of products above, we add twice.  */
-        zk += 2 * X[lim] * Y[lim];
+        zk += 2 * X[lim] * (mantissa_store_t) Y[lim];
 
       for (i = 1, j = k - 1; i < j; i++, j--)
-	zk += (X[i] + X[j]) * (Y[i] + Y[j]);
+	zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]);
 
       zk -= diag[k - 1];
 
-      u = (zk + CUTTER) - CUTTER;
-      if (u > zk)
-	u -= RADIX;
-      Z[k--] = zk - u;
-      zk = u * RADIXI;
+      DIV_RADIX (zk, Z[k]);
+      k--;
     }
   Z[k] = zk;
 
@@ -774,7 +760,7 @@ SECTION
 __sqr (const mp_no *x, mp_no *y, int p)
 {
   long i, j, k, ip;
-  double u, yk;
+  mantissa_store_t yk;
 
   /* Is z=0?  */
   if (__glibc_unlikely (X[0] == ZERO))
@@ -798,11 +784,11 @@ __sqr (const mp_no *x, mp_no *y, int p)
 
   while (k > p)
     {
-      double yk2 = 0.0;
+      mantissa_store_t yk2 = 0;
       long lim = k / 2;
 
       if (k % 2 == 0)
-	yk += X[lim] * X[lim];
+	yk += X[lim] * (mantissa_store_t) X[lim];
 
       /* In __mul, this loop (and the one within the next while loop) run
          between a range to calculate the mantissa as follows:
@@ -814,36 +800,30 @@ __sqr (const mp_no *x, mp_no *y, int p)
 	 result.  For cases where the range size is even, the mid-point needs
 	 to be added separately (above).  */
       for (i = k - p, j = p; i < j; i++, j--)
-	yk2 += X[i] * X[j];
+	yk2 += X[i] * (mantissa_store_t) X[j];
 
       yk += 2.0 * yk2;
 
-      u = (yk + CUTTER) - CUTTER;
-      if (u > yk)
-	u -= RADIX;
-      Y[k--] = yk - u;
-      yk = u * RADIXI;
+      DIV_RADIX (yk, Y[k]);
+      k--;
     }
 
   while (k > 1)
     {
-      double yk2 = 0.0;
+      mantissa_store_t yk2 = 0;
       long lim = k / 2;
 
       if (k % 2 == 0)
-	yk += X[lim] * X[lim];
+	yk += X[lim] * (mantissa_store_t) X[lim];
 
       /* Likewise for this loop.  */
       for (i = 1, j = k - 1; i < j; i++, j--)
-	yk2 += X[i] * X[j];
+	yk2 += X[i] * (mantissa_store_t) X[j];
 
       yk += 2.0 * yk2;
 
-      u = (yk + CUTTER) - CUTTER;
-      if (u > yk)
-	u -= RADIX;
-      Y[k--] = yk - u;
-      yk = u * RADIXI;
+      DIV_RADIX (yk, Y[k]);
+      k--;
     }
   Y[k] = yk;