This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: [PATCH][ppc] Use generic mpa.c code for everything except __sqrand __mul
- From: Siddhesh Poyarekar <siddhesh at redhat dot com>
- To: libc-alpha at sourceware dot org
- Cc: rsa at us dot ibm dot com
- Date: Thu, 28 Feb 2013 11:33:58 +0530
- Subject: Re: [PATCH][ppc] Use generic mpa.c code for everything except __sqrand __mul
- References: <20130221104222.GE7508@spoyarek.pnq.redhat.com>
Ping!
On Thu, Feb 21, 2013 at 04:12:22PM +0530, Siddhesh Poyarekar wrote:
> 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);
> - }
> -}