1 /* ellpe.c
2 *
3 * Complete elliptic integral of the second kind
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double m1, y, ellpe();
10 *
11 * y = ellpe( m1 );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Approximates the integral
18 *
19 *
20 * pi/2
21 * -
22 * | | 2
23 * E(m) = | sqrt( 1 - m sin t ) dt
24 * | |
25 * -
26 * 0
27 *
28 * Where m = 1 - m1, using the approximation
29 *
30 * P(x) - x log x Q(x).
31 *
32 * Though there are no singularities, the argument m1 is used
33 * rather than m for compatibility with ellpk().
34 *
35 * E(1) = 1; E(0) = pi/2.
36 *
37 *
38 * ACCURACY:
39 *
40 * Relative error:
41 * arithmetic domain # trials peak rms
42 * DEC 0, 1 13000 3.1e-17 9.4e-18
43 * IEEE 0, 1 10000 2.1e-16 7.3e-17
44 *
45 *
46 * ERROR MESSAGES:
47 *
48 * message condition value returned
49 * ellpe domain x<0, x>1 0.0
50 *
51 */
52
53 /* ellpe.c */
54
55 /* Elliptic integral of second kind */
56
57 /*
58 Cephes Math Library, Release 2.1: February, 1989
59 Copyright 1984, 1987, 1989 by Stephen L. Moshier
60 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
61 */
62
63 #include "mconf.h"
64 #include "cephes.h"
65
66 #ifdef UNK
67 static double P[] = {
68 1.53552577301013293365E-4,
69 2.50888492163602060990E-3,
70 8.68786816565889628429E-3,
71 1.07350949056076193403E-2,
72 7.77395492516787092951E-3,
73 7.58395289413514708519E-3,
74 1.15688436810574127319E-2,
75 2.18317996015557253103E-2,
76 5.68051945617860553470E-2,
77 4.43147180560990850618E-1,
78 1.00000000000000000299E0
79 };
80 static double Q[] = {
81 3.27954898576485872656E-5,
82 1.00962792679356715133E-3,
83 6.50609489976927491433E-3,
84 1.68862163993311317300E-2,
85 2.61769742454493659583E-2,
86 3.34833904888224918614E-2,
87 4.27180926518931511717E-2,
88 5.85936634471101055642E-2,
89 9.37499997197644278445E-2,
90 2.49999999999888314361E-1
91 };
92 #endif
93
94 #ifdef DEC
95 static unsigned short P[] = {
96 0035041,0001364,0141572,0117555,
97 0036044,0066032,0130027,0033404,
98 0036416,0053617,0064456,0102632,
99 0036457,0161100,0061177,0122612,
100 0036376,0136251,0012403,0124162,
101 0036370,0101316,0151715,0131613,
102 0036475,0105477,0050317,0133272,
103 0036662,0154232,0024645,0171552,
104 0037150,0126220,0047054,0030064,
105 0037742,0162057,0167645,0165612,
106 0040200,0000000,0000000,0000000
107 };
108 static unsigned short Q[] = {
109 0034411,0106743,0115771,0055462,
110 0035604,0052575,0155171,0045540,
111 0036325,0030424,0064332,0167756,
112 0036612,0052366,0063006,0115175,
113 0036726,0070430,0004533,0124654,
114 0037011,0022741,0030675,0030711,
115 0037056,0174452,0127062,0132122,
116 0037157,0177750,0142041,0072523,
117 0037277,0177777,0173137,0002627,
118 0037577,0177777,0177777,0101101
119 };
120 #endif
121
122 #ifdef IBMPC
123 static unsigned short P[] = {
124 0x53ee,0x986f,0x205e,0x3f24,
125 0xe6e0,0x5602,0x8d83,0x3f64,
126 0xd0b3,0xed25,0xcaf1,0x3f81,
127 0xf4b1,0x0c4f,0xfc48,0x3f85,
128 0x750e,0x22a0,0xd795,0x3f7f,
129 0xb671,0xda79,0x1059,0x3f7f,
130 0xf6d7,0xea19,0xb167,0x3f87,
131 0xbe6d,0x4534,0x5b13,0x3f96,
132 0x8607,0x09c5,0x1592,0x3fad,
133 0xbd71,0xfdf4,0x5c85,0x3fdc,
134 0x0000,0x0000,0x0000,0x3ff0
135 };
136 static unsigned short Q[] = {
137 0x2b66,0x737f,0x31bc,0x3f01,
138 0x296c,0xbb4f,0x8aaf,0x3f50,
139 0x5dfe,0x8d1b,0xa622,0x3f7a,
140 0xd350,0xccc0,0x4a9e,0x3f91,
141 0x7535,0x012b,0xce23,0x3f9a,
142 0xa639,0x2637,0x24bc,0x3fa1,
143 0x568a,0x55c6,0xdf25,0x3fa5,
144 0x2eaa,0x1884,0xfffd,0x3fad,
145 0xe0b3,0xfecb,0xffff,0x3fb7,
146 0xf048,0xffff,0xffff,0x3fcf
147 };
148 #endif
149
150 #ifdef MIEEE
151 static unsigned short P[] = {
152 0x3f24,0x205e,0x986f,0x53ee,
153 0x3f64,0x8d83,0x5602,0xe6e0,
154 0x3f81,0xcaf1,0xed25,0xd0b3,
155 0x3f85,0xfc48,0x0c4f,0xf4b1,
156 0x3f7f,0xd795,0x22a0,0x750e,
157 0x3f7f,0x1059,0xda79,0xb671,
158 0x3f87,0xb167,0xea19,0xf6d7,
159 0x3f96,0x5b13,0x4534,0xbe6d,
160 0x3fad,0x1592,0x09c5,0x8607,
161 0x3fdc,0x5c85,0xfdf4,0xbd71,
162 0x3ff0,0x0000,0x0000,0x0000
163 };
164 static unsigned short Q[] = {
165 0x3f01,0x31bc,0x737f,0x2b66,
166 0x3f50,0x8aaf,0xbb4f,0x296c,
167 0x3f7a,0xa622,0x8d1b,0x5dfe,
168 0x3f91,0x4a9e,0xccc0,0xd350,
169 0x3f9a,0xce23,0x012b,0x7535,
170 0x3fa1,0x24bc,0x2637,0xa639,
171 0x3fa5,0xdf25,0x55c6,0x568a,
172 0x3fad,0xfffd,0x1884,0x2eaa,
173 0x3fb7,0xffff,0xfecb,0xe0b3,
174 0x3fcf,0xffff,0xffff,0xf048
175 };
176 #endif
177
ellpe(x)178 double ellpe(x)
179 double x;
180 {
181
182 if( (x <= 0.0) || (x > 1.0) )
183 {
184 if( x == 0.0 )
185 return( 1.0 );
186 mtherr( "ellpe", DOMAIN );
187 return( 0.0 );
188 }
189 return( polevl(x,P,10) - log(x) * (x * polevl(x,Q,9)) );
190 }
191