1 files changed, 173 insertions, 0 deletions
diff --git a/libm/ldouble/ellpel.c b/libm/ldouble/ellpel.c
new file mode 100644
index 000000000..6965db066
--- /dev/null
+++ b/libm/ldouble/ellpel.c
@@ -0,0 +1,173 @@
+/*							ellpel.c
+ *
+ *	Complete elliptic integral of the second kind
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double m1, y, ellpel();
+ *
+ * y = ellpel( m1 );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Approximates the integral
+ *
+ *
+ *            pi/2
+ *             -
+ *            | |                 2
+ * E(m)  =    |    sqrt( 1 - m sin t ) dt
+ *          | |    
+ *           -
+ *            0
+ *
+ * Where m = 1 - m1, using the approximation
+ *
+ *      P(x)  -  x log x Q(x).
+ *
+ * Though there are no singularities, the argument m1 is used
+ * rather than m for compatibility with ellpk().
+ *
+ * E(1) = 1; E(0) = pi/2.
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE       0, 1       10000       1.1e-19     3.5e-20
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * ellpel domain     x<0, x>1            0.0
+ *
+ */
+
+/*							ellpe.c		*/
+
+/* Elliptic integral of second kind */
+
+/*
+Cephes Math Library, Release 2.3:  October, 1995
+Copyright 1984, 1987, 1989, 1995 by Stephen L. Moshier
+*/
+
+#include <math.h>
+
+#if UNK
+static long double P[12] = {
+ 3.198937812032341294902E-5L,
+ 7.742523238588775116241E-4L,
+ 4.140384701571542000550E-3L,
+ 7.963509564694454269086E-3L,
+ 7.280911706839967541799E-3L,
+ 5.044067167184043853799E-3L,
+ 5.076832243257395296304E-3L,
+ 7.155775630578315248130E-3L,
+ 1.154485760526450950611E-2L,
+ 2.183137319801117971860E-2L,
+ 5.680519271556930583433E-2L,
+ 4.431471805599467050354E-1L,
+};
+static long double Q[12] = {
+ 6.393938134301205485085E-6L,
+ 2.741404591220851603273E-4L,
+ 2.480876752984331133799E-3L,
+ 8.770638497964078750003E-3L,
+ 1.676835725889463343319E-2L,
+ 2.281970801531577700830E-2L,
+ 2.767367465121309044166E-2L,
+ 3.364167778770018154356E-2L,
+ 4.272453406734691973083E-2L,
+ 5.859374951483909267451E-2L,
+ 9.374999999923942267270E-2L,
+ 2.499999999999998643587E-1L,
+};
+#endif
+#if IBMPC
+static short P[] = {
+0x7a78,0x5a02,0x554d,0x862c,0x3ff0, XPD
+0x34db,0xa965,0x31a3,0xcaf7,0x3ff4, XPD
+0xca6c,0x6c00,0x1071,0x87ac,0x3ff7, XPD
+0x4cdb,0x125d,0x6149,0x8279,0x3ff8, XPD
+0xadbd,0x3d8f,0xb6d5,0xee94,0x3ff7, XPD
+0x8189,0xcd0e,0xb3c2,0xa548,0x3ff7, XPD
+0x32b5,0xdd64,0x8e39,0xa65b,0x3ff7, XPD
+0x0344,0xc9db,0xff27,0xea7a,0x3ff7, XPD
+0xba2d,0x806a,0xa476,0xbd26,0x3ff8, XPD
+0xc3e0,0x30fa,0xb53d,0xb2d7,0x3ff9, XPD
+0x23b8,0x4d33,0x8fcf,0xe8ac,0x3ffa, XPD
+0xbc79,0xa39f,0x2fef,0xe2e4,0x3ffd, XPD
+};
+static short Q[] = {
+0x89f1,0xe234,0x82a6,0xd68b,0x3fed, XPD
+0x202a,0x96b3,0x8273,0x8fba,0x3ff3, XPD
+0xc183,0xfc45,0x3484,0xa296,0x3ff6, XPD
+0x683e,0xe201,0xb960,0x8fb2,0x3ff8, XPD
+0x721a,0x1b6a,0xcb41,0x895d,0x3ff9, XPD
+0x4eee,0x295f,0x6574,0xbaf0,0x3ff9, XPD
+0x3ade,0xc98f,0xe6f2,0xe2b3,0x3ff9, XPD
+0xd470,0x1784,0xdb1e,0x89cb,0x3ffa, XPD
+0xa649,0xe5c1,0xebc8,0xaeff,0x3ffa, XPD
+0x84c0,0xa8f5,0xffde,0xefff,0x3ffa, XPD
+0x5506,0xf94f,0xffff,0xbfff,0x3ffb, XPD
+0xd8e7,0xffff,0xffff,0xffff,0x3ffc, XPD
+};
+#endif
+#if MIEEE
+static long P[36] = {
+0x3ff00000,0x862c554d,0x5a027a78,
+0x3ff40000,0xcaf731a3,0xa96534db,
+0x3ff70000,0x87ac1071,0x6c00ca6c,
+0x3ff80000,0x82796149,0x125d4cdb,
+0x3ff70000,0xee94b6d5,0x3d8fadbd,
+0x3ff70000,0xa548b3c2,0xcd0e8189,
+0x3ff70000,0xa65b8e39,0xdd6432b5,
+0x3ff70000,0xea7aff27,0xc9db0344,
+0x3ff80000,0xbd26a476,0x806aba2d,
+0x3ff90000,0xb2d7b53d,0x30fac3e0,
+0x3ffa0000,0xe8ac8fcf,0x4d3323b8,
+0x3ffd0000,0xe2e42fef,0xa39fbc79,
+};
+static long Q[36] = {
+0x3fed0000,0xd68b82a6,0xe23489f1,
+0x3ff30000,0x8fba8273,0x96b3202a,
+0x3ff60000,0xa2963484,0xfc45c183,
+0x3ff80000,0x8fb2b960,0xe201683e,
+0x3ff90000,0x895dcb41,0x1b6a721a,
+0x3ff90000,0xbaf06574,0x295f4eee,
+0x3ff90000,0xe2b3e6f2,0xc98f3ade,
+0x3ffa0000,0x89cbdb1e,0x1784d470,
+0x3ffa0000,0xaeffebc8,0xe5c1a649,
+0x3ffa0000,0xefffffde,0xa8f584c0,
+0x3ffb0000,0xbfffffff,0xf94f5506,
+0x3ffc0000,0xffffffff,0xffffd8e7,
+};
+#endif
+
+#ifdef ANSIPROT
+extern long double polevll ( long double, void *, int );
+extern long double logl ( long double );
+#else
+long double polevll(), logl();
+#endif
+
+long double ellpel(x)
+long double x;
+{
+
+if( (x <= 0.0L) || (x > 1.0L) )
+	{
+	if( x == 0.0L )
+		return( 1.0L );
+	mtherr( "ellpel", DOMAIN );
+	return( 0.0L );
+	}
+return( 1.0L + x * polevll(x,P,11) - logl(x) * (x * polevll(x,Q,11)) );
+}