123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141 |
- /* ellpe.c
- *
- * Complete elliptic integral of the second kind
- *
- *
- *
- * SYNOPSIS:
- *
- * double m1, y, ellpe();
- *
- * y = ellpe( 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
- * DEC 0, 1 13000 3.1e-17 9.4e-18
- * IEEE 0, 1 10000 2.1e-16 7.3e-17
- *
- *
- * ERROR MESSAGES:
- *
- * message condition value returned
- * ellpe domain x<0, x>1 0.0
- *
- */
- /* ellpe.c */
- /* Elliptic integral of second kind */
- /*
- Cephes Math Library, Release 2.8: June, 2000
- Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
- */
- #include "mconf.h"
- #ifdef UNK
- static double P[] = {1.53552577301013293365E-4, 2.50888492163602060990E-3,
- 8.68786816565889628429E-3, 1.07350949056076193403E-2,
- 7.77395492516787092951E-3, 7.58395289413514708519E-3,
- 1.15688436810574127319E-2, 2.18317996015557253103E-2,
- 5.68051945617860553470E-2, 4.43147180560990850618E-1,
- 1.00000000000000000299E0};
- static double Q[] = {3.27954898576485872656E-5, 1.00962792679356715133E-3,
- 6.50609489976927491433E-3, 1.68862163993311317300E-2,
- 2.61769742454493659583E-2, 3.34833904888224918614E-2,
- 4.27180926518931511717E-2, 5.85936634471101055642E-2,
- 9.37499997197644278445E-2, 2.49999999999888314361E-1};
- #endif
- #ifdef DEC
- static unsigned short P[] = {
- 0035041, 0001364, 0141572, 0117555, 0036044, 0066032, 0130027, 0033404,
- 0036416, 0053617, 0064456, 0102632, 0036457, 0161100, 0061177, 0122612,
- 0036376, 0136251, 0012403, 0124162, 0036370, 0101316, 0151715, 0131613,
- 0036475, 0105477, 0050317, 0133272, 0036662, 0154232, 0024645, 0171552,
- 0037150, 0126220, 0047054, 0030064, 0037742, 0162057, 0167645, 0165612,
- 0040200, 0000000, 0000000, 0000000};
- static unsigned short Q[] = {
- 0034411, 0106743, 0115771, 0055462, 0035604, 0052575, 0155171, 0045540,
- 0036325, 0030424, 0064332, 0167756, 0036612, 0052366, 0063006, 0115175,
- 0036726, 0070430, 0004533, 0124654, 0037011, 0022741, 0030675, 0030711,
- 0037056, 0174452, 0127062, 0132122, 0037157, 0177750, 0142041, 0072523,
- 0037277, 0177777, 0173137, 0002627, 0037577, 0177777, 0177777, 0101101};
- #endif
- #ifdef IBMPC
- static unsigned short P[] = {
- 0x53ee, 0x986f, 0x205e, 0x3f24, 0xe6e0, 0x5602, 0x8d83, 0x3f64, 0xd0b3,
- 0xed25, 0xcaf1, 0x3f81, 0xf4b1, 0x0c4f, 0xfc48, 0x3f85, 0x750e, 0x22a0,
- 0xd795, 0x3f7f, 0xb671, 0xda79, 0x1059, 0x3f7f, 0xf6d7, 0xea19, 0xb167,
- 0x3f87, 0xbe6d, 0x4534, 0x5b13, 0x3f96, 0x8607, 0x09c5, 0x1592, 0x3fad,
- 0xbd71, 0xfdf4, 0x5c85, 0x3fdc, 0x0000, 0x0000, 0x0000, 0x3ff0};
- static unsigned short Q[] = {
- 0x2b66, 0x737f, 0x31bc, 0x3f01, 0x296c, 0xbb4f, 0x8aaf, 0x3f50,
- 0x5dfe, 0x8d1b, 0xa622, 0x3f7a, 0xd350, 0xccc0, 0x4a9e, 0x3f91,
- 0x7535, 0x012b, 0xce23, 0x3f9a, 0xa639, 0x2637, 0x24bc, 0x3fa1,
- 0x568a, 0x55c6, 0xdf25, 0x3fa5, 0x2eaa, 0x1884, 0xfffd, 0x3fad,
- 0xe0b3, 0xfecb, 0xffff, 0x3fb7, 0xf048, 0xffff, 0xffff, 0x3fcf};
- #endif
- #ifdef MIEEE
- static unsigned short P[] = {
- 0x3f24, 0x205e, 0x986f, 0x53ee, 0x3f64, 0x8d83, 0x5602, 0xe6e0, 0x3f81,
- 0xcaf1, 0xed25, 0xd0b3, 0x3f85, 0xfc48, 0x0c4f, 0xf4b1, 0x3f7f, 0xd795,
- 0x22a0, 0x750e, 0x3f7f, 0x1059, 0xda79, 0xb671, 0x3f87, 0xb167, 0xea19,
- 0xf6d7, 0x3f96, 0x5b13, 0x4534, 0xbe6d, 0x3fad, 0x1592, 0x09c5, 0x8607,
- 0x3fdc, 0x5c85, 0xfdf4, 0xbd71, 0x3ff0, 0x0000, 0x0000, 0x0000};
- static unsigned short Q[] = {
- 0x3f01, 0x31bc, 0x737f, 0x2b66, 0x3f50, 0x8aaf, 0xbb4f, 0x296c,
- 0x3f7a, 0xa622, 0x8d1b, 0x5dfe, 0x3f91, 0x4a9e, 0xccc0, 0xd350,
- 0x3f9a, 0xce23, 0x012b, 0x7535, 0x3fa1, 0x24bc, 0x2637, 0xa639,
- 0x3fa5, 0xdf25, 0x55c6, 0x568a, 0x3fad, 0xfffd, 0x1884, 0x2eaa,
- 0x3fb7, 0xffff, 0xfecb, 0xe0b3, 0x3fcf, 0xffff, 0xffff, 0xf048};
- #endif
- #ifdef ANSIPROT
- extern double polevl(double, void *, int);
- extern double log(double);
- #else
- double polevl(), log();
- #endif
- double ellpe(x) double x;
- {
- if ((x <= 0.0) || (x > 1.0)) {
- if (x == 0.0)
- return (1.0);
- mtherr("ellpe", DOMAIN);
- return (0.0);
- }
- return (polevl(x, P, 10) - log(x) * (x * polevl(x, Q, 9)));
- }
|