693 Opensource replacement of sunwlibm

Reviewed by: Igor Kozhukhov ikozhukhov@gmail.com
Reviewed by: Keith M Wesolowski <keith.wesolowski@joyent.com>
Reviewed by: Richard Lowe <richlowe@richlowe.net>
Approved by: Dan McDonald <danmcd@omniti.com>

1 files changed, 196 insertions, 0 deletions

diff --git a/usr/src/lib/libm/common/R/atanf.c b/usr/src/lib/libm/common/R/atanf.c
new file mode 100644
index 0000000000..b68c48cca8
--- /dev/null
+++ b/usr/src/lib/libm/common/R/atanf.c
@@ -0,0 +1,196 @@
+/*
+ * CDDL HEADER START
+ *
+ * The contents of this file are subject to the terms of the
+ * Common Development and Distribution License (the "License").
+ * You may not use this file except in compliance with the License.
+ *
+ * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
+ * or http://www.opensolaris.org/os/licensing.
+ * See the License for the specific language governing permissions
+ * and limitations under the License.
+ *
+ * When distributing Covered Code, include this CDDL HEADER in each
+ * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
+ * If applicable, add the following below this CDDL HEADER, with the
+ * fields enclosed by brackets "[]" replaced with your own identifying
+ * information: Portions Copyright [yyyy] [name of copyright owner]
+ *
+ * CDDL HEADER END
+ */
+
+/*
+ * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
+ */
+/*
+ * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
+ * Use is subject to license terms.
+ */
+
+#pragma weak atanf = __atanf
+
+/* INDENT OFF */
+/*
+ * float atanf(float x);
+ * Table look-up algorithm
+ * By K.C. Ng, March 9, 1989
+ *
+ * Algorithm.
+ *
+ * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)).
+ * We use poly1(x) to approximate atan(x) for x in [0,1/8] with
+ * error (relative)
+ * 	|(atan(x)-poly1(x))/x|<= 2^-115.94 	long double
+ * 	|(atan(x)-poly1(x))/x|<= 2^-58.85	double
+ * 	|(atan(x)-poly1(x))/x|<= 2^-25.53 	float
+ * and use poly2(x) to approximate atan(x) for x in [0,1/65] with
+ * error (absolute)
+ *	|atan(x)-poly2(x)|<= 2^-122.15		long double
+ *	|atan(x)-poly2(x)|<= 2^-64.79		double
+ *	|atan(x)-poly2(x)|<= 2^-35.36		float
+ * and use poly3(x) to approximate atan(x) for x in [1/8,7/16] with
+ * error (relative, on for single precision)
+ * 	|(atan(x)-poly1(x))/x|<= 2^-25.53 	float
+ *
+ * Here poly1-3 are odd polynomial with the following form:
+ *		x + x^3*(a1+x^2*(a2+...))
+ *
+ * (0). Purge off Inf and NaN and 0
+ * (1). Reduce x to positive by atan(x) = -atan(-x).
+ * (2). For x <= 1/8, use
+ *	(2.1) if x < 2^(-prec/2-2), atan(x) =  x  with inexact
+ *	(2.2) Otherwise
+ *		atan(x) = poly1(x)
+ * (3). For x >= 8 then
+ *	(3.1) if x >= 2^(prec+2),   atan(x) = atan(inf) - pio2lo
+ *	(3.2) if x >= 2^(prec/3+2), atan(x) = atan(inf) - 1/x
+ *	(3.3) if x >  65,           atan(x) = atan(inf) - poly2(1/x)
+ *	(3.4) Otherwise,	    atan(x) = atan(inf) - poly1(1/x)
+ *
+ * (4). Now x is in (0.125, 8)
+ *      Find y that match x to 4.5 bit after binary (easy).
+ *	If iy is the high word of y, then
+ *		single : j = (iy - 0x3e000000) >> 19
+ *		(single is modified to (iy-0x3f000000)>>19)
+ *		double : j = (iy - 0x3fc00000) >> 16
+ *		quad   : j = (iy - 0x3ffc0000) >> 12
+ *
+ *	Let s = (x-y)/(1+x*y). Then
+ *	atan(x) = atan(y) + poly1(s)
+ *		= _TBL_r_atan_hi[j] + (_TBL_r_atan_lo[j] + poly2(s) )
+ *
+ *	Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
+ *
+ */
+
+#include "libm.h"
+
+extern const float _TBL_r_atan_hi[], _TBL_r_atan_lo[];
+static const float
+	big	=   1.0e37F,
+	one	=   1.0F,
+	p1	=  -3.333185951111688247225368498733544672172e-0001F,
+	p2	=   1.969352894213455405211341983203180636021e-0001F,
+	q1	=  -3.332921964095646819563419704110132937456e-0001F,
+	a1	=  -3.333323465223893614063523351509338934592e-0001F,
+	a2	=   1.999425625935277805494082274808174062403e-0001F,
+	a3	=  -1.417547090509737780085769846290301788559e-0001F,
+	a4	=   1.016250813871991983097273733227432685084e-0001F,
+	a5	=  -5.137023693688358515753093811791755221805e-0002F,
+	pio2hi	=   1.570796371e+0000F,
+	pio2lo	=  -4.371139000e-0008F;
+/* INDENT ON */
+
+float
+atanf(float xx) {
+	float x, y, z, r, p, s;
+	volatile double dummy;
+	int ix, iy, sign, j;
+
+	x = xx;
+	ix = *(int *) &x;
+	sign = ix & 0x80000000;
+	ix ^= sign;
+
+	/* for |x| < 1/8 */
+	if (ix < 0x3e000000) {
+		if (ix < 0x38800000) {	/* if |x| < 2**(-prec/2-2) */
+			dummy = big + x;	/* get inexact flag if x != 0 */
+#ifdef lint
+			dummy = dummy;
+#endif
+			return (x);
+		}
+		z = x * x;
+		if (ix < 0x3c000000) {	/* if |x| < 2**(-prec/4-1) */
+			x = x + (x * z) * p1;
+			return (x);
+		} else {
+			x = x + (x * z) * (p1 + z * p2);
+			return (x);
+		}
+	}
+
+	/* for |x| >= 8.0 */
+	if (ix >= 0x41000000) {
+		*(int *) &x = ix;
+		if (ix < 0x42820000) {	/* x <  65 */
+			r = one / x;
+			z = r * r;
+			y = r * (one + z * (p1 + z * p2));	/* poly1 */
+			y -= pio2lo;
+		} else if (ix < 0x44800000) {	/* x <  2**(prec/3+2) */
+			r = one / x;
+			z = r * r;
+			y = r * (one + z * q1);	/* poly2 */
+			y -= pio2lo;
+		} else if (ix < 0x4c800000) {	/* x <  2**(prec+2) */
+			y = one / x - pio2lo;
+		} else if (ix < 0x7f800000) {	/* x <  inf */
+			y = -pio2lo;
+		} else {		/* x is inf or NaN */
+			if (ix > 0x7f800000) {
+				return (x * x);	/* - -> * for Cheetah */
+			}
+			y = -pio2lo;
+		}
+
+		if (sign == 0)
+			x = pio2hi - y;
+		else
+			x = y - pio2hi;
+		return (x);
+	}
+
+
+	/* now x is between 1/8 and 8 */
+	if (ix < 0x3f000000) {	/* between 1/8 and 1/2 */
+		z = x * x;
+		x = x + (x * z) * (a1 + z * (a2 + z * (a3 + z * (a4 +
+			z * a5))));
+		return (x);
+	}
+	*(int *) &x = ix;
+	iy = (ix + 0x00040000) & 0x7ff80000;
+	*(int *) &y = iy;
+	j = (iy - 0x3f000000) >> 19;
+
+	if (ix == iy)
+		p = x - y;	/* p=0.0 */
+	else {
+		if (sign == 0)
+			s = (x - y) / (one + x * y);
+		else
+			s = (y - x) / (one + x * y);
+		z = s * s;
+		p = s * (one + z * q1);
+	}
+	if (sign == 0) {
+		r = p + _TBL_r_atan_lo[j];
+		x = r + _TBL_r_atan_hi[j];
+	} else {
+		r = p - _TBL_r_atan_lo[j];
+		x = r - _TBL_r_atan_hi[j];
+	}
+	return (x);
+}