uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 308 insertions, 0 deletions

diff --git a/libm/double/sindg.c b/libm/double/sindg.c
new file mode 100644
index 000000000..8057ab68d
--- /dev/null
+++ b/libm/double/sindg.c
@@ -0,0 +1,308 @@
+/*							sindg.c
+ *
+ *	Circular sine of angle in degrees
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, sindg();
+ *
+ * y = sindg( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Range reduction is into intervals of 45 degrees.
+ *
+ * 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 P(x**2).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain      # trials      peak         rms
+ *    DEC       +-1000        3100      3.3e-17      9.0e-18
+ *    IEEE      +-1000       30000      2.3e-16      5.6e-17
+ * 
+ * ERROR MESSAGES:
+ *
+ *   message           condition        value returned
+ * sindg total loss   x > 8.0e14 (DEC)      0.0
+ *                    x > 1.0e14 (IEEE)
+ *
+ */
+/*							cosdg.c
+ *
+ *	Circular cosine of angle in degrees
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, cosdg();
+ *
+ * y = cosdg( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Range reduction is into intervals of 45 degrees.
+ *
+ * Two polynomial approximating functions are employed.
+ * Between 0 and pi/4 the cosine is approximated by
+ *      1  -  x**2 P(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
+ *    DEC      +-1000         3400       3.5e-17     9.1e-18
+ *    IEEE     +-1000        30000       2.1e-16     5.7e-17
+ *  See also sin().
+ *
+ */
+
+/* Cephes Math Library Release 2.0:  April, 1987
+ * Copyright 1985, 1987 by Stephen L. Moshier
+ * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */
+
+#include <math.h>
+
+#ifdef UNK
+static double sincof[] = {
+ 1.58962301572218447952E-10,
+-2.50507477628503540135E-8,
+ 2.75573136213856773549E-6,
+-1.98412698295895384658E-4,
+ 8.33333333332211858862E-3,
+-1.66666666666666307295E-1
+};
+static double coscof[] = {
+ 1.13678171382044553091E-11,
+-2.08758833757683644217E-9,
+ 2.75573155429816611547E-7,
+-2.48015872936186303776E-5,
+ 1.38888888888806666760E-3,
+-4.16666666666666348141E-2,
+ 4.99999999999999999798E-1
+};
+static double PI180 = 1.74532925199432957692E-2; /* pi/180 */
+static double lossth = 1.0e14;
+#endif
+
+#ifdef DEC
+static unsigned short sincof[] = {
+0030056,0143750,0177170,0073013,
+0131727,0027455,0044510,0132205,
+0033470,0167432,0131752,0042263,
+0135120,0006400,0146776,0174027,
+0036410,0104210,0104207,0137202,
+0137452,0125252,0125252,0125103
+};
+static unsigned short coscof[] = {
+0027107,0176030,0153315,0110312,
+0131017,0072476,0007450,0123243,
+0032623,0171174,0070066,0146445,
+0134320,0006400,0147355,0163313,
+0035666,0005540,0133012,0165067,
+0137052,0125252,0125252,0125206,
+0040000,0000000,0000000,0000000
+};
+static unsigned short P1[] = {0036616,0175065,0011224,0164711};
+#define PI180 *(double *)P1
+static double lossth = 8.0e14;
+#endif
+
+#ifdef IBMPC
+static unsigned short sincof[] = {
+0x0ec1,0x1fcf,0xd8fd,0x3de5,
+0x1691,0xa929,0xe5e5,0xbe5a,
+0x4896,0x567d,0x1de3,0x3ec7,
+0xdf03,0x19bf,0x01a0,0xbf2a,
+0xf7d0,0x1110,0x1111,0x3f81,
+0x5548,0x5555,0x5555,0xbfc5
+};
+static unsigned short coscof[] = {
+0xb219,0x1ad9,0xff83,0x3da8,
+0x14d4,0xc1e5,0xeea7,0xbe21,
+0xd9a5,0x8e06,0x7e4f,0x3e92,
+0xbcd9,0x19dd,0x01a0,0xbefa,
+0x5d47,0x16c1,0xc16c,0x3f56,
+0x5551,0x5555,0x5555,0xbfa5,
+0x0000,0x0000,0x0000,0x3fe0
+};
+
+static unsigned short P1[] = {0x9d39,0xa252,0xdf46,0x3f91};
+#define PI180 *(double *)P1
+static double lossth = 1.0e14;
+#endif
+
+#ifdef MIEEE
+static unsigned short sincof[] = {
+0x3de5,0xd8fd,0x1fcf,0x0ec1,
+0xbe5a,0xe5e5,0xa929,0x1691,
+0x3ec7,0x1de3,0x567d,0x4896,
+0xbf2a,0x01a0,0x19bf,0xdf03,
+0x3f81,0x1111,0x1110,0xf7d0,
+0xbfc5,0x5555,0x5555,0x5548
+};
+static unsigned short coscof[] = {
+0x3da8,0xff83,0x1ad9,0xb219,
+0xbe21,0xeea7,0xc1e5,0x14d4,
+0x3e92,0x7e4f,0x8e06,0xd9a5,
+0xbefa,0x01a0,0x19dd,0xbcd9,
+0x3f56,0xc16c,0x16c1,0x5d47,
+0xbfa5,0x5555,0x5555,0x5551,
+0x3fe0,0x0000,0x0000,0x0000
+};
+
+static unsigned short P1[] = {
+0x3f91,0xdf46,0xa252,0x9d39
+};
+#define PI180 *(double *)P1
+static double lossth = 1.0e14;
+#endif
+
+#ifdef ANSIPROT
+extern double polevl ( double, void *, int );
+extern double floor ( double );
+extern double ldexp ( double, int );
+#else
+double polevl(), floor(), ldexp();
+#endif
+extern double PIO4;
+
+double sindg(x)
+double x;
+{
+double y, z, zz;
+int j, sign;
+
+/* make argument positive but save the sign */
+sign = 1;
+if( x < 0 )
+	{
+	x = -x;
+	sign = -1;
+	}
+
+if( x > lossth )
+	{
+	mtherr( "sindg", TLOSS );
+	return(0.0);
+	}
+
+y = floor( x/45.0 ); /* 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;
+	}
+
+z = x - y * 45.0; /* x mod 45 degrees */
+z *= PI180;	/* multiply by pi/180 to convert to radians */
+zz = z * z;
+
+if( (j==1) || (j==2) )
+	{
+	y = 1.0 - zz * polevl( zz, coscof, 6 );
+	}
+else
+	{
+	y = z  +  z * (zz * polevl( zz, sincof, 5 ));
+	}
+
+if(sign < 0)
+	y = -y;
+
+return(y);
+}
+
+
+
+
+
+double cosdg(x)
+double x;
+{
+double y, z, zz;
+int j, sign;
+
+/* make argument positive */
+sign = 1;
+if( x < 0 )
+	x = -x;
+
+if( x > lossth )
+	{
+	mtherr( "cosdg", TLOSS );
+	return(0.0);
+	}
+
+y = floor( x/45.0 );
+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 */
+j = z;
+if( j & 1 )	/* map zeros to origin */
+	{
+	j += 1;
+	y += 1.0;
+	}
+j = j & 07;
+if( j > 3)
+	{
+	j -=4;
+	sign = -sign;
+	}
+
+if( j > 1 )
+	sign = -sign;
+
+z = x - y * 45.0; /* x mod 45 degrees */
+z *= PI180;	/* multiply by pi/180 to convert to radians */
+
+zz = z * z;
+
+if( (j==1) || (j==2) )
+	{
+	y = z  +  z * (zz * polevl( zz, sincof, 5 ));
+	}
+else
+	{
+	y = 1.0 - zz * polevl( zz, coscof, 6 );
+	}
+
+if(sign < 0)
+	y = -y;
+
+return(y);
+}