diff options
| -rw-r--r-- | ChangeLog | 7 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpa.c | 4 | ||||
| -rw-r--r-- | sysdeps/powerpc/powerpc32/power4/fpu/mpa.c | 632 | ||||
| -rw-r--r-- | sysdeps/powerpc/powerpc64/power4/fpu/mpa.c | 632 |
4 files changed, 19 insertions, 1256 deletions
diff --git a/ChangeLog b/ChangeLog index dd781886c4..f10a14684b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,12 @@ 2013-03-07 Siddhesh Poyarekar <siddhesh@redhat.com> + * 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. + [BZ #12723] * posix/Makefile (tests): Add tst-pathconf. * posix/tst-pathconf.c: New test case. 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); - } -} |