summaryrefslogtreecommitdiff
path: root/libm/float/ellpkf.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/float/ellpkf.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/float/ellpkf.c')
-rw-r--r--libm/float/ellpkf.c128
1 files changed, 128 insertions, 0 deletions
diff --git a/libm/float/ellpkf.c b/libm/float/ellpkf.c
new file mode 100644
index 000000000..2cc13d90a
--- /dev/null
+++ b/libm/float/ellpkf.c
@@ -0,0 +1,128 @@
+/* ellpkf.c
+ *
+ * Complete elliptic integral of the first kind
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float m1, y, ellpkf();
+ *
+ * y = ellpkf( m1 );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Approximates the integral
+ *
+ *
+ *
+ * pi/2
+ * -
+ * | |
+ * | dt
+ * K(m) = | ------------------
+ * | 2
+ * | | sqrt( 1 - m sin t )
+ * -
+ * 0
+ *
+ * where m = 1 - m1, using the approximation
+ *
+ * P(x) - log x Q(x).
+ *
+ * The argument m1 is used rather than m so that the logarithmic
+ * singularity at m = 1 will be shifted to the origin; this
+ * preserves maximum accuracy.
+ *
+ * K(0) = pi/2.
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0,1 30000 1.3e-7 3.4e-8
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * ellpkf domain x<0, x>1 0.0
+ *
+ */
+
+/* ellpk.c */
+
+
+/*
+Cephes Math Library, Release 2.0: April, 1987
+Copyright 1984, 1987 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+#include <math.h>
+
+static float P[] =
+{
+ 1.37982864606273237150E-4,
+ 2.28025724005875567385E-3,
+ 7.97404013220415179367E-3,
+ 9.85821379021226008714E-3,
+ 6.87489687449949877925E-3,
+ 6.18901033637687613229E-3,
+ 8.79078273952743772254E-3,
+ 1.49380448916805252718E-2,
+ 3.08851465246711995998E-2,
+ 9.65735902811690126535E-2,
+ 1.38629436111989062502E0
+};
+
+static float Q[] =
+{
+ 2.94078955048598507511E-5,
+ 9.14184723865917226571E-4,
+ 5.94058303753167793257E-3,
+ 1.54850516649762399335E-2,
+ 2.39089602715924892727E-2,
+ 3.01204715227604046988E-2,
+ 3.73774314173823228969E-2,
+ 4.88280347570998239232E-2,
+ 7.03124996963957469739E-2,
+ 1.24999999999870820058E-1,
+ 4.99999999999999999821E-1
+};
+static float C1 = 1.3862943611198906188E0; /* log(4) */
+
+extern float MACHEPF, MAXNUMF;
+
+float polevlf(float, float *, int);
+float p1evlf(float, float *, int);
+float logf(float);
+float ellpkf(float xx)
+{
+float x;
+
+x = xx;
+if( (x < 0.0) || (x > 1.0) )
+ {
+ mtherr( "ellpkf", DOMAIN );
+ return( 0.0 );
+ }
+
+if( x > MACHEPF )
+ {
+ return( polevlf(x,P,10) - logf(x) * polevlf(x,Q,10) );
+ }
+else
+ {
+ if( x == 0.0 )
+ {
+ mtherr( "ellpkf", SING );
+ return( MAXNUMF );
+ }
+ else
+ {
+ return( C1 - 0.5 * logf(x) );
+ }
+ }
+}