From fe72622ebeb251c2dca1aabf9528bcfb8e5a7a9e Mon Sep 17 00:00:00 2001
From: Bruce Evans <bde@FreeBSD.org>
Date: Wed, 5 Jul 2006 22:59:33 +0000
Subject: [PATCH] Fixed tanh(-0.0) on ia64 and optimizeed tanh(x) for 2**-55 <=
 |x| < 2**-28 as a side effect, by merging with the float precision version of
 tanh() and the double precision version of sinh().

For tiny x, tanh(x) ~= x, and we used the expression x*(one+x) to
return this value (x) and set the inexact flag iff x != 0.  This
doesn't work on ia64 since gcc -O does the dubious optimization
x*(one+x) = x+x*x so as to use fma, so the sign of -0.0 was lost.

Instead, handle tiny x in the same as sinh(), although this is imperfect:
- return x directly and set the inexact flag in a less efficient way.
- increased the threshold for non-tinyness from 2**-55 to 2**-28 so that
  many more cases are optimized than are pessimized.

Updated some comments and fixed bugs in others (ranges for half-open
intervals mostly had the open end backwards, and there were nearby style
bugs).
---
 lib/msun/src/s_tanh.c | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/lib/msun/src/s_tanh.c b/lib/msun/src/s_tanh.c
index e3bbf47e27dd..a05ad30e4173 100644
--- a/lib/msun/src/s_tanh.c
+++ b/lib/msun/src/s_tanh.c
@@ -24,14 +24,14 @@ static char rcsid[] = "$FreeBSD$";
  *				       x    -x
  *				      e  + e
  *	1. reduce x to non-negative by tanh(-x) = -tanh(x).
- *	2.  0      <= x <= 2**-55 : tanh(x) := x*(one+x)
+ *	2.  0      <= x <  2**-28 : tanh(x) := x with inexact if x != 0
  *					        -t
- *	    2**-55 <  x <=  1     : tanh(x) := -----; t = expm1(-2x)
+ *	    2**-28 <= x <  1      : tanh(x) := -----; t = expm1(-2x)
  *					       t + 2
  *						     2
- *	    1      <= x <=  22.0  : tanh(x) := 1-  ----- ; t=expm1(2x)
+ *	    1      <= x <  22     : tanh(x) := 1 - -----; t = expm1(2x)
  *						   t + 2
- *	    22.0   <  x <= INF    : tanh(x) := 1.
+ *	    22     <= x <= INF    : tanh(x) := 1.
  *
  * Special cases:
  *	tanh(NaN) is NaN;
@@ -41,7 +41,7 @@ static char rcsid[] = "$FreeBSD$";
 #include "math.h"
 #include "math_private.h"
 
-static const double one=1.0, two=2.0, tiny = 1.0e-300;
+static const double one = 1.0, two = 2.0, tiny = 1.0e-300, huge = 1.0e300;
 
 double
 tanh(double x)
@@ -49,7 +49,6 @@ tanh(double x)
 	double t,z;
 	int32_t jx,ix;
 
-    /* High word of |x|. */
 	GET_HIGH_WORD(jx,x);
 	ix = jx&0x7fffffff;
 
@@ -61,8 +60,9 @@ tanh(double x)
 
     /* |x| < 22 */
 	if (ix < 0x40360000) {		/* |x|<22 */
-	    if (ix<0x3c800000) 		/* |x|<2**-55 */
-		return x*(one+x);    	/* tanh(small) = small */
+	    if (ix<0x3e300000) {	/* |x|<2**-28 */
+		if(huge+x>one) return x; /* tanh(tiny) = tiny with inexact */
+	    }
 	    if (ix>=0x3ff00000) {	/* |x|>=1  */
 		t = expm1(two*fabs(x));
 		z = one - two/(t+two);
@@ -70,9 +70,9 @@ tanh(double x)
 	        t = expm1(-two*fabs(x));
 	        z= -t/(t+two);
 	    }
-    /* |x| > 22, return +-1 */
+    /* |x| >= 22, return +-1 */
 	} else {
-	    z = one - tiny;		/* raised inexact flag */
+	    z = one - tiny;		/* raise inexact flag */
 	}
 	return (jx>=0)? z: -z;
 }