1 /* ellpk.c
2 *
3 * Complete elliptic integral of the first kind
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double m1, y, ellpk();
10 *
11 * y = ellpk( m1 );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Approximates the integral
18 *
19 *
20 *
21 * pi/2
22 * -
23 * | |
24 * | dt
25 * K(m) = | ------------------
26 * | 2
27 * | | sqrt( 1 - m sin t )
28 * -
29 * 0
30 *
31 * where m = 1 - m1, using the approximation
32 *
33 * P(x) - log x Q(x).
34 *
35 * The argument m1 is used rather than m so that the logarithmic
36 * singularity at m = 1 will be shifted to the origin; this
37 * preserves maximum accuracy.
38 *
39 * K(0) = pi/2.
40 *
41 * ACCURACY:
42 *
43 * Relative error:
44 * arithmetic domain # trials peak rms
45 * DEC 0,1 16000 3.5e-17 1.1e-17
46 * IEEE 0,1 30000 2.5e-16 6.8e-17
47 *
48 * ERROR MESSAGES:
49 *
50 * message condition value returned
51 * ellpk domain x<0, x>1 0.0
52 *
53 */
54
55 /* ellpk.c */
56
57
58 /*
59 Cephes Math Library, Release 2.0: April, 1987
60 Copyright 1984, 1987 by Stephen L. Moshier
61 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
62 */
63
64 #include "mconf.h"
65 #include "cephes.h"
66
67
68 #ifdef DEC
69 static unsigned short P[] =
70 {
71 0035020,0127576,0040430,0051544,
72 0036025,0070136,0042703,0153716,
73 0036402,0122614,0062555,0077777,
74 0036441,0102130,0072334,0025172,
75 0036341,0043320,0117242,0172076,
76 0036312,0146456,0077242,0154141,
77 0036420,0003467,0013727,0035407,
78 0036564,0137263,0110651,0020237,
79 0036775,0001330,0144056,0020305,
80 0037305,0144137,0157521,0141734,
81 0040261,0071027,0173721,0147572
82 };
83 static unsigned short Q[] =
84 {
85 0034366,0130371,0103453,0077633,
86 0035557,0122745,0173515,0113016,
87 0036302,0124470,0167304,0074473,
88 0036575,0132403,0117226,0117576,
89 0036703,0156271,0047124,0147733,
90 0036766,0137465,0002053,0157312,
91 0037031,0014423,0154274,0176515,
92 0037107,0177747,0143216,0016145,
93 0037217,0177777,0172621,0074000,
94 0037377,0177777,0177776,0156435,
95 0040000,0000000,0000000,0000000
96 };
97 static union us2d_t ac1 = {{0040261,0071027,0173721,0147572}};
98 #define C1 ac1.d
99 #endif
100
101 #ifdef IBMPC
102 static unsigned short P[] =
103 {
104 0x0a6d,0xc823,0x15ef,0x3f22,
105 0x7afa,0xc8b8,0xae0b,0x3f62,
106 0xb000,0x8cad,0x54b1,0x3f80,
107 0x854f,0x0e9b,0x308b,0x3f84,
108 0x5e88,0x13d4,0x28da,0x3f7c,
109 0x5b0c,0xcfd4,0x59a5,0x3f79,
110 0xe761,0xe2fa,0x00e6,0x3f82,
111 0x2414,0x7235,0x97d6,0x3f8e,
112 0xc419,0x1905,0xa05b,0x3f9f,
113 0x387c,0xfbea,0xb90b,0x3fb8,
114 0x39ef,0xfefa,0x2e42,0x3ff6
115 };
116 static unsigned short Q[] =
117 {
118 0x6ff3,0x30e5,0xd61f,0x3efe,
119 0xb2c2,0xbee9,0xf4bc,0x3f4d,
120 0x8f27,0x1dd8,0x5527,0x3f78,
121 0xd3f0,0x73d2,0xb6a0,0x3f8f,
122 0x99fb,0x29ca,0x7b97,0x3f98,
123 0x7bd9,0xa085,0xd7e6,0x3f9e,
124 0x9faa,0x7b17,0x2322,0x3fa3,
125 0xc38d,0xf8d1,0xfffc,0x3fa8,
126 0x2f00,0xfeb2,0xffff,0x3fb1,
127 0xdba4,0xffff,0xffff,0x3fbf,
128 0x0000,0x0000,0x0000,0x3fe0
129 };
130 static union us2d_t ac1 = {{0x39ef,0xfefa,0x2e42,0x3ff6}};
131 #define C1 ac1.d
132 #endif
133
134 #ifdef MIEEE
135 static unsigned short P[] =
136 {
137 0x3f22,0x15ef,0xc823,0x0a6d,
138 0x3f62,0xae0b,0xc8b8,0x7afa,
139 0x3f80,0x54b1,0x8cad,0xb000,
140 0x3f84,0x308b,0x0e9b,0x854f,
141 0x3f7c,0x28da,0x13d4,0x5e88,
142 0x3f79,0x59a5,0xcfd4,0x5b0c,
143 0x3f82,0x00e6,0xe2fa,0xe761,
144 0x3f8e,0x97d6,0x7235,0x2414,
145 0x3f9f,0xa05b,0x1905,0xc419,
146 0x3fb8,0xb90b,0xfbea,0x387c,
147 0x3ff6,0x2e42,0xfefa,0x39ef
148 };
149 static unsigned short Q[] =
150 {
151 0x3efe,0xd61f,0x30e5,0x6ff3,
152 0x3f4d,0xf4bc,0xbee9,0xb2c2,
153 0x3f78,0x5527,0x1dd8,0x8f27,
154 0x3f8f,0xb6a0,0x73d2,0xd3f0,
155 0x3f98,0x7b97,0x29ca,0x99fb,
156 0x3f9e,0xd7e6,0xa085,0x7bd9,
157 0x3fa3,0x2322,0x7b17,0x9faa,
158 0x3fa8,0xfffc,0xf8d1,0xc38d,
159 0x3fb1,0xffff,0xfeb2,0x2f00,
160 0x3fbf,0xffff,0xffff,0xdba4,
161 0x3fe0,0x0000,0x0000,0x0000
162 };
163 static union us2d_t ac1 = { {0x3ff6,0x2e42,0xfefa,0x39ef} };
164 #define C1 ac1.d
165 #endif
166
167 #ifdef UNK
168 static double P[] =
169 {
170 1.37982864606273237150E-4,
171 2.28025724005875567385E-3,
172 7.97404013220415179367E-3,
173 9.85821379021226008714E-3,
174 6.87489687449949877925E-3,
175 6.18901033637687613229E-3,
176 8.79078273952743772254E-3,
177 1.49380448916805252718E-2,
178 3.08851465246711995998E-2,
179 9.65735902811690126535E-2,
180 1.38629436111989062502E0
181 };
182
183 static double Q[] =
184 {
185 2.94078955048598507511E-5,
186 9.14184723865917226571E-4,
187 5.94058303753167793257E-3,
188 1.54850516649762399335E-2,
189 2.39089602715924892727E-2,
190 3.01204715227604046988E-2,
191 3.73774314173823228969E-2,
192 4.88280347570998239232E-2,
193 7.03124996963957469739E-2,
194 1.24999999999870820058E-1,
195 4.99999999999999999821E-1
196 };
197 static double C1 = 1.3862943611198906188E0; /* log(4) */
198 #endif
199
200 extern double MACHEP, MAXNUM;
201
ellpk(x)202 double ellpk(x)
203 double x;
204 {
205
206 if( (x < 0.0) || (x > 1.0) )
207 {
208 mtherr( "ellpk", DOMAIN );
209 return( 0.0 );
210 }
211
212 if( x > MACHEP )
213 {
214 return( polevl(x,P,10) - log(x) * polevl(x,Q,10) );
215 }
216 else
217 {
218 if( x == 0.0 )
219 {
220 mtherr( "ellpk", SING );
221 return( MAXNUM );
222 }
223 else
224 {
225 return( C1 - 0.5 * log(x) );
226 }
227 }
228 }
229