uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 283 insertions, 0 deletions

diff --git a/libm/float/sinf.c b/libm/float/sinf.c
new file mode 100644
index 000000000..2f1bb45b8
--- /dev/null
+++ b/libm/float/sinf.c
@@ -0,0 +1,283 @@
+/*							sinf.c
+ *
+ *	Circular sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, sinf();
+ *
+ * y = sinf( 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
+ *    IEEE    -4096,+4096   100,000      1.2e-7     3.0e-8
+ *    IEEE    -8192,+8192   100,000      3.0e-7     3.0e-8
+ * 
+ * ERROR MESSAGES:
+ *
+ *   message           condition        value returned
+ * sin total loss      x > 2^24              0.0
+ *
+ * Partial loss of accuracy begins to occur at x = 2^13
+ * = 8192. Results may be meaningless for x >= 2^24
+ * The routine as implemented flags a TLOSS error
+ * for x >= 2^24 and returns 0.0.
+ */
+
+/*							cosf.c
+ *
+ *	Circular cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, cosf();
+ *
+ * y = cosf( 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    -8192,+8192   100,000      3.0e-7     3.0e-8
+ */
+
+/*
+Cephes Math Library Release 2.2:  June, 1992
+Copyright 1985, 1987, 1988, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+
+/* Single precision circular sine
+ * test interval: [-pi/4, +pi/4]
+ * trials: 10000
+ * peak relative error: 6.8e-8
+ * rms relative error: 2.6e-8
+ */
+#include <math.h>
+
+
+static float FOPI = 1.27323954473516;
+
+extern float PIO4F;
+/* Note, these constants are for a 32-bit significand: */
+/*
+static float DP1 =  0.7853851318359375;
+static float DP2 =  1.30315311253070831298828125e-5;
+static float DP3 =  3.03855025325309630e-11;
+static float lossth = 65536.;
+*/
+
+/* These are for a 24-bit significand: */
+static float DP1 = 0.78515625;
+static float DP2 = 2.4187564849853515625e-4;
+static float DP3 = 3.77489497744594108e-8;
+static float lossth = 8192.;
+static float T24M1 = 16777215.;
+
+static float sincof[] = {
+-1.9515295891E-4,
+ 8.3321608736E-3,
+-1.6666654611E-1
+};
+static float coscof[] = {
+ 2.443315711809948E-005,
+-1.388731625493765E-003,
+ 4.166664568298827E-002
+};
+
+float sinf( float xx )
+{
+float *p;
+float x, y, z;
+register unsigned long j;
+register int sign;
+
+sign = 1;
+x = xx;
+if( xx < 0 )
+	{
+	sign = -1;
+	x = -xx;
+	}
+if( x > T24M1 )
+	{
+	mtherr( "sinf", TLOSS );
+	return(0.0);
+	}
+j = FOPI * x; /* integer part of x/(PI/4) */
+y = j;
+/* map zeros to origin */
+if( j & 1 )
+	{
+	j += 1;
+	y += 1.0;
+	}
+j &= 7; /* octant modulo 360 degrees */
+/* reflect in x axis */
+if( j > 3)
+	{
+	sign = -sign;
+	j -= 4;
+	}
+
+if( x > lossth )
+	{
+	mtherr( "sinf", PLOSS );
+	x = x - y * PIO4F;
+	}
+else
+	{
+/* Extended precision modular arithmetic */
+	x = ((x - y * DP1) - y * DP2) - y * DP3;
+	}
+/*einits();*/
+z = x * x;
+if( (j==1) || (j==2) )
+	{
+/* measured relative error in +/- pi/4 is 7.8e-8 */
+/*
+	y = ((  2.443315711809948E-005 * z
+	  - 1.388731625493765E-003) * z
+	  + 4.166664568298827E-002) * z * z;
+*/
+	p = coscof;
+	y = *p++;
+	y = y * z + *p++;
+	y = y * z + *p++;
+	y *= z * z;
+	y -= 0.5 * z;
+	y += 1.0;
+	}
+else
+	{
+/* Theoretical relative error = 3.8e-9 in [-pi/4, +pi/4] */
+/*
+	y = ((-1.9515295891E-4 * z
+	     + 8.3321608736E-3) * z
+	     - 1.6666654611E-1) * z * x;
+	y += x;
+*/
+	p = sincof;
+	y = *p++;
+	y = y * z + *p++;
+	y = y * z + *p++;
+	y *= z * x;
+	y += x;
+	}
+/*einitd();*/
+if(sign < 0)
+	y = -y;
+return( y);
+}
+
+
+/* Single precision circular cosine
+ * test interval: [-pi/4, +pi/4]
+ * trials: 10000
+ * peak relative error: 8.3e-8
+ * rms relative error: 2.2e-8
+ */
+
+float cosf( float xx )
+{
+float x, y, z;
+int j, sign;
+
+/* make argument positive */
+sign = 1;
+x = xx;
+if( x < 0 )
+	x = -x;
+
+if( x > T24M1 )
+	{
+	mtherr( "cosf", TLOSS );
+	return(0.0);
+	}
+
+j = FOPI * x; /* integer part of x/PIO4 */
+y = j;
+/* integer and fractional part modulo one octant */
+if( j & 1 )	/* map zeros to origin */
+	{
+	j += 1;
+	y += 1.0;
+	}
+j &= 7;
+if( j > 3)
+	{
+	j -=4;
+	sign = -sign;
+	}
+
+if( j > 1 )
+	sign = -sign;
+
+if( x > lossth )
+	{
+	mtherr( "cosf", PLOSS );
+	x = x - y * PIO4F;
+	}
+else
+/* Extended precision modular arithmetic */
+	x = ((x - y * DP1) - y * DP2) - y * DP3;
+
+z = x * x;
+
+if( (j==1) || (j==2) )
+	{
+	y = (((-1.9515295891E-4 * z
+	     + 8.3321608736E-3) * z
+	     - 1.6666654611E-1) * z * x)
+	     + x;
+	}
+else
+	{
+	y = ((  2.443315711809948E-005 * z
+	  - 1.388731625493765E-003) * z
+	  + 4.166664568298827E-002) * z * z;
+	y -= 0.5 * z;
+	y += 1.0;
+	}
+if(sign < 0)
+	y = -y;
+return( y );
+}
+