1 files changed, 113 insertions, 0 deletions
diff --git a/libm/float/ellikf.c b/libm/float/ellikf.c
new file mode 100644
index 000000000..8ec890926
--- /dev/null
+++ b/libm/float/ellikf.c
@@ -0,0 +1,113 @@
+/*							ellikf.c
+ *
+ *	Incomplete elliptic integral of the first kind
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float phi, m, y, ellikf();
+ *
+ * y = ellikf( phi, m );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Approximates the integral
+ *
+ *
+ *
+ *                phi
+ *                 -
+ *                | |
+ *                |           dt
+ * F(phi\m)  =    |    ------------------
+ *                |                   2
+ *              | |    sqrt( 1 - m sin t )
+ *               -
+ *                0
+ *
+ * of amplitude phi and modulus m, using the arithmetic -
+ * geometric mean algorithm.
+ *
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Tested at random points with phi in [0, 2] and m in
+ * [0, 1].
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0,2         10000       2.9e-7      5.8e-8
+ *
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.2:  July, 1992
+Copyright 1984, 1987, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+/*	Incomplete elliptic integral of first kind	*/
+
+#include <math.h>
+extern float PIF, PIO2F, MACHEPF;
+
+#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
+
+#ifdef ANSIC
+float sqrtf(float), logf(float), sinf(float), tanf(float), atanf(float);
+#else
+float sqrtf(), logf(), sinf(), tanf(), atanf();
+#endif
+
+
+float ellikf( float phia, float ma )
+{
+float phi, m, a, b, c, temp;
+float t;
+int d, mod, sign;
+
+phi = phia;
+m = ma;
+if( m == 0.0 )
+	return( phi );
+if( phi < 0.0 )
+	{
+	phi = -phi;
+	sign = -1;
+	}
+else
+	sign = 0;
+a = 1.0;
+b = 1.0 - m;
+if( b == 0.0 )
+	return(  logf(  tanf( 0.5*(PIO2F + phi) )  )   );
+b = sqrtf(b);
+c = sqrtf(m);
+d = 1;
+t = tanf( phi );
+mod = (phi + PIO2F)/PIF;
+
+while( fabsf(c/a) > MACHEPF )
+	{
+	temp = b/a;
+	phi = phi + atanf(t*temp) + mod * PIF;
+	mod = (phi + PIO2F)/PIF;
+	t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
+	c = ( a - b )/2.0;
+	temp = sqrtf( a * b );
+	a = ( a + b )/2.0;
+	b = temp;
+	d += d;
+	}
+
+temp = (atanf(t) + mod * PIF)/(d * a);
+if( sign < 0 )
+	temp = -temp;
+return( temp );
+}