Revert configurable mantissa patch

Reverts d22ca8cdfb98001d03772ef264b244930d439b3f

since it is severely broken on 32-bit.

3 files changed, 60 insertions, 92 deletions

diff --git a/sysdeps/ieee754/dbl-64/mpa-arch.h b/sysdeps/ieee754/dbl-64/mpa-arch.h
deleted file mode 100644
index 7de9d51ae2..0000000000
--- a/sysdeps/ieee754/dbl-64/mpa-arch.h
+++ /dev/null
@@ -1,47 +0,0 @@
-/* Overridable constants and operations.
-   Copyright (C) 2013 Free Software Foundation, Inc.
-
-   This program is free software; you can redistribute it and/or modify
-   it under the terms of the GNU Lesser General Public License as published by
-   the Free Software Foundation; either version 2.1 of the License, or
-   (at your option) any later version.
-
-   This program is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-   GNU Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public License
-   along with this program; if not, see <http://www.gnu.org/licenses/>.  */
-
-#include <stdint.h>
-
-typedef long mantissa_t;
-typedef int64_t mantissa_store_t;
-
-#define TWOPOW(i) (1L << i)
-
-#define RADIX_EXP 24
-#define  RADIX TWOPOW (RADIX_EXP)		/* 2^24    */
-
-/* Divide D by RADIX and put the remainder in R.  D must be a non-negative
-   integral value.  */
-#define DIV_RADIX(d, r) \
-  ({									      \
-    r = d & (RADIX - 1);						      \
-    d >>= RADIX_EXP;							      \
-  })
-
-/* Put the integer component of a double X in R and retain the fraction in
-   X.  This is used in extracting mantissa digits for MP_NO by using the
-   integer portion of the current value of the number as the current mantissa
-   digit and then scaling by RADIX to get the next mantissa digit in the same
-   manner.  */
-#define INTEGER_OF(x, i) \
-  ({									      \
-    i = (mantissa_t) x;							      \
-    x -= i;								      \
-  })
-
-/* Align IN down to F.  The code assumes that F is a power of two.  */
-#define ALIGN_DOWN_TO(in, f) ((in) & -(f))
diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
index 860e859ae8..0766476544 100644
--- a/sysdeps/ieee754/dbl-64/mpa.c
+++ b/sysdeps/ieee754/dbl-64/mpa.c
@@ -125,8 +125,7 @@ norm (const mp_no *x, double *y, int p)
 {
 #define R RADIXI
   long i;
-  double c;
-  mantissa_t a, u, v, z[5];
+  double a, c, u, v, z[5];
   if (p < 5)
     {
       if (p == 1)
@@ -148,14 +147,17 @@ norm (const mp_no *x, double *y, int p)
 
       for (i = 2; i < 5; i++)
 	{
-	  mantissa_t d, r;
-	  d = X[i] * a;
-	  DIV_RADIX (d, r);
-	  z[i] = r;
-	  z[i - 1] += d;
+	  z[i] = X[i] * a;
+	  u = (z[i] + CUTTER) - CUTTER;
+	  if (u > z[i])
+	    u -= RADIX;
+	  z[i] -= u;
+	  z[i - 1] += u * RADIXI;
 	}
 
-      u = ALIGN_DOWN_TO (z[3], TWO19);
+      u = (z[3] + TWO71) - TWO71;
+      if (u > z[3])
+	u -= TWO19;
       v = z[3] - u;
 
       if (v == TWO18)
@@ -198,8 +200,7 @@ denorm (const mp_no *x, double *y, int p)
 {
   long i, k;
   long p2 = p;
-  double c;
-  mantissa_t u, z[5];
+  double c, u, z[5];
 
 #define R RADIXI
   if (EX < -44 || (EX == -44 && X[1] < TWO5))
@@ -279,7 +280,9 @@ denorm (const mp_no *x, double *y, int p)
       z[3] = X[k];
     }
 
-  u = ALIGN_DOWN_TO (z[3], TWO5);
+  u = (z[3] + TWO57) - TWO57;
+  if (u > z[3])
+    u -= TWO5;
 
   if (u == z[3])
     {
@@ -327,6 +330,7 @@ __dbl_mp (double x, mp_no *y, int p)
 {
   long i, n;
   long p2 = p;
+  double u;
 
   /* Sign.  */
   if (x == ZERO)
@@ -352,7 +356,11 @@ __dbl_mp (double x, mp_no *y, int p)
   n = MIN (p2, 4);
   for (i = 1; i <= n; i++)
     {
-      INTEGER_OF (x, Y[i]);
+      u = (x + TWO52) - TWO52;
+      if (u > x)
+	u -= ONE;
+      Y[i] = u;
+      x -= u;
       x *= RADIX;
     }
   for (; i <= p2; i++)
@@ -369,7 +377,7 @@ add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
   long i, j, k;
   long p2 = p;
-  mantissa_t zk;
+  double zk;
 
   EZ = EX;
 
@@ -437,7 +445,7 @@ sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
   long i, j, k;
   long p2 = p;
-  mantissa_t zk;
+  double zk;
 
   EZ = EX;
   i = p2;
@@ -613,9 +621,9 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
   long i, j, k, ip, ip2;
   long p2 = p;
-  mantissa_store_t zk;
+  double u, zk;
   const mp_no *a;
-  mantissa_store_t *diag;
+  double *diag;
 
   /* Is z=0?  */
   if (__glibc_unlikely (X[0] * Y[0] == ZERO))
@@ -672,8 +680,8 @@ __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 (mantissa_store_t));
-  mantissa_store_t d = ZERO;
+  diag = alloca (k * sizeof (double));
+  double d = ZERO;
   for (i = 1; i <= ip; i++)
     {
       d += X[i] * Y[i];
@@ -696,8 +704,11 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
 
       zk -= diag[k - 1];
 
-      DIV_RADIX (zk, Z[k]);
-      k--;
+      u = (zk + CUTTER) - CUTTER;
+      if (u > zk)
+	u -= RADIX;
+      Z[k--] = zk - u;
+      zk = u * RADIXI;
     }
 
   /* The real deal.  Mantissa digit Z[k] is the sum of all X[i] * Y[j] where i
@@ -727,8 +738,11 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
 
       zk -= diag[k - 1];
 
-      DIV_RADIX (zk, Z[k]);
-      k--;
+      u = (zk + CUTTER) - CUTTER;
+      if (u > zk)
+	u -= RADIX;
+      Z[k--] = zk - u;
+      zk = u * RADIXI;
     }
   Z[k] = zk;
 
@@ -760,7 +774,7 @@ SECTION
 __sqr (const mp_no *x, mp_no *y, int p)
 {
   long i, j, k, ip;
-  mantissa_store_t yk;
+  double u, yk;
 
   /* Is z=0?  */
   if (__glibc_unlikely (X[0] == ZERO))
@@ -784,7 +798,7 @@ __sqr (const mp_no *x, mp_no *y, int p)
 
   while (k > p)
     {
-      mantissa_store_t yk2 = 0;
+      double yk2 = 0.0;
       long lim = k / 2;
 
       if (k % 2 == 0)
@@ -804,13 +818,16 @@ __sqr (const mp_no *x, mp_no *y, int p)
 
       yk += 2.0 * yk2;
 
-      DIV_RADIX (yk, Y[k]);
-      k--;
+      u = (yk + CUTTER) - CUTTER;
+      if (u > yk)
+	u -= RADIX;
+      Y[k--] = yk - u;
+      yk = u * RADIXI;
     }
 
   while (k > 1)
     {
-      mantissa_store_t yk2 = 0;
+      double yk2 = 0.0;
       long lim = k / 2;
 
       if (k % 2 == 0)
@@ -822,8 +839,11 @@ __sqr (const mp_no *x, mp_no *y, int p)
 
       yk += 2.0 * yk2;
 
-      DIV_RADIX (yk, Y[k]);
-      k--;
+      u = (yk + CUTTER) - CUTTER;
+      if (u > yk)
+	u -= RADIX;
+      Y[k--] = yk - u;
+      yk = u * RADIXI;
     }
   Y[k] = yk;
 
diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h
index 54044a0586..168b334ed0 100644
--- a/sysdeps/ieee754/dbl-64/mpa.h
+++ b/sysdeps/ieee754/dbl-64/mpa.h
@@ -35,7 +35,6 @@
 /* Common types and definition                                          */
 /************************************************************************/
 
-#include <mpa-arch.h>
 
 /* The mp_no structure holds the details of a multi-precision floating point
    number.
@@ -62,7 +61,7 @@
 typedef struct
 {
   int e;
-  mantissa_t d[40];
+  double d[40];
 } mp_no;
 
 typedef union
@@ -83,13 +82,9 @@ extern const mp_no mptwo;
 
 #define ABS(x)   ((x) <  0  ? -(x) : (x))
 
-#ifndef RADIXI
-# define  RADIXI    0x1.0p-24		/* 2^-24   */
-#endif
-
-#ifndef TWO52
-# define  TWO52     0x1.0p52		/* 2^52    */
-#endif
+#define  RADIX     0x1.0p24		/* 2^24    */
+#define  RADIXI    0x1.0p-24		/* 2^-24   */
+#define  CUTTER    0x1.0p76		/* 2^76    */
 
 #define  ZERO      0.0			/* 0       */
 #define  MZERO     -0.0			/* 0 with the sign bit set */
@@ -97,13 +92,13 @@ extern const mp_no mptwo;
 #define  MONE      -1.0			/* -1      */
 #define  TWO       2.0			/*  2      */
 
-#define  TWO5      TWOPOW (5)		/* 2^5     */
-#define  TWO8      TWOPOW (8)		/* 2^52    */
-#define  TWO10     TWOPOW (10)		/* 2^10    */
-#define  TWO18     TWOPOW (18)		/* 2^18    */
-#define  TWO19     TWOPOW (19)		/* 2^19    */
-#define  TWO23     TWOPOW (23)		/* 2^23    */
-
+#define  TWO5      0x1.0p5		/* 2^5     */
+#define  TWO8      0x1.0p8		/* 2^52    */
+#define  TWO10     0x1.0p10		/* 2^10    */
+#define  TWO18     0x1.0p18		/* 2^18    */
+#define  TWO19     0x1.0p19		/* 2^19    */
+#define  TWO23     0x1.0p23		/* 2^23    */
+#define  TWO52     0x1.0p52		/* 2^52    */
 #define  TWO57     0x1.0p57		/* 2^57    */
 #define  TWO71     0x1.0p71		/* 2^71    */
 #define  TWOM1032  0x1.0p-1032		/* 2^-1032 */