Fixed some magic numbers.

The threshold for not being tiny was too small.  Use the usual 2**-12
threshold.  As for sinhf, use a different method (now the same as for
sinhf) to set the inexact flag for tiny nonzero x so that the larger
threshold works, although this method is imperfect.  As for sinhf,
this change is not just an optimization, since the general code that
we fell into has accuracy problems even for tiny x.  On amd64, avoiding
it fixes tanhf on 2*13495596 args with errors of between 1 and 1.3
ulps and thus reduces the total number of args with errors of >= 1 ulp
from 37533748 to 5271278; the maximum error is unchanged at 2.2 ulps.

The magic number 22 is log(DBL_MAX)/2 plus slop.  This is bogus for
float precision.  Use 9 (log(FLT_MAX)/2 plus less slop than for
double precision).  Unlike for coshf and tanhf, this is just an
optimization, and MAX isn't misspelled EPSILON in the commit log.

I started testing with nonstandard rounding modes, and verified that
the chosen thresholds work for all modes modulo problems not related
to thresholds.  The best thresholds are not very dependent on the mode,
at least for tanhf.

1 files changed, 8 insertions, 8 deletions

diff --git a/lib/msun/src/s_tanhf.c b/lib/msun/src/s_tanhf.c
index 74cb2db..8f062c8 100644
--- a/lib/msun/src/s_tanhf.c
+++ b/lib/msun/src/s_tanhf.c
@@ -20,8 +20,7 @@ static char rcsid[] = "$FreeBSD$";
 #include "math.h"
 #include "math_private.h"
 
-static const float one=1.0, two=2.0, tiny = 1.0e-30;
-
+static const float one=1.0, two=2.0, tiny = 1.0e-30, huge = 1.0e30;
 float
 tanhf(float x)
 {
@@ -37,10 +36,11 @@ tanhf(float x)
 	    else       return one/x-one;    /* tanh(NaN) = NaN */
 	}
 
-    /* |x| < 22 */
-	if (ix < 0x41b00000) {		/* |x|<22 */
-	    if (ix<0x24000000) 		/* |x|<2**-55 */
-		return x*(one+x);    	/* tanh(small) = small */
+    /* |x| < 9 */
+	if (ix < 0x41100000) {		/* |x|<9 */
+	    if (ix<0x39800000) {	/* |x|<2**-12 */
+		if(huge+x>one) return x; /* tanh(tiny) = tiny with inexact */
+	    }
 	    if (ix>=0x3f800000) {	/* |x|>=1  */
 		t = expm1f(two*fabsf(x));
 		z = one - two/(t+two);
@@ -48,9 +48,9 @@ tanhf(float x)
 	        t = expm1f(-two*fabsf(x));
 	        z= -t/(t+two);
 	    }
-    /* |x| > 22, return +-1 */
+    /* |x| >= 9, return +-1 */
 	} else {
-	    z = one - tiny;		/* raised inexact flag */
+	    z = one - tiny;		/* raise inexact flag */
 	}
 	return (jx>=0)? z: -z;
 }