diff options
author | Ulrich Drepper <drepper@gmail.com> | 2011-10-12 11:27:51 -0400 |
---|---|---|
committer | Ulrich Drepper <drepper@gmail.com> | 2011-10-12 11:27:51 -0400 |
commit | 0ac5ae2335292908f39031b1ea9fe8edce433c0f (patch) | |
tree | f9d26c8abc0de39d18d4c13e70f6022cdc6b461f /sysdeps/ieee754/flt-32 | |
parent | a843a204a3e8a0dd53584dad3668771abaec84ac (diff) | |
download | glibc-0ac5ae2335292908f39031b1ea9fe8edce433c0f.tar.gz glibc-0ac5ae2335292908f39031b1ea9fe8edce433c0f.tar.xz glibc-0ac5ae2335292908f39031b1ea9fe8edce433c0f.zip |
Optimize libm
libm is now somewhat integrated with gcc's -ffinite-math-only option and lots of the wrapper functions have been optimized.
Diffstat (limited to 'sysdeps/ieee754/flt-32')
22 files changed, 377 insertions, 749 deletions
diff --git a/sysdeps/ieee754/flt-32/e_acosf.c b/sysdeps/ieee754/flt-32/e_acosf.c index 0d85c4210d..a258e2f918 100644 --- a/sysdeps/ieee754/flt-32/e_acosf.c +++ b/sysdeps/ieee754/flt-32/e_acosf.c @@ -8,23 +8,15 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_acosf.c,v 1.5 1995/05/12 04:57:16 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ -static const float -#else -static float -#endif +static const float one = 1.0000000000e+00, /* 0x3F800000 */ pi = 3.1415925026e+00, /* 0x40490fda */ pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */ @@ -40,12 +32,8 @@ qS2 = 2.0209457874e+00, /* 0x4001572d */ qS3 = -6.8828397989e-01, /* 0xbf303361 */ qS4 = 7.7038154006e-02; /* 0x3d9dc62e */ -#ifdef __STDC__ - float __ieee754_acosf(float x) -#else - float __ieee754_acosf(x) - float x; -#endif +float +__ieee754_acosf(float x) { float z,p,q,r,w,s,c,df; int32_t hx,ix; @@ -87,3 +75,4 @@ qS4 = 7.7038154006e-02; /* 0x3d9dc62e */ return (float)2.0*(df+w); } } +strong_alias (__ieee754_acosf, __acosf_finite) diff --git a/sysdeps/ieee754/flt-32/e_acoshf.c b/sysdeps/ieee754/flt-32/e_acoshf.c index c607f72117..db8f6ec462 100644 --- a/sysdeps/ieee754/flt-32/e_acoshf.c +++ b/sysdeps/ieee754/flt-32/e_acoshf.c @@ -8,7 +8,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -21,9 +21,9 @@ static char rcsid[] = "$NetBSD: e_acoshf.c,v 1.5 1995/05/12 04:57:20 jtc Exp $"; #include "math_private.h" #ifdef __STDC__ -static const float +static const float #else -static float +static float #endif one = 1.0, ln2 = 6.9314718246e-01; /* 0x3f317218 */ @@ -34,7 +34,7 @@ ln2 = 6.9314718246e-01; /* 0x3f317218 */ float __ieee754_acoshf(x) float x; #endif -{ +{ float t; int32_t hx; GET_FLOAT_WORD(hx,x); @@ -42,8 +42,8 @@ ln2 = 6.9314718246e-01; /* 0x3f317218 */ return (x-x)/(x-x); } else if(hx >=0x4d800000) { /* x > 2**28 */ if(hx >=0x7f800000) { /* x is inf of NaN */ - return x+x; - } else + return x+x; + } else return __ieee754_logf(x)+ln2; /* acosh(huge)=log(2x) */ } else if (hx==0x3f800000) { return 0.0; /* acosh(1) = 0 */ @@ -55,3 +55,4 @@ ln2 = 6.9314718246e-01; /* 0x3f317218 */ return __log1pf(t+__sqrtf((float)2.0*t+t*t)); } } +strong_alias (__ieee754_acoshf, __acoshf_finite) diff --git a/sysdeps/ieee754/flt-32/e_asinf.c b/sysdeps/ieee754/flt-32/e_asinf.c index b0c835c83c..7296ba3e13 100644 --- a/sysdeps/ieee754/flt-32/e_asinf.c +++ b/sysdeps/ieee754/flt-32/e_asinf.c @@ -14,11 +14,11 @@ */ /* - Modifications for single precision expansion are + Modifications for single precision expansion are Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> - and are incorporated herein by permission of the author. The author + and are incorporated herein by permission of the author. The author reserves the right to distribute this material elsewhere under different - copying permissions. These modifications are distributed here under + copying permissions. These modifications are distributed here under the following terms: This library is free software; you can redistribute it and/or @@ -93,7 +93,7 @@ p4 = 4.216630880e-2f; t = w*0.5f; p = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)))); s = __ieee754_sqrtf(t); - if(ix>=0x3F79999A) { /* if |x| > 0.975 */ + if(ix>=0x3F79999A) { /* if |x| > 0.975 */ t = pio2_hi-(2.0f*(s+s*p)-pio2_lo); } else { int32_t iw; @@ -108,3 +108,4 @@ p4 = 4.216630880e-2f; } if(hx>0) return t; else return -t; } +strong_alias (__ieee754_asinf, __asinf_finite) diff --git a/sysdeps/ieee754/flt-32/e_atan2f.c b/sysdeps/ieee754/flt-32/e_atan2f.c index c0cafb16b8..abbde88bdf 100644 --- a/sysdeps/ieee754/flt-32/e_atan2f.c +++ b/sysdeps/ieee754/flt-32/e_atan2f.c @@ -13,18 +13,10 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_atan2f.c,v 1.4 1995/05/10 20:44:53 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float -#else -static float -#endif tiny = 1.0e-30, zero = 0.0, pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */ @@ -32,12 +24,8 @@ pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */ pi = 3.1415927410e+00, /* 0x40490fdb */ pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */ -#ifdef __STDC__ - float __ieee754_atan2f(float y, float x) -#else - float __ieee754_atan2f(y,x) - float y,x; -#endif +float +__ieee754_atan2f (float y, float x) { float z; int32_t k,m,hx,hy,ix,iy; @@ -56,7 +44,7 @@ pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */ if(iy==0) { switch(m) { case 0: - case 1: return y; /* atan(+-0,+anything)=+-0 */ + case 1: return y; /* atan(+-0,+anything)=+-0 */ case 2: return pi+tiny;/* atan(+0,-anything) = pi */ case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ } @@ -87,19 +75,20 @@ pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */ /* compute y/x */ k = (iy-ix)>>23; - if(k > 60) z=pi_o_2+(float)0.5*pi_lo; /* |y/x| > 2**60 */ - else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */ + if(k > 60) z=pi_o_2+(float)0.5*pi_lo; /* |y/x| > 2**60 */ + else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */ else z=__atanf(fabsf(y/x)); /* safe to do y/x */ switch (m) { case 0: return z ; /* atan(+,+) */ case 1: { - u_int32_t zh; + u_int32_t zh; GET_FLOAT_WORD(zh,z); SET_FLOAT_WORD(z,zh ^ 0x80000000); } return z ; /* atan(-,+) */ case 2: return pi-(z-pi_lo);/* atan(+,-) */ default: /* case 3 */ - return (z-pi_lo)-pi;/* atan(-,-) */ + return (z-pi_lo)-pi;/* atan(-,-) */ } } +strong_alias (__ieee754_atan2f, __atan2f_finite) diff --git a/sysdeps/ieee754/flt-32/e_atanhf.c b/sysdeps/ieee754/flt-32/e_atanhf.c index f26a15bbc1..ddd18ab300 100644 --- a/sysdeps/ieee754/flt-32/e_atanhf.c +++ b/sysdeps/ieee754/flt-32/e_atanhf.c @@ -1,58 +1,70 @@ -/* e_atanhf.c -- float version of e_atanh.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ +/* Copyright (C) 2011 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ + The GNU C Library 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. + + The GNU C Library 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 the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + + +/* __ieee754_atanh(x) + Method : + 1.Reduced x to positive by atanh(-x) = -atanh(x) + 2.For x>=0.5 + 1 2x x + atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) + 2 1 - x 1 - x + + For x<0.5 + atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_atanhf.c,v 1.4 1995/05/10 20:44:56 jtc Exp $"; -#endif + Special cases: + atanh(x) is NaN if |x| > 1 with signal; + atanh(NaN) is that NaN with no signal; + atanh(+-1) is +-INF with signal. + */ + +#include <inttypes.h> #include "math.h" #include "math_private.h" -#ifdef __STDC__ -static const float one = 1.0, huge = 1e30; -#else -static float one = 1.0, huge = 1e30; -#endif - -#ifdef __STDC__ -static const float zero = 0.0; -#else -static float zero = 0.0; -#endif - -#ifdef __STDC__ - float __ieee754_atanhf(float x) -#else - float __ieee754_atanhf(x) - float x; -#endif +static const float huge = 1e30; + +float +__ieee754_atanhf (float x) { - float t; - int32_t hx,ix; - GET_FLOAT_WORD(hx,x); - ix = hx&0x7fffffff; - if (ix>0x3f800000) /* |x|>1 */ - return (x-x)/(x-x); - if(ix==0x3f800000) - return x/zero; - if(ix<0x31800000&&(huge+x)>zero) return x; /* x<2**-28 */ - SET_FLOAT_WORD(x,ix); - if(ix<0x3f000000) { /* x < 0.5 */ - t = x+x; - t = (float)0.5*__log1pf(t+t*x/(one-x)); - } else - t = (float)0.5*__log1pf((x+x)/(one-x)); - if(hx>=0) return t; else return -t; + float xa = fabsf (x); + float t; + if (xa < 0.5f) + { + if (__builtin_expect (xa < 0x1.0p-28f, 0) && (huge + x) > 0.0f) + return x; + + t = xa + xa; + t = 0.5f * __log1pf (t + t * xa / (1.0f - xa)); + } + else if (__builtin_expect (xa < 1.0f, 1)) + t = 0.5f * __log1pf ((xa + xa) / (1.0f - xa)); + else + { + if (xa > 1.0f) + return (x - x) / (x - x); + + return x / 0.0f; + } + + return __copysignf (t, x); } +strong_alias (__ieee754_atanhf, __atanhf_finite) diff --git a/sysdeps/ieee754/flt-32/e_coshf.c b/sysdeps/ieee754/flt-32/e_coshf.c index 223fbeea20..1887639a6b 100644 --- a/sysdeps/ieee754/flt-32/e_coshf.c +++ b/sysdeps/ieee754/flt-32/e_coshf.c @@ -1,5 +1,6 @@ /* e_coshf.c -- float version of e_cosh.c. * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + * Optimizations by Ulrich Drepper <drepper@gmail.com>, 2011 */ /* @@ -13,26 +14,14 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_coshf.c,v 1.6 1996/04/08 15:43:41 phil Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float huge = 1.0e30; static const float one = 1.0, half=0.5; -#else -static float one = 1.0, half=0.5, huge = 1.0e30; -#endif -#ifdef __STDC__ - float __ieee754_coshf(float x) -#else - float __ieee754_coshf(x) - float x; -#endif +float +__ieee754_coshf (float x) { float t,w; int32_t ix; @@ -40,19 +29,17 @@ static float one = 1.0, half=0.5, huge = 1.0e30; GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; - /* x is INF or NaN */ - if(ix>=0x7f800000) return x*x; - - /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ - if(ix<0x3eb17218) { - t = __expm1f(fabsf(x)); - w = one+t; - if (ix<0x24000000) return w; /* cosh(tiny) = 1 */ - return one+(t*t)/(w+w); - } - - /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ + /* |x| in [0,22] */ if (ix < 0x41b00000) { + /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ + if(ix<0x3eb17218) { + t = __expm1f(fabsf(x)); + w = one+t; + if (ix<0x24000000) return w; /* cosh(tiny) = 1 */ + return one+(t*t)/(w+w); + } + + /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ t = __ieee754_expf(fabsf(x)); return half*t+half/t; } @@ -67,6 +54,10 @@ static float one = 1.0, half=0.5, huge = 1.0e30; return t*w; } + /* x is INF or NaN */ + if(ix>=0x7f800000) return x*x; + /* |x| > overflowthresold, cosh(x) overflow */ return huge*huge; } +strong_alias (__ieee754_coshf, __coshf_finite) diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c index 194222a0cd..0703cea403 100644 --- a/sysdeps/ieee754/flt-32/e_exp2f.c +++ b/sysdeps/ieee754/flt-32/e_exp2f.c @@ -1,5 +1,6 @@ /* Single-precision floating point 2^x. - Copyright (C) 1997,1998,2000,2001,2005,2006 Free Software Foundation, Inc. + Copyright (C) 1997,1998,2000,2001,2005,2006,2011 + Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Geoffrey Keating <geoffk@ozemail.com.au> @@ -126,3 +127,4 @@ __ieee754_exp2f (float x) /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ return TWO127*x; } +strong_alias (__ieee754_exp2f, __exp2f_finite) diff --git a/sysdeps/ieee754/flt-32/e_fmodf.c b/sysdeps/ieee754/flt-32/e_fmodf.c index 47b312392c..e82a9ceab7 100644 --- a/sysdeps/ieee754/flt-32/e_fmodf.c +++ b/sysdeps/ieee754/flt-32/e_fmodf.c @@ -8,16 +8,12 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_fmodf.c,v 1.4 1995/05/10 20:45:10 jtc Exp $"; -#endif - -/* +/* * __ieee754_fmodf(x,y) * Return x mod y in exact arithmetic * Method: shift and subtract @@ -26,18 +22,10 @@ static char rcsid[] = "$NetBSD: e_fmodf.c,v 1.4 1995/05/10 20:45:10 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float one = 1.0, Zero[] = {0.0, -0.0,}; -#else -static float one = 1.0, Zero[] = {0.0, -0.0,}; -#endif -#ifdef __STDC__ - float __ieee754_fmodf(float x, float y) -#else - float __ieee754_fmodf(x,y) - float x,y ; -#endif +float +__ieee754_fmodf (float x, float y) { int32_t n,hx,hy,hz,ix,iy,sx,i; @@ -66,13 +54,13 @@ static float one = 1.0, Zero[] = {0.0, -0.0,}; } else iy = (hy>>23)-127; /* set up {hx,lx}, {hy,ly} and align y to x */ - if(ix >= -126) + if(ix >= -126) hx = 0x00800000|(0x007fffff&hx); else { /* subnormal x, shift x to normal */ n = -126-ix; hx = hx<<n; } - if(iy >= -126) + if(iy >= -126) hy = 0x00800000|(0x007fffff&hy); else { /* subnormal y, shift y to normal */ n = -126-iy; @@ -85,17 +73,17 @@ static float one = 1.0, Zero[] = {0.0, -0.0,}; hz=hx-hy; if(hz<0){hx = hx+hx;} else { - if(hz==0) /* return sign(x)*0 */ + if(hz==0) /* return sign(x)*0 */ return Zero[(u_int32_t)sx>>31]; - hx = hz+hz; + hx = hz+hz; } } hz=hx-hy; if(hz>=0) {hx=hz;} /* convert back to floating value and restore the sign */ - if(hx==0) /* return sign(x)*0 */ - return Zero[(u_int32_t)sx>>31]; + if(hx==0) /* return sign(x)*0 */ + return Zero[(u_int32_t)sx>>31]; while(hx<0x00800000) { /* normalize x */ hx = hx+hx; iy -= 1; @@ -111,3 +99,4 @@ static float one = 1.0, Zero[] = {0.0, -0.0,}; } return x; /* exact output */ } +strong_alias (__ieee754_fmodf, __fmodf_finite) diff --git a/sysdeps/ieee754/flt-32/e_gammaf_r.c b/sysdeps/ieee754/flt-32/e_gammaf_r.c index 926c84f368..aeeddf1a96 100644 --- a/sysdeps/ieee754/flt-32/e_gammaf_r.c +++ b/sysdeps/ieee754/flt-32/e_gammaf_r.c @@ -1,5 +1,5 @@ /* Implementation of gamma function according to ISO C. - Copyright (C) 1997, 1999, 2001, 2004 Free Software Foundation, Inc. + Copyright (C) 1997, 1999, 2001, 2004, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -32,19 +32,20 @@ __ieee754_gammaf_r (float x, int *signgamp) GET_FLOAT_WORD (hx, x); - if ((hx & 0x7fffffff) == 0) + if (__builtin_expect ((hx & 0x7fffffff) == 0, 0)) { /* Return value for x == 0 is Inf with divide by zero exception. */ *signgamp = 0; return 1.0 / x; } - if (hx < 0 && (u_int32_t) hx < 0xff800000 && __rintf (x) == x) + if (__builtin_expect (hx < 0, 0) + && (u_int32_t) hx < 0xff800000 && __rintf (x) == x) { /* Return value for integer x < 0 is NaN with invalid exception. */ *signgamp = 0; return (x - x) / (x - x); } - if (hx == 0xff800000) + if (__builtin_expect (hx == 0xff800000, 0)) { /* x == -Inf. According to ISO this is NaN. */ *signgamp = 0; @@ -54,3 +55,4 @@ __ieee754_gammaf_r (float x, int *signgamp) /* XXX FIXME. */ return __ieee754_expf (__ieee754_lgammaf_r (x, signgamp)); } +strong_alias (__ieee754_gammaf_r, __gammaf_r_finite) diff --git a/sysdeps/ieee754/flt-32/e_hypotf.c b/sysdeps/ieee754/flt-32/e_hypotf.c index a8e1a52d2b..7ec8ae66bc 100644 --- a/sysdeps/ieee754/flt-32/e_hypotf.c +++ b/sysdeps/ieee754/flt-32/e_hypotf.c @@ -13,19 +13,11 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_hypotf.c,v 1.5 1995/05/12 04:57:30 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ - float __ieee754_hypotf(float x, float y) -#else - float __ieee754_hypotf(x,y) - float x, y; -#endif +float +__ieee754_hypotf(float x, float y) { float a,b,t1,t2,y1,y2,w; int32_t j,k,ha,hb; @@ -39,7 +31,7 @@ static char rcsid[] = "$NetBSD: e_hypotf.c,v 1.5 1995/05/12 04:57:30 jtc Exp $"; SET_FLOAT_WORD(b,hb); /* b <- |b| */ if((ha-hb)>0xf000000) {return a+b;} /* x/y > 2**30 */ k=0; - if(ha > 0x58800000) { /* a>2**50 */ + if(__builtin_expect(ha > 0x58800000, 0)) { /* a>2**50 */ if(ha >= 0x7f800000) { /* Inf or NaN */ w = a+b; /* for sNaN */ if(ha == 0x7f800000) w = a; @@ -51,15 +43,15 @@ static char rcsid[] = "$NetBSD: e_hypotf.c,v 1.5 1995/05/12 04:57:30 jtc Exp $"; SET_FLOAT_WORD(a,ha); SET_FLOAT_WORD(b,hb); } - if(hb < 0x26800000) { /* b < 2**-50 */ + if(__builtin_expect(hb < 0x26800000, 0)) { /* b < 2**-50 */ if(hb <= 0x007fffff) { /* subnormal b or 0 */ - if(hb==0) return a; + if(hb==0) return a; SET_FLOAT_WORD(t1,0x7e800000); /* t1=2^126 */ b *= t1; a *= t1; k -= 126; } else { /* scale a and b by 2^60 */ - ha += 0x1e000000; /* a *= 2^60 */ + ha += 0x1e000000; /* a *= 2^60 */ hb += 0x1e000000; /* b *= 2^60 */ k -= 60; SET_FLOAT_WORD(a,ha); @@ -85,3 +77,4 @@ static char rcsid[] = "$NetBSD: e_hypotf.c,v 1.5 1995/05/12 04:57:30 jtc Exp $"; return t1*w; } else return w; } +strong_alias (__ieee754_hypotf, __hypotf_finite) diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c index 8c499e614e..d2da43f929 100644 --- a/sysdeps/ieee754/flt-32/e_j0f.c +++ b/sysdeps/ieee754/flt-32/e_j0f.c @@ -13,29 +13,17 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_j0f.c,v 1.4 1995/05/10 20:45:25 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static float pzerof(float), qzerof(float); -#else -static float pzerof(), qzerof(); -#endif -#ifdef __STDC__ static const float -#else -static float -#endif -huge = 1e30, +huge = 1e30, one = 1.0, invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ tpi = 6.3661974669e-01, /* 0x3f22f983 */ - /* R0/S0 on [0, 2.00] */ + /* R0/S0 on [0, 2.00] */ R02 = 1.5625000000e-02, /* 0x3c800000 */ R03 = -1.8997929874e-04, /* 0xb947352e */ R04 = 1.8295404516e-06, /* 0x35f58e88 */ @@ -45,18 +33,10 @@ S02 = 1.1692678527e-04, /* 0x38f53697 */ S03 = 5.1354652442e-07, /* 0x3509daa6 */ S04 = 1.1661400734e-09; /* 0x30a045e8 */ -#ifdef __STDC__ static const float zero = 0.0; -#else -static float zero = 0.0; -#endif -#ifdef __STDC__ - float __ieee754_j0f(float x) -#else - float __ieee754_j0f(x) - float x; -#endif +float +__ieee754_j0f(float x) { float z, s,c,ss,cc,r,u,v; int32_t hx,ix; @@ -72,7 +52,7 @@ static float zero = 0.0; if(ix<0x7f000000) { /* make sure x+x not overflow */ z = -__cosf(x+x); if ((s*c)<zero) cc = z/ss; - else ss = z/cc; + else ss = z/cc; } /* * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) @@ -87,8 +67,8 @@ static float zero = 0.0; } if(ix<0x39000000) { /* |x| < 2**-13 */ if(huge+x>one) { /* raise inexact if x != 0 */ - if(ix<0x32000000) return one; /* |x|<2**-27 */ - else return one - (float)0.25*x*x; + if(ix<0x32000000) return one; /* |x|<2**-27 */ + else return one - (float)0.25*x*x; } } z = x*x; @@ -101,12 +81,9 @@ static float zero = 0.0; return((one+u)*(one-u)+z*(r/s)); } } +strong_alias (__ieee754_j0f, __j0f_finite) -#ifdef __STDC__ static const float -#else -static float -#endif u00 = -7.3804296553e-02, /* 0xbd9726b5 */ u01 = 1.7666645348e-01, /* 0x3e34e80d */ u02 = -1.3818567619e-02, /* 0xbc626746 */ @@ -119,52 +96,48 @@ v02 = 7.6006865129e-05, /* 0x389f65e0 */ v03 = 2.5915085189e-07, /* 0x348b216c */ v04 = 4.4111031494e-10; /* 0x2ff280c2 */ -#ifdef __STDC__ - float __ieee754_y0f(float x) -#else - float __ieee754_y0f(x) - float x; -#endif +float +__ieee754_y0f(float x) { float z, s,c,ss,cc,u,v; int32_t hx,ix; GET_FLOAT_WORD(hx,x); - ix = 0x7fffffff&hx; + ix = 0x7fffffff&hx; /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */ if(ix>=0x7f800000) return one/(x+x*x); - if(ix==0) return -HUGE_VALF+x; /* -inf and overflow exception. */ - if(hx<0) return zero/(zero*x); - if(ix >= 0x40000000) { /* |x| >= 2.0 */ - /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) - * where x0 = x-pi/4 - * Better formula: - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) - * = 1/sqrt(2) * (sin(x) + cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ + if(ix==0) return -HUGE_VALF+x; /* -inf and overflow exception. */ + if(hx<0) return zero/(zero*x); + if(ix >= 0x40000000) { /* |x| >= 2.0 */ + /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) + * where x0 = x-pi/4 + * Better formula: + * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) + * = 1/sqrt(2) * (sin(x) + cos(x)) + * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) + * = 1/sqrt(2) * (sin(x) - cos(x)) + * To avoid cancellation, use + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + * to compute the worse one. + */ __sincosf (x, &s, &c); - ss = s-c; - cc = s+c; + ss = s-c; + cc = s+c; /* * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) */ - if(ix<0x7f000000) { /* make sure x+x not overflow */ - z = -__cosf(x+x); - if ((s*c)<zero) cc = z/ss; - else ss = z/cc; - } - if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x); - else { - u = pzerof(x); v = qzerof(x); - z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x); - } - return z; + if(ix<0x7f000000) { /* make sure x+x not overflow */ + z = -__cosf(x+x); + if ((s*c)<zero) cc = z/ss; + else ss = z/cc; + } + if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x); + else { + u = pzerof(x); v = qzerof(x); + z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x); + } + return z; } if(ix<=0x32000000) { /* x < 2**-27 */ return(u00 + tpi*__ieee754_logf(x)); @@ -174,21 +147,18 @@ v04 = 4.4111031494e-10; /* 0x2ff280c2 */ v = one+z*(v01+z*(v02+z*(v03+z*v04))); return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x))); } +strong_alias (__ieee754_y0f, __y0f_finite) /* The asymptotic expansions of pzero is * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. * For x >= 2, We approximate pzero by - * pzero(x) = 1 + (R/S) + * pzero(x) = 1 + (R/S) * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10 - * S = 1 + pS0*s^2 + ... + pS4*s^10 + * S = 1 + pS0*s^2 + ... + pS4*s^10 * and * | pzero(x)-1-R/S | <= 2 ** ( -60.26) */ -#ifdef __STDC__ static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#else -static float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#endif 0.0000000000e+00, /* 0x00000000 */ -7.0312500000e-02, /* 0xbd900000 */ -8.0816707611e+00, /* 0xc1014e86 */ @@ -196,22 +166,14 @@ static float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -2.4852163086e+03, /* 0xc51b5376 */ -5.2530439453e+03, /* 0xc5a4285a */ }; -#ifdef __STDC__ static const float pS8[5] = { -#else -static float pS8[5] = { -#endif 1.1653436279e+02, /* 0x42e91198 */ 3.8337448730e+03, /* 0x456f9beb */ 4.0597855469e+04, /* 0x471e95db */ 1.1675296875e+05, /* 0x47e4087c */ 4.7627726562e+04, /* 0x473a0bba */ }; -#ifdef __STDC__ static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#else -static float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#endif -1.1412546255e-11, /* 0xad48c58a */ -7.0312492549e-02, /* 0xbd8fffff */ -4.1596107483e+00, /* 0xc0851b88 */ @@ -219,11 +181,7 @@ static float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -3.3123129272e+02, /* 0xc3a59d9b */ -3.4643338013e+02, /* 0xc3ad3779 */ }; -#ifdef __STDC__ static const float pS5[5] = { -#else -static float pS5[5] = { -#endif 6.0753936768e+01, /* 0x42730408 */ 1.0512523193e+03, /* 0x44836813 */ 5.9789707031e+03, /* 0x45bad7c4 */ @@ -231,11 +189,7 @@ static float pS5[5] = { 2.4060581055e+03, /* 0x451660ee */ }; -#ifdef __STDC__ static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#else -static float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#endif -2.5470459075e-09, /* 0xb12f081b */ -7.0311963558e-02, /* 0xbd8fffb8 */ -2.4090321064e+00, /* 0xc01a2d95 */ @@ -243,11 +197,7 @@ static float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -5.8079170227e+01, /* 0xc2685112 */ -3.1447946548e+01, /* 0xc1fb9565 */ }; -#ifdef __STDC__ static const float pS3[5] = { -#else -static float pS3[5] = { -#endif 3.5856033325e+01, /* 0x420f6c94 */ 3.6151397705e+02, /* 0x43b4c1ca */ 1.1936077881e+03, /* 0x44953373 */ @@ -255,11 +205,7 @@ static float pS3[5] = { 1.7358093262e+02, /* 0x432d94b8 */ }; -#ifdef __STDC__ static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#else -static float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#endif -8.8753431271e-08, /* 0xb3be98b7 */ -7.0303097367e-02, /* 0xbd8ffb12 */ -1.4507384300e+00, /* 0xbfb9b1cc */ @@ -267,11 +213,7 @@ static float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -1.1193166733e+01, /* 0xc1331736 */ -3.2336456776e+00, /* 0xc04ef40d */ }; -#ifdef __STDC__ static const float pS2[5] = { -#else -static float pS2[5] = { -#endif 2.2220300674e+01, /* 0x41b1c32d */ 1.3620678711e+02, /* 0x430834f0 */ 2.7047027588e+02, /* 0x43873c32 */ @@ -279,18 +221,10 @@ static float pS2[5] = { 1.4657617569e+01, /* 0x416a859a */ }; -#ifdef __STDC__ - static float pzerof(float x) -#else - static float pzerof(x) - float x; -#endif +static float +pzerof(float x) { -#ifdef __STDC__ const float *p,*q; -#else - float *p,*q; -#endif float z,r,s; int32_t ix; GET_FLOAT_WORD(ix,x); @@ -309,17 +243,13 @@ static float pS2[5] = { /* For x >= 8, the asymptotic expansions of qzero is * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. * We approximate pzero by - * qzero(x) = s*(-1.25 + (R/S)) + * qzero(x) = s*(-1.25 + (R/S)) * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10 - * S = 1 + qS0*s^2 + ... + qS5*s^12 + * S = 1 + qS0*s^2 + ... + qS5*s^12 * and * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) */ -#ifdef __STDC__ static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#else -static float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#endif 0.0000000000e+00, /* 0x00000000 */ 7.3242187500e-02, /* 0x3d960000 */ 1.1768206596e+01, /* 0x413c4a93 */ @@ -327,11 +257,7 @@ static float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 8.8591972656e+03, /* 0x460a6cca */ 3.7014625000e+04, /* 0x471096a0 */ }; -#ifdef __STDC__ static const float qS8[6] = { -#else -static float qS8[6] = { -#endif 1.6377603149e+02, /* 0x4323c6aa */ 8.0983447266e+03, /* 0x45fd12c2 */ 1.4253829688e+05, /* 0x480b3293 */ @@ -340,11 +266,7 @@ static float qS8[6] = { -3.4389928125e+05, /* 0xc8a7eb69 */ }; -#ifdef __STDC__ static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#else -static float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#endif 1.8408595828e-11, /* 0x2da1ec79 */ 7.3242180049e-02, /* 0x3d95ffff */ 5.8356351852e+00, /* 0x40babd86 */ @@ -352,11 +274,7 @@ static float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 1.0272437744e+03, /* 0x448067cd */ 1.9899779053e+03, /* 0x44f8bf4b */ }; -#ifdef __STDC__ static const float qS5[6] = { -#else -static float qS5[6] = { -#endif 8.2776611328e+01, /* 0x42a58da0 */ 2.0778142090e+03, /* 0x4501dd07 */ 1.8847289062e+04, /* 0x46933e94 */ @@ -365,11 +283,7 @@ static float qS5[6] = { -5.3543427734e+03, /* 0xc5a752be */ }; -#ifdef __STDC__ static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#else -static float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#endif 4.3774099900e-09, /* 0x3196681b */ 7.3241114616e-02, /* 0x3d95ff70 */ 3.3442313671e+00, /* 0x405607e3 */ @@ -377,11 +291,7 @@ static float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 1.7080809021e+02, /* 0x432acedf */ 1.6673394775e+02, /* 0x4326bbe4 */ }; -#ifdef __STDC__ static const float qS3[6] = { -#else -static float qS3[6] = { -#endif 4.8758872986e+01, /* 0x42430916 */ 7.0968920898e+02, /* 0x44316c1c */ 3.7041481934e+03, /* 0x4567825f */ @@ -390,11 +300,7 @@ static float qS3[6] = { -1.4924745178e+02, /* 0xc3153f59 */ }; -#ifdef __STDC__ static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#else -static float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#endif 1.5044444979e-07, /* 0x342189db */ 7.3223426938e-02, /* 0x3d95f62a */ 1.9981917143e+00, /* 0x3fffc4bf */ @@ -402,11 +308,7 @@ static float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 3.1666231155e+01, /* 0x41fd5471 */ 1.6252708435e+01, /* 0x4182058c */ }; -#ifdef __STDC__ static const float qS2[6] = { -#else -static float qS2[6] = { -#endif 3.0365585327e+01, /* 0x41f2ecb8 */ 2.6934811401e+02, /* 0x4386ac8f */ 8.4478375244e+02, /* 0x44533229 */ @@ -415,18 +317,10 @@ static float qS2[6] = { -5.3109550476e+00, /* 0xc0a9f358 */ }; -#ifdef __STDC__ - static float qzerof(float x) -#else - static float qzerof(x) - float x; -#endif +static float +qzerof(float x) { -#ifdef __STDC__ const float *p,*q; -#else - float *p,*q; -#endif float s,r,z; int32_t ix; GET_FLOAT_WORD(ix,x); diff --git a/sysdeps/ieee754/flt-32/e_j1f.c b/sysdeps/ieee754/flt-32/e_j1f.c index 71bb2515af..bb335a7403 100644 --- a/sysdeps/ieee754/flt-32/e_j1f.c +++ b/sysdeps/ieee754/flt-32/e_j1f.c @@ -13,24 +13,12 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_j1f.c,v 1.4 1995/05/10 20:45:31 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static float ponef(float), qonef(float); -#else -static float ponef(), qonef(); -#endif -#ifdef __STDC__ static const float -#else -static float -#endif huge = 1e30, one = 1.0, invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ @@ -46,25 +34,17 @@ s03 = 1.1771846857e-06, /* 0x359dffc2 */ s04 = 5.0463624390e-09, /* 0x31ad6446 */ s05 = 1.2354227016e-11; /* 0x2d59567e */ -#ifdef __STDC__ static const float zero = 0.0; -#else -static float zero = 0.0; -#endif -#ifdef __STDC__ - float __ieee754_j1f(float x) -#else - float __ieee754_j1f(x) - float x; -#endif +float +__ieee754_j1f(float x) { float z, s,c,ss,cc,r,u,v,y; int32_t hx,ix; GET_FLOAT_WORD(hx,x); ix = hx&0x7fffffff; - if(ix>=0x7f800000) return one/x; + if(__builtin_expect(ix>=0x7f800000, 0)) return one/x; y = fabsf(x); if(ix >= 0x40000000) { /* |x| >= 2.0 */ __sincosf (y, &s, &c); @@ -73,7 +53,7 @@ static float zero = 0.0; if(ix<0x7f000000) { /* make sure y+y not overflow */ z = __cosf(y+y); if ((s*c)>zero) cc = z/ss; - else ss = z/cc; + else ss = z/cc; } /* * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) @@ -85,9 +65,9 @@ static float zero = 0.0; z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrtf(y); } if(hx<0) return -z; - else return z; + else return z; } - if(ix<0x32000000) { /* |x|<2**-27 */ + if(__builtin_expect(ix<0x32000000, 0)) { /* |x|<2**-27 */ if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */ } z = x*x; @@ -96,23 +76,16 @@ static float zero = 0.0; r *= x; return(x*(float)0.5+r/s); } +strong_alias (__ieee754_j1f, __j1f_finite) -#ifdef __STDC__ static const float U0[5] = { -#else -static float U0[5] = { -#endif -1.9605709612e-01, /* 0xbe48c331 */ 5.0443872809e-02, /* 0x3d4e9e3c */ -1.9125689287e-03, /* 0xbafaaf2a */ 2.3525259166e-05, /* 0x37c5581c */ -9.1909917899e-08, /* 0xb3c56003 */ }; -#ifdef __STDC__ static const float V0[5] = { -#else -static float V0[5] = { -#endif 1.9916731864e-02, /* 0x3ca3286a */ 2.0255257550e-04, /* 0x3954644b */ 1.3560879779e-06, /* 0x35b602d4 */ @@ -120,73 +93,67 @@ static float V0[5] = { 1.6655924903e-11, /* 0x2d9281cf */ }; -#ifdef __STDC__ - float __ieee754_y1f(float x) -#else - float __ieee754_y1f(x) - float x; -#endif +float +__ieee754_y1f(float x) { float z, s,c,ss,cc,u,v; int32_t hx,ix; GET_FLOAT_WORD(hx,x); - ix = 0x7fffffff&hx; + ix = 0x7fffffff&hx; /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if(ix>=0x7f800000) return one/(x+x*x); - if(ix==0) return -HUGE_VALF+x; /* -inf and overflow exception. */ - if(hx<0) return zero/(zero*x); - if(ix >= 0x40000000) { /* |x| >= 2.0 */ + if(__builtin_expect(ix>=0x7f800000, 0)) return one/(x+x*x); + if(__builtin_expect(ix==0, 0)) + return -HUGE_VALF+x; /* -inf and overflow exception. */ + if(__builtin_expect(hx<0, 0)) return zero/(zero*x); + if(ix >= 0x40000000) { /* |x| >= 2.0 */ __sincosf (x, &s, &c); - ss = -s-c; - cc = s-c; - if(ix<0x7f000000) { /* make sure x+x not overflow */ - z = __cosf(x+x); - if ((s*c)>zero) cc = z/ss; - else ss = z/cc; - } - /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) - * where x0 = x-3pi/4 - * Better formula: - * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = -1/sqrt(2) * (cos(x) + sin(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ - if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x); - else { - u = ponef(x); v = qonef(x); - z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x); - } - return z; - } - if(ix<=0x24800000) { /* x < 2**-54 */ - return(-tpi/x); - } - z = x*x; - u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); - v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); - return(x*(u/v) + tpi*(__ieee754_j1f(x)*__ieee754_logf(x)-one/x)); + ss = -s-c; + cc = s-c; + if(ix<0x7f000000) { /* make sure x+x not overflow */ + z = __cosf(x+x); + if ((s*c)>zero) cc = z/ss; + else ss = z/cc; + } + /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) + * where x0 = x-3pi/4 + * Better formula: + * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) + * = 1/sqrt(2) * (sin(x) - cos(x)) + * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) + * = -1/sqrt(2) * (cos(x) + sin(x)) + * To avoid cancellation, use + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + * to compute the worse one. + */ + if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x); + else { + u = ponef(x); v = qonef(x); + z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x); + } + return z; + } + if(__builtin_expect(ix<=0x24800000, 0)) { /* x < 2**-54 */ + return(-tpi/x); + } + z = x*x; + u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); + v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); + return(x*(u/v) + tpi*(__ieee754_j1f(x)*__ieee754_logf(x)-one/x)); } +strong_alias (__ieee754_y1f, __y1f_finite) /* For x >= 8, the asymptotic expansions of pone is * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. * We approximate pone by - * pone(x) = 1 + (R/S) + * pone(x) = 1 + (R/S) * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 - * S = 1 + ps0*s^2 + ... + ps4*s^10 + * S = 1 + ps0*s^2 + ... + ps4*s^10 * and * | pone(x)-1-R/S | <= 2 ** ( -60.06) */ -#ifdef __STDC__ static const float pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#else -static float pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#endif 0.0000000000e+00, /* 0x00000000 */ 1.1718750000e-01, /* 0x3df00000 */ 1.3239480972e+01, /* 0x4153d4ea */ @@ -194,11 +161,7 @@ static float pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 3.8747453613e+03, /* 0x45722bed */ 7.9144794922e+03, /* 0x45f753d6 */ }; -#ifdef __STDC__ static const float ps8[5] = { -#else -static float ps8[5] = { -#endif 1.1420736694e+02, /* 0x42e46a2c */ 3.6509309082e+03, /* 0x45642ee5 */ 3.6956207031e+04, /* 0x47105c35 */ @@ -206,11 +169,7 @@ static float ps8[5] = { 3.0804271484e+04, /* 0x46f0a88b */ }; -#ifdef __STDC__ static const float pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#else -static float pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#endif 1.3199052094e-11, /* 0x2d68333f */ 1.1718749255e-01, /* 0x3defffff */ 6.8027510643e+00, /* 0x40d9b023 */ @@ -218,11 +177,7 @@ static float pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 5.1763616943e+02, /* 0x440168b7 */ 5.2871520996e+02, /* 0x44042dc6 */ }; -#ifdef __STDC__ static const float ps5[5] = { -#else -static float ps5[5] = { -#endif 5.9280597687e+01, /* 0x426d1f55 */ 9.9140142822e+02, /* 0x4477d9b1 */ 5.3532670898e+03, /* 0x45a74a23 */ @@ -230,11 +185,7 @@ static float ps5[5] = { 1.5040468750e+03, /* 0x44bc0180 */ }; -#ifdef __STDC__ static const float pr3[6] = { -#else -static float pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#endif 3.0250391081e-09, /* 0x314fe10d */ 1.1718686670e-01, /* 0x3defffab */ 3.9329774380e+00, /* 0x407bb5e7 */ @@ -242,11 +193,7 @@ static float pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 9.1055007935e+01, /* 0x42b61c2a */ 4.8559066772e+01, /* 0x42423c7c */ }; -#ifdef __STDC__ static const float ps3[5] = { -#else -static float ps3[5] = { -#endif 3.4791309357e+01, /* 0x420b2a4d */ 3.3676245117e+02, /* 0x43a86198 */ 1.0468714600e+03, /* 0x4482dbe3 */ @@ -254,11 +201,7 @@ static float ps3[5] = { 1.0378793335e+02, /* 0x42cf936c */ }; -#ifdef __STDC__ static const float pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#else -static float pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#endif 1.0771083225e-07, /* 0x33e74ea8 */ 1.1717621982e-01, /* 0x3deffa16 */ 2.3685150146e+00, /* 0x401795c0 */ @@ -266,11 +209,7 @@ static float pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 1.7693971634e+01, /* 0x418d8d41 */ 5.0735230446e+00, /* 0x40a25a4d */ }; -#ifdef __STDC__ static const float ps2[5] = { -#else -static float ps2[5] = { -#endif 2.1436485291e+01, /* 0x41ab7dec */ 1.2529022980e+02, /* 0x42fa9499 */ 2.3227647400e+02, /* 0x436846c7 */ @@ -278,48 +217,36 @@ static float ps2[5] = { 8.3646392822e+00, /* 0x4105d590 */ }; -#ifdef __STDC__ - static float ponef(float x) -#else - static float ponef(x) - float x; -#endif +static float +ponef(float x) { -#ifdef __STDC__ const float *p,*q; -#else - float *p,*q; -#endif float z,r,s; - int32_t ix; + int32_t ix; GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; - if(ix>=0x41000000) {p = pr8; q= ps8;} - else if(ix>=0x40f71c58){p = pr5; q= ps5;} - else if(ix>=0x4036db68){p = pr3; q= ps3;} - else if(ix>=0x40000000){p = pr2; q= ps2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return one+ r/s; + if(ix>=0x41000000) {p = pr8; q= ps8;} + else if(ix>=0x40f71c58){p = pr5; q= ps5;} + else if(ix>=0x4036db68){p = pr3; q= ps3;} + else if(ix>=0x40000000){p = pr2; q= ps2;} + z = one/(x*x); + r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); + s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); + return one+ r/s; } /* For x >= 8, the asymptotic expansions of qone is * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. * We approximate pone by - * qone(x) = s*(0.375 + (R/S)) + * qone(x) = s*(0.375 + (R/S)) * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 - * S = 1 + qs1*s^2 + ... + qs6*s^12 + * S = 1 + qs1*s^2 + ... + qs6*s^12 * and * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) */ -#ifdef __STDC__ static const float qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#else -static float qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#endif 0.0000000000e+00, /* 0x00000000 */ -1.0253906250e-01, /* 0xbdd20000 */ -1.6271753311e+01, /* 0xc1822c8d */ @@ -327,11 +254,7 @@ static float qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -1.1849806641e+04, /* 0xc639273a */ -4.8438511719e+04, /* 0xc73d3683 */ }; -#ifdef __STDC__ static const float qs8[6] = { -#else -static float qs8[6] = { -#endif 1.6139537048e+02, /* 0x43216537 */ 7.8253862305e+03, /* 0x45f48b17 */ 1.3387534375e+05, /* 0x4802bcd6 */ @@ -340,11 +263,7 @@ static float qs8[6] = { -2.9449025000e+05, /* 0xc88fcb48 */ }; -#ifdef __STDC__ static const float qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#else -static float qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#endif -2.0897993405e-11, /* 0xadb7d219 */ -1.0253904760e-01, /* 0xbdd1fffe */ -8.0564479828e+00, /* 0xc100e736 */ @@ -352,11 +271,7 @@ static float qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -1.3731937256e+03, /* 0xc4aba633 */ -2.6124443359e+03, /* 0xc523471c */ }; -#ifdef __STDC__ static const float qs5[6] = { -#else -static float qs5[6] = { -#endif 8.1276550293e+01, /* 0x42a28d98 */ 1.9917987061e+03, /* 0x44f8f98f */ 1.7468484375e+04, /* 0x468878f8 */ @@ -365,11 +280,7 @@ static float qs5[6] = { -4.7191835938e+03, /* 0xc5937978 */ }; -#ifdef __STDC__ static const float qr3[6] = { -#else -static float qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#endif -5.0783124372e-09, /* 0xb1ae7d4f */ -1.0253783315e-01, /* 0xbdd1ff5b */ -4.6101160049e+00, /* 0xc0938612 */ @@ -377,11 +288,7 @@ static float qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -2.2824453735e+02, /* 0xc3643e9a */ -2.1921012878e+02, /* 0xc35b35cb */ }; -#ifdef __STDC__ static const float qs3[6] = { -#else -static float qs3[6] = { -#endif 4.7665153503e+01, /* 0x423ea91e */ 6.7386511230e+02, /* 0x4428775e */ 3.3801528320e+03, /* 0x45534272 */ @@ -390,11 +297,7 @@ static float qs3[6] = { -1.3520118713e+02, /* 0xc3073381 */ }; -#ifdef __STDC__ static const float qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#else -static float qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#endif -1.7838172539e-07, /* 0xb43f8932 */ -1.0251704603e-01, /* 0xbdd1f475 */ -2.7522056103e+00, /* 0xc0302423 */ @@ -402,11 +305,7 @@ static float qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -4.2325313568e+01, /* 0xc2294d1f */ -2.1371921539e+01, /* 0xc1aaf9b2 */ }; -#ifdef __STDC__ static const float qs2[6] = { -#else -static float qs2[6] = { -#endif 2.9533363342e+01, /* 0x41ec4454 */ 2.5298155212e+02, /* 0x437cfb47 */ 7.5750280762e+02, /* 0x443d602e */ @@ -415,18 +314,10 @@ static float qs2[6] = { -4.9594988823e+00, /* 0xc09eb437 */ }; -#ifdef __STDC__ - static float qonef(float x) -#else - static float qonef(x) - float x; -#endif +static float +qonef(float x) { -#ifdef __STDC__ const float *p,*q; -#else - float *p,*q; -#endif float s,r,z; int32_t ix; GET_FLOAT_WORD(ix,x); diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c index dd3d551a39..1e55485b02 100644 --- a/sysdeps/ieee754/flt-32/e_jnf.c +++ b/sysdeps/ieee754/flt-32/e_jnf.c @@ -13,33 +13,17 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_jnf.c,v 1.5 1995/05/10 20:45:37 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float -#else -static float -#endif two = 2.0000000000e+00, /* 0x40000000 */ one = 1.0000000000e+00; /* 0x3F800000 */ -#ifdef __STDC__ static const float zero = 0.0000000000e+00; -#else -static float zero = 0.0000000000e+00; -#endif -#ifdef __STDC__ - float __ieee754_jnf(int n, float x) -#else - float __ieee754_jnf(n,x) - int n; float x; -#endif +float +__ieee754_jnf(int n, float x) { int32_t i,hx,ix, sgn; float a, b, temp, di; @@ -51,7 +35,7 @@ static float zero = 0.0000000000e+00; GET_FLOAT_WORD(hx,x); ix = 0x7fffffff&hx; /* if J(n,NaN) is NaN */ - if(ix>0x7f800000) return x+x; + if(__builtin_expect(ix>0x7f800000, 0)) return x+x; if(n<0){ n = -n; x = -x; @@ -61,7 +45,7 @@ static float zero = 0.0000000000e+00; if(n==1) return(__ieee754_j1f(x)); sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ x = fabsf(x); - if(ix==0||ix>=0x7f800000) /* if x is 0 or inf */ + if(__builtin_expect(ix==0||ix>=0x7f800000, 0)) /* if x is 0 or inf */ b = zero; else if((float)n<=x) { /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ @@ -106,7 +90,7 @@ static float zero = 0.0000000000e+00; * 1 * w - ----------------- * 1 - * w+h - --------- + * w+h - --------- * w+2h - ... * * To determine how many terms needed, let @@ -144,18 +128,18 @@ static float zero = 0.0000000000e+00; tmp = tmp*__ieee754_logf(fabsf(v*tmp)); if(tmp<(float)8.8721679688e+01) { for(i=n-1,di=(float)(i+i);i>0;i--){ - temp = b; + temp = b; b *= di; b = b/x - a; - a = temp; + a = temp; di -= two; } } else { for(i=n-1,di=(float)(i+i);i>0;i--){ - temp = b; + temp = b; b *= di; b = b/x - a; - a = temp; + a = temp; di -= two; /* scale b to avoid spurious overflow */ if(b>(float)1e10) { @@ -179,13 +163,10 @@ static float zero = 0.0000000000e+00; } if(sgn==1) return -b; else return b; } +strong_alias (__ieee754_jnf, __jnf_finite) -#ifdef __STDC__ - float __ieee754_ynf(int n, float x) -#else - float __ieee754_ynf(n,x) - int n; float x; -#endif +float +__ieee754_ynf(int n, float x) { int32_t i,hx,ix; u_int32_t ib; @@ -195,9 +176,10 @@ static float zero = 0.0000000000e+00; GET_FLOAT_WORD(hx,x); ix = 0x7fffffff&hx; /* if Y(n,NaN) is NaN */ - if(ix>0x7f800000) return x+x; - if(ix==0) return -HUGE_VALF+x; /* -inf and overflow exception. */ - if(hx<0) return zero/(zero*x); + if(__builtin_expect(ix>0x7f800000, 0)) return x+x; + if(__builtin_expect(ix==0, 0)) + return -HUGE_VALF+x; /* -inf and overflow exception. */ + if(__builtin_expect(hx<0, 0)) return zero/(zero*x); sign = 1; if(n<0){ n = -n; @@ -205,7 +187,7 @@ static float zero = 0.0000000000e+00; } if(n==0) return(__ieee754_y0f(x)); if(n==1) return(sign*__ieee754_y1f(x)); - if(ix==0x7f800000) return zero; + if(__builtin_expect(ix==0x7f800000, 0)) return zero; a = __ieee754_y0f(x); b = __ieee754_y1f(x); @@ -219,3 +201,4 @@ static float zero = 0.0000000000e+00; } if(sign>0) return b; else return -b; } +strong_alias (__ieee754_ynf, __ynf_finite) diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c index 0ed2610085..cbee9db8fb 100644 --- a/sysdeps/ieee754/flt-32/e_lgammaf_r.c +++ b/sysdeps/ieee754/flt-32/e_lgammaf_r.c @@ -13,18 +13,10 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_lgammaf_r.c,v 1.3 1995/05/10 20:45:47 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float -#else -static float -#endif two23= 8.3886080000e+06, /* 0x4b000000 */ half= 5.0000000000e-01, /* 0x3f000000 */ one = 1.0000000000e+00, /* 0x3f800000 */ @@ -92,18 +84,10 @@ w4 = -5.9518753551e-04, /* 0xba1c065c */ w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ -#ifdef __STDC__ static const float zero= 0.0000000000e+00; -#else -static float zero= 0.0000000000e+00; -#endif -#ifdef __STDC__ - static float sin_pif(float x) -#else - static float sin_pif(x) - float x; -#endif +static float +sin_pif(float x) { float y,z; int n,ix; @@ -124,16 +108,16 @@ static float zero= 0.0000000000e+00; y = (float)2.0*(y - __floorf(y)); /* y = |x| mod 2.0 */ n = (int) (y*(float)4.0); } else { - if(ix>=0x4b800000) { - y = zero; n = 0; /* y must be even */ - } else { - if(ix<0x4b000000) z = y+two23; /* exact */ + if(ix>=0x4b800000) { + y = zero; n = 0; /* y must be even */ + } else { + if(ix<0x4b000000) z = y+two23; /* exact */ GET_FLOAT_WORD(n,z); n &= 1; - y = n; - n<<= 2; - } - } + y = n; + n<<= 2; + } + } switch (n) { case 0: y = __kernel_sinf(pi*y,zero,0); break; case 1: @@ -148,12 +132,8 @@ static float zero= 0.0000000000e+00; } -#ifdef __STDC__ - float __ieee754_lgammaf_r(float x, int *signgamp) -#else - float __ieee754_lgammaf_r(x,signgamp) - float x; int *signgamp; -#endif +float +__ieee754_lgammaf_r(float x, int *signgamp) { float t,y,z,nadj,p,p1,p2,p3,q,r,w; int i,hx,ix; @@ -163,21 +143,22 @@ static float zero= 0.0000000000e+00; /* purge off +-inf, NaN, +-0, and negative arguments */ *signgamp = 1; ix = hx&0x7fffffff; - if(ix>=0x7f800000) return x*x; - if(ix==0) + if(__builtin_expect(ix>=0x7f800000, 0)) return x*x; + if(__builtin_expect(ix==0, 0)) { if (hx < 0) *signgamp = -1; return one/fabsf(x); } - if(ix<0x1c800000) { /* |x|<2**-70, return -log(|x|) */ + if(__builtin_expect(ix<0x1c800000, 0)) { + /* |x|<2**-70, return -log(|x|) */ if(hx<0) { - *signgamp = -1; - return -__ieee754_logf(-x); + *signgamp = -1; + return -__ieee754_logf(-x); } else return -__ieee754_logf(x); } if(hx<0) { - if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ + if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ return x/zero; t = sin_pif(x); if(t==zero) return one/fabsf(t); /* -integer */ @@ -190,15 +171,15 @@ static float zero= 0.0000000000e+00; if (ix==0x3f800000||ix==0x40000000) r = 0; /* for x < 2.0 */ else if(ix<0x40000000) { - if(ix<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ + if(ix<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ r = -__ieee754_logf(x); if(ix>=0x3f3b4a20) {y = one-x; i= 0;} else if(ix>=0x3e6d3308) {y= x-(tc-one); i=1;} - else {y = x; i=2;} + else {y = x; i=2;} } else { - r = zero; - if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */ - else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */ + r = zero; + if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */ + else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */ else {y=x-one;i=2;} } switch(i) { @@ -222,7 +203,7 @@ static float zero= 0.0000000000e+00; r += (-(float)0.5*y + p1/p2); } } - else if(ix<0x41000000) { /* x < 8.0 */ + else if(ix<0x41000000) { /* x < 8.0 */ i = (int)x; t = zero; y = x-(float)i; @@ -251,3 +232,4 @@ static float zero= 0.0000000000e+00; if(hx<0) r = nadj - r; return r; } +strong_alias (__ieee754_lgammaf_r, __lgammaf_r_finite) diff --git a/sysdeps/ieee754/flt-32/e_log10f.c b/sysdeps/ieee754/flt-32/e_log10f.c index cea3d9156b..72dcea6250 100644 --- a/sysdeps/ieee754/flt-32/e_log10f.c +++ b/sysdeps/ieee754/flt-32/e_log10f.c @@ -13,55 +13,41 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_log10f.c,v 1.5 1995/05/10 20:45:53 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float -#else -static float -#endif two25 = 3.3554432000e+07, /* 0x4c000000 */ ivln10 = 4.3429449201e-01, /* 0x3ede5bd9 */ log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ log10_2lo = 7.9034151668e-07; /* 0x355427db */ -#ifdef __STDC__ static const float zero = 0.0; -#else -static float zero = 0.0; -#endif -#ifdef __STDC__ - float __ieee754_log10f(float x) -#else - float __ieee754_log10f(x) - float x; -#endif +float +__ieee754_log10f(float x) { float y,z; int32_t i,k,hx; GET_FLOAT_WORD(hx,x); - k=0; - if (hx < 0x00800000) { /* x < 2**-126 */ - if ((hx&0x7fffffff)==0) - return -two25/(x-x); /* log(+-0)=-inf */ - if (hx<0) return (x-x)/(x-x); /* log(-#) = NaN */ - k -= 25; x *= two25; /* subnormal number, scale up x */ + k=0; + if (hx < 0x00800000) { /* x < 2**-126 */ + if (__builtin_expect((hx&0x7fffffff)==0, 0)) + return -two25/(x-x); /* log(+-0)=-inf */ + if (__builtin_expect(hx<0, 0)) + return (x-x)/(x-x); /* log(-#) = NaN */ + k -= 25; x *= two25; /* subnormal number, scale up x */ GET_FLOAT_WORD(hx,x); - } - if (hx >= 0x7f800000) return x+x; + } + if (__builtin_expect(hx >= 0x7f800000, 0)) return x+x; k += (hx>>23)-127; i = ((u_int32_t)k&0x80000000)>>31; - hx = (hx&0x007fffff)|((0x7f-i)<<23); - y = (float)(k+i); + hx = (hx&0x007fffff)|((0x7f-i)<<23); + y = (float)(k+i); SET_FLOAT_WORD(x,hx); z = y*log10_2lo + ivln10*__ieee754_logf(x); return z+y*log10_2hi; } +strong_alias (__ieee754_log10f, __log10f_finite) diff --git a/sysdeps/ieee754/flt-32/e_log2f.c b/sysdeps/ieee754/flt-32/e_log2f.c index af3c6eadfe..7453214516 100644 --- a/sysdeps/ieee754/flt-32/e_log2f.c +++ b/sysdeps/ieee754/flt-32/e_log2f.c @@ -18,11 +18,7 @@ #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float -#else -static float -#endif ln2 = 0.69314718055994530942, two25 = 3.355443200e+07, /* 0x4c000000 */ Lg1 = 6.6666668653e-01, /* 3F2AAAAB */ @@ -33,18 +29,10 @@ Lg5 = 1.8183572590e-01, /* 3E3A3325 */ Lg6 = 1.5313838422e-01, /* 3E1CD04F */ Lg7 = 1.4798198640e-01; /* 3E178897 */ -#ifdef __STDC__ static const float zero = 0.0; -#else -static float zero = 0.0; -#endif -#ifdef __STDC__ - float __ieee754_log2f(float x) -#else - float __ieee754_log2f(x) - float x; -#endif +float +__ieee754_log2f(float x) { float hfsq,f,s,z,R,w,t1,t2,dk; int32_t k,ix,i,j; @@ -53,13 +41,14 @@ static float zero = 0.0; k=0; if (ix < 0x00800000) { /* x < 2**-126 */ - if ((ix&0x7fffffff)==0) + if (__builtin_expect((ix&0x7fffffff)==0, 0)) return -two25/(x-x); /* log(+-0)=-inf */ - if (ix<0) return (x-x)/(x-x); /* log(-#) = NaN */ + if (__builtin_expect(ix<0, 0)) + return (x-x)/(x-x); /* log(-#) = NaN */ k -= 25; x *= two25; /* subnormal number, scale up x */ GET_FLOAT_WORD(ix,x); } - if (ix >= 0x7f800000) return x+x; + if (__builtin_expect(ix >= 0x7f800000, 0)) return x+x; k += (ix>>23)-127; ix &= 0x007fffff; i = (ix+(0x95f64<<3))&0x800000; @@ -72,7 +61,7 @@ static float zero = 0.0; R = f*f*((float)0.5-(float)0.33333333333333333*f); return dk-(R-f)/ln2; } - s = f/((float)2.0+f); + s = f/((float)2.0+f); z = s*s; i = ix-(0x6147a<<3); w = z*z; @@ -88,3 +77,4 @@ static float zero = 0.0; return dk-((s*(f-R))-f)/ln2; } } +strong_alias (__ieee754_log2f, __log2f_finite) diff --git a/sysdeps/ieee754/flt-32/e_logf.c b/sysdeps/ieee754/flt-32/e_logf.c index de8f869df4..b870b31957 100644 --- a/sysdeps/ieee754/flt-32/e_logf.c +++ b/sysdeps/ieee754/flt-32/e_logf.c @@ -13,18 +13,10 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_logf.c,v 1.4 1995/05/10 20:45:54 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float -#else -static float -#endif ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ two25 = 3.355443200e+07, /* 0x4c000000 */ @@ -36,18 +28,10 @@ Lg5 = 1.8183572590e-01, /* 3E3A3325 */ Lg6 = 1.5313838422e-01, /* 3E1CD04F */ Lg7 = 1.4798198640e-01; /* 3E178897 */ -#ifdef __STDC__ static const float zero = 0.0; -#else -static float zero = 0.0; -#endif -#ifdef __STDC__ - float __ieee754_logf(float x) -#else - float __ieee754_logf(x) - float x; -#endif +float +__ieee754_logf(float x) { float hfsq,f,s,z,R,w,t1,t2,dk; int32_t k,ix,i,j; @@ -56,13 +40,14 @@ static float zero = 0.0; k=0; if (ix < 0x00800000) { /* x < 2**-126 */ - if ((ix&0x7fffffff)==0) + if (__builtin_expect((ix&0x7fffffff)==0, 0)) return -two25/(x-x); /* log(+-0)=-inf */ - if (ix<0) return (x-x)/(x-x); /* log(-#) = NaN */ + if (__builtin_expect(ix<0, 0)) + return (x-x)/(x-x); /* log(-#) = NaN */ k -= 25; x *= two25; /* subnormal number, scale up x */ GET_FLOAT_WORD(ix,x); } - if (ix >= 0x7f800000) return x+x; + if (__builtin_expect(ix >= 0x7f800000, 0)) return x+x; k += (ix>>23)-127; ix &= 0x007fffff; i = (ix+(0x95f64<<3))&0x800000; @@ -76,9 +61,9 @@ static float zero = 0.0; } R = f*f*((float)0.5-(float)0.33333333333333333*f); if(k==0) return f-R; else {dk=(float)k; - return dk*ln2_hi-((R-dk*ln2_lo)-f);} + return dk*ln2_hi-((R-dk*ln2_lo)-f);} } - s = f/((float)2.0+f); + s = f/((float)2.0+f); dk = (float)k; z = s*s; i = ix-(0x6147a<<3); @@ -97,3 +82,4 @@ static float zero = 0.0; return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); } } +strong_alias (__ieee754_logf, __logf_finite) diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c index 9f520804a4..460055204f 100644 --- a/sysdeps/ieee754/flt-32/e_powf.c +++ b/sysdeps/ieee754/flt-32/e_powf.c @@ -13,20 +13,12 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_powf.c,v 1.7 1996/04/08 15:43:44 phil Exp $"; -#endif - #include "math.h" #include "math_private.h" static const float huge = 1.0e+30, tiny = 1.0e-30; -#ifdef __STDC__ static const float -#else -static float -#endif bp[] = {1.0, 1.5,}, dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */ dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */ @@ -57,12 +49,8 @@ ivln2 = 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */ ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ -#ifdef __STDC__ - float __ieee754_powf(float x, float y) -#else - float __ieee754_powf(x,y) - float x, y; -#endif +float +__ieee754_powf(float x, float y) { float z,ax,z_h,z_l,p_h,p_l; float y1,t1,t2,r,s,t,u,v,w; @@ -81,8 +69,8 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ if(x == -1.0 && isinf(y)) return one; /* +-NaN return x+y */ - if(ix > 0x7f800000 || - iy > 0x7f800000) + if(__builtin_expect(ix > 0x7f800000 || + iy > 0x7f800000, 0)) return x+y; /* determine if y is an odd int when x < 0 @@ -101,26 +89,26 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ } /* special value of y */ - if (iy==0x7f800000) { /* y is +-inf */ + if (__builtin_expect(iy==0x7f800000, 0)) { /* y is +-inf */ if (ix==0x3f800000) - return y - y; /* inf**+-1 is NaN */ + return y - y; /* inf**+-1 is NaN */ else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */ - return (hy>=0)? y: zero; + return (hy>=0)? y: zero; else /* (|x|<1)**-,+inf = inf,0 */ - return (hy<0)?-y: zero; + return (hy<0)?-y: zero; } if(iy==0x3f800000) { /* y is +-1 */ if(hy<0) return one/x; else return x; } if(hy==0x40000000) return x*x; /* y is 2 */ if(hy==0x3f000000) { /* y is 0.5 */ - if(hx>=0) /* x >= +0 */ + if(__builtin_expect(hx>=0, 1)) /* x >= +0 */ return __ieee754_sqrtf(x); } ax = fabsf(x); /* special value of x */ - if(ix==0x7f800000||ix==0||ix==0x3f800000){ + if(__builtin_expect(ix==0x7f800000||ix==0||ix==0x3f800000, 0)){ z = ax; /*x is +-0,+-inf,+-1*/ if(hy<0) z = one/z; /* z = (1/|x|) */ if(hx<0) { @@ -133,10 +121,11 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ } /* (x<0)**(non-int) is NaN */ - if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x); + if(__builtin_expect(((((u_int32_t)hx>>31)-1)|yisint)==0, 0)) + return (x-x)/(x-x); /* |y| is huge */ - if(iy>0x4d000000) { /* if |y| > 2**27 */ + if(__builtin_expect(iy>0x4d000000, 0)) { /* if |y| > 2**27 */ /* over/underflow if x is not close to one */ if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny; if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny; @@ -214,14 +203,14 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ p_h = y1*t1; z = p_l+p_h; GET_FLOAT_WORD(j,z); - if (j>0x43000000) /* if z > 128 */ + if (__builtin_expect(j>0x43000000, 0)) /* if z > 128 */ return s*huge*huge; /* overflow */ - else if (j==0x43000000) { /* if z == 128 */ + else if (__builtin_expect(j==0x43000000, 0)) { /* if z == 128 */ if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ } - else if ((j&0x7fffffff)>0x43160000) /* z <= -150 */ + else if (__builtin_expect((j&0x7fffffff)>0x43160000, 0))/* z <= -150 */ return s*tiny*tiny; /* underflow */ - else if ((u_int32_t) j==0xc3160000){ /* z == -150 */ + else if (__builtin_expect((u_int32_t) j==0xc3160000, 0)){/* z == -150*/ if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ } /* @@ -255,3 +244,4 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ else SET_FLOAT_WORD(z,j); return s*z; } +strong_alias (__ieee754_powf, __powf_finite) diff --git a/sysdeps/ieee754/flt-32/e_remainderf.c b/sysdeps/ieee754/flt-32/e_remainderf.c index 90d0d366d4..aaf15dfe54 100644 --- a/sysdeps/ieee754/flt-32/e_remainderf.c +++ b/sysdeps/ieee754/flt-32/e_remainderf.c @@ -8,31 +8,19 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_remainderf.c,v 1.4 1995/05/10 20:46:08 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float zero = 0.0; -#else -static float zero = 0.0; -#endif -#ifdef __STDC__ - float __ieee754_remainderf(float x, float p) -#else - float __ieee754_remainderf(x,p) - float x,p; -#endif +float +__ieee754_remainderf(float x, float p) { int32_t hx,hp; u_int32_t sx; @@ -45,7 +33,7 @@ static float zero = 0.0; hx &= 0x7fffffff; /* purge off exception values */ - if(hp==0) return (x*p)/(x*p); /* p = 0 */ + if(hp==0) return (x*p)/(x*p); /* p = 0 */ if((hx>=0x7f800000)|| /* x not finite */ ((hp>0x7f800000))) /* p is NaN */ return (x*p)/(x*p); @@ -71,3 +59,4 @@ static float zero = 0.0; SET_FLOAT_WORD(x,hx^sx); return x; } +strong_alias (__ieee754_remainderf, __remainderf_finite) diff --git a/sysdeps/ieee754/flt-32/e_sinhf.c b/sysdeps/ieee754/flt-32/e_sinhf.c index 045f6f121e..5813963488 100644 --- a/sysdeps/ieee754/flt-32/e_sinhf.c +++ b/sysdeps/ieee754/flt-32/e_sinhf.c @@ -8,31 +8,19 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_sinhf.c,v 1.4 1995/05/10 20:46:15 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float one = 1.0, shuge = 1.0e37; -#else -static float one = 1.0, shuge = 1.0e37; -#endif -#ifdef __STDC__ - float __ieee754_sinhf(float x) -#else - float __ieee754_sinhf(x) - float x; -#endif -{ +float +__ieee754_sinhf(float x) +{ float t,w,h; int32_t ix,jx; @@ -40,13 +28,13 @@ static float one = 1.0, shuge = 1.0e37; ix = jx&0x7fffffff; /* x is INF or NaN */ - if(ix>=0x7f800000) return x+x; + if(__builtin_expect(ix>=0x7f800000, 0)) return x+x; h = 0.5; if (jx<0) h = -h; /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ if (ix < 0x41b00000) { /* |x|<22 */ - if (ix<0x31800000) /* |x|<2**-28 */ + if (__builtin_expect(ix<0x31800000, 0)) /* |x|<2**-28 */ if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ t = __expm1f(fabsf(x)); if(ix<0x3f800000) return h*((float)2.0*t-t*t/(t+one)); @@ -66,3 +54,4 @@ static float one = 1.0, shuge = 1.0e37; /* |x| > overflowthresold, sinh(x) overflow */ return x*shuge; } +strong_alias (__ieee754_sinhf, __sinhf_finite) diff --git a/sysdeps/ieee754/flt-32/e_sqrtf.c b/sysdeps/ieee754/flt-32/e_sqrtf.c index 7648ef4bca..6d6688c52c 100644 --- a/sysdeps/ieee754/flt-32/e_sqrtf.c +++ b/sysdeps/ieee754/flt-32/e_sqrtf.c @@ -8,43 +8,31 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_sqrtf.c,v 1.4 1995/05/10 20:46:19 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float one = 1.0, tiny=1.0e-30; -#else -static float one = 1.0, tiny=1.0e-30; -#endif -#ifdef __STDC__ - float __ieee754_sqrtf(float x) -#else - float __ieee754_sqrtf(x) - float x; -#endif +float +__ieee754_sqrtf(float x) { float z; - int32_t sign = (int)0x80000000; + int32_t sign = (int)0x80000000; int32_t ix,s,q,m,t,i; u_int32_t r; GET_FLOAT_WORD(ix,x); /* take care of Inf and NaN */ - if((ix&0x7f800000)==0x7f800000) { + if((ix&0x7f800000)==0x7f800000) { return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf sqrt(-inf)=sNaN */ - } + } /* take care of zero */ if(ix<=0) { if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */ @@ -69,12 +57,12 @@ static float one = 1.0, tiny=1.0e-30; r = 0x01000000; /* r = moving bit from right to left */ while(r!=0) { - t = s+r; - if(t<=ix) { - s = t+r; - ix -= t; - q += r; - } + t = s+r; + if(t<=ix) { + s = t+r; + ix -= t; + q += r; + } ix += ix; r>>=1; } @@ -83,7 +71,7 @@ static float one = 1.0, tiny=1.0e-30; if(ix!=0) { z = one-tiny; /* trigger inexact flag */ if (z>=one) { - z = one+tiny; + z = one+tiny; if (z>one) q += 2; else @@ -95,3 +83,4 @@ static float one = 1.0, tiny=1.0e-30; SET_FLOAT_WORD(z,ix); return z; } +strong_alias (__ieee754_sqrtf, __sqrtf_finite) diff --git a/sysdeps/ieee754/flt-32/s_asinhf.c b/sysdeps/ieee754/flt-32/s_asinhf.c index fac256d37a..aa46f90004 100644 --- a/sysdeps/ieee754/flt-32/s_asinhf.c +++ b/sysdeps/ieee754/flt-32/s_asinhf.c @@ -13,46 +13,36 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_asinhf.c,v 1.5 1995/05/12 04:57:39 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float -#else -static float -#endif one = 1.0000000000e+00, /* 0x3F800000 */ ln2 = 6.9314718246e-01, /* 0x3f317218 */ huge= 1.0000000000e+30; -#ifdef __STDC__ - float __asinhf(float x) -#else - float __asinhf(x) - float x; -#endif +float +__asinhf(float x) { - float t,w; + float w; int32_t hx,ix; GET_FLOAT_WORD(hx,x); ix = hx&0x7fffffff; - if(ix>=0x7f800000) return x+x; /* x is inf or NaN */ - if(ix< 0x38000000) { /* |x|<2**-14 */ + if(__builtin_expect(ix< 0x38000000, 0)) { /* |x|<2**-14 */ if(huge+x>one) return x; /* return x inexact except 0 */ } - if(ix>0x47000000) { /* |x| > 2**14 */ + if(__builtin_expect(ix>0x47000000, 0)) { /* |x| > 2**14 */ + if(ix>=0x7f800000) return x+x; /* x is inf or NaN */ w = __ieee754_logf(fabsf(x))+ln2; - } else if (ix>0x40000000) { /* 2**14 > |x| > 2.0 */ - t = fabsf(x); - w = __ieee754_logf((float)2.0*t+one/(__ieee754_sqrtf(x*x+one)+t)); - } else { /* 2.0 > |x| > 2**-14 */ - t = x*x; - w =__log1pf(fabsf(x)+t/(one+__ieee754_sqrtf(one+t))); + } else { + float xa = fabsf(x); + if (ix>0x40000000) { /* 2**14 > |x| > 2.0 */ + w = __ieee754_logf(2.0f*xa+one/(__ieee754_sqrtf(xa*xa+one)+xa)); + } else { /* 2.0 > |x| > 2**-14 */ + float t = xa*xa; + w =__log1pf(xa+t/(one+__ieee754_sqrtf(one+t))); + } } - if(hx>0) return w; else return -w; + return __copysignf(w, x); } weak_alias (__asinhf, asinhf) |