uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 171 insertions, 0 deletions

diff --git a/libm/double/ellpj.c b/libm/double/ellpj.c
new file mode 100644
index 000000000..327fc56e8
--- /dev/null
+++ b/libm/double/ellpj.c
@@ -0,0 +1,171 @@
+/*							ellpj.c
+ *
+ *	Jacobian Elliptic Functions
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double u, m, sn, cn, dn, phi;
+ * int ellpj();
+ *
+ * ellpj( u, m, _&sn, _&cn, _&dn, _&phi );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ *
+ * Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
+ * and dn(u|m) of parameter m between 0 and 1, and real
+ * argument u.
+ *
+ * These functions are periodic, with quarter-period on the
+ * real axis equal to the complete elliptic integral
+ * ellpk(1.0-m).
+ *
+ * Relation to incomplete elliptic integral:
+ * If u = ellik(phi,m), then sn(u|m) = sin(phi),
+ * and cn(u|m) = cos(phi).  Phi is called the amplitude of u.
+ *
+ * Computation is by means of the arithmetic-geometric mean
+ * algorithm, except when m is within 1e-9 of 0 or 1.  In the
+ * latter case with m close to 1, the approximation applies
+ * only for phi < pi/2.
+ *
+ * ACCURACY:
+ *
+ * Tested at random points with u between 0 and 10, m between
+ * 0 and 1.
+ *
+ *            Absolute error (* = relative error):
+ * arithmetic   function   # trials      peak         rms
+ *    DEC       sn           1800       4.5e-16     8.7e-17
+ *    IEEE      phi         10000       9.2e-16*    1.4e-16*
+ *    IEEE      sn          50000       4.1e-15     4.6e-16
+ *    IEEE      cn          40000       3.6e-15     4.4e-16
+ *    IEEE      dn          10000       1.3e-12     1.8e-14
+ *
+ *  Peak error observed in consistency check using addition
+ * theorem for sn(u+v) was 4e-16 (absolute).  Also tested by
+ * the above relation to the incomplete elliptic integral.
+ * Accuracy deteriorates when u is large.
+ *
+ */
+
+/*							ellpj.c		*/
+
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*/
+
+#include <math.h>
+#ifdef ANSIPROT
+extern double sqrt ( double );
+extern double fabs ( double );
+extern double sin ( double );
+extern double cos ( double );
+extern double asin ( double );
+extern double tanh ( double );
+extern double sinh ( double );
+extern double cosh ( double );
+extern double atan ( double );
+extern double exp ( double );
+#else
+double sqrt(), fabs(), sin(), cos(), asin(), tanh();
+double sinh(), cosh(), atan(), exp();
+#endif
+extern double PIO2, MACHEP;
+
+int ellpj( u, m, sn, cn, dn, ph )
+double u, m;
+double *sn, *cn, *dn, *ph;
+{
+double ai, b, phi, t, twon;
+double a[9], c[9];
+int i;
+
+
+/* Check for special cases */
+
+if( m < 0.0 || m > 1.0 )
+	{
+	mtherr( "ellpj", DOMAIN );
+	*sn = 0.0;
+	*cn = 0.0;
+	*ph = 0.0;
+	*dn = 0.0;
+	return(-1);
+	}
+if( m < 1.0e-9 )
+	{
+	t = sin(u);
+	b = cos(u);
+	ai = 0.25 * m * (u - t*b);
+	*sn = t - ai*b;
+	*cn = b + ai*t;
+	*ph = u - ai;
+	*dn = 1.0 - 0.5*m*t*t;
+	return(0);
+	}
+
+if( m >= 0.9999999999 )
+	{
+	ai = 0.25 * (1.0-m);
+	b = cosh(u);
+	t = tanh(u);
+	phi = 1.0/b;
+	twon = b * sinh(u);
+	*sn = t + ai * (twon - u)/(b*b);
+	*ph = 2.0*atan(exp(u)) - PIO2 + ai*(twon - u)/b;
+	ai *= t * phi;
+	*cn = phi - ai * (twon - u);
+	*dn = phi + ai * (twon + u);
+	return(0);
+	}
+
+
+/*	A. G. M. scale		*/
+a[0] = 1.0;
+b = sqrt(1.0 - m);
+c[0] = sqrt(m);
+twon = 1.0;
+i = 0;
+
+while( fabs(c[i]/a[i]) > MACHEP )
+	{
+	if( i > 7 )
+		{
+		mtherr( "ellpj", OVERFLOW );
+		goto done;
+		}
+	ai = a[i];
+	++i;
+	c[i] = ( ai - b )/2.0;
+	t = sqrt( ai * b );
+	a[i] = ( ai + b )/2.0;
+	b = t;
+	twon *= 2.0;
+	}
+
+done:
+
+/* backward recurrence */
+phi = twon * a[i] * u;
+do
+	{
+	t = c[i] * sin(phi) / a[i];
+	b = phi;
+	phi = (asin(t) + phi)/2.0;
+	}
+while( --i );
+
+*sn = sin(phi);
+t = cos(phi);
+*cn = t;
+*dn = t/cos(phi-b);
+*ph = phi;
+return(0);
+}