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

Fixed some magic numbers.

Browse Source 
The threshold for not being tiny was too small.  Use the usual 2**-12
threshold.  This change is not just an optimization, since the general
code that we fell into has accuracy problems even for tiny x.  Avoiding
it fixes 2*1366 args with errors of more than 1 ulp, with a maximum
error of 1.167 ulps.

The magic number 22 is log(DBL_EPSILON)/2 plus slop.  This is bogus
for float precision.  Use 9 (~log(FLT_EPSILON)/2 plus less slop than
for double precision).  The code for handling the interval
[2**-28, 9_was_22] has accuracy problems even for [9, 22], so this
change happens to fix errors of more than 1 ulp in about 2*17000
cases.  It leaves such errors in about 2*1074000 cases, with a max
error of 1.242 ulps.

The threshold for switching from returning exp(x)/2 to returning
exp(x/2)^2/2 was a little smaller than necessary.  As for coshf(),
This was not quite harmless since the exp(x/2)^2/2 case is inaccurate,
and fixing it avoids accuracy problems in 2*6 cases, leaving problems
in 2*19997 cases.

Fixed naming errors in pseudo-code in comments.
This commit is contained in:
Bruce Evans 2005-11-13 00:41:46 +00:00
parent c24b7984fc
commit 25efbfb212
1 changed files with 6 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
lib/msun/src/e_sinhf.c
View File

@ -36,19 +36,19 @@ __ieee754_sinhf(float 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 */
/* |x| in [0,9], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x41100000) { /* |x|<9 */
if (ix<0x39800000) /* |x|<2**-12 */
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));
return h*(t+t/(t+one));
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x42b17180) return h*__ieee754_expf(fabsf(x));
/* |x| in [9, logf(maxfloat)] return 0.5*exp(|x|) */
if (ix < 0x42b17217) return h*__ieee754_expf(fabsf(x));
/* |x| in [log(maxdouble), overflowthresold] */
/* |x| in [logf(maxfloat), overflowthresold] */
if (ix<=0x42b2d4fc) {
w = __ieee754_expf((float)0.5*fabsf(x));
t = h*w;