Module Name: src Committed By: christos Date: Sat Aug 27 10:01:08 UTC 2016
Modified Files: src/lib/libm/src: e_pow.c e_powf.c Log Message: sync with FreeBSD To generate a diff of this commit: cvs rdiff -u -r1.16 -r1.17 src/lib/libm/src/e_pow.c cvs rdiff -u -r1.15 -r1.16 src/lib/libm/src/e_powf.c Please note that diffs are not public domain; they are subject to the copyright notices on the relevant files.
Modified files: Index: src/lib/libm/src/e_pow.c diff -u src/lib/libm/src/e_pow.c:1.16 src/lib/libm/src/e_pow.c:1.17 --- src/lib/libm/src/e_pow.c:1.16 Fri Apr 23 15:17:07 2010 +++ src/lib/libm/src/e_pow.c Sat Aug 27 06:01:08 2016 @@ -1,9 +1,8 @@ -/* @(#)e_pow.c 5.1 93/09/24 */ +/* @(#)e_pow.c 1.5 04/04/22 SMI */ /* * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * Copyright (C) 2004 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. @@ -12,7 +11,7 @@ #include <sys/cdefs.h> #if defined(LIBM_SCCS) && !defined(lint) -__RCSID("$NetBSD: e_pow.c,v 1.16 2010/04/23 19:17:07 drochner Exp $"); +__RCSID("$NetBSD: e_pow.c,v 1.17 2016/08/27 10:01:08 christos Exp $"); #endif /* __ieee754_pow(x,y) return x**y @@ -29,13 +28,13 @@ __RCSID("$NetBSD: e_pow.c,v 1.16 2010/04 * Special cases: * 1. (anything) ** 0 is 1 * 2. (anything) ** 1 is itself - * 3. (anything) ** NAN is NAN + * 3. (anything) ** NAN is NAN except 1 ** NAN = 1 * 4. NAN ** (anything except 0) is NAN * 5. +-(|x| > 1) ** +INF is +INF * 6. +-(|x| > 1) ** -INF is +0 * 7. +-(|x| < 1) ** +INF is +0 * 8. +-(|x| < 1) ** -INF is +INF - * 9. +-1 ** +-INF is NAN + * 9. +-1 ** +-INF is 1 * 10. +0 ** (+anything except 0, NAN) is +0 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 * 12. +0 ** (-anything except 0, NAN) is +INF @@ -113,10 +112,13 @@ __ieee754_pow(double x, double y) /* y==zero: x**0 = 1 */ if((iy|ly)==0) return one; - /* +-NaN return x+y */ + /* x==1: 1**y = 1, even if y is NaN */ + if (hx==0x3ff00000 && lx == 0) return one; + + /* y!=zero: result is NaN if either arg is NaN */ if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) || iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) - return x+y; + return (x+0.0)+(y+0.0); /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer @@ -142,7 +144,7 @@ __ieee754_pow(double x, double y) if(ly==0) { if (iy==0x7ff00000) { /* y is +-inf */ if(((ix-0x3ff00000)|lx)==0) - return y - y; /* inf**+-1 is NaN */ + return one; /* (-1)**+-inf is 1 */ else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */ return (hy>=0)? y: zero; else /* (|x|<1)**-,+inf = inf,0 */ @@ -174,7 +176,11 @@ __ieee754_pow(double x, double y) } } + /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be n = (hx>>31)+1; + but ANSI C says a right shift of a signed negative quantity is + implementation defined. */ + n = ((u_int32_t)hx>>31)-1; /* (x<0)**(non-int) is NaN */ if((n|yisint)==0) return (x-x)/(x-x); Index: src/lib/libm/src/e_powf.c diff -u src/lib/libm/src/e_powf.c:1.15 src/lib/libm/src/e_powf.c:1.16 --- src/lib/libm/src/e_powf.c:1.15 Fri Apr 23 15:17:07 2010 +++ src/lib/libm/src/e_powf.c Sat Aug 27 06:01:08 2016 @@ -15,15 +15,13 @@ #include <sys/cdefs.h> #if defined(LIBM_SCCS) && !defined(lint) -__RCSID("$NetBSD: e_powf.c,v 1.15 2010/04/23 19:17:07 drochner Exp $"); +__RCSID("$NetBSD: e_powf.c,v 1.16 2016/08/27 10:01:08 christos Exp $"); #endif #include "namespace.h" #include "math.h" #include "math_private.h" -static const float huge = 1.0e+30, tiny = 1.0e-30; - static const float bp[] = {1.0, 1.5,}, dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */ @@ -32,6 +30,8 @@ zero = 0.0, one = 1.0, two = 2.0, two24 = 16777216.0, /* 0x4b800000 */ +huge = 1.0e30, +tiny = 1.0e-30, /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ L1 = 6.0000002384e-01, /* 0x3f19999a */ L2 = 4.2857143283e-01, /* 0x3edb6db7 */ @@ -49,8 +49,8 @@ lg2_h = 6.93145752e-01, /* 0x3f317200 lg2_l = 1.42860654e-06, /* 0x35bfbe8c */ ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */ cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */ -cp_h = 9.6179199219e-01, /* 0x3f763800 =head of cp */ -cp_l = 4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */ +cp_h = 9.6191406250e-01, /* 0x3f764000 =12b cp */ +cp_l = -1.1736857402e-04, /* 0xb8f623c6 =tail of cp_h */ 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*/ @@ -59,7 +59,7 @@ float __ieee754_powf(float x, float y) { float z,ax,z_h,z_l,p_h,p_l; - float yy1,t1,t2,r,s,t,u,v,w; + float yy1,t1,t2,r,s,sn,t,u,v,w; int32_t i,j,k,yisint,n; int32_t hx,hy,ix,iy,is; @@ -70,10 +70,13 @@ __ieee754_powf(float x, float y) /* y==zero: x**0 = 1 */ if(iy==0) return one; - /* +-NaN return x+y */ + /* x==1: 1**y = 1, even if y is NaN */ + if (hx==0x3f800000) return one; + + /* y!=zero: result is NaN if either arg is NaN */ if(ix > 0x7f800000 || iy > 0x7f800000) - return x+y; + return (x+0.0F)+(y+0.0F); /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer @@ -84,8 +87,8 @@ __ieee754_powf(float x, float y) if(hx<0) { if(iy>=0x4b800000) yisint = 2; /* even integer y */ else if(iy>=0x3f800000) { - k = (iy>>23)-0x7f; /* exponent */ - j = iy>>(23-k); + k = ((uint32_t)iy>>23)-0x7f; /* exponent */ + j = (uint32_t)iy>>(23-k); if((j<<(23-k))==iy) yisint = 2-(j&1); } } @@ -93,7 +96,7 @@ __ieee754_powf(float x, float y) /* special value of y */ if (iy==0x7f800000) { /* y is +-inf */ if (ix==0x3f800000) - return y - y; /* inf**+-1 is NaN */ + return one; /* (-1)**+-inf is NaN */ else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */ return (hy>=0)? y: zero; else /* (|x|<1)**-,+inf = inf,0 */ @@ -122,17 +125,22 @@ __ieee754_powf(float x, float y) return z; } + n = ((u_int32_t)hx>>31)-1; + /* (x<0)**(non-int) is NaN */ - if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x); + if((n|yisint)==0) return (x-x)/(x-x); + + sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */ + if((n|(yisint-1))==0) sn = -one;/* (-ve)**(odd int) */ /* |y| is huge */ if(iy>0x4d000000) { /* 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; + if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny; + if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny; /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ - t = ax-one; /* t has 20 trailing zeros */ + t = ax-1; /* t has 20 trailing zeros */ w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25)); u = ivln2_h*t; /* ivln2_h has 16 sig. bits */ v = t*ivln2_l-w*ivln2; @@ -146,7 +154,7 @@ __ieee754_powf(float x, float y) /* take care subnormal number */ if(ix<0x00800000) {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); } - n += ((ix)>>23)-0x7f; + n += (((uint32_t)ix)>>23)-0x7f; j = ix&0x007fffff; /* determine interval */ ix = j|0x3f800000; /* normalize ix */ @@ -163,7 +171,8 @@ __ieee754_powf(float x, float y) GET_FLOAT_WORD(is,s_h); SET_FLOAT_WORD(s_h,is&0xfffff000); /* t_h=ax+bp[k] High */ - SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x0040000+(k<<21)); + is = (((uint32_t)ix>>1)&0xfffff000)|0x20000000; + SET_FLOAT_WORD(t_h,is+0x00400000+(k<<21)); t_l = ax - (t_h-bp[k]); s_l = v*((u-s_h*t_h)-s_h*t_l); /* compute log(ax) */ @@ -193,10 +202,6 @@ __ieee754_powf(float x, float y) t2 = z_l-(((t1-t)-dp_h[k])-z_h); } - s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ - if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0) - s = -one; /* (-ve)**(odd int) */ - /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */ GET_FLOAT_WORD(is,y); SET_FLOAT_WORD(yy1,is&0xfffff000); @@ -205,32 +210,32 @@ __ieee754_powf(float x, float y) z = p_l+p_h; GET_FLOAT_WORD(j,z); if (j>0x43000000) /* if z > 128 */ - return s*huge*huge; /* overflow */ + return sn*huge*huge; /* overflow */ else if (j==0x43000000) { /* if z == 128 */ - if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ - } - else if ((uint32_t)j==0xc3160000){ /* z == -150 */ - if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ + if(p_l+ovt>z-p_h) return sn*huge*huge; /* overflow */ } else if ((j&0x7fffffff)>0x43160000) /* z <= -150 */ - return s*tiny*tiny; /* underflow */ + return sn*tiny*tiny; /* underflow */ + else if (j==(int)0xc3160000){ /* z == -150 */ + if(p_l<=z-p_h) return sn*tiny*tiny; /* underflow */ + } /* * compute 2**(p_h+p_l) */ i = j&0x7fffffff; - k = (i>>23)-0x7f; + k = ((uint32_t)i>>23)-0x7f; n = 0; if(i>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */ n = j+(0x00800000>>(k+1)); - k = ((n&0x7fffffff)>>23)-0x7f; /* new k for n */ + k = (((uint32_t)n&0x7fffffff)>>23)-0x7f; /* new k for n */ SET_FLOAT_WORD(t,n&~(0x007fffff>>k)); - n = ((n&0x007fffff)|0x00800000)>>(23-k); + n = (((uint32_t)n&0x007fffff)|0x00800000)>>(23-k); if(j<0) n = -n; p_h -= t; } t = p_l+p_h; GET_FLOAT_WORD(is,t); - SET_FLOAT_WORD(t,is&0xfffff000); + SET_FLOAT_WORD(t,is&0xffff8000); u = t*lg2_h; v = (p_l-(t-p_h))*lg2+t*lg2_l; z = u+v; @@ -243,5 +248,5 @@ __ieee754_powf(float x, float y) j += (n<<23); if((j>>23)<=0) z = scalbnf(z,n); /* subnormal output */ else SET_FLOAT_WORD(z,j); - return s*z; + return sn*z; }