summaryrefslogtreecommitdiff
path: root/libm/float/ellikf.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/float/ellikf.c')
-rw-r--r--libm/float/ellikf.c113
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 );
+}