Port a bugfix from FDLIBM 5.3. The bug really only applies to tan()

and not tanf() because float type can't represent numbers large enough
to trigger the problem.  However, there seems to be a precedent that
the float versions of the fdlibm routines should mirror their double
counterparts.

Also update to the FDLIBM 5.3 license.

Obtained from:	FDLIBM
Reviewed by:	exhaustive comparison

1 files changed, 24 insertions, 8 deletions

diff --git a/lib/msun/src/k_tanf.c b/lib/msun/src/k_tanf.c
index 35bbe61..41c522c 100644
--- a/lib/msun/src/k_tanf.c
+++ b/lib/msun/src/k_tanf.c
@@ -4,9 +4,8 @@
 
 /*
  * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
  * software is freely granted, provided that this notice
  * is preserved.
@@ -46,12 +45,29 @@ __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((ix|(iy+1))==0) return one/fabsf(x);
-		else return (iy==1)? x: -one/x;
-	    }
-	    }
+	if(ix<0x31800000) {			/* x < 2**-28 */
+		if ((int) x == 0) {		/* generate inexact */
+			if ((ix | (iy + 1)) == 0)
+				return one / fabsf(x);
+			else {
+				if (iy == 1)
+					return x;
+				else {	/* compute -1 / (x+y) carefully */
+					double 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;