summaryrefslogtreecommitdiff
path: root/libm/ldouble/ellikl.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/ldouble/ellikl.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/ldouble/ellikl.c')
-rw-r--r--libm/ldouble/ellikl.c148
1 files changed, 148 insertions, 0 deletions
diff --git a/libm/ldouble/ellikl.c b/libm/ldouble/ellikl.c
new file mode 100644
index 000000000..4eeffe0f5
--- /dev/null
+++ b/libm/ldouble/ellikl.c
@@ -0,0 +1,148 @@
+/* ellikl.c
+ *
+ * Incomplete elliptic integral of the first kind
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double phi, m, y, ellikl();
+ *
+ * y = ellikl( 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 30000 3.6e-18 4.1e-19
+ *
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.3: November, 1995
+Copyright 1984, 1987, 1995 by Stephen L. Moshier
+*/
+
+/* Incomplete elliptic integral of first kind */
+
+#include <math.h>
+#ifdef ANSIPROT
+extern long double sqrtl ( long double );
+extern long double fabsl ( long double );
+extern long double logl ( long double );
+extern long double tanl ( long double );
+extern long double atanl ( long double );
+extern long double floorl ( long double );
+extern long double ellpkl ( long double );
+long double ellikl ( long double, long double );
+#else
+long double sqrtl(), fabsl(), logl(), tanl(), atanl(), floorl(), ellpkl();
+long double ellikl();
+#endif
+extern long double PIL, PIO2L, MACHEPL, MAXNUML;
+
+long double ellikl( phi, m )
+long double phi, m;
+{
+long double a, b, c, e, temp, t, K;
+int d, mod, sign, npio2;
+
+if( m == 0.0L )
+ return( phi );
+a = 1.0L - m;
+if( a == 0.0L )
+ {
+ if( fabsl(phi) >= PIO2L )
+ {
+ mtherr( "ellikl", SING );
+ return( MAXNUML );
+ }
+ return( logl( tanl( 0.5L*(PIO2L + phi) ) ) );
+ }
+npio2 = floorl( phi/PIO2L );
+if( npio2 & 1 )
+ npio2 += 1;
+if( npio2 )
+ {
+ K = ellpkl( a );
+ phi = phi - npio2 * PIO2L;
+ }
+else
+ K = 0.0L;
+if( phi < 0.0L )
+ {
+ phi = -phi;
+ sign = -1;
+ }
+else
+ sign = 0;
+b = sqrtl(a);
+t = tanl( phi );
+if( fabsl(t) > 10.0L )
+ {
+ /* Transform the amplitude */
+ e = 1.0L/(b*t);
+ /* ... but avoid multiple recursions. */
+ if( fabsl(e) < 10.0L )
+ {
+ e = atanl(e);
+ if( npio2 == 0 )
+ K = ellpkl( a );
+ temp = K - ellikl( e, m );
+ goto done;
+ }
+ }
+a = 1.0L;
+c = sqrtl(m);
+d = 1;
+mod = 0;
+
+while( fabsl(c/a) > MACHEPL )
+ {
+ temp = b/a;
+ phi = phi + atanl(t*temp) + mod * PIL;
+ mod = (phi + PIO2L)/PIL;
+ t = t * ( 1.0L + temp )/( 1.0L - temp * t * t );
+ c = 0.5L * ( a - b );
+ temp = sqrtl( a * b );
+ a = 0.5L * ( a + b );
+ b = temp;
+ d += d;
+ }
+
+temp = (atanl(t) + mod * PIL)/(d * a);
+
+done:
+if( sign < 0 )
+ temp = -temp;
+temp += npio2 * K;
+return( temp );
+}