1 files changed, 148 insertions, 0 deletions
diff --git a/libm/double/ellik.c b/libm/double/ellik.c
new file mode 100644
index 000000000..1c9053676
--- /dev/null
+++ b/libm/double/ellik.c
@@ -0,0 +1,148 @@
+/*							ellik.c
+ *
+ *	Incomplete elliptic integral of the first kind
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double phi, m, y, ellik();
+ *
+ * y = ellik( 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 m in [0, 1] and phi as indicated.
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     -10,10       200000      7.4e-16     1.0e-16
+ *
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*/
+
+/*	Incomplete elliptic integral of first kind	*/
+
+#include <math.h>
+#ifdef ANSIPROT
+extern double sqrt ( double );
+extern double fabs ( double );
+extern double log ( double );
+extern double tan ( double );
+extern double atan ( double );
+extern double floor ( double );
+extern double ellpk ( double );
+double ellik ( double, double );
+#else
+double sqrt(), fabs(), log(), tan(), atan(), floor(), ellpk();
+double ellik();
+#endif
+extern double PI, PIO2, MACHEP, MAXNUM;
+
+double ellik( phi, m )
+double phi, m;
+{
+double a, b, c, e, temp, t, K;
+int d, mod, sign, npio2;
+
+if( m == 0.0 )
+	return( phi );
+a = 1.0 - m;
+if( a == 0.0 )
+	{
+	if( fabs(phi) >= PIO2 )
+		{
+		mtherr( "ellik", SING );
+		return( MAXNUM );
+		}
+	return(  log(  tan( (PIO2 + phi)/2.0 )  )   );
+	}
+npio2 = floor( phi/PIO2 );
+if( npio2 & 1 )
+	npio2 += 1;
+if( npio2 )
+	{
+	K = ellpk( a );
+	phi = phi - npio2 * PIO2;
+	}
+else
+	K = 0.0;
+if( phi < 0.0 )
+	{
+	phi = -phi;
+	sign = -1;
+	}
+else
+	sign = 0;
+b = sqrt(a);
+t = tan( phi );
+if( fabs(t) > 10.0 )
+	{
+	/* Transform the amplitude */
+	e = 1.0/(b*t);
+	/* ... but avoid multiple recursions.  */
+	if( fabs(e) < 10.0 )
+		{
+		e = atan(e);
+		if( npio2 == 0 )
+			K = ellpk( a );
+		temp = K - ellik( e, m );
+		goto done;
+		}
+	}
+a = 1.0;
+c = sqrt(m);
+d = 1;
+mod = 0;
+
+while( fabs(c/a) > MACHEP )
+	{
+	temp = b/a;
+	phi = phi + atan(t*temp) + mod * PI;
+	mod = (phi + PIO2)/PI;
+	t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
+	c = ( a - b )/2.0;
+	temp = sqrt( a * b );
+	a = ( a + b )/2.0;
+	b = temp;
+	d += d;
+	}
+
+temp = (atan(t) + mod * PI)/(d * a);
+
+done:
+if( sign < 0 )
+	temp = -temp;
+temp += npio2 * K;
+return( temp );
+}