diff options
author | bde <bde@FreeBSD.org> | 2005-11-13 00:41:46 +0000 |
---|---|---|
committer | bde <bde@FreeBSD.org> | 2005-11-13 00:41:46 +0000 |
commit | f63f109c0bfa1b21d8d5360a6bcf0e0f2e86b7ee (patch) | |
tree | 4798349d4b1e2b7dc1da97fcd15a875421f24073 | |
parent | 3f7e4f15381a1bfc788cd99f2dcc2cb8a8f5b460 (diff) | |
download | FreeBSD-src-f63f109c0bfa1b21d8d5360a6bcf0e0f2e86b7ee.zip FreeBSD-src-f63f109c0bfa1b21d8d5360a6bcf0e0f2e86b7ee.tar.gz |
Fixed some magic numbers.
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.
-rw-r--r-- | lib/msun/src/e_sinhf.c | 12 |
1 files changed, 6 insertions, 6 deletions
diff --git a/lib/msun/src/e_sinhf.c b/lib/msun/src/e_sinhf.c index 7e022b0..50388fa 100644 --- a/lib/msun/src/e_sinhf.c +++ b/lib/msun/src/e_sinhf.c @@ -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; |