Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/elliptic.cs @ 3839

Last change on this file since 3839 was 3839, checked in by mkommend, 14 years ago

implemented first version of LR (ticket #1012)

File size: 15.4 KB
Line 
1/*************************************************************************
2Cephes Math Library Release 2.8:  June, 2000
3Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
4
5Contributors:
6    * Sergey Bochkanov (ALGLIB project). Translation from C to
7      pseudocode.
8
9See subroutines comments for additional copyrights.
10
11>>> SOURCE LICENSE >>>
12This program is free software; you can redistribute it and/or modify
13it under the terms of the GNU General Public License as published by
14the Free Software Foundation (www.fsf.org); either version 2 of the
15License, or (at your option) any later version.
16
17This program is distributed in the hope that it will be useful,
18but WITHOUT ANY WARRANTY; without even the implied warranty of
19MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20GNU General Public License for more details.
21
22A copy of the GNU General Public License is available at
23http://www.fsf.org/licensing/licenses
24
25>>> END OF LICENSE >>>
26*************************************************************************/
27
28using System;
29
30namespace 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}
Note: See TracBrowser for help on using the repository browser.