This is the mail archive of the
glibc-cvs@sourceware.org
mailing list for the glibc project.
GNU C Library master sources branch master updated. glibc-2.17-356-g82a9811
- From: siddhesh at sourceware dot org
- To: glibc-cvs at sourceware dot org
- Date: 7 Mar 2013 07:38:53 -0000
- Subject: GNU C Library master sources branch master updated. glibc-2.17-356-g82a9811
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GNU C Library master sources".
The branch, master has been updated
via 82a9811d29c00161c7c8ea7f70e2cc30988e192e (commit)
from adbb8027be47b3295367019b2f45863ea3d6c727 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
http://sources.redhat.com/git/gitweb.cgi?p=glibc.git;a=commitdiff;h=82a9811d29c00161c7c8ea7f70e2cc30988e192e
commit 82a9811d29c00161c7c8ea7f70e2cc30988e192e
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date: Thu Mar 7 12:23:29 2013 +0530
Use generic mpa.c code for everything except __mul and __sqr
diff --git a/ChangeLog b/ChangeLog
index dd78188..f10a146 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 8fc2626..0766476 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);
- }
-}
-----------------------------------------------------------------------
Summary of changes:
ChangeLog | 7 +
sysdeps/ieee754/dbl-64/mpa.c | 4 +
sysdeps/powerpc/powerpc32/power4/fpu/mpa.c | 632 +---------------------------
sysdeps/powerpc/powerpc64/power4/fpu/mpa.c | 632 +---------------------------
4 files changed, 19 insertions(+), 1256 deletions(-)
hooks/post-receive
--
GNU C Library master sources