aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--ChangeLog7
-rw-r--r--sysdeps/ieee754/dbl-64/mpa-arch.h47
-rw-r--r--sysdeps/ieee754/dbl-64/mpa.c78
-rw-r--r--sysdeps/ieee754/dbl-64/mpa.h27
-rw-r--r--sysdeps/powerpc/power4/fpu/mpa-arch.h56
5 files changed, 67 insertions, 148 deletions
diff --git a/ChangeLog b/ChangeLog
index fe92b82479..dc84119138 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+2013-03-15 Siddhesh Poyarekar <siddhesh@redhat.com>
+
+ * sysdeps/ieee754/dbl-64/mpa-arch.h: Remove.
+ * sysdeps/ieee754/dbl-64/mpa.c: Revert last change.
+ * sysdeps/ieee754/dbl-64/mpa.h: Revert last change.
+ * sysdeps/powerpc/power4/fpu/mpa-arch.h: Remove.
+
2013-03-15 Adhemerval Zanella <azanella@linux.vnet.ibm.com>
* sysdeps/unix/sysv/linux/powerpc/bits/libc-vdso.h (VDSO_IFUNC_RET): Add
diff --git a/sysdeps/ieee754/dbl-64/mpa-arch.h b/sysdeps/ieee754/dbl-64/mpa-arch.h
deleted file mode 100644
index 7de9d51ae2..0000000000
--- a/sysdeps/ieee754/dbl-64/mpa-arch.h
+++ /dev/null
@@ -1,47 +0,0 @@
-/* Overridable constants and operations.
- Copyright (C) 2013 Free Software Foundation, Inc.
-
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation; either version 2.1 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU Lesser General Public License for more details.
-
- 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/>. */
-
-#include <stdint.h>
-
-typedef long mantissa_t;
-typedef int64_t mantissa_store_t;
-
-#define TWOPOW(i) (1L << i)
-
-#define RADIX_EXP 24
-#define RADIX TWOPOW (RADIX_EXP) /* 2^24 */
-
-/* Divide D by RADIX and put the remainder in R. D must be a non-negative
- integral value. */
-#define DIV_RADIX(d, r) \
- ({ \
- r = d & (RADIX - 1); \
- d >>= RADIX_EXP; \
- })
-
-/* Put the integer component of a double X in R and retain the fraction in
- X. This is used in extracting mantissa digits for MP_NO by using the
- integer portion of the current value of the number as the current mantissa
- digit and then scaling by RADIX to get the next mantissa digit in the same
- manner. */
-#define INTEGER_OF(x, i) \
- ({ \
- i = (mantissa_t) x; \
- x -= i; \
- })
-
-/* Align IN down to F. The code assumes that F is a power of two. */
-#define ALIGN_DOWN_TO(in, f) ((in) & -(f))
diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
index 860e859ae8..0766476544 100644
--- a/sysdeps/ieee754/dbl-64/mpa.c
+++ b/sysdeps/ieee754/dbl-64/mpa.c
@@ -125,8 +125,7 @@ norm (const mp_no *x, double *y, int p)
{
#define R RADIXI
long i;
- double c;
- mantissa_t a, u, v, z[5];
+ double a, c, u, v, z[5];
if (p < 5)
{
if (p == 1)
@@ -148,14 +147,17 @@ norm (const mp_no *x, double *y, int p)
for (i = 2; i < 5; i++)
{
- mantissa_t d, r;
- d = X[i] * a;
- DIV_RADIX (d, r);
- z[i] = r;
- z[i - 1] += d;
+ 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 = ALIGN_DOWN_TO (z[3], TWO19);
+ u = (z[3] + TWO71) - TWO71;
+ if (u > z[3])
+ u -= TWO19;
v = z[3] - u;
if (v == TWO18)
@@ -198,8 +200,7 @@ denorm (const mp_no *x, double *y, int p)
{
long i, k;
long p2 = p;
- double c;
- mantissa_t u, z[5];
+ double c, u, z[5];
#define R RADIXI
if (EX < -44 || (EX == -44 && X[1] < TWO5))
@@ -279,7 +280,9 @@ denorm (const mp_no *x, double *y, int p)
z[3] = X[k];
}
- u = ALIGN_DOWN_TO (z[3], TWO5);
+ u = (z[3] + TWO57) - TWO57;
+ if (u > z[3])
+ u -= TWO5;
if (u == z[3])
{
@@ -327,6 +330,7 @@ __dbl_mp (double x, mp_no *y, int p)
{
long i, n;
long p2 = p;
+ double u;
/* Sign. */
if (x == ZERO)
@@ -352,7 +356,11 @@ __dbl_mp (double x, mp_no *y, int p)
n = MIN (p2, 4);
for (i = 1; i <= n; i++)
{
- INTEGER_OF (x, Y[i]);
+ u = (x + TWO52) - TWO52;
+ if (u > x)
+ u -= ONE;
+ Y[i] = u;
+ x -= u;
x *= RADIX;
}
for (; i <= p2; i++)
@@ -369,7 +377,7 @@ add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
{
long i, j, k;
long p2 = p;
- mantissa_t zk;
+ double zk;
EZ = EX;
@@ -437,7 +445,7 @@ sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
{
long i, j, k;
long p2 = p;
- mantissa_t zk;
+ double zk;
EZ = EX;
i = p2;
@@ -613,9 +621,9 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
{
long i, j, k, ip, ip2;
long p2 = p;
- mantissa_store_t zk;
+ double u, zk;
const mp_no *a;
- mantissa_store_t *diag;
+ double *diag;
/* Is z=0? */
if (__glibc_unlikely (X[0] * Y[0] == ZERO))
@@ -672,8 +680,8 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
/* Precompute sums of diagonal elements so that we can directly use them
later. See the next comment to know we why need them. */
- diag = alloca (k * sizeof (mantissa_store_t));
- mantissa_store_t d = ZERO;
+ diag = alloca (k * sizeof (double));
+ double d = ZERO;
for (i = 1; i <= ip; i++)
{
d += X[i] * Y[i];
@@ -696,8 +704,11 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
zk -= diag[k - 1];
- DIV_RADIX (zk, Z[k]);
- k--;
+ u = (zk + CUTTER) - CUTTER;
+ if (u > zk)
+ u -= RADIX;
+ Z[k--] = zk - u;
+ zk = u * RADIXI;
}
/* The real deal. Mantissa digit Z[k] is the sum of all X[i] * Y[j] where i
@@ -727,8 +738,11 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
zk -= diag[k - 1];
- DIV_RADIX (zk, Z[k]);
- k--;
+ u = (zk + CUTTER) - CUTTER;
+ if (u > zk)
+ u -= RADIX;
+ Z[k--] = zk - u;
+ zk = u * RADIXI;
}
Z[k] = zk;
@@ -760,7 +774,7 @@ SECTION
__sqr (const mp_no *x, mp_no *y, int p)
{
long i, j, k, ip;
- mantissa_store_t yk;
+ double u, yk;
/* Is z=0? */
if (__glibc_unlikely (X[0] == ZERO))
@@ -784,7 +798,7 @@ __sqr (const mp_no *x, mp_no *y, int p)
while (k > p)
{
- mantissa_store_t yk2 = 0;
+ double yk2 = 0.0;
long lim = k / 2;
if (k % 2 == 0)
@@ -804,13 +818,16 @@ __sqr (const mp_no *x, mp_no *y, int p)
yk += 2.0 * yk2;
- DIV_RADIX (yk, Y[k]);
- k--;
+ u = (yk + CUTTER) - CUTTER;
+ if (u > yk)
+ u -= RADIX;
+ Y[k--] = yk - u;
+ yk = u * RADIXI;
}
while (k > 1)
{
- mantissa_store_t yk2 = 0;
+ double yk2 = 0.0;
long lim = k / 2;
if (k % 2 == 0)
@@ -822,8 +839,11 @@ __sqr (const mp_no *x, mp_no *y, int p)
yk += 2.0 * yk2;
- DIV_RADIX (yk, Y[k]);
- k--;
+ u = (yk + CUTTER) - CUTTER;
+ if (u > yk)
+ u -= RADIX;
+ Y[k--] = yk - u;
+ yk = u * RADIXI;
}
Y[k] = yk;
diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h
index 54044a0586..168b334ed0 100644
--- a/sysdeps/ieee754/dbl-64/mpa.h
+++ b/sysdeps/ieee754/dbl-64/mpa.h
@@ -35,7 +35,6 @@
/* Common types and definition */
/************************************************************************/
-#include <mpa-arch.h>
/* The mp_no structure holds the details of a multi-precision floating point
number.
@@ -62,7 +61,7 @@
typedef struct
{
int e;
- mantissa_t d[40];
+ double d[40];
} mp_no;
typedef union
@@ -83,13 +82,9 @@ extern const mp_no mptwo;
#define ABS(x) ((x) < 0 ? -(x) : (x))
-#ifndef RADIXI
-# define RADIXI 0x1.0p-24 /* 2^-24 */
-#endif
-
-#ifndef TWO52
-# define TWO52 0x1.0p52 /* 2^52 */
-#endif
+#define RADIX 0x1.0p24 /* 2^24 */
+#define RADIXI 0x1.0p-24 /* 2^-24 */
+#define CUTTER 0x1.0p76 /* 2^76 */
#define ZERO 0.0 /* 0 */
#define MZERO -0.0 /* 0 with the sign bit set */
@@ -97,13 +92,13 @@ extern const mp_no mptwo;
#define MONE -1.0 /* -1 */
#define TWO 2.0 /* 2 */
-#define TWO5 TWOPOW (5) /* 2^5 */
-#define TWO8 TWOPOW (8) /* 2^52 */
-#define TWO10 TWOPOW (10) /* 2^10 */
-#define TWO18 TWOPOW (18) /* 2^18 */
-#define TWO19 TWOPOW (19) /* 2^19 */
-#define TWO23 TWOPOW (23) /* 2^23 */
-
+#define TWO5 0x1.0p5 /* 2^5 */
+#define TWO8 0x1.0p8 /* 2^52 */
+#define TWO10 0x1.0p10 /* 2^10 */
+#define TWO18 0x1.0p18 /* 2^18 */
+#define TWO19 0x1.0p19 /* 2^19 */
+#define TWO23 0x1.0p23 /* 2^23 */
+#define TWO52 0x1.0p52 /* 2^52 */
#define TWO57 0x1.0p57 /* 2^57 */
#define TWO71 0x1.0p71 /* 2^71 */
#define TWOM1032 0x1.0p-1032 /* 2^-1032 */
diff --git a/sysdeps/powerpc/power4/fpu/mpa-arch.h b/sysdeps/powerpc/power4/fpu/mpa-arch.h
deleted file mode 100644
index 9007c9d0ca..0000000000
--- a/sysdeps/powerpc/power4/fpu/mpa-arch.h
+++ /dev/null
@@ -1,56 +0,0 @@
-/* Overridable constants and operations.
- Copyright (C) 2013 Free Software Foundation, Inc.
-
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation; either version 2.1 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU Lesser General Public License for more details.
-
- 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/>. */
-
-typedef double mantissa_t;
-typedef double mantissa_store_t;
-
-#define TWOPOW(i) (0x1.0p##i)
-
-#define RADIX TWOPOW (24) /* 2^24 */
-#define CUTTER TWOPOW (76) /* 2^76 */
-#define RADIXI 0x1.0p-24 /* 2^-24 */
-#define TWO52 TWOPOW (52) /* 2^52 */
-
-/* Divide D by RADIX and put the remainder in R. */
-#define DIV_RADIX(d,r) \
- ({ \
- double u = ((d) + CUTTER) - CUTTER; \
- if (u > (d)) \
- u -= RADIX; \
- r = (d) - u; \
- (d) = u * RADIXI; \
- })
-
-/* Put the integer component of a double X in R and retain the fraction in
- X. */
-#define INTEGER_OF(x, r) \
- ({ \
- double u = ((x) + TWO52) - TWO52; \
- if (u > (x)) \
- u -= ONE; \
- (r) = u; \
- (x) -= u; \
- })
-
-/* Align IN down to a multiple of F, where F is a power of two. */
-#define ALIGN_DOWN_TO(in, f) \
- ({ \
- double factor = f * TWO52; \
- double u = (in + factor) - factor; \
- if (u > in) \
- u -= f; \
- u; \
- })
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-07 20:55:22 +0000