Quelle  ellpel.c

/* 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 "mconf.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)) );
}