From 152a4c4166c18974fe1fc9d9a8535d78956aca26 Mon Sep 17 00:00:00 2001
From: das <das@FreeBSD.org>
Date: Wed, 2 Jun 2004 04:39:44 +0000
Subject: 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
---
 lib/msun/src/k_tanf.c | 32 ++++++++++++++++++++++++--------
 1 file changed, 24 insertions(+), 8 deletions(-)

(limited to 'lib/msun/src')

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;
-- 
cgit v1.1