From patchwork Mon Apr 27 10:34:22 2015 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Wilco Dijkstra X-Patchwork-Id: 6458 Received: (qmail 72240 invoked by alias); 27 Apr 2015 10:34:48 -0000 Mailing-List: contact libc-alpha-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: libc-alpha-owner@sourceware.org Delivered-To: mailing list libc-alpha@sourceware.org Received: (qmail 72228 invoked by uid 89); 27 Apr 2015 10:34:46 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-0.6 required=5.0 tests=AWL, BAYES_50, SPF_PASS autolearn=ham version=3.3.2 X-HELO: eu-smtp-delivery-143.mimecast.com From: "Wilco Dijkstra" To: Subject: [PATCH] Replace ABS macros with fabs Date: Mon, 27 Apr 2015 11:34:22 +0100 Message-ID: <000301d080d5$ba014860$2e03d920$@com> MIME-Version: 1.0 X-MC-Unique: ZqBpWpygTqKAzjs38xOhtA-1 GLIBC contains many uses of ABS macros which GCC is not able to optimize into fabs due to strict IEEE mode being used. On AArch64 GCC fails to use CSEL so several math functions contain many unnecessary and unpredictable branches. This patch removes the various ABS macros and replaces uses with fabs (or in one case abs) which is more efficient on all targets. OK for commit? 2015-04-27 Wilco Dijkstra * stdio-common/printf_fp.c (___printf_fp): Use abs. * stdlib/gmp-impl.h (ABS): Remove define. (ABSIZ): Remove. * sysdeps/ieee754/dbl-64/branred.c (__branred): Use fabs. * sysdeps/ieee754/dbl-64/dla.h (EADD): Use fabs. (ESUB): Use fabs. (ADD2): Use fabs. (SUB2): Use fabs. (ADD2A): Use fabs. (SUB2A): Use fabs. * sysdeps/ieee754/dbl-64/e_asin.c (__ieee754_asin): Use fabs. * sysdeps/ieee754/dbl-64/e_log.c (__ieee754_log): Use fabs. * sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Use fabs. (log1): Use fabs. (my_log2): Use fabs. * sysdeps/ieee754/dbl-64/e_remainder.c (__ieee754_remainder): Use fabs. * sysdeps/ieee754/dbl-64/mpa.h (ABS): Remove define. * sysdeps/ieee754/dbl-64/mpatan.c (__mpatan): Use fabs. * sysdeps/ieee754/dbl-64/mydefs.h (ABS): Remove define. * sysdeps/ieee754/dbl-64/s_sin.c (__sin): Use fabs. (__cos): Use fabs. (slow): Use fabs. (slow2): Use fabs. (sloww): Use fabs. (sloww1): Use fabs. (sloww2): Use fabs. (bslow1): Use fabs. (bslow2): Use fabs. (cslow2): Use fabs. (csloww): Use fabs. (csloww1): Use fabs. (csloww2): Use fabs. * sysdeps/ieee754/dbl-64/sincos32.c (__mpranred): Use fabs. --- stdio-common/printf_fp.c | 2 +- stdlib/gmp-impl.h | 2 - sysdeps/ieee754/dbl-64/branred.c | 3 +- sysdeps/ieee754/dbl-64/dla.h | 26 ++++++------ sysdeps/ieee754/dbl-64/e_asin.c | 29 ++++++------- sysdeps/ieee754/dbl-64/e_log.c | 3 +- sysdeps/ieee754/dbl-64/e_pow.c | 12 +++--- sysdeps/ieee754/dbl-64/e_remainder.c | 15 +++---- sysdeps/ieee754/dbl-64/mpa.h | 2 - sysdeps/ieee754/dbl-64/mpatan.c | 3 +- sysdeps/ieee754/dbl-64/mydefs.h | 1 - sysdeps/ieee754/dbl-64/s_sin.c | 79 ++++++++++++++++++------------------ sysdeps/ieee754/dbl-64/sincos32.c | 3 +- 13 files changed, 92 insertions(+), 88 deletions(-) diff --git a/stdio-common/printf_fp.c b/stdio-common/printf_fp.c index f9ea379..575842b 100644 --- a/stdio-common/printf_fp.c +++ b/stdio-common/printf_fp.c @@ -449,7 +449,7 @@ ___printf_fp (FILE *fp, efficient to use variables of the fixed maximum size but because this would be really big it could lead to memory problems. */ { - mp_size_t bignum_size = ((ABS (p.exponent) + BITS_PER_MP_LIMB - 1) + mp_size_t bignum_size = ((abs (p.exponent) + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB + (LDBL_MANT_DIG / BITS_PER_MP_LIMB > 2 ? 8 : 4)) * sizeof (mp_limb_t); diff --git a/stdlib/gmp-impl.h b/stdlib/gmp-impl.h index 314a4b3..2cbb21b 100644 --- a/stdlib/gmp-impl.h +++ b/stdlib/gmp-impl.h @@ -64,7 +64,6 @@ along with the GNU MP Library; see the file COPYING.LIB. If not, see #define inline /* Empty */ #endif -#define ABS(x) (x >= 0 ? x : -x) #ifndef MIN #define MIN(l,o) ((l) < (o) ? (l) : (o)) #endif @@ -74,7 +73,6 @@ along with the GNU MP Library; see the file COPYING.LIB. If not, see /* Field access macros. */ #define SIZ(x) ((x)->_mp_size) -#define ABSIZ(x) ABS (SIZ (x)) #define PTR(x) ((x)->_mp_d) #define EXP(x) ((x)->_mp_exp) #define PREC(x) ((x)->_mp_prec) diff --git a/sysdeps/ieee754/dbl-64/branred.c b/sysdeps/ieee754/dbl-64/branred.c index 331cde2..7757a9c 100644 --- a/sysdeps/ieee754/dbl-64/branred.c +++ b/sysdeps/ieee754/dbl-64/branred.c @@ -35,6 +35,7 @@ #include "endian.h" #include "mydefs.h" #include "branred.h" +#include #include #ifndef SECTION @@ -123,7 +124,7 @@ __branred(double x, double *a, double *aa) sum=sum1+sum2; b=b1+b2; - bb = (ABS(b1)>ABS(b2))? (b1-b)+b2 : (b2-b)+b1; + bb = (fabs(b1)>fabs(b2))? (b1-b)+b2 : (b2-b)+b1; if (b > 0.5) {b-=1.0; sum+=1.0;} else if (b < -0.5) diff --git a/sysdeps/ieee754/dbl-64/dla.h b/sysdeps/ieee754/dbl-64/dla.h index a4d0b8b..641644e 100644 --- a/sysdeps/ieee754/dbl-64/dla.h +++ b/sysdeps/ieee754/dbl-64/dla.h @@ -17,6 +17,8 @@ * along with this program; if not, see . */ +#include + /***********************************************************************/ /*MODULE_NAME: dla.h */ /* */ @@ -44,7 +46,7 @@ /* z+zz = x+y exactly. */ #define EADD(x,y,z,zz) \ - z=(x)+(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x)); + z=(x)+(y); zz=(fabs(x)>fabs(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x)); /* Exact subtraction of two single-length floating point numbers, Dekker. */ @@ -52,7 +54,7 @@ /* z+zz = x-y exactly. */ #define ESUB(x,y,z,zz) \ - z=(x)-(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z))); + z=(x)-(y); zz=(fabs(x)>fabs(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z))); /* Exact multiplication of two single-length floating point numbers, */ @@ -94,7 +96,7 @@ /* storage variables of type double. */ #define ADD2(x, xx, y, yy, z, zz, r, s) \ - r = (x) + (y); s = (ABS (x) > ABS (y)) ? \ + r = (x) + (y); s = (fabs (x) > fabs (y)) ? \ (((((x) - r) + (y)) + (yy)) + (xx)) : \ (((((y) - r) + (x)) + (xx)) + (yy)); \ z = r + s; zz = (r - z) + s; @@ -107,7 +109,7 @@ /* storage variables of type double. */ #define SUB2(x, xx, y, yy, z, zz, r, s) \ - r = (x) - (y); s = (ABS (x) > ABS (y)) ? \ + r = (x) - (y); s = (fabs (x) > fabs (y)) ? \ (((((x) - r) - (y)) - (yy)) + (xx)) : \ ((((x) - ((y) + r)) + (xx)) - (yy)); \ z = r + s; zz = (r - z) + s; @@ -144,16 +146,16 @@ #define ADD2A(x, xx, y, yy, z, zz, r, rr, s, ss, u, uu, w) \ r = (x) + (y); \ - if (ABS (x) > ABS (y)) { rr = ((x) - r) + (y); s = (rr + (yy)) + (xx); } \ + if (fabs (x) > fabs (y)) { rr = ((x) - r) + (y); s = (rr + (yy)) + (xx); } \ else { rr = ((y) - r) + (x); s = (rr + (xx)) + (yy); } \ if (rr != 0.0) { \ z = r + s; zz = (r - z) + s; } \ else { \ - ss = (ABS (xx) > ABS (yy)) ? (((xx) - s) + (yy)) : (((yy) - s) + (xx));\ + ss = (fabs (xx) > fabs (yy)) ? (((xx) - s) + (yy)) : (((yy) - s) + (xx));\ u = r + s; \ - uu = (ABS (r) > ABS (s)) ? ((r - u) + s) : ((s - u) + r); \ + uu = (fabs (r) > fabs (s)) ? ((r - u) + s) : ((s - u) + r); \ w = uu + ss; z = u + w; \ - zz = (ABS (u) > ABS (w)) ? ((u - z) + w) : ((w - z) + u); } + zz = (fabs (u) > fabs (w)) ? ((u - z) + w) : ((w - z) + u); } /* Double-length subtraction, slower but more accurate than SUB2. */ @@ -165,13 +167,13 @@ #define SUB2A(x, xx, y, yy, z, zz, r, rr, s, ss, u, uu, w) \ r = (x) - (y); \ - if (ABS (x) > ABS (y)) { rr = ((x) - r) - (y); s = (rr - (yy)) + (xx); } \ + if (fabs (x) > fabs (y)) { rr = ((x) - r) - (y); s = (rr - (yy)) + (xx); } \ else { rr = (x) - ((y) + r); s = (rr + (xx)) - (yy); } \ if (rr != 0.0) { \ z = r + s; zz = (r - z) + s; } \ else { \ - ss = (ABS (xx) > ABS (yy)) ? (((xx) - s) - (yy)) : ((xx) - ((yy) + s)); \ + ss = (fabs (xx) > fabs (yy)) ? (((xx) - s) - (yy)) : ((xx) - ((yy) + s)); \ u = r + s; \ - uu = (ABS (r) > ABS (s)) ? ((r - u) + s) : ((s - u) + r); \ + uu = (fabs (r) > fabs (s)) ? ((r - u) + s) : ((s - u) + r); \ w = uu + ss; z = u + w; \ - zz = (ABS (u) > ABS (w)) ? ((u - z) + w) : ((w - z) + u); } + zz = (fabs (u) > fabs (w)) ? ((u - z) + w) : ((w - z) + u); } diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c index e90f47c..8ed7054 100644 --- a/sysdeps/ieee754/dbl-64/e_asin.c +++ b/sysdeps/ieee754/dbl-64/e_asin.c @@ -39,6 +39,7 @@ #include "powtwo.tbl" #include "MathLib.h" #include "uasncs.h" +#include #include #ifndef SECTION @@ -94,9 +95,9 @@ __ieee754_asin(double x){ __doasin(x,0,w); if (w[0]==(w[0]+1.00000001*w[1])) return w[0]; else { - y=ABS(x); - res=ABS(w[0]); - res1=ABS(w[0]+1.1*w[1]); + y=fabs(x); + res=fabs(w[0]); + res1=fabs(w[0]+1.1*w[1]); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -125,11 +126,11 @@ __ieee754_asin(double x){ res1=res+1.1*cor; z=0.5*(res1-res); __dubsin(res,z,w); - z=(w[0]-ABS(x))+w[1]; + z=(w[0]-fabs(x))+w[1]; if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); else { - y=ABS(x); + y=fabs(x); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -158,11 +159,11 @@ __ieee754_asin(double x){ res1=res+1.1*cor; z=0.5*(res1-res); __dubsin(res,z,w); - z=(w[0]-ABS(x))+w[1]; + z=(w[0]-fabs(x))+w[1]; if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); else { - y=ABS(x); + y=fabs(x); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -193,11 +194,11 @@ __ieee754_asin(double x){ y=hp0.x-res; z=((hp0.x-y)-res)+(hp1.x-z); __dubcos(y,z,w); - z=(w[0]-ABS(x))+w[1]; + z=(w[0]-fabs(x))+w[1]; if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); else { - y=ABS(x); + y=fabs(x); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -231,11 +232,11 @@ __ieee754_asin(double x){ z=y+hp1.x; y=(y-z)+hp1.x; __dubcos(z,y,w); - z=(w[0]-ABS(x))+w[1]; + z=(w[0]-fabs(x))+w[1]; if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); else { - y=ABS(x); + y=fabs(x); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -270,11 +271,11 @@ __ieee754_asin(double x){ z=y+hp1.x; y=(y-z)+hp1.x; __dubcos(z,y,w); - z=(w[0]-ABS(x))+w[1]; + z=(w[0]-fabs(x))+w[1]; if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); else { - y=ABS(x); + y=fabs(x); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -308,7 +309,7 @@ __ieee754_asin(double x){ cor = (res1-res)+cor; if (res==(res+1.0000001*cor)) return (m>0)?res:-res; else { - y=ABS(x); + y=fabs(x); res1=res+1.1*cor; return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c index 7a498e8..a9117fb 100644 --- a/sysdeps/ieee754/dbl-64/e_log.c +++ b/sysdeps/ieee754/dbl-64/e_log.c @@ -38,6 +38,7 @@ #include #include "mpa.h" #include "MathLib.h" +#include #include #include @@ -93,7 +94,7 @@ __ieee754_log (double x) /* Regular values of x */ w = x - 1; - if (__glibc_likely (ABS (w) > U03)) + if (__glibc_likely (fabs (w) > U03)) goto case_03; /* log (1) is +0 in all rounding modes. */ diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index ba932f4..8a1f72f 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.c @@ -97,7 +97,7 @@ __ieee754_pow (double x, double y) /* Avoid internal underflow for tiny y. The exact value of y does not matter if |y| <= 2**-64. */ - if (ABS (y) < 0x1p-64) + if (fabs (y) < 0x1p-64) y = y < 0 ? -0x1p-64 : 0x1p-64; z = log1 (x, &aa, &error); /* x^y =e^(y log (X)) */ t = y * CN; @@ -110,7 +110,7 @@ __ieee754_pow (double x, double y) aa = y2 * a1 + y * a2; a1 = a + aa; a2 = (a - a1) + aa; - error = error * ABS (y); + error = error * fabs (y); t = __exp1 (a1, a2, 1.9e16 * error); /* return -10 or 0 if wasn't computed exactly */ retval = (t > 0) ? t : power1 (x, y); } @@ -127,7 +127,7 @@ __ieee754_pow (double x, double y) if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0) || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000) /* NaN */ return y; - if (ABS (y) > 1.0e20) + if (fabs (y) > 1.0e20) return (y > 0) ? 0 : 1.0 / 0.0; k = checkint (y); if (k == -1) @@ -232,7 +232,7 @@ power1 (double x, double y) aa = ((y1 * a1 - a) + y1 * a2 + y2 * a1) + y2 * a2 + aa * y; a1 = a + aa; a2 = (a - a1) + aa; - error = error * ABS (y); + error = error * fabs (y); t = __exp1 (a1, a2, 1.9e16 * error); return (t >= 0) ? t : __slowpow (x, y, z); } @@ -292,7 +292,7 @@ log1 (double x, double *delta, double *error) * (r7 + t * r8))))) - 0.5 * t2 * (t + t1)); res = e1 + e2; - *error = 1.0e-21 * ABS (t); + *error = 1.0e-21 * fabs (t); *delta = (e1 - res) + e2; return res; } /* |x-1| < 1.5*2**-10 */ @@ -398,7 +398,7 @@ my_log2 (double x, double *delta, double *error) e2 = ((((t - e1) + z) + zz) + t * t * t * (ss3 + t * (s4 + t * (s5 + t * (s6 + t * (s7 + t * s8)))))); res = e1 + e2; - *error = 1.0e-25 * ABS (t); + *error = 1.0e-25 * fabs (t); *delta = (e1 - res) + e2; return res; } diff --git a/sysdeps/ieee754/dbl-64/e_remainder.c b/sysdeps/ieee754/dbl-64/e_remainder.c index ba479b4..7f3a02b 100644 --- a/sysdeps/ieee754/dbl-64/e_remainder.c +++ b/sysdeps/ieee754/dbl-64/e_remainder.c @@ -33,6 +33,7 @@ #include "mydefs.h" #include "urem.h" #include "MathLib.h" +#include #include /**************************************************************************/ @@ -66,7 +67,7 @@ __ieee754_remainder (double x, double y) return (xx != 0) ? xx : ((x > 0) ? ZERO.x : nZERO.x); else { - if (ABS (xx) > 0.5 * t.x) + if (fabs (xx) > 0.5 * t.x) return (z > d) ? xx - t.x : xx + t.x; else return xx; @@ -98,10 +99,10 @@ __ieee754_remainder (double x, double y) z = u.x * r.x; d = (z + big.x) - big.x; u.x = (u.x - d * w.x) - d * ww.x; - if (ABS (u.x) < 0.5 * t.x) + if (fabs (u.x) < 0.5 * t.x) return (u.x != 0) ? u.x : ((x > 0) ? ZERO.x : nZERO.x); else - if (ABS (u.x) > 0.5 * t.x) + if (fabs (u.x) > 0.5 * t.x) return (d > z) ? u.x + t.x : u.x - t.x; else { @@ -114,7 +115,7 @@ __ieee754_remainder (double x, double y) { if (kx < 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0)) { - y = ABS (y) * t128.x; + y = fabs (y) * t128.x; z = __ieee754_remainder (x, y) * t128.x; z = __ieee754_remainder (z, y) * tm128.x; return z; @@ -124,10 +125,10 @@ __ieee754_remainder (double x, double y) if ((kx & 0x7ff00000) == 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0)) { - y = ABS (y); + y = fabs (y); z = 2.0 * __ieee754_remainder (0.5 * x, y); - d = ABS (z); - if (d <= ABS (d - y)) + d = fabs (z); + if (d <= fabs (d - y)) return z; else return (z > 0) ? z - y : z + y; diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h index 9e87ea6..47dd6c4 100644 --- a/sysdeps/ieee754/dbl-64/mpa.h +++ b/sysdeps/ieee754/dbl-64/mpa.h @@ -81,8 +81,6 @@ extern const mp_no __mptwo; #define EY y->e #define EZ z->e -#define ABS(x) ((x) < 0 ? -(x) : (x)) - #ifndef RADIXI # define RADIXI 0x1.0p-24 /* 2^-24 */ #endif diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c index 004ed7a..d9ac3d6 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.c +++ b/sysdeps/ieee754/dbl-64/mpatan.c @@ -32,6 +32,7 @@ #include "endian.h" #include "mpa.h" +#include #ifndef SECTION # define SECTION @@ -65,7 +66,7 @@ __mpatan (mp_no *x, mp_no *y, int p) else { __mp_dbl (x, &dx, p); - dx = ABS (dx); + dx = fabs (dx); for (m = 6; m > 0; m--) { if (dx > __atan_xm[m].d) diff --git a/sysdeps/ieee754/dbl-64/mydefs.h b/sysdeps/ieee754/dbl-64/mydefs.h index f24a89c..3158514 100644 --- a/sysdeps/ieee754/dbl-64/mydefs.h +++ b/sysdeps/ieee754/dbl-64/mydefs.h @@ -30,7 +30,6 @@ typedef int int4; typedef union { int4 i[2]; double x; } mynumber; -#define ABS(x) (((x) > 0) ? (x) : -(x)) #define max(x, y) (((y) > (x)) ? (y) : (x)) #define min(x, y) (((y) < (x)) ? (y) : (x)) #endif diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index 4c38fa0..ea89ad5 100644 --- a/sysdeps/ieee754/dbl-64/s_sin.c +++ b/sysdeps/ieee754/dbl-64/s_sin.c @@ -52,6 +52,7 @@ #include "mydefs.h" #include "usncs.h" #include "MathLib.h" +#include #include #include @@ -355,7 +356,7 @@ __sin (double x) da = xn * mp3; a = y - da; da = (y - a) - da; - eps = ABS (x) * 1.2e-30; + eps = fabs (x) * 1.2e-30; switch (n) { /* quarter of unit circle */ @@ -530,7 +531,7 @@ __cos (double x) else if (k < 0x3feb6000) { /* 2^-27 < |x| < 0.855469 */ - y = ABS (x); + y = fabs (x); u.x = big + y; y = y - (u.x - big); res = do_cos (u, y, &cor); @@ -539,7 +540,7 @@ __cos (double x) else if (k < 0x400368fd) { /* 0.855469 <|x|<2.426265 */ ; - y = hp0 - ABS (x); + y = hp0 - fabs (x); a = y + hp1; da = (y - a) + hp1; xx = a * a; @@ -582,7 +583,7 @@ __cos (double x) da = xn * mp3; a = y - da; da = (y - a) - da; - eps = ABS (x) * 1.2e-30; + eps = fabs (x) * 1.2e-30; switch (n) { @@ -741,7 +742,7 @@ slow (double x) return res; else { - __dubsin (ABS (x), 0, w); + __dubsin (fabs (x), 0, w); if (w[0] == w[0] + 1.000000001 * w[1]) return (x > 0) ? w[0] : -w[0]; else @@ -760,7 +761,7 @@ slow1 (double x) { mynumber u; double w[2], y, cor, res; - y = ABS (x); + y = fabs (x); u.x = big + y; y = y - (u.x - big); res = do_sin_slow (u, y, 0, 0, &cor); @@ -768,7 +769,7 @@ slow1 (double x) return (x > 0) ? res : -res; else { - __dubsin (ABS (x), 0, w); + __dubsin (fabs (x), 0, w); if (w[0] == w[0] + 1.000000005 * w[1]) return (x > 0) ? w[0] : -w[0]; else @@ -787,7 +788,7 @@ slow2 (double x) mynumber u; double w[2], y, y1, y2, cor, res, del; - y = ABS (x); + y = fabs (x); y = hp0 - y; if (y >= 0) { @@ -806,7 +807,7 @@ slow2 (double x) return (x > 0) ? res : -res; else { - y = ABS (x) - hp0; + y = fabs (x) - hp0; y1 = y - hp1; y2 = (y - y1) - hp1; __docos (y1, y2, w); @@ -834,9 +835,9 @@ sloww (double x, double dx, double orig) int4 n; res = TAYLOR_SLOW (x, dx, cor); if (cor > 0) - cor = 1.0005 * cor + ABS (orig) * 3.1e-30; + cor = 1.0005 * cor + fabs (orig) * 3.1e-30; else - cor = 1.0005 * cor - ABS (orig) * 3.1e-30; + cor = 1.0005 * cor - fabs (orig) * 3.1e-30; if (res == res + cor) return res; @@ -844,9 +845,9 @@ sloww (double x, double dx, double orig) { (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); if (w[1] > 0) - cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-30; + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; else - cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-30; + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; if (w[0] == w[0] + cor) return (x > 0) ? w[0] : -w[0]; @@ -870,9 +871,9 @@ sloww (double x, double dx, double orig) } (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); if (w[1] > 0) - cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-40; + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; else - cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-40; + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; if (w[0] == w[0] + cor) return (a > 0) ? w[0] : -w[0]; @@ -898,7 +899,7 @@ sloww1 (double x, double dx, double orig, int m) u.x = big + x; y = x - (u.x - big); - res = do_sin_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); + res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); if (res == res + cor) return (m > 0) ? res : -res; @@ -907,9 +908,9 @@ sloww1 (double x, double dx, double orig, int m) __dubsin (x, dx, w); if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); if (w[0] == w[0] + cor) return (m > 0) ? w[0] : -w[0]; @@ -934,7 +935,7 @@ sloww2 (double x, double dx, double orig, int n) u.x = big + x; y = x - (u.x - big); - res = do_cos_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); + res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); if (res == res + cor) return (n & 2) ? -res : res; @@ -943,9 +944,9 @@ sloww2 (double x, double dx, double orig, int n) __docos (x, dx, w); if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); if (w[0] == w[0] + cor) return (n & 2) ? -w[0] : w[0]; @@ -1000,7 +1001,7 @@ bsloww1 (double x, double dx, double orig, int n) mynumber u; double w[2], y, cor, res; - y = ABS (x); + y = fabs (x); u.x = big + y; y = y - (u.x - big); dx = (x > 0) ? dx : -dx; @@ -1009,7 +1010,7 @@ bsloww1 (double x, double dx, double orig, int n) return (x > 0) ? res : -res; else { - __dubsin (ABS (x), dx, w); + __dubsin (fabs (x), dx, w); if (w[1] > 0) cor = 1.000000005 * w[1] + 1.1e-24; @@ -1037,7 +1038,7 @@ bsloww2 (double x, double dx, double orig, int n) mynumber u; double w[2], y, cor, res; - y = ABS (x); + y = fabs (x); u.x = big + y; y = y - (u.x - big); dx = (x > 0) ? dx : -dx; @@ -1046,7 +1047,7 @@ bsloww2 (double x, double dx, double orig, int n) return (n & 2) ? -res : res; else { - __docos (ABS (x), dx, w); + __docos (fabs (x), dx, w); if (w[1] > 0) cor = 1.000000005 * w[1] + 1.1e-24; @@ -1072,7 +1073,7 @@ cslow2 (double x) mynumber u; double w[2], y, cor, res; - y = ABS (x); + y = fabs (x); u.x = big + y; y = y - (u.x - big); res = do_cos_slow (u, y, 0, 0, &cor); @@ -1080,7 +1081,7 @@ cslow2 (double x) return res; else { - y = ABS (x); + y = fabs (x); __docos (y, 0, w); if (w[0] == w[0] + 1.000000005 * w[1]) return w[0]; @@ -1109,9 +1110,9 @@ csloww (double x, double dx, double orig) res = TAYLOR_SLOW (x, dx, cor); if (cor > 0) - cor = 1.0005 * cor + ABS (orig) * 3.1e-30; + cor = 1.0005 * cor + fabs (orig) * 3.1e-30; else - cor = 1.0005 * cor - ABS (orig) * 3.1e-30; + cor = 1.0005 * cor - fabs (orig) * 3.1e-30; if (res == res + cor) return res; @@ -1120,9 +1121,9 @@ csloww (double x, double dx, double orig) (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); if (w[1] > 0) - cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-30; + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; else - cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-30; + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; if (w[0] == w[0] + cor) return (x > 0) ? w[0] : -w[0]; @@ -1147,9 +1148,9 @@ csloww (double x, double dx, double orig) (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); if (w[1] > 0) - cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-40; + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; else - cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-40; + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; if (w[0] == w[0] + cor) return (a > 0) ? w[0] : -w[0]; @@ -1175,7 +1176,7 @@ csloww1 (double x, double dx, double orig, int m) u.x = big + x; y = x - (u.x - big); - res = do_sin_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); + res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); if (res == res + cor) return (m > 0) ? res : -res; @@ -1183,9 +1184,9 @@ csloww1 (double x, double dx, double orig, int m) { __dubsin (x, dx, w); if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); if (w[0] == w[0] + cor) return (m > 0) ? w[0] : -w[0]; else @@ -1210,7 +1211,7 @@ csloww2 (double x, double dx, double orig, int n) u.x = big + x; y = x - (u.x - big); - res = do_cos_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); + res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); if (res == res + cor) return (n) ? -res : res; @@ -1218,9 +1219,9 @@ csloww2 (double x, double dx, double orig, int n) { __docos (x, dx, w); if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); if (w[0] == w[0] + cor) return (n) ? -w[0] : w[0]; else diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c index 01cdd30..d942a1f 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ b/sysdeps/ieee754/dbl-64/sincos32.c @@ -42,6 +42,7 @@ #include "endian.h" #include "mpa.h" #include "sincos32.h" +#include #include #include @@ -318,7 +319,7 @@ __mpranred (double x, mp_no *y, int p) int i, k, n; mp_no a, b, c; - if (ABS (x) < 2.8e14) + if (fabs (x) < 2.8e14) { t = (x * hpinv.d + toint.d); xn = t - toint.d;