1 files changed, 125 insertions, 0 deletions
diff --git a/private/fp32/tran/tan.c b/private/fp32/tran/tan.c
new file mode 100644
index 000000000..2100595a7
--- /dev/null
+++ b/private/fp32/tran/tan.c
@@ -0,0 +1,125 @@
+/***
+*tan.c - tangent
+*
+*	Copyright (c) 1991-1991, Microsoft Corporation.	All rights reserved.
+*
+*Purpose:
+*
+*Revision History:
+*   8-15-91	GDP	written
+*  12-30-91	GDP	support IEEE exceptions
+*  03-11-91	GDP	use 66 significant bits for representing pi
+*			support FP_TLOSS
+*******************************************************************************/
+#include <math.h>
+#include <trans.h>
+
+/* constants */
+static double const TWO_OVER_PI	= 0.63661977236758134308;
+static double const EPS  = 1.05367121277235079465e-8; /* 2^(-53/2) */
+static double const YMAX = 2.98156826864790199324e8; /* 2^(53/2)*PI/2 */
+
+//
+// The sum of C1 and C2 is a representation of PI/2 with 66 bits in the
+// significand (same as x87). (PI/2 = 2 * 0.c90fdaa2 2168c234 c h)
+//
+
+static _dbl _C1	 = {SET_DBL (0x3ff921fb, 0x54400000)};
+static _dbl _C2	 = {SET_DBL (0x3dd0b461, 0x1a600000)};
+#define C1  (_C1.dbl)
+#define C2  (_C2.dbl)
+
+/* constants for the rational approximation */
+/* p0 = 1.0  is not used (avoid mult by 1) */
+static double const p1 = -0.13338350006421960681e+0;
+static double const p2 =  0.34248878235890589960e-2;
+static double const p3 = -0.17861707342254426711e-4;
+static double const q0 =  0.10000000000000000000e+1;
+static double const q1 = -0.46671683339755294240e+0;
+static double const q2 =  0.25663832289440112864e-1;
+static double const q3 = -0.31181531907010027307e-3;
+static double const q4 =  0.49819433993786512270e-6;
+
+
+#define Q(g)   ((((q4 * (g) + q3) * (g) + q2) * (g) + q1) * (g) + q0)
+#define P(g,f)	(((p3 * (g) + p2) * (g) + p1) * (g) * (f) + (f))
+
+#define ISODD(i) ((i)&0x1)
+
+
+/***
+*double tan(double x) - tangent
+*
+*Purpose:
+*   Compute the tangent of a number.
+*   The algorithm (reduction / rational approximation) is
+*   taken from Cody & Waite.
+*
+*Entry:
+*
+*Exit:
+*
+*Exceptions:
+*   P, I
+*   no exception if x is denormal: return x
+*******************************************************************************/
+double tan(double x)
+{
+    unsigned int savedcw;
+    unsigned long n;
+    double xn,xnum,xden;
+    double f,g,result;
+
+    /* save user fp control word */
+    savedcw = _maskfp();
+
+    if (IS_D_SPECIAL(x)){
+	switch(_sptype(x)) {
+	case T_PINF:
+	case T_NINF:
+	    return _except1(FP_I,OP_TAN,x,QNAN_TAN1,savedcw);
+	case T_QNAN:
+	    return _handle_qnan1(OP_TAN, x, savedcw);
+	default: //T_SNAN
+	    return _except1(FP_I,OP_TAN,x,_s2qnan(x),savedcw);
+	}
+    }
+
+    if (x == 0.0)
+	RETURN(savedcw, x);
+
+    if (ABS(x) > YMAX) {
+
+	// The argument is too large to produce a meaningful result,
+	// so this is treated as an invalid operation.
+	// We also set the (extra) FP_TLOSS flag for matherr
+	// support
+
+	return _except1(FP_TLOSS | FP_I,OP_TAN,x,QNAN_TAN2,savedcw);
+    }
+
+    xn = _frnd(x * TWO_OVER_PI);
+    n = (unsigned long) xn;
+
+
+    /* assume there is a guard digit for addition */
+    f = (x - xn * C1) - xn * C2;
+    if (ABS(f) < EPS) {
+	xnum = f;
+	xden = 1;
+    }
+    else {
+	g = f*f;
+	xnum = P(g,f);
+	xden = Q(g);
+    }
+
+    if (ISODD(n)) {
+	xnum = -xnum;
+	result = xden/xnum;
+    }
+    else
+	result = xnum/xden;
+
+    RETURN_INEXACT1(OP_TAN,x,result,savedcw);
+}