From 1e3803454e5ff517609c96166fcfaf966369920f Mon Sep 17 00:00:00 2001 From: Siddhesh Poyarekar Date: Fri, 15 Mar 2013 23:11:50 +0530 Subject: Revert configurable mantissa patch Reverts d22ca8cdfb98001d03772ef264b244930d439b3f since it is severely broken on 32-bit. --- ChangeLog | 7 ++++ sysdeps/ieee754/dbl-64/mpa-arch.h | 47 --------------------- sysdeps/ieee754/dbl-64/mpa.c | 78 ++++++++++++++++++++++------------- sysdeps/ieee754/dbl-64/mpa.h | 27 +++++------- sysdeps/powerpc/power4/fpu/mpa-arch.h | 56 ------------------------- 5 files changed, 67 insertions(+), 148 deletions(-) delete mode 100644 sysdeps/ieee754/dbl-64/mpa-arch.h delete mode 100644 sysdeps/powerpc/power4/fpu/mpa-arch.h diff --git a/ChangeLog b/ChangeLog index fe92b82479..dc84119138 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,10 @@ +2013-03-15 Siddhesh Poyarekar + + * 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 * 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 . */ - -#include - -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 /* 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 . */ - -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; \ - }) -- cgit 1.4.1