d/freebsd-skq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fixed the sign of the result in some overflow and underflow cases (ones

Browse Source 
where the exponent is an odd integer and the base is negative).

Obtained from:	fdlibm-5.3

Sun finally released a new version of fdlibm just a coupe of weeks
ago.  It only fixes 3 bugs (this one, another one in pow() that we
already have (rev.1.9), and one in tan().  I've learned too much about
powf() lately, so this fix was easy to merge.  The patch is not verbatim,
because our base version has many differences for portability and I
didn't like global renaming of an unrelated variable to keep it separate
from the sign variable.  This patch uses a new variable named sn for
the sign.
This commit is contained in:
bde 2004-06-01 19:28:38 +00:00
parent 719aa077cb
commit 152787c7f1
1 changed files with 18 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
lib/msun/src/e_pow.c
View File

@ -99,7 +99,7 @@ double
__ieee754_pow(double x, double y)
{
double z,ax,z_h,z_l,p_h,p_l;
double y1,t1,t2,r,s,t,u,v,w;
double y1,t1,t2,r,s,sn,t,u,v,w;
int32_t i,j,k,yisint,n;
int32_t hx,hy,ix,iy;
u_int32_t lx,ly;
@ -172,12 +172,17 @@ __ieee754_pow(double x, double y)
}
}
/* (x<0)**(non-int) is NaN */
/* CYGNUS LOCAL: This used to be
if((((hx>>31)+1)|yisint)==0) return (x-x)/(x-x);
/* 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. */
if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
n = ((u_int32_t)hx>>31)-1;
/* (x<0)**(non-int) is NaN */
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>0x41e00000) { /* if |y| > 2**31 */
@ -186,8 +191,8 @@ __ieee754_pow(double x, double y)
if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
}
/* over/underflow if x is not close to one */
if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
if(ix<0x3fefffff) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
if(ix>0x3ff00000) 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-1; /* t has 20 trailing zeros */
@ -247,10 +252,6 @@ __ieee754_pow(double x, double 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 y1+y2 and compute (y1+y2)*(t1+t2) */
y1 = y;
SET_LOW_WORD(y1,0);
@ -260,15 +261,15 @@ __ieee754_pow(double x, double y)
EXTRACT_WORDS(j,i,z);
if (j>=0x40900000) { /* z >= 1024 */
if(((j-0x40900000)|i)!=0) /* if z > 1024 */
return s*huge*huge; /* overflow */
return sn*huge*huge; /* overflow */
else {
if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
if(p_l+ovt>z-p_h) return sn*huge*huge; /* overflow */
}
} else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
return s*tiny*tiny; /* underflow */
if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
return sn*tiny*tiny; /* underflow */
else {
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
if(p_l<=z-p_h) return sn*tiny*tiny; /* underflow */
}
}
/*
@ -300,5 +301,5 @@ __ieee754_pow(double x, double y)
j += (n<<20);
if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
else SET_HIGH_WORD(z,j);
return s*z;
return sn*z;
}