uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 376 insertions, 0 deletions

diff --git a/libm/ldouble/atanl.c b/libm/ldouble/atanl.c
new file mode 100644
index 000000000..9e6d9af3c
--- /dev/null
+++ b/libm/ldouble/atanl.c
@@ -0,0 +1,376 @@
+/*							atanl.c
+ *
+ *	Inverse circular tangent, long double precision
+ *      (arctangent)
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, atanl();
+ *
+ * y = atanl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns radian angle between -pi/2 and +pi/2 whose tangent
+ * is x.
+ *
+ * Range reduction is from four intervals into the interval
+ * from zero to  tan( pi/8 ).  The approximant uses a rational
+ * function of degree 3/4 of the form x + x**3 P(x)/Q(x).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10, 10    150000       1.3e-19     3.0e-20
+ *
+ */
+/*							atan2l()
+ *
+ *	Quadrant correct inverse circular tangent,
+ *	long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, z, atan2l();
+ *
+ * z = atan2l( y, x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns radian angle whose tangent is y/x.
+ * Define compile time symbol ANSIC = 1 for ANSI standard,
+ * range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range
+ * 0 to 2PI, args (x,y).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10, 10     60000       1.7e-19     3.2e-20
+ * See atan.c.
+ *
+ */
+
+/*							atan.c */
+
+
+/*
+Cephes Math Library Release 2.7:  May, 1998
+Copyright 1984, 1990, 1998 by Stephen L. Moshier
+*/
+
+
+#include <math.h>
+
+#ifdef UNK
+static long double P[] = {
+-8.6863818178092187535440E-1L,
+-1.4683508633175792446076E1L,
+-6.3976888655834347413154E1L,
+-9.9988763777265819915721E1L,
+-5.0894116899623603312185E1L,
+};
+static long double Q[] = {
+/* 1.00000000000000000000E0L,*/
+ 2.2981886733594175366172E1L,
+ 1.4399096122250781605352E2L,
+ 3.6144079386152023162701E2L,
+ 3.9157570175111990631099E2L,
+ 1.5268235069887081006606E2L,
+};
+
+/* tan( 3*pi/8 ) */
+static long double T3P8 =  2.41421356237309504880169L;
+
+/* tan( pi/8 ) */
+static long double TP8 =  4.1421356237309504880169e-1L;
+#endif
+
+
+#ifdef IBMPC
+static unsigned short P[] = {
+0x8ece,0xce53,0x1266,0xde5f,0xbffe, XPD
+0x07e6,0xa061,0xa6bf,0xeaef,0xc002, XPD
+0x53ee,0xf291,0x557f,0xffe8,0xc004, XPD
+0xf9d6,0xeda6,0x3f3e,0xc7fa,0xc005, XPD
+0xb6c3,0x6abc,0x9361,0xcb93,0xc004, XPD
+};
+static unsigned short Q[] = {
+/*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
+0x54d4,0x894e,0xe76e,0xb7da,0x4003, XPD
+0x76b9,0x7a46,0xafa2,0x8ffd,0x4006, XPD
+0xe3a9,0xe9c0,0x6bee,0xb4b8,0x4007, XPD
+0xabc1,0x50a7,0xb098,0xc3c9,0x4007, XPD
+0x891c,0x100d,0xae89,0x98ae,0x4006, XPD
+};
+
+/* tan( 3*pi/8 ) = 2.41421356237309504880 */
+static unsigned short T3P8A[] = {0x3242,0xfcef,0x7999,0x9a82,0x4000, XPD};
+#define T3P8 *(long double *)T3P8A
+
+/* tan( pi/8 ) = 0.41421356237309504880 */
+static unsigned short TP8A[] = {0x9211,0xe779,0xcccf,0xd413,0x3ffd, XPD};
+#define TP8 *(long double *)TP8A
+#endif
+
+#ifdef MIEEE
+static unsigned long P[] = {
+0xbffe0000,0xde5f1266,0xce538ece,
+0xc0020000,0xeaefa6bf,0xa06107e6,
+0xc0040000,0xffe8557f,0xf29153ee,
+0xc0050000,0xc7fa3f3e,0xeda6f9d6,
+0xc0040000,0xcb939361,0x6abcb6c3,
+};
+static unsigned long Q[] = {
+/*0x3fff0000,0x80000000,0x00000000,*/
+0x40030000,0xb7dae76e,0x894e54d4,
+0x40060000,0x8ffdafa2,0x7a4676b9,
+0x40070000,0xb4b86bee,0xe9c0e3a9,
+0x40070000,0xc3c9b098,0x50a7abc1,
+0x40060000,0x98aeae89,0x100d891c,
+};
+
+/* tan( 3*pi/8 ) = 2.41421356237309504880 */
+static long T3P8A[] = {0x40000000,0x9a827999,0xfcef3242};
+#define T3P8 *(long double *)T3P8A
+
+/* tan( pi/8 ) = 0.41421356237309504880 */
+static long TP8A[] = {0x3ffd0000,0xd413cccf,0xe7799211};
+#define TP8 *(long double *)TP8A
+#endif
+
+#ifdef ANSIPROT
+extern long double polevll ( long double, void *, int );
+extern long double p1evll ( long double, void *, int );
+extern long double fabsl ( long double );
+extern int signbitl ( long double );
+extern int isnanl ( long double );
+long double atanl ( long double );
+#else
+long double polevll(), p1evll(), fabsl(), signbitl(), isnanl();
+long double atanl();
+#endif
+#ifdef INFINITIES
+extern long double INFINITYL;
+#endif
+#ifdef NANS
+extern long double NANL;
+#endif
+#ifdef MINUSZERO
+extern long double NEGZEROL;
+#endif
+
+long double atanl(x)
+long double x;
+{
+extern long double PIO2L, PIO4L;
+long double y, z;
+short sign;
+
+#ifdef MINUSZERO
+if( x == 0.0L )
+	return(x);
+#endif
+#ifdef INFINITIES
+if( x == INFINITYL )
+	return( PIO2L );
+if( x == -INFINITYL )
+	return( -PIO2L );
+#endif
+/* make argument positive and save the sign */
+sign = 1;
+if( x < 0.0L )
+	{
+	sign = -1;
+	x = -x;
+	}
+
+/* range reduction */
+if( x > T3P8 )
+	{
+	y = PIO2L;
+	x = -( 1.0L/x );
+	}
+
+else if( x > TP8 )
+	{
+	y = PIO4L;
+	x = (x-1.0L)/(x+1.0L);
+	}
+else
+	y = 0.0L;
+
+/* rational form in x**2 */
+z = x * x;
+y = y + ( polevll( z, P, 4 ) / p1evll( z, Q, 5 ) ) * z * x + x;
+
+if( sign < 0 )
+	y = -y;
+
+return(y);
+}
+
+/*							atan2	*/
+
+
+extern long double PIL, PIO2L, MAXNUML;
+
+#if ANSIC
+long double atan2l( y, x )
+#else
+long double atan2l( x, y )
+#endif
+long double x, y;
+{
+long double z, w;
+short code;
+
+code = 0;
+
+if( x < 0.0L )
+	code = 2;
+if( y < 0.0L )
+	code |= 1;
+
+#ifdef NANS
+if( isnanl(x) )
+	return(x);
+if( isnanl(y) )
+	return(y);
+#endif
+#ifdef MINUSZERO
+if( y == 0.0L )
+	{
+	if( signbitl(y) )
+		{
+		if( x > 0.0L )
+			z = y;
+		else if( x < 0.0L )
+			z = -PIL;
+		else
+			{
+			if( signbitl(x) )
+				z = -PIL;
+			else
+				z = y;
+			}
+		}
+	else /* y is +0 */
+		{
+		if( x == 0.0L )
+			{
+			if( signbitl(x) )
+				z = PIL;
+			else
+				z = 0.0L;
+			}
+		else if( x > 0.0L )
+			z = 0.0L;
+		else
+			z = PIL;
+		}
+	return z;
+	}
+if( x == 0.0L )
+	{
+	if( y > 0.0L )
+		z = PIO2L;
+	else
+		z = -PIO2L;
+	return z;
+	}
+#endif /* MINUSZERO */
+#ifdef INFINITIES
+if( x == INFINITYL )
+	{
+	if( y == INFINITYL )
+		z = 0.25L * PIL;
+	else if( y == -INFINITYL )
+		z = -0.25L * PIL;
+	else if( y < 0.0L )
+		z = NEGZEROL;
+	else
+		z = 0.0L;
+	return z;
+	}
+if( x == -INFINITYL )
+	{
+	if( y == INFINITYL )
+		z = 0.75L * PIL;
+	else if( y == -INFINITYL )
+		z = -0.75L * PIL;
+	else if( y >= 0.0L )
+		z = PIL;
+	else
+		z = -PIL;
+	return z;
+	}
+if( y == INFINITYL )
+	return( PIO2L );
+if( y == -INFINITYL )
+	return( -PIO2L );
+#endif /* INFINITIES */
+
+#ifdef INFINITIES
+if( x == 0.0L )
+#else
+if( fabsl(x) <= (fabsl(y) / MAXNUML) )
+#endif
+	{
+	if( code & 1 )
+		{
+#if ANSIC
+		return( -PIO2L );
+#else
+		return( 3.0L*PIO2L );
+#endif
+		}
+	if( y == 0.0L )
+		return( 0.0L );
+	return( PIO2L );
+	}
+
+if( y == 0.0L )
+	{
+	if( code & 2 )
+		return( PIL );
+	return( 0.0L );
+	}
+
+
+switch( code )
+	{
+	default:
+#if ANSIC
+	case 0:
+	case 1: w = 0.0L; break;
+	case 2: w = PIL; break;
+	case 3: w = -PIL; break;
+#else
+	case 0: w = 0.0L; break;
+	case 1: w = 2.0L * PIL; break;
+	case 2:
+	case 3: w = PIL; break;
+#endif
+	}
+
+z = w + atanl( y/x );
+#ifdef MINUSZERO
+if( z == 0.0L && y < 0.0L )
+	z = NEGZEROL;
+#endif
+return( z );
+}