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: Wed, 6 Mar 2013 15:42:31 +0530
- Subject: Re: [PATCH][ppc] Use generic mpa.c code for everything except __sqrand __mul
- References: <20130221104222.GE7508@spoyarek.pnq.redhat.com><20130228060358.GD28148@spoyarek.pnq.redhat.com>
Ping!
On Thu, Feb 28, 2013 at 11:33:58AM +0530, Siddhesh Poyarekar wrote:
> 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);
> > - }
> > -}