Use generic mpa.c code for everything except __mul and __sqr

3 files changed, 12 insertions, 1256 deletions

diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
index 8fc2626f76..0766476544 100644
--- a/sysdeps/ieee754/dbl-64/mpa.c
+++ b/sysdeps/ieee754/dbl-64/mpa.c
@@ -611,6 +611,7 @@ __sub (const mp_no *x, const mp_no *y, mp_no *z, int p)
     }
 }
 
+#ifndef NO__MUL
 /* Multiply *X and *Y and store result in *Z.  X and Y may overlap but not X
    and Z or Y and Z.  For P in [1, 2, 3], the exact result is truncated to P
    digits.  In case P > 3 the error is bounded by 1.001 ULP.  */
@@ -761,7 +762,9 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
   EZ = e;
   Z[0] = X[0] * Y[0];
 }
+#endif
 
+#ifndef NO__SQR
 /* Square *X and store result in *Y.  X and Y may not overlap.  For P in
    [1, 2, 3], the exact result is truncated to P digits.  In case P > 3 the
    error is bounded by 1.001 ULP.  This is a faster special case of
@@ -862,6 +865,7 @@ __sqr (const mp_no *x, mp_no *y, int p)
 
   EY = e;
 }
+#endif
 
 /* Invert *X and store in *Y.  Relative error bound:
    - For P = 2: 1.001 * R ^ (1 - P)
diff --git a/sysdeps/powerpc/powerpc32/power4/fpu/mpa.c b/sysdeps/powerpc/powerpc32/power4/fpu/mpa.c
index b22664772e..ef7ada749a 100644
--- a/sysdeps/powerpc/powerpc32/power4/fpu/mpa.c
+++ b/sysdeps/powerpc/powerpc32/power4/fpu/mpa.c
@@ -17,582 +17,12 @@
  * 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/>.
  */
-/************************************************************************/
-/*  MODULE_NAME: mpa.c                                                  */
-/*                                                                      */
-/*  FUNCTIONS:                                                          */
-/*               mcr                                                    */
-/*               acr                                                    */
-/*               cpy                                                    */
-/*               norm                                                   */
-/*               denorm                                                 */
-/*               mp_dbl                                                 */
-/*               dbl_mp                                                 */
-/*               add_magnitudes                                         */
-/*               sub_magnitudes                                         */
-/*               add                                                    */
-/*               sub                                                    */
-/*               mul                                                    */
-/*               inv                                                    */
-/*               dvd                                                    */
-/*                                                                      */
-/* Arithmetic functions for multiple precision numbers.                 */
-/* Relative errors are bounded                                          */
-/************************************************************************/
 
+/* Define __mul and __sqr and use the rest from generic code.  */
+#define NO__MUL
+#define NO__SQR
 
