This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[PATCH][ppc] Use generic mpa.c code for everything except __sqr and__mul


Hi,

Here's a patch that applies on top of all the patches I've posted so
far to make the powerpc mpa.c use much of the generic mpa.c,
overriding only the __mul and __sqr functions since the powerpc
implementation of those functions currrently need to be different.  I
have compared the generated code using objdump to verify that the
generated code does not change in powerpc due to this.  The only thing
that changes is the order of function definition in the binary between
__dvd and __sqr.

The intent of doing this is to ensure that only functions that are
really necessary to be powerpc-only remain duplicated and everything
else is common code.

Siddhesh

	* sysdeps/ieee754/dbl-64/mpa.c [!NO__MUL]: Define __mul.
	[!NO__SQR]: Define __sqr.
	* sysdeps/powerpc/powerpc32/power4/fpu/mpa.c: define NO__MUL
	and NO__SQR.  Remove all code except __mul and __sqr.  Include
	sysdeps/ieee754/dbl-64/mpa.c.
	* sysdeps/powerpc/powerpc64/power4/fpu/mpa.c: Likewise.

diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
index cad227a..c246c8c 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 b226647..ef7ada7 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 b226647..ef7ada7 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);
-    }
-}


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]