1 | /*************************************************************************
|
---|
2 | Cephes Math Library Release 2.8: June, 2000
|
---|
3 | Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
|
---|
4 |
|
---|
5 | Contributors:
|
---|
6 | * Sergey Bochkanov (ALGLIB project). Translation from C to
|
---|
7 | pseudocode.
|
---|
8 |
|
---|
9 | See subroutines comments for additional copyrights.
|
---|
10 |
|
---|
11 | >>> SOURCE LICENSE >>>
|
---|
12 | This program is free software; you can redistribute it and/or modify
|
---|
13 | it under the terms of the GNU General Public License as published by
|
---|
14 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
15 | License, or (at your option) any later version.
|
---|
16 |
|
---|
17 | This program is distributed in the hope that it will be useful,
|
---|
18 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
19 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
20 | GNU General Public License for more details.
|
---|
21 |
|
---|
22 | A copy of the GNU General Public License is available at
|
---|
23 | http://www.fsf.org/licensing/licenses
|
---|
24 |
|
---|
25 | >>> END OF LICENSE >>>
|
---|
26 | *************************************************************************/
|
---|
27 |
|
---|
28 | using System;
|
---|
29 |
|
---|
30 | namespace alglib
|
---|
31 | {
|
---|
32 | public class elliptic
|
---|
33 | {
|
---|
34 | /*************************************************************************
|
---|
35 | Complete elliptic integral of the first kind
|
---|
36 |
|
---|
37 | Approximates the integral
|
---|
38 |
|
---|
39 |
|
---|
40 |
|
---|
41 | pi/2
|
---|
42 | -
|
---|
43 | | |
|
---|
44 | | dt
|
---|
45 | K(m) = | ------------------
|
---|
46 | | 2
|
---|
47 | | | sqrt( 1 - m sin t )
|
---|
48 | -
|
---|
49 | 0
|
---|
50 |
|
---|
51 | using the approximation
|
---|
52 |
|
---|
53 | P(x) - log x Q(x).
|
---|
54 |
|
---|
55 | ACCURACY:
|
---|
56 |
|
---|
57 | Relative error:
|
---|
58 | arithmetic domain # trials peak rms
|
---|
59 | IEEE 0,1 30000 2.5e-16 6.8e-17
|
---|
60 |
|
---|
61 | Cephes Math Library, Release 2.8: June, 2000
|
---|
62 | Copyright 1984, 1987, 2000 by Stephen L. Moshier
|
---|
63 | *************************************************************************/
|
---|
64 | public static double ellipticintegralk(double m)
|
---|
65 | {
|
---|
66 | double result = 0;
|
---|
67 |
|
---|
68 | result = ellipticintegralkhighprecision(1.0-m);
|
---|
69 | return result;
|
---|
70 | }
|
---|
71 |
|
---|
72 |
|
---|
73 | /*************************************************************************
|
---|
74 | Complete elliptic integral of the first kind
|
---|
75 |
|
---|
76 | Approximates the integral
|
---|
77 |
|
---|
78 |
|
---|
79 |
|
---|
80 | pi/2
|
---|
81 | -
|
---|
82 | | |
|
---|
83 | | dt
|
---|
84 | K(m) = | ------------------
|
---|
85 | | 2
|
---|
86 | | | sqrt( 1 - m sin t )
|
---|
87 | -
|
---|
88 | 0
|
---|
89 |
|
---|
90 | where m = 1 - m1, using the approximation
|
---|
91 |
|
---|
92 | P(x) - log x Q(x).
|
---|
93 |
|
---|
94 | The argument m1 is used rather than m so that the logarithmic
|
---|
95 | singularity at m = 1 will be shifted to the origin; this
|
---|
96 | preserves maximum accuracy.
|
---|
97 |
|
---|
98 | K(0) = pi/2.
|
---|
99 |
|
---|
100 | ACCURACY:
|
---|
101 |
|
---|
102 | Relative error:
|
---|
103 | arithmetic domain # trials peak rms
|
---|
104 | IEEE 0,1 30000 2.5e-16 6.8e-17
|
---|
105 |
|
---|
106 | Àëãîðèòì âçÿò èç áèáëèîòåêè Cephes
|
---|
107 | *************************************************************************/
|
---|
108 | public static double ellipticintegralkhighprecision(double m1)
|
---|
109 | {
|
---|
110 | double result = 0;
|
---|
111 | double p = 0;
|
---|
112 | double q = 0;
|
---|
113 |
|
---|
114 | if( (double)(m1)<=(double)(AP.Math.MachineEpsilon) )
|
---|
115 | {
|
---|
116 | result = 1.3862943611198906188E0-0.5*Math.Log(m1);
|
---|
117 | }
|
---|
118 | else
|
---|
119 | {
|
---|
120 | p = 1.37982864606273237150E-4;
|
---|
121 | p = p*m1+2.28025724005875567385E-3;
|
---|
122 | p = p*m1+7.97404013220415179367E-3;
|
---|
123 | p = p*m1+9.85821379021226008714E-3;
|
---|
124 | p = p*m1+6.87489687449949877925E-3;
|
---|
125 | p = p*m1+6.18901033637687613229E-3;
|
---|
126 | p = p*m1+8.79078273952743772254E-3;
|
---|
127 | p = p*m1+1.49380448916805252718E-2;
|
---|
128 | p = p*m1+3.08851465246711995998E-2;
|
---|
129 | p = p*m1+9.65735902811690126535E-2;
|
---|
130 | p = p*m1+1.38629436111989062502E0;
|
---|
131 | q = 2.94078955048598507511E-5;
|
---|
132 | q = q*m1+9.14184723865917226571E-4;
|
---|
133 | q = q*m1+5.94058303753167793257E-3;
|
---|
134 | q = q*m1+1.54850516649762399335E-2;
|
---|
135 | q = q*m1+2.39089602715924892727E-2;
|
---|
136 | q = q*m1+3.01204715227604046988E-2;
|
---|
137 | q = q*m1+3.73774314173823228969E-2;
|
---|
138 | q = q*m1+4.88280347570998239232E-2;
|
---|
139 | q = q*m1+7.03124996963957469739E-2;
|
---|
140 | q = q*m1+1.24999999999870820058E-1;
|
---|
141 | q = q*m1+4.99999999999999999821E-1;
|
---|
142 | result = p-q*Math.Log(m1);
|
---|
143 | }
|
---|
144 | return result;
|
---|
145 | }
|
---|
146 |
|
---|
147 |
|
---|
148 | /*************************************************************************
|
---|
149 | Incomplete elliptic integral of the first kind F(phi|m)
|
---|
150 |
|
---|
151 | Approximates the integral
|
---|
152 |
|
---|
153 |
|
---|
154 |
|
---|
155 | phi
|
---|
156 | -
|
---|
157 | | |
|
---|
158 | | dt
|
---|
159 | F(phi_\m) = | ------------------
|
---|
160 | | 2
|
---|
161 | | | sqrt( 1 - m sin t )
|
---|
162 | -
|
---|
163 | 0
|
---|
164 |
|
---|
165 | of amplitude phi and modulus m, using the arithmetic -
|
---|
166 | geometric mean algorithm.
|
---|
167 |
|
---|
168 |
|
---|
169 |
|
---|
170 |
|
---|
171 | ACCURACY:
|
---|
172 |
|
---|
173 | Tested at random points with m in [0, 1] and phi as indicated.
|
---|
174 |
|
---|
175 | Relative error:
|
---|
176 | arithmetic domain # trials peak rms
|
---|
177 | IEEE -10,10 200000 7.4e-16 1.0e-16
|
---|
178 |
|
---|
179 | Cephes Math Library Release 2.8: June, 2000
|
---|
180 | Copyright 1984, 1987, 2000 by Stephen L. Moshier
|
---|
181 | *************************************************************************/
|
---|
182 | public static double incompleteellipticintegralk(double phi,
|
---|
183 | double m)
|
---|
184 | {
|
---|
185 | double result = 0;
|
---|
186 | double a = 0;
|
---|
187 | double b = 0;
|
---|
188 | double c = 0;
|
---|
189 | double e = 0;
|
---|
190 | double temp = 0;
|
---|
191 | double pio2 = 0;
|
---|
192 | double t = 0;
|
---|
193 | double k = 0;
|
---|
194 | int d = 0;
|
---|
195 | int md = 0;
|
---|
196 | int s = 0;
|
---|
197 | int npio2 = 0;
|
---|
198 |
|
---|
199 | pio2 = 1.57079632679489661923;
|
---|
200 | if( (double)(m)==(double)(0) )
|
---|
201 | {
|
---|
202 | result = phi;
|
---|
203 | return result;
|
---|
204 | }
|
---|
205 | a = 1-m;
|
---|
206 | if( (double)(a)==(double)(0) )
|
---|
207 | {
|
---|
208 | result = Math.Log(Math.Tan(0.5*(pio2+phi)));
|
---|
209 | return result;
|
---|
210 | }
|
---|
211 | npio2 = (int)Math.Floor(phi/pio2);
|
---|
212 | if( npio2%2!=0 )
|
---|
213 | {
|
---|
214 | npio2 = npio2+1;
|
---|
215 | }
|
---|
216 | if( npio2!=0 )
|
---|
217 | {
|
---|
218 | k = ellipticintegralk(1-a);
|
---|
219 | phi = phi-npio2*pio2;
|
---|
220 | }
|
---|
221 | else
|
---|
222 | {
|
---|
223 | k = 0;
|
---|
224 | }
|
---|
225 | if( (double)(phi)<(double)(0) )
|
---|
226 | {
|
---|
227 | phi = -phi;
|
---|
228 | s = -1;
|
---|
229 | }
|
---|
230 | else
|
---|
231 | {
|
---|
232 | s = 0;
|
---|
233 | }
|
---|
234 | b = Math.Sqrt(a);
|
---|
235 | t = Math.Tan(phi);
|
---|
236 | if( (double)(Math.Abs(t))>(double)(10) )
|
---|
237 | {
|
---|
238 | e = 1.0/(b*t);
|
---|
239 | if( (double)(Math.Abs(e))<(double)(10) )
|
---|
240 | {
|
---|
241 | e = Math.Atan(e);
|
---|
242 | if( npio2==0 )
|
---|
243 | {
|
---|
244 | k = ellipticintegralk(1-a);
|
---|
245 | }
|
---|
246 | temp = k-incompleteellipticintegralk(e, m);
|
---|
247 | if( s<0 )
|
---|
248 | {
|
---|
249 | temp = -temp;
|
---|
250 | }
|
---|
251 | result = temp+npio2*k;
|
---|
252 | return result;
|
---|
253 | }
|
---|
254 | }
|
---|
255 | a = 1.0;
|
---|
256 | c = Math.Sqrt(m);
|
---|
257 | d = 1;
|
---|
258 | md = 0;
|
---|
259 | while( (double)(Math.Abs(c/a))>(double)(AP.Math.MachineEpsilon) )
|
---|
260 | {
|
---|
261 | temp = b/a;
|
---|
262 | phi = phi+Math.Atan(t*temp)+md*Math.PI;
|
---|
263 | md = (int)((phi+pio2)/Math.PI);
|
---|
264 | t = t*(1.0+temp)/(1.0-temp*t*t);
|
---|
265 | c = 0.5*(a-b);
|
---|
266 | temp = Math.Sqrt(a*b);
|
---|
267 | a = 0.5*(a+b);
|
---|
268 | b = temp;
|
---|
269 | d = d+d;
|
---|
270 | }
|
---|
271 | temp = (Math.Atan(t)+md*Math.PI)/(d*a);
|
---|
272 | if( s<0 )
|
---|
273 | {
|
---|
274 | temp = -temp;
|
---|
275 | }
|
---|
276 | result = temp+npio2*k;
|
---|
277 | return result;
|
---|
278 | }
|
---|
279 |
|
---|
280 |
|
---|
281 | /*************************************************************************
|
---|
282 | Complete elliptic integral of the second kind
|
---|
283 |
|
---|
284 | Approximates the integral
|
---|
285 |
|
---|
286 |
|
---|
287 | pi/2
|
---|
288 | -
|
---|
289 | | | 2
|
---|
290 | E(m) = | sqrt( 1 - m sin t ) dt
|
---|
291 | | |
|
---|
292 | -
|
---|
293 | 0
|
---|
294 |
|
---|
295 | using the approximation
|
---|
296 |
|
---|
297 | P(x) - x log x Q(x).
|
---|
298 |
|
---|
299 | ACCURACY:
|
---|
300 |
|
---|
301 | Relative error:
|
---|
302 | arithmetic domain # trials peak rms
|
---|
303 | IEEE 0, 1 10000 2.1e-16 7.3e-17
|
---|
304 |
|
---|
305 | Cephes Math Library, Release 2.8: June, 2000
|
---|
306 | Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
|
---|
307 | *************************************************************************/
|
---|
308 | public static double ellipticintegrale(double m)
|
---|
309 | {
|
---|
310 | double result = 0;
|
---|
311 | double p = 0;
|
---|
312 | double q = 0;
|
---|
313 |
|
---|
314 | System.Diagnostics.Debug.Assert((double)(m)>=(double)(0) & (double)(m)<=(double)(1), "Domain error in EllipticIntegralE: m<0 or m>1");
|
---|
315 | m = 1-m;
|
---|
316 | if( (double)(m)==(double)(0) )
|
---|
317 | {
|
---|
318 | result = 1;
|
---|
319 | return result;
|
---|
320 | }
|
---|
321 | p = 1.53552577301013293365E-4;
|
---|
322 | p = p*m+2.50888492163602060990E-3;
|
---|
323 | p = p*m+8.68786816565889628429E-3;
|
---|
324 | p = p*m+1.07350949056076193403E-2;
|
---|
325 | p = p*m+7.77395492516787092951E-3;
|
---|
326 | p = p*m+7.58395289413514708519E-3;
|
---|
327 | p = p*m+1.15688436810574127319E-2;
|
---|
328 | p = p*m+2.18317996015557253103E-2;
|
---|
329 | p = p*m+5.68051945617860553470E-2;
|
---|
330 | p = p*m+4.43147180560990850618E-1;
|
---|
331 | p = p*m+1.00000000000000000299E0;
|
---|
332 | q = 3.27954898576485872656E-5;
|
---|
333 | q = q*m+1.00962792679356715133E-3;
|
---|
334 | q = q*m+6.50609489976927491433E-3;
|
---|
335 | q = q*m+1.68862163993311317300E-2;
|
---|
336 | q = q*m+2.61769742454493659583E-2;
|
---|
337 | q = q*m+3.34833904888224918614E-2;
|
---|
338 | q = q*m+4.27180926518931511717E-2;
|
---|
339 | q = q*m+5.85936634471101055642E-2;
|
---|
340 | q = q*m+9.37499997197644278445E-2;
|
---|
341 | q = q*m+2.49999999999888314361E-1;
|
---|
342 | result = p-q*m*Math.Log(m);
|
---|
343 | return result;
|
---|
344 | }
|
---|
345 |
|
---|
346 |
|
---|
347 | /*************************************************************************
|
---|
348 | Incomplete elliptic integral of the second kind
|
---|
349 |
|
---|
350 | Approximates the integral
|
---|
351 |
|
---|
352 |
|
---|
353 | phi
|
---|
354 | -
|
---|
355 | | |
|
---|
356 | | 2
|
---|
357 | E(phi_\m) = | sqrt( 1 - m sin t ) dt
|
---|
358 | |
|
---|
359 | | |
|
---|
360 | -
|
---|
361 | 0
|
---|
362 |
|
---|
363 | of amplitude phi and modulus m, using the arithmetic -
|
---|
364 | geometric mean algorithm.
|
---|
365 |
|
---|
366 | ACCURACY:
|
---|
367 |
|
---|
368 | Tested at random arguments with phi in [-10, 10] and m in
|
---|
369 | [0, 1].
|
---|
370 | Relative error:
|
---|
371 | arithmetic domain # trials peak rms
|
---|
372 | IEEE -10,10 150000 3.3e-15 1.4e-16
|
---|
373 |
|
---|
374 | Cephes Math Library Release 2.8: June, 2000
|
---|
375 | Copyright 1984, 1987, 1993, 2000 by Stephen L. Moshier
|
---|
376 | *************************************************************************/
|
---|
377 | public static double incompleteellipticintegrale(double phi,
|
---|
378 | double m)
|
---|
379 | {
|
---|
380 | double result = 0;
|
---|
381 | double pio2 = 0;
|
---|
382 | double a = 0;
|
---|
383 | double b = 0;
|
---|
384 | double c = 0;
|
---|
385 | double e = 0;
|
---|
386 | double temp = 0;
|
---|
387 | double lphi = 0;
|
---|
388 | double t = 0;
|
---|
389 | double ebig = 0;
|
---|
390 | int d = 0;
|
---|
391 | int md = 0;
|
---|
392 | int npio2 = 0;
|
---|
393 | int s = 0;
|
---|
394 |
|
---|
395 | pio2 = 1.57079632679489661923;
|
---|
396 | if( (double)(m)==(double)(0) )
|
---|
397 | {
|
---|
398 | result = phi;
|
---|
399 | return result;
|
---|
400 | }
|
---|
401 | lphi = phi;
|
---|
402 | npio2 = (int)Math.Floor(lphi/pio2);
|
---|
403 | if( npio2%2!=0 )
|
---|
404 | {
|
---|
405 | npio2 = npio2+1;
|
---|
406 | }
|
---|
407 | lphi = lphi-npio2*pio2;
|
---|
408 | if( (double)(lphi)<(double)(0) )
|
---|
409 | {
|
---|
410 | lphi = -lphi;
|
---|
411 | s = -1;
|
---|
412 | }
|
---|
413 | else
|
---|
414 | {
|
---|
415 | s = 1;
|
---|
416 | }
|
---|
417 | a = 1.0-m;
|
---|
418 | ebig = ellipticintegrale(m);
|
---|
419 | if( (double)(a)==(double)(0) )
|
---|
420 | {
|
---|
421 | temp = Math.Sin(lphi);
|
---|
422 | if( s<0 )
|
---|
423 | {
|
---|
424 | temp = -temp;
|
---|
425 | }
|
---|
426 | result = temp+npio2*ebig;
|
---|
427 | return result;
|
---|
428 | }
|
---|
429 | t = Math.Tan(lphi);
|
---|
430 | b = Math.Sqrt(a);
|
---|
431 |
|
---|
432 | //
|
---|
433 | // Thanks to Brian Fitzgerald <fitzgb@mml0.meche.rpi.edu>
|
---|
434 | // for pointing out an instability near odd multiples of pi/2
|
---|
435 | //
|
---|
436 | if( (double)(Math.Abs(t))>(double)(10) )
|
---|
437 | {
|
---|
438 |
|
---|
439 | //
|
---|
440 | // Transform the amplitude
|
---|
441 | //
|
---|
442 | e = 1.0/(b*t);
|
---|
443 |
|
---|
444 | //
|
---|
445 | // ... but avoid multiple recursions.
|
---|
446 | //
|
---|
447 | if( (double)(Math.Abs(e))<(double)(10) )
|
---|
448 | {
|
---|
449 | e = Math.Atan(e);
|
---|
450 | temp = ebig+m*Math.Sin(lphi)*Math.Sin(e)-incompleteellipticintegrale(e, m);
|
---|
451 | if( s<0 )
|
---|
452 | {
|
---|
453 | temp = -temp;
|
---|
454 | }
|
---|
455 | result = temp+npio2*ebig;
|
---|
456 | return result;
|
---|
457 | }
|
---|
458 | }
|
---|
459 | c = Math.Sqrt(m);
|
---|
460 | a = 1.0;
|
---|
461 | d = 1;
|
---|
462 | e = 0.0;
|
---|
463 | md = 0;
|
---|
464 | while( (double)(Math.Abs(c/a))>(double)(AP.Math.MachineEpsilon) )
|
---|
465 | {
|
---|
466 | temp = b/a;
|
---|
467 | lphi = lphi+Math.Atan(t*temp)+md*Math.PI;
|
---|
468 | md = (int)((lphi+pio2)/Math.PI);
|
---|
469 | t = t*(1.0+temp)/(1.0-temp*t*t);
|
---|
470 | c = 0.5*(a-b);
|
---|
471 | temp = Math.Sqrt(a*b);
|
---|
472 | a = 0.5*(a+b);
|
---|
473 | b = temp;
|
---|
474 | d = d+d;
|
---|
475 | e = e+c*Math.Sin(lphi);
|
---|
476 | }
|
---|
477 | temp = ebig/ellipticintegralk(m);
|
---|
478 | temp = temp*((Math.Atan(t)+md*Math.PI)/(d*a));
|
---|
479 | temp = temp+e;
|
---|
480 | if( s<0 )
|
---|
481 | {
|
---|
482 | temp = -temp;
|
---|
483 | }
|
---|
484 | result = temp+npio2*ebig;
|
---|
485 | return result;
|
---|
486 | return result;
|
---|
487 | }
|
---|
488 | }
|
---|
489 | }
|
---|