uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 387 insertions, 0 deletions

diff --git a/libm/double/sin.c b/libm/double/sin.c
new file mode 100644
index 000000000..24746d79d
--- /dev/null
+++ b/libm/double/sin.c
@@ -0,0 +1,387 @@
+/*							sin.c
+ *
+ *	Circular sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, sin();
+ *
+ * y = sin( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Range reduction is into intervals of pi/4.  The reduction
+ * error is nearly eliminated by contriving an extended precision
+ * modular arithmetic.
+ *
+ * Two polynomial approximating functions are employed.
+ * Between 0 and pi/4 the sine is approximated by
+ *      x  +  x**3 P(x**2).
+ * Between pi/4 and pi/2 the cosine is represented as
+ *      1  -  x**2 Q(x**2).
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain      # trials      peak         rms
+ *    DEC       0, 10       150000       3.0e-17     7.8e-18
+ *    IEEE -1.07e9,+1.07e9  130000       2.1e-16     5.4e-17
+ * 
+ * ERROR MESSAGES:
+ *
+ *   message           condition        value returned
+ * sin total loss   x > 1.073741824e9      0.0
+ *
+ * Partial loss of accuracy begins to occur at x = 2**30
+ * = 1.074e9.  The loss is not gradual, but jumps suddenly to
+ * about 1 part in 10e7.  Results may be meaningless for
+ * x > 2**49 = 5.6e14.  The routine as implemented flags a
+ * TLOSS error for x > 2**30 and returns 0.0.
+ */
+/*							cos.c
+ *
+ *	Circular cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, cos();
+ *
+ * y = cos( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Range reduction is into intervals of pi/4.  The reduction
+ * error is nearly eliminated by contriving an extended precision
+ * modular arithmetic.
+ *
+ * Two polynomial approximating functions are employed.
+ * Between 0 and pi/4 the cosine is approximated by
+ *      1  -  x**2 Q(x**2).
+ * Between pi/4 and pi/2 the sine is represented as
+ *      x  +  x**3 P(x**2).
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain      # trials      peak         rms
+ *    IEEE -1.07e9,+1.07e9  130000       2.1e-16     5.4e-17
+ *    DEC        0,+1.07e9   17000       3.0e-17     7.2e-18
+ */
+
+/*							sin.c	*/
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1985, 1995, 2000 by Stephen L. Moshier
+*/
+
+#include <math.h>
+
+#ifdef UNK
+static double sincof[] = {
+ 1.58962301576546568060E-10,
+-2.50507477628578072866E-8,
+ 2.75573136213857245213E-6,
+-1.98412698295895385996E-4,
+ 8.33333333332211858878E-3,
+-1.66666666666666307295E-1,
+};
+static double coscof[6] = {
+-1.13585365213876817300E-11,
+ 2.08757008419747316778E-9,
+-2.75573141792967388112E-7,
+ 2.48015872888517045348E-5,
+-1.38888888888730564116E-3,
+ 4.16666666666665929218E-2,
+};
+static double DP1 =   7.85398125648498535156E-1;
+static double DP2 =   3.77489470793079817668E-8;
+static double DP3 =   2.69515142907905952645E-15;
+/* static double lossth = 1.073741824e9; */
+#endif
+
+#ifdef DEC
+static unsigned short sincof[] = {
+0030056,0143750,0177214,0163153,
+0131727,0027455,0044510,0175352,
+0033470,0167432,0131752,0042414,
+0135120,0006400,0146776,0174027,
+0036410,0104210,0104207,0137202,
+0137452,0125252,0125252,0125103,
+};
+static unsigned short coscof[24] = {
+0127107,0151115,0002060,0152325,
+0031017,0072353,0155161,0174053,
+0132623,0171173,0172542,0057056,
+0034320,0006400,0147102,0023652,
+0135666,0005540,0133012,0076213,
+0037052,0125252,0125252,0125126,
+};
+/*  7.853981629014015197753906250000E-1 */
+static unsigned short P1[] = {0040111,0007732,0120000,0000000,};
+/*  4.960467869796758577649598009884E-10 */
+static unsigned short P2[] = {0030410,0055060,0100000,0000000,};
+/*  2.860594363054915898381331279295E-18 */
+static unsigned short P3[] = {0021523,0011431,0105056,0001560,};
+#define DP1 *(double *)P1
+#define DP2 *(double *)P2
+#define DP3 *(double *)P3
+#endif
+
+#ifdef IBMPC
+static unsigned short sincof[] = {
+0x9ccd,0x1fd1,0xd8fd,0x3de5,
+0x1f5d,0xa929,0xe5e5,0xbe5a,
+0x48a1,0x567d,0x1de3,0x3ec7,
+0xdf03,0x19bf,0x01a0,0xbf2a,
+0xf7d0,0x1110,0x1111,0x3f81,
+0x5548,0x5555,0x5555,0xbfc5,
+};
+static unsigned short coscof[24] = {
+0x1a9b,0xa086,0xfa49,0xbda8,
+0x3f05,0x7b4e,0xee9d,0x3e21,
+0x4bc6,0x7eac,0x7e4f,0xbe92,
+0x44f5,0x19c8,0x01a0,0x3efa,
+0x4f91,0x16c1,0xc16c,0xbf56,
+0x554b,0x5555,0x5555,0x3fa5,
+};
+/*
+  7.85398125648498535156E-1,
+  3.77489470793079817668E-8,
+  2.69515142907905952645E-15,
+*/
+static unsigned short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
+static unsigned short P2[] = {0x0000,0x0000,0x442d,0x3e64};
+static unsigned short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
+#define DP1 *(double *)P1
+#define DP2 *(double *)P2
+#define DP3 *(double *)P3
+#endif
+
+#ifdef MIEEE
+static unsigned short sincof[] = {
+0x3de5,0xd8fd,0x1fd1,0x9ccd,
+0xbe5a,0xe5e5,0xa929,0x1f5d,
+0x3ec7,0x1de3,0x567d,0x48a1,
+0xbf2a,0x01a0,0x19bf,0xdf03,
+0x3f81,0x1111,0x1110,0xf7d0,
+0xbfc5,0x5555,0x5555,0x5548,
+};
+static unsigned short coscof[24] = {
+0xbda8,0xfa49,0xa086,0x1a9b,
+0x3e21,0xee9d,0x7b4e,0x3f05,
+0xbe92,0x7e4f,0x7eac,0x4bc6,
+0x3efa,0x01a0,0x19c8,0x44f5,
+0xbf56,0xc16c,0x16c1,0x4f91,
+0x3fa5,0x5555,0x5555,0x554b,
+};
+static unsigned short P1[] = {0x3fe9,0x21fb,0x4000,0x0000};
+static unsigned short P2[] = {0x3e64,0x442d,0x0000,0x0000};
+static unsigned short P3[] = {0x3ce8,0x4698,0x98cc,0x5170};
+#define DP1 *(double *)P1
+#define DP2 *(double *)P2
+#define DP3 *(double *)P3
+#endif
+
+#ifdef ANSIPROT
+extern double polevl ( double, void *, int );
+extern double p1evl ( double, void *, int );
+extern double floor ( double );
+extern double ldexp ( double, int );
+extern int isnan ( double );
+extern int isfinite ( double );
+#else
+double polevl(), floor(), ldexp();
+int isnan(), isfinite();
+#endif
+extern double PIO4;
+static double lossth = 1.073741824e9;
+#ifdef NANS
+extern double NAN;
+#endif
+#ifdef INFINITIES
+extern double INFINITY;
+#endif
+
+
+double sin(x)
+double x;
+{
+double y, z, zz;
+int j, sign;
+
+#ifdef MINUSZERO
+if( x == 0.0 )
+	return(x);
+#endif
+#ifdef NANS
+if( isnan(x) )
+	return(x);
+if( !isfinite(x) )
+	{
+	mtherr( "sin", DOMAIN );
+	return(NAN);
+	}
+#endif
+/* make argument positive but save the sign */
+sign = 1;
+if( x < 0 )
+	{
+	x = -x;
+	sign = -1;
+	}
+
+if( x > lossth )
+	{
+	mtherr( "sin", TLOSS );
+	return(0.0);
+	}
+
+y = floor( x/PIO4 ); /* integer part of x/PIO4 */
+
+/* strip high bits of integer part to prevent integer overflow */
+z = ldexp( y, -4 );
+z = floor(z);           /* integer part of y/8 */
+z = y - ldexp( z, 4 );  /* y - 16 * (y/16) */
+
+j = z; /* convert to integer for tests on the phase angle */
+/* map zeros to origin */
+if( j & 1 )
+	{
+	j += 1;
+	y += 1.0;
+	}
+j = j & 07; /* octant modulo 360 degrees */
+/* reflect in x axis */
+if( j > 3)
+	{
+	sign = -sign;
+	j -= 4;
+	}
+
+/* Extended precision modular arithmetic */
+z = ((x - y * DP1) - y * DP2) - y * DP3;
+
+zz = z * z;
+
+if( (j==1) || (j==2) )
+	{
+	y = 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
+	}
+else
+	{
+/*	y = z  +  z * (zz * polevl( zz, sincof, 5 ));*/
+	y = z  +  z * z * z * polevl( zz, sincof, 5 );
+	}
+
+if(sign < 0)
+	y = -y;
+
+return(y);
+}
+
+
+
+
+
+double cos(x)
+double x;
+{
+double y, z, zz;
+long i;
+int j, sign;
+
+#ifdef NANS
+if( isnan(x) )
+	return(x);
+if( !isfinite(x) )
+	{
+	mtherr( "cos", DOMAIN );
+	return(NAN);
+	}
+#endif
+
+/* make argument positive */
+sign = 1;
+if( x < 0 )
+	x = -x;
+
+if( x > lossth )
+	{
+	mtherr( "cos", TLOSS );
+	return(0.0);
+	}
+
+y = floor( x/PIO4 );
+z = ldexp( y, -4 );
+z = floor(z);		/* integer part of y/8 */
+z = y - ldexp( z, 4 );  /* y - 16 * (y/16) */
+
+/* integer and fractional part modulo one octant */
+i = z;
+if( i & 1 )	/* map zeros to origin */
+	{
+	i += 1;
+	y += 1.0;
+	}
+j = i & 07;
+if( j > 3)
+	{
+	j -=4;
+	sign = -sign;
+	}
+
+if( j > 1 )
+	sign = -sign;
+
+/* Extended precision modular arithmetic */
+z = ((x - y * DP1) - y * DP2) - y * DP3;
+
+zz = z * z;
+
+if( (j==1) || (j==2) )
+	{
+/*	y = z  +  z * (zz * polevl( zz, sincof, 5 ));*/
+	y = z  +  z * z * z * polevl( zz, sincof, 5 );
+	}
+else
+	{
+	y = 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
+	}
+
+if(sign < 0)
+	y = -y;
+
+return(y);
+}
+
+
+
+
+
+/* Degrees, minutes, seconds to radians: */
+
+/* 1 arc second, in radians = 4.8481368110953599358991410e-5 */
+#ifdef DEC
+static unsigned short P648[] = {034513,054170,0176773,0116043,};
+#define P64800 *(double *)P648
+#else
+static double P64800 = 4.8481368110953599358991410e-5;
+#endif
+
+double radian(d,m,s)
+double d,m,s;
+{
+
+return( ((d*60.0 + m)*60.0 + s)*P64800 );
+}