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.

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;