-#include "endian.h"
-#include "mpa.h"
-#include <sys/param.h>
-
-const mp_no mpone = {1, {1.0, 1.0}};
-const mp_no mptwo = {1, {1.0, 2.0}};
-
-/* Compare mantissa of two multiple precision numbers regardless of the sign
-   and exponent of the numbers.  */
-static int
-mcr (const mp_no *x, const mp_no *y, int p)
-{
-  long i;
-  long p2 = p;
-  for (i = 1; i <= p2; i++)
-    {
-      if (X[i] == Y[i])
-	continue;
-      else if (X[i] > Y[i])
-	return 1;
-      else
-	return -1;
-    }
-  return 0;
-}
-
-/* Compare the absolute values of two multiple precision numbers.  */
-int
-__acr (const mp_no *x, const mp_no *y, int p)
-{
-  long i;
-
-  if (X[0] == ZERO)
-    {
-      if (Y[0] == ZERO)
-	i = 0;
-      else
-	i = -1;
-    }
-  else if (Y[0] == ZERO)
-    i = 1;
-  else
-    {
-      if (EX > EY)
-	i = 1;
-      else if (EX < EY)
-	i = -1;
-      else
-	i = mcr (x, y, p);
-    }
-
-  return i;
-}
-
-/* Copy multiple precision number X into Y.  They could be the same
-   number.  */
-void
-__cpy (const mp_no *x, mp_no *y, int p)
-{
-  long i;
-
-  EY = EX;
-  for (i = 0; i <= p; i++)
-    Y[i] = X[i];
-}
-
-/* Convert a multiple precision number *X into a double precision
-   number *Y, normalized case  (|x| >= 2**(-1022))).  */
-static void
-norm (const mp_no *x, double *y, int p)
-{
-#define R RADIXI
-  long i;
-  double a, c, u, v, z[5];
-  if (p < 5)
-    {
-      if (p == 1)
-	c = X[1];
-      else if (p == 2)
-	c = X[1] + R * X[2];
-      else if (p == 3)
-	c = X[1] + R * (X[2] + R * X[3]);
-      else if (p == 4)
-	c = (X[1] + R * X[2]) + R * R * (X[3] + R * X[4]);
-    }
-  else
-    {
-      for (a = ONE, z[1] = X[1]; z[1] < TWO23;)
-	{
-	  a *= TWO;
-	  z[1] *= TWO;
-	}
-
-      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;
-	}
-
-      u = (z[3] + TWO71) - TWO71;
-      if (u > z[3])
-	u -= TWO19;
-      v = z[3] - u;
-
-      if (v == TWO18)
-	{
-	  if (z[4] == ZERO)
-	    {
-	      for (i = 5; i <= p; i++)
-		{
-		  if (X[i] == ZERO)
-		    continue;
-		  else
-		    {
-		      z[3] += ONE;
-		      break;
-		    }
-		}
-	    }
-	  else
-	    z[3] += ONE;
-	}
-
-      c = (z[1] + R * (z[2] + R * z[3])) / a;
-    }
-
-  c *= X[0];
-
-  for (i = 1; i < EX; i++)
-    c *= RADIX;
-  for (i = 1; i > EX; i--)
-    c *= RADIXI;
-
-  *y = c;
-#undef R
-}
-
-/* Convert a multiple precision number *X into a double precision
-   number *Y, Denormal case  (|x| < 2**(-1022))).  */
-static void
-denorm (const mp_no *x, double *y, int p)
-{
-  long i, k;
-  long p2 = p;
-  double c, u, z[5];
-
-#define R RADIXI
-  if (EX < -44 || (EX == -44 && X[1] < TWO5))
-    {
-      *y = ZERO;
-      return;
-    }
-
-  if (p2 == 1)
-    {
-      if (EX == -42)
-	{
-	  z[1] = X[1] + TWO10;
-	  z[2] = ZERO;
-	  z[3] = ZERO;
-	  k = 3;
-	}
-      else if (EX == -43)
-	{
-	  z[1] = TWO10;
-	  z[2] = X[1];
-	  z[3] = ZERO;
-	  k = 2;
-	}
-      else
-	{
-	  z[1] = TWO10;
-	  z[2] = ZERO;
-	  z[3] = X[1];
-	  k = 1;
-	}
-    }
-  else if (p2 == 2)
-    {
-      if (EX == -42)
-	{
-	  z[1] = X[1] + TWO10;
-	  z[2] = X[2];
-	  z[3] = ZERO;
-	  k = 3;
-	}
-      else if (EX == -43)
-	{
-	  z[1] = TWO10;
-	  z[2] = X[1];
-	  z[3] = X[2];
-	  k = 2;
-	}
-      else
-	{
-	  z[1] = TWO10;
-	  z[2] = ZERO;
-	  z[3] = X[1];
-	  k = 1;
-	}
-    }
-  else
-    {
-      if (EX == -42)
-	{
-	  z[1] = X[1] + TWO10;
-	  z[2] = X[2];
-	  k = 3;
-	}
-      else if (EX == -43)
-	{
-	  z[1] = TWO10;
-	  z[2] = X[1];
-	  k = 2;
-	}
-      else
-	{
-	  z[1] = TWO10;
-	  z[2] = ZERO;
-	  k = 1;
-	}
-      z[3] = X[k];
-    }
-
-  u = (z[3] + TWO57) - TWO57;
-  if (u > z[3])
-    u -= TWO5;
-
-  if (u == z[3])
-    {
-      for (i = k + 1; i <= p2; i++)
-	{
-	  if (X[i] == ZERO)
-	    continue;
-	  else
-	    {
-	      z[3] += ONE;
-	      break;
-	    }
-	}
-    }
-
-  c = X[0] * ((z[1] + R * (z[2] + R * z[3])) - TWO10);
-
-  *y = c * TWOM1032;
-#undef R
-}
-
-/* Convert multiple precision number *X into double precision number *Y.  The
-   result is correctly rounded to the nearest/even.  */
-void
-__mp_dbl (const mp_no *x, double *y, int p)
-{
-  if (X[0] == ZERO)
-    {
-      *y = ZERO;
-      return;
-    }
-
-  if (__glibc_likely (EX > -42 || (EX == -42 && X[1] >= TWO10)))
-    norm (x, y, p);
-  else
-    denorm (x, y, p);
-}
-
-/* Get the multiple precision equivalent of X into *Y.  If the precision is too
-   small, the result is truncated.  */
-void
-__dbl_mp (double x, mp_no *y, int p)
-{
-  long i, n;
-  long p2 = p;
-  double u;
-
-  /* Sign.  */
-  if (x == ZERO)
-    {
-      Y[0] = ZERO;
-      return;
-    }
-  else if (x > ZERO)
-    Y[0] = ONE;
-  else
-    {
-      Y[0] = MONE;
-      x = -x;
-    }
-
-  /* Exponent.  */
-  for (EY = ONE; x >= RADIX; EY += ONE)
-    x *= RADIXI;
-  for (; x < ONE; EY -= ONE)
-    x *= RADIX;
-
-  /* Digits.  */
-  n = MIN (p2, 4);
-  for (i = 1; i <= n; i++)
-    {
-      u = (x + TWO52) - TWO52;
-      if (u > x)
-	u -= ONE;
-      Y[i] = u;
-      x -= u;
-      x *= RADIX;
-    }
-  for (; i <= p2; i++)
-    Y[i] = ZERO;
-}
-
-/* Add magnitudes of *X and *Y assuming that abs (*X) >= abs (*Y) > 0.  The
-   sign of the sum *Z is not changed.  X and Y may overlap but not X and Z or
-   Y and Z.  No guard digit is used.  The result equals the exact sum,
-   truncated.  */
-static void
-add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  long i, j, k;
-  long p2 = p;
-  double zk;
-
-  EZ = EX;
-
-  i = p2;
-  j = p2 + EY - EX;
-  k = p2 + 1;
-
-  if (__glibc_unlikely (j < 1))
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  zk = ZERO;
-
-  for (; j > 0; i--, j--)
-    {
-      zk += X[i] + Y[j];
-      if (zk >= RADIX)
-	{
-	  Z[k--] = zk - RADIX;
-	  zk = ONE;
-	}
-      else
-        {
-	  Z[k--] = zk;
-	  zk = ZERO;
-	}
-    }
-
-  for (; i > 0; i--)
-    {
-      zk += X[i];
-      if (zk >= RADIX)
-	{
-	  Z[k--] = zk - RADIX;
-	  zk = ONE;
-	}
-      else
-        {
-	  Z[k--] = zk;
-	  zk = ZERO;
-	}
-    }
-
-  if (zk == ZERO)
-    {
-      for (i = 1; i <= p2; i++)
-	Z[i] = Z[i + 1];
-    }
-  else
-    {
-      Z[1] = zk;
-      EZ += ONE;
-    }
-}
-
-/* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0.
-   The sign of the difference *Z is not changed.  X and Y may overlap but not X
-   and Z or Y and Z.  One guard digit is used.  The error is less than one
-   ULP.  */
-static void
-sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  long i, j, k;
-  long p2 = p;
-  double zk;
-
-  EZ = EX;
-  i = p2;
-  j = p2 + EY - EX;
-  k = p2;
-
-  /* Y is too small compared to X, copy X over to the result.  */
-  if (__glibc_unlikely (j < 1))
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  /* The relevant least significant digit in Y is non-zero, so we factor it in
-     to enhance accuracy.  */
-  if (j < p2 && Y[j + 1] > ZERO)
-    {
-      Z[k + 1] = RADIX - Y[j + 1];
-      zk = MONE;
-    }
-  else
-    zk = Z[k + 1] = ZERO;
-
-  /* Subtract and borrow.  */
-  for (; j > 0; i--, j--)
-    {
-      zk += (X[i] - Y[j]);
-      if (zk < ZERO)
-	{
-	  Z[k--] = zk + RADIX;
-	  zk = MONE;
-	}
-      else
-        {
-	  Z[k--] = zk;
-	  zk = ZERO;
-	}
-    }
-
-  /* We're done with digits from Y, so it's just digits in X.  */
-  for (; i > 0; i--)
-    {
-      zk += X[i];
-      if (zk < ZERO)
-	{
-	  Z[k--] = zk + RADIX;
-	  zk = MONE;
-	}
-      else
-        {
-	  Z[k--] = zk;
-	  zk = ZERO;
-	}
-    }
-
-  /* Normalize.  */
-  for (i = 1; Z[i] == ZERO; i++);
-  EZ = EZ - i + 1;
-  for (k = 1; i <= p2 + 1;)
-    Z[k++] = Z[i++];
-  for (; k <= p2;)
-    Z[k++] = ZERO;
-}
-
-/* Add *X and *Y and store the result in *Z.  X and Y may overlap, but not X
-   and Z or Y and Z.  One guard digit is used.  The error is less than one
-   ULP.  */
-void
-__add (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  int n;
-
-  if (X[0] == ZERO)
-    {
-      __cpy (y, z, p);
-      return;
-    }
-  else if (Y[0] == ZERO)
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  if (X[0] == Y[0])
-    {
-      if (__acr (x, y, p) > 0)
-	{
-	  add_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else
-	{
-	  add_magnitudes (y, x, z, p);
-	  Z[0] = Y[0];
-	}
-    }
-  else
-    {
-      if ((n = __acr (x, y, p)) == 1)
-	{
-	  sub_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else if (n == -1)
-	{
-	  sub_magnitudes (y, x, z, p);
-	  Z[0] = Y[0];
-	}
-      else
-	Z[0] = ZERO;
-    }
-}
-
-/* Subtract *Y from *X and return the result in *Z.  X and Y may overlap but
-   not X and Z or Y and Z.  One guard digit is used.  The error is less than
-   one ULP.  */
-void
-__sub (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  int n;
-
-  if (X[0] == ZERO)
-    {
-      __cpy (y, z, p);
-      Z[0] = -Z[0];
-      return;
-    }
-  else if (Y[0] == ZERO)
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  if (X[0] != Y[0])
-    {
-      if (__acr (x, y, p) > 0)
-	{
-	  add_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else
-	{
-	  add_magnitudes (y, x, z, p);
-	  Z[0] = -Y[0];
-	}
-    }
-  else
-    {
-      if ((n = __acr (x, y, p)) == 1)
-	{
-	  sub_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else if (n == -1)
-	{
-	  sub_magnitudes (y, x, z, p);
-	  Z[0] = -Y[0];
-	}
-      else
-	Z[0] = ZERO;
-    }
-}
+#include <sysdeps/ieee754/dbl-64/mpa.c>
 
 /* Multiply *X and *Y and store result in *Z.  X and Y may overlap but not X
    and Z or Y and Z.  For P in [1, 2, 3], the exact result is truncated to P
@@ -781,57 +211,3 @@ __sqr (const mp_no *x, mp_no *y, int p)
       EY--;
     }
 }
-
-/* Invert *X and store in *Y.  Relative error bound:
-   - For P = 2: 1.001 * R ^ (1 - P)
-   - For P = 3: 1.063 * R ^ (1 - P)
-   - For P > 3: 2.001 * R ^ (1 - P)
-
-   *X = 0 is not permissible.  */
-static void
-__inv (const mp_no *x, mp_no *y, int p)
-{
-  long i;
-  double t;
-  mp_no z, w;
-  static const int np1[] =
-    { 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
-    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
-  };
-
-  __cpy (x, &z, p);
-  z.e = 0;
-  __mp_dbl (&z, &t, p);
-  t = ONE / t;
-  __dbl_mp (t, y, p);
-  EY -= EX;
-
-  for (i = 0; i < np1[p]; i++)
-    {
-      __cpy (y, &w, p);
-      __mul (x, &w, y, p);
-      __sub (&mptwo, y, &z, p);
-      __mul (&w, &z, y, p);
-    }
-}
-
-/* Divide *X by *Y and store result in *Z.  X and Y may overlap but not X and Z
-   or Y and Z.  Relative error bound:
-   - For P = 2: 2.001 * R ^ (1 - P)
-   - For P = 3: 2.063 * R ^ (1 - P)
-   - For P > 3: 3.001 * R ^ (1 - P)
-
-   *X = 0 is not permissible.  */
-void
-__dvd (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  mp_no w;
-
-  if (X[0] == ZERO)
-    Z[0] = ZERO;
-  else
-    {
-      __inv (y, &w, p);
-      __mul (x, &w, z, p);
-    }
-}
diff --git a/sysdeps/powerpc/powerpc64/power4/fpu/mpa.c b/sysdeps/powerpc/powerpc64/power4/fpu/mpa.c
index b22664772e..ef7ada749a 100644
--- a/sysdeps/powerpc/powerpc64/power4/fpu/mpa.c
+++ b/sysdeps/powerpc/powerpc64/power4/fpu/mpa.c
@@ -17,582 +17,12 @@
  * 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/>.
  */
-/************************************************************************/
-/*  MODULE_NAME: mpa.c                                                  */
-/*                                                                      */
-/*  FUNCTIONS:                                                          */
-/*               mcr                                                    */
-/*               acr                                                    */
-/*               cpy                                                    */
-/*               norm                                                   */
-/*               denorm                                                 */
-/*               mp_dbl                                                 */
-/*               dbl_mp                                                 */
-/*               add_magnitudes                                         */
-/*               sub_magnitudes                                         */
-/*               add                                                    */
-/*               sub                                                    */
-/*               mul                                                    */
-/*               inv                                                    */
-/*               dvd                                                    */
-/*                                                                      */
-/* Arithmetic functions for multiple precision numbers.                 */
-/* Relative errors are bounded                                          */
-/************************************************************************/
 
+/* Define __mul and __sqr and use the rest from generic code.  */
+#define NO__MUL
+#define NO__SQR
 
-#include "endian.h"
-#include "mpa.h"
-#include <sys/param.h>
-
-const mp_no mpone = {1, {1.0, 1.0}};
-const mp_no mptwo = {1, {1.0, 2.0}};
-
-/* Compare mantissa of two multiple precision numbers regardless of the sign
-   and exponent of the numbers.  */
-static int
-mcr (const mp_no *x, const mp_no *y, int p)
-{
-  long i;
-  long p2 = p;
-  for (i = 1; i <= p2; i++)
-    {
-      if (X[i] == Y[i])
-	continue;
-      else if (X[i] > Y[i])
-	return 1;
-      else
-	return -1;
-    }
-  return 0;
-}
-
-/* Compare the absolute values of two multiple precision numbers.  */
-int
-__acr (const mp_no *x, const mp_no *y, int p)
-{
-  long i;
-
-  if (X[0] == ZERO)
-    {
-      if (Y[0] == ZERO)
-	i = 0;
-      else
-	i = -1;
-    }
-  else if (Y[0] == ZERO)
-    i = 1;
-  else
-    {
-      if (EX > EY)
-	i = 1;
-      else if (EX < EY)
-	i = -1;
-      else
-	i = mcr (x, y, p);
-    }
-
-  return i;
-}
-
-/* Copy multiple precision number X into Y.  They could be the same
-   number.  */
-void
-__cpy (const mp_no *x, mp_no *y, int p)
-{
-  long i;
-
-  EY = EX;
-  for (i = 0; i <= p; i++)
-    Y[i] = X[i];
-}
-
-/* Convert a multiple precision number *X into a double precision
-   number *Y, normalized case  (|x| >= 2**(-1022))).  */
-static void
-norm (const mp_no *x, double *y, int p)
-{
-#define R RADIXI
-  long i;
-  double a, c, u, v, z[5];
-  if (p < 5)
-    {
-      if (p == 1)
-	c = X[1];
-      else if (p == 2)
-	c = X[1] + R * X[2];
-      else if (p == 3)
-	c = X[1] + R * (X[2] + R * X[3]);
-      else if (p == 4)
-	c = (X[1] + R * X[2]) + R * R * (X[3] + R * X[4]);
-    }
-  else
-    {
-      for (a = ONE, z[1] = X[1]; z[1] < TWO23;)
-	{
-	  a *= TWO;
-	  z[1] *= TWO;
-	}
-
-      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;
-	}
-
-      u = (z[3] + TWO71) - TWO71;
-      if (u > z[3])
-	u -= TWO19;
-      v = z[3] - u;
-
-      if (v == TWO18)
-	{
-	  if (z[4] == ZERO)
-	    {
-	      for (i = 5; i <= p; i++)
-		{
-		  if (X[i] == ZERO)
-		    continue;
-		  else
-		    {
-		      z[3] += ONE;
-		      break;
-		    }
-		}
-	    }
-	  else
-	    z[3] += ONE;
-	}
-
-      c = (z[1] + R * (z[2] + R * z[3])) / a;
-    }
-
-  c *= X[0];
-
-  for (i = 1; i < EX; i++)
-    c *= RADIX;
-  for (i = 1; i > EX; i--)
-    c *= RADIXI;
-
-  *y = c;
-#undef R
-}
-
-/* Convert a multiple precision number *X into a double precision
-   number *Y, Denormal case  (|x| < 2**(-1022))).  */
-static void
-denorm (const mp_no *x, double *y, int p)
-{
-  long i, k;
-  long p2 = p;
-  double c, u, z[5];
-
-#define R RADIXI
-  if (EX < -44 || (EX == -44 && X[1] < TWO5))
-    {
-      *y = ZERO;
-      return;
-    }
-
-  if (p2 == 1)
-    {
-      if (EX == -42)
-	{
-	  z[1] = X[1] + TWO10;
-	  z[2] = ZERO;
-	  z[3] = ZERO;
-	  k = 3;
-	}
-      else if (EX == -43)
-	{
-	  z[1] = TWO10;
-	  z[2] = X[1];
-	  z[3] = ZERO;
-	  k = 2;
-	}
-      else
-	{
-	  z[1] = TWO10;
-	  z[2] = ZERO;
-	  z[3] = X[1];
-	  k = 1;
-	}
-    }
-  else if (p2 == 2)
-    {
-      if (EX == -42)
-	{
-	  z[1] = X[1] + TWO10;
-	  z[2] = X[2];
-	  z[3] = ZERO;
-	  k = 3;
-	}
-      else if (EX == -43)
-	{
-	  z[1] = TWO10;
-	  z[2] = X[1];
-	  z[3] = X[2];
-	  k = 2;
-	}
-      else
-	{
-	  z[1] = TWO10;
-	  z[2] = ZERO;
-	  z[3] = X[1];
-	  k = 1;
-	}
-    }
-  else
-    {
-      if (EX == -42)
-	{
-	  z[1] = X[1] + TWO10;
-	  z[2] = X[2];
-	  k = 3;
-	}
-      else if (EX == -43)
-	{
-	  z[1] = TWO10;
-	  z[2] = X[1];
-	  k = 2;
-	}
-      else
-	{
-	  z[1] = TWO10;
-	  z[2] = ZERO;
-	  k = 1;
-	}
-      z[3] = X[k];
-    }
-
-  u = (z[3] + TWO57) - TWO57;
-  if (u > z[3])
-    u -= TWO5;
-
-  if (u == z[3])
-    {
-      for (i = k + 1; i <= p2; i++)
-	{
-	  if (X[i] == ZERO)
-	    continue;
-	  else
-	    {
-	      z[3] += ONE;
-	      break;
-	    }
-	}
-    }
-
-  c = X[0] * ((z[1] + R * (z[2] + R * z[3])) - TWO10);
-
-  *y = c * TWOM1032;
-#undef R
-}
-
-/* Convert multiple precision number *X into double precision number *Y.  The
-   result is correctly rounded to the nearest/even.  */
-void
-__mp_dbl (const mp_no *x, double *y, int p)
-{
-  if (X[0] == ZERO)
-    {
-      *y = ZERO;
-      return;
-    }
-
-  if (__glibc_likely (EX > -42 || (EX == -42 && X[1] >= TWO10)))
-    norm (x, y, p);
-  else
-    denorm (x, y, p);
-}
-
-/* Get the multiple precision equivalent of X into *Y.  If the precision is too
-   small, the result is truncated.  */
-void
-__dbl_mp (double x, mp_no *y, int p)
-{
-  long i, n;
-  long p2 = p;
-  double u;
-
-  /* Sign.  */
-  if (x == ZERO)
-    {
-      Y[0] = ZERO;
-      return;
-    }
-  else if (x > ZERO)
-    Y[0] = ONE;
-  else
-    {
-      Y[0] = MONE;
-      x = -x;
-    }
-
-  /* Exponent.  */
-  for (EY = ONE; x >= RADIX; EY += ONE)
-    x *= RADIXI;
-  for (; x < ONE; EY -= ONE)
-    x *= RADIX;
-
-  /* Digits.  */
-  n = MIN (p2, 4);
-  for (i = 1; i <= n; i++)
-    {
-      u = (x + TWO52) - TWO52;
-      if (u > x)
-	u -= ONE;
-      Y[i] = u;
-      x -= u;
-      x *= RADIX;
-    }
-  for (; i <= p2; i++)
-    Y[i] = ZERO;
-}
-
-/* Add magnitudes of *X and *Y assuming that abs (*X) >= abs (*Y) > 0.  The
-   sign of the sum *Z is not changed.  X and Y may overlap but not X and Z or
-   Y and Z.  No guard digit is used.  The result equals the exact sum,
-   truncated.  */
-static void
-add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  long i, j, k;
-  long p2 = p;
-  double zk;
-
-  EZ = EX;
-
-  i = p2;
-  j = p2 + EY - EX;
-  k = p2 + 1;
-
-  if (__glibc_unlikely (j < 1))
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  zk = ZERO;
-
-  for (; j > 0; i--, j--)
-    {
-      zk += X[i] + Y[j];
-      if (zk >= RADIX)
-	{
-	  Z[k--] = zk - RADIX;
-	  zk = ONE;
-	}
-      else
-        {
-	  Z[k--] = zk;
-	  zk = ZERO;
-	}
-    }
-
-  for (; i > 0; i--)
-    {
-      zk += X[i];
-      if (zk >= RADIX)
-	{
-	  Z[k--] = zk - RADIX;
-	  zk = ONE;
-	}
-      else
-        {
-	  Z[k--] = zk;
-	  zk = ZERO;
-	}
-    }
-
-  if (zk == ZERO)
-    {
-      for (i = 1; i <= p2; i++)
-	Z[i] = Z[i + 1];
-    }
-  else
-    {
-      Z[1] = zk;
-      EZ += ONE;
-    }
-}
-
-/* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0.
-   The sign of the difference *Z is not changed.  X and Y may overlap but not X
-   and Z or Y and Z.  One guard digit is used.  The error is less than one
-   ULP.  */
-static void
-sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  long i, j, k;
-  long p2 = p;
-  double zk;
-
-  EZ = EX;
-  i = p2;
-  j = p2 + EY - EX;
-  k = p2;
-
-  /* Y is too small compared to X, copy X over to the result.  */
-  if (__glibc_unlikely (j < 1))
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  /* The relevant least significant digit in Y is non-zero, so we factor it in
-     to enhance accuracy.  */
-  if (j < p2 && Y[j + 1] > ZERO)
-    {
-      Z[k + 1] = RADIX - Y[j + 1];
-      zk = MONE;
-    }
-  else
-    zk = Z[k + 1] = ZERO;
-
-  /* Subtract and borrow.  */
-  for (; j > 0; i--, j--)
-    {
-      zk += (X[i] - Y[j]);
-      if (zk < ZERO)
-	{
-	  Z[k--] = zk + RADIX;
-	  zk = MONE;
-	}
-      else
-        {
-	  Z[k--] = zk;
-	  zk = ZERO;
-	}
-    }
-
-  /* We're done with digits from Y, so it's just digits in X.  */
-  for (; i > 0; i--)
-    {
-      zk += X[i];
-      if (zk < ZERO)
-	{
-	  Z[k--] = zk + RADIX;
-	  zk = MONE;
-	}
-      else
-        {
-	  Z[k--] = zk;
-	  zk = ZERO;
-	}
-    }
-
-  /* Normalize.  */
-  for (i = 1; Z[i] == ZERO; i++);
-  EZ = EZ - i + 1;
-  for (k = 1; i <= p2 + 1;)
-    Z[k++] = Z[i++];
-  for (; k <= p2;)
-    Z[k++] = ZERO;
-}
-
-/* Add *X and *Y and store the result in *Z.  X and Y may overlap, but not X
-   and Z or Y and Z.  One guard digit is used.  The error is less than one
-   ULP.  */
-void
-__add (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  int n;
-
-  if (X[0] == ZERO)
-    {
-      __cpy (y, z, p);
-      return;
-    }
-  else if (Y[0] == ZERO)
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  if (X[0] == Y[0])
-    {
-      if (__acr (x, y, p) > 0)
-	{
-	  add_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else
-	{
-	  add_magnitudes (y, x, z, p);
-	  Z[0] = Y[0];
-	}
-    }
-  else
-    {
-      if ((n = __acr (x, y, p)) == 1)
-	{
-	  sub_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else if (n == -1)
-	{
-	  sub_magnitudes (y, x, z, p);
-	  Z[0] = Y[0];
-	}
-      else
-	Z[0] = ZERO;
-    }
-}
-
-/* Subtract *Y from *X and return the result in *Z.  X and Y may overlap but
-   not X and Z or Y and Z.  One guard digit is used.  The error is less than
-   one ULP.  */
-void
-__sub (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  int n;
-
-  if (X[0] == ZERO)
-    {
-      __cpy (y, z, p);
-      Z[0] = -Z[0];
-      return;
-    }
-  else if (Y[0] == ZERO)
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  if (X[0] != Y[0])
-    {
-      if (__acr (x, y, p) > 0)
-	{
-	  add_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else
-	{
-	  add_magnitudes (y, x, z, p);
-	  Z[0] = -Y[0];
-	}
-    }
-  else
-    {
-      if ((n = __acr (x, y, p)) == 1)
-	{
-	  sub_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else if (n == -1)
-	{
-	  sub_magnitudes (y, x, z, p);
-	  Z[0] = -Y[0];
-	}
-      else
-	Z[0] = ZERO;
-    }
-}
+#include <sysdeps/ieee754/dbl-64/mpa.c>
 
 /* Multiply *X and *Y and store result in *Z.  X and Y may overlap but not X
    and Z or Y and Z.  For P in [1, 2, 3], the exact result is truncated to P
@@ -781,57 +211,3 @@ __sqr (const mp_no *x, mp_no *y, int p)
       EY--;
     }
 }
-
-/* Invert *X and store in *Y.  Relative error bound:
-   - For P = 2: 1.001 * R ^ (1 - P)
-   - For P = 3: 1.063 * R ^ (1 - P)
-   - For P > 3: 2.001 * R ^ (1 - P)
-
-   *X = 0 is not permissible.  */
-static void
-__inv (const mp_no *x, mp_no *y, int p)
-{
-  long i;
-  double t;
-  mp_no z, w;
-  static const int np1[] =
-    { 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
-    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
-  };
-
-  __cpy (x, &z, p);
-  z.e = 0;
-  __mp_dbl (&z, &t, p);
-  t = ONE / t;
-  __dbl_mp (t, y, p);
-  EY -= EX;
-
-  for (i = 0; i < np1[p]; i++)
-    {
-      __cpy (y, &w, p);
-      __mul (x, &w, y, p);
-      __sub (&mptwo, y, &z, p);
-      __mul (&w, &z, y, p);
-    }
-}
-
-/* Divide *X by *Y and store result in *Z.  X and Y may overlap but not X and Z
-   or Y and Z.  Relative error bound:
-   - For P = 2: 2.001 * R ^ (1 - P)
-   - For P = 3: 2.063 * R ^ (1 - P)
-   - For P > 3: 3.001 * R ^ (1 - P)
-
-   *X = 0 is not permissible.  */
-void
-__dvd (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  mp_no w;
-
-  if (X[0] == ZERO)
-    Z[0] = ZERO;
-  else
-    {
-      __inv (y, &w, p);
-      __mul (x, &w, z, p);
-    }
-}