Moved the optimization for tiny x from __kernel_tan[f](x) to tan[f](x)

so that it can be faster for tiny x and avoided for reduced x.

This improves things a little differently than for cosine and sine.
We still need to reclassify x in the "kernel" functions, but we get
an extra optimization for tiny x, and an overall optimization since
tiny reduced x rarely happens.  We also get optimizations for space
and style.  A large block of poorly duplicated code to fix a special
case is no longer needed.  This supersedes the fixes in k_sin.c revs
1.9 and 1.11 and k_sinf.c 1.8 and 1.10.

Fixed wrong constant for the cutoff for "tiny" in tanf().  It was
2**-28, but should be almost the same as the cutoff in sinf() (2**-12).
The incorrect cutoff protected us from the bugs fixed in k_sinf.c 1.8
and 1.10, except 4 cases of reduced args passed the cutoff and needed
special handling in theory although not in practice.  Now we essentially
use a cutoff of 0 for the case of reduced args, so we now have 0 special
args instead of 4.

This change makes no difference to the results for sinf() (since it
only changes the algorithm for the 4 special args and the results for
those happen not to change), but it changes lots of results for sin().
Exhaustive testing is impossible for sin(), but exhaustive testing
for sinf() (relative to a version with the old algorithm and a fixed
cutoff) shows that the changes in the error are either reductions or
from 0.5-epsilon ulps to 0.5+epsilon ulps.  The new method just uses
some extra terms in approximations so it tends to give more accurate
results, and there are apparently no problems from having extra
accuracy.  On amd64 with -O1, on all float args the error range in ulps
is reduced from (0.500, 0.665] to [0.335, 0.500) in 24168 cases and
increased from 0.500-epsilon to 0.500+epsilon in 24 cases.  Non-
exhaustive testing by ucbtest shows no differences.

4 files changed, 14 insertions, 46 deletions

diff --git a/lib/msun/src/k_tan.c b/lib/msun/src/k_tan.c
index 15145a6..13737d0 100644
--- a/lib/msun/src/k_tan.c
+++ b/lib/msun/src/k_tan.c
@@ -16,14 +16,16 @@ static char rcsid[] = "$FreeBSD$";
 #endif
 
 /* __kernel_tan( x, y, k )
- * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
+ * kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
  * Input x is assumed to be bounded by ~pi/4 in magnitude.
  * Input y is the tail of x.
  * Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned.
  *
  * Algorithm
  *	1. Since tan(-x) = -tan(x), we need only to consider positive x.
- *	2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
+ *	2. Callers must return tan(-0) = -0 without calling here since our
+ *	   odd polynomial is not evaluated in a way that preserves -0.
+ *	   Callers may do the optimization tan(x) ~ x for tiny x.
  *	3. tan(x) is approximated by a odd polynomial of degree 27 on
  *	   [0,0.67434]
  *		  	         3             27
@@ -81,27 +83,6 @@ __kernel_tan(double x, double y, int iy) {
 
 	GET_HIGH_WORD(hx,x);
 	ix = hx & 0x7fffffff;			/* high word of |x| */
-	if (ix < 0x3e300000) {			/* x < 2**-28 */
-		if ((int) x == 0) {		/* generate inexact */
-			u_int32_t low;
-			GET_LOW_WORD(low,x);
-			{
-				if (iy == 1)
-					return x;
-				else {	/* compute -1 / (x+y) carefully */
-					double a, t;
-
-					z = w = x + y;
-					SET_LOW_WORD(z, 0);
-					v = y - (z - x);
-					t = a = -one / w;
-					SET_LOW_WORD(t, 0);
-					s = one + t * z;
-					return t + a * (s + t * v);
-				}
-			}
-		}
-	}
 	if (ix >= 0x3FE59428) {	/* |x| >= 0.6744 */
 		if (hx < 0) {
 			x = -x;
diff --git a/lib/msun/src/k_tanf.c b/lib/msun/src/k_tanf.c
index e344787..8e19309 100644
--- a/lib/msun/src/k_tanf.c
+++ b/lib/msun/src/k_tanf.c
@@ -45,27 +45,6 @@ __kernel_tanf(float x, float y, int iy)
 	int32_t ix,hx;
 	GET_FLOAT_WORD(hx,x);
 	ix = hx&0x7fffffff;	/* high word of |x| */
-	if(ix<0x31800000) {			/* x < 2**-28 */
-		if ((int) x == 0) {		/* generate inexact */
-			{
-				if (iy == 1)
-					return x;
-				else {	/* compute -1 / (x+y) carefully */
-					float a, t;
-
-					z = w = x + y;
-					GET_FLOAT_WORD(ix, z);
-					SET_FLOAT_WORD(z, ix & 0xfffff000);
-					v = y - (z - x);
-					t = a = -one / w;
-					GET_FLOAT_WORD(ix, t);
-					SET_FLOAT_WORD(t, ix & 0xfffff000);
-					s = one + t * z;
-					return t + a * (s + t * v);
-				}
-			}
-		}
-	}
 	if(ix>=0x3f2ca140) { 			/* |x|>=0.6744 */
 	    if(hx<0) {x = -x; y = -y;}
 	    z = pio4-x;
diff --git a/lib/msun/src/s_tan.c b/lib/msun/src/s_tan.c
index 5094d30..d707ecd 100644
--- a/lib/msun/src/s_tan.c
+++ b/lib/msun/src/s_tan.c
@@ -58,7 +58,11 @@ tan(double x)
 
     /* |x| ~< pi/4 */
 	ix &= 0x7fffffff;
-	if(ix <= 0x3fe921fb) return __kernel_tan(x,z,1);
+	if(ix <= 0x3fe921fb) {
+	    if(ix<0x3e300000)			/* x < 2**-28 */
+		if((int)x==0) return x;		/* generate inexact */
+	    return __kernel_tan(x,z,1);
+	}
 
     /* tan(Inf or NaN) is NaN */
 	else if (ix>=0x7ff00000) return x-x;		/* NaN */
diff --git a/lib/msun/src/s_tanf.c b/lib/msun/src/s_tanf.c
index c81caaf..81c84e4 100644
--- a/lib/msun/src/s_tanf.c
+++ b/lib/msun/src/s_tanf.c
@@ -30,7 +30,11 @@ tanf(float x)
 
     /* |x| ~< pi/4 */
 	ix &= 0x7fffffff;
-	if(ix <= 0x3f490fda) return __kernel_tanf(x,z,1);
+	if(ix <= 0x3f490fda) {
+	    if(ix<0x39800000)			/* |x| < 2**-12 */
+		if(((int)x)==0) return x;	/* generate inexact */
+	    return __kernel_tanf(x,z,1);
+	}
 
     /* tan(Inf or NaN) is NaN */
 	else if (ix>=0x7f800000) return x-x;		/* NaN */