Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/expintegrals.cs @ 3932

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

implemented first version of LR (ticket #1012)
File size: 15.0 KB
Line 
1/*************************************************************************
2Cephes Math Library Release 2.8:  June, 2000
3Copyright 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 expintegrals
33    {
34        /*************************************************************************
35        Exponential integral Ei(x)
36
37                      x
38                       -     t
39                      | |   e
40           Ei(x) =   -|-   ---  dt .
41                    | |     t
42                     -
43                    -inf
44
45        Not defined for x <= 0.
46        See also expn.c.
47
48
49
50        ACCURACY:
51
52                             Relative error:
53        arithmetic   domain     # trials      peak         rms
54           IEEE       0,100       50000      8.6e-16     1.3e-16
55
56        Cephes Math Library Release 2.8:  May, 1999
57        Copyright 1999 by Stephen L. Moshier
58        *************************************************************************/
59        public static double exponentialintegralei(double x)
60        {
61            double result = 0;
62            double eul = 0;
63            double f = 0;
64            double f1 = 0;
65            double f2 = 0;
66            double w = 0;
67
68            eul = 0.5772156649015328606065;
69            if( (double)(x)<=(double)(0) )
70            {
71                result = 0;
72                return result;
73            }
74            if( (double)(x)<(double)(2) )
75            {
76                f1 = -5.350447357812542947283;
77                f1 = f1*x+218.5049168816613393830;
78                f1 = f1*x-4176.572384826693777058;
79                f1 = f1*x+55411.76756393557601232;
80                f1 = f1*x-331338.1331178144034309;
81                f1 = f1*x+1592627.163384945414220;
82                f2 = 1.000000000000000000000;
83                f2 = f2*x-52.50547959112862969197;
84                f2 = f2*x+1259.616186786790571525;
85                f2 = f2*x-17565.49581973534652631;
86                f2 = f2*x+149306.2117002725991967;
87                f2 = f2*x-729494.9239640527645655;
88                f2 = f2*x+1592627.163384945429726;
89                f = f1/f2;
90                result = eul+Math.Log(x)+x*f;
91                return result;
92            }
93            if( (double)(x)<(double)(4) )
94            {
95                w = 1/x;
96                f1 = 1.981808503259689673238E-2;
97                f1 = f1*w-1.271645625984917501326;
98                f1 = f1*w-2.088160335681228318920;
99                f1 = f1*w+2.755544509187936721172;
100                f1 = f1*w-4.409507048701600257171E-1;
101                f1 = f1*w+4.665623805935891391017E-2;
102                f1 = f1*w-1.545042679673485262580E-3;
103                f1 = f1*w+7.059980605299617478514E-5;
104                f2 = 1.000000000000000000000;
105                f2 = f2*w+1.476498670914921440652;
106                f2 = f2*w+5.629177174822436244827E-1;
107                f2 = f2*w+1.699017897879307263248E-1;
108                f2 = f2*w+2.291647179034212017463E-2;
109                f2 = f2*w+4.450150439728752875043E-3;
110                f2 = f2*w+1.727439612206521482874E-4;
111                f2 = f2*w+3.953167195549672482304E-5;
112                f = f1/f2;
113                result = Math.Exp(x)*w*(1+w*f);
114                return result;
115            }
116            if( (double)(x)<(double)(8) )
117            {
118                w = 1/x;
119                f1 = -1.373215375871208729803;
120                f1 = f1*w-7.084559133740838761406E-1;
121                f1 = f1*w+1.580806855547941010501;
122                f1 = f1*w-2.601500427425622944234E-1;
123                f1 = f1*w+2.994674694113713763365E-2;
124                f1 = f1*w-1.038086040188744005513E-3;
125                f1 = f1*w+4.371064420753005429514E-5;
126                f1 = f1*w+2.141783679522602903795E-6;
127                f2 = 1.000000000000000000000;
128                f2 = f2*w+8.585231423622028380768E-1;
129                f2 = f2*w+4.483285822873995129957E-1;
130                f2 = f2*w+7.687932158124475434091E-2;
131                f2 = f2*w+2.449868241021887685904E-2;
132                f2 = f2*w+8.832165941927796567926E-4;
133                f2 = f2*w+4.590952299511353531215E-4;
134                f2 = f2*w+-4.729848351866523044863E-6;
135                f2 = f2*w+2.665195537390710170105E-6;
136                f = f1/f2;
137                result = Math.Exp(x)*w*(1+w*f);
138                return result;
139            }
140            if( (double)(x)<(double)(16) )
141            {
142                w = 1/x;
143                f1 = -2.106934601691916512584;
144                f1 = f1*w+1.732733869664688041885;
145                f1 = f1*w-2.423619178935841904839E-1;
146                f1 = f1*w+2.322724180937565842585E-2;
147                f1 = f1*w+2.372880440493179832059E-4;
148                f1 = f1*w-8.343219561192552752335E-5;
149                f1 = f1*w+1.363408795605250394881E-5;
150                f1 = f1*w-3.655412321999253963714E-7;
151                f1 = f1*w+1.464941733975961318456E-8;
152                f1 = f1*w+6.176407863710360207074E-10;
153                f2 = 1.000000000000000000000;
154                f2 = f2*w-2.298062239901678075778E-1;
155                f2 = f2*w+1.105077041474037862347E-1;
156                f2 = f2*w-1.566542966630792353556E-2;
157                f2 = f2*w+2.761106850817352773874E-3;
158                f2 = f2*w-2.089148012284048449115E-4;
159                f2 = f2*w+1.708528938807675304186E-5;
160                f2 = f2*w-4.459311796356686423199E-7;
161                f2 = f2*w+1.394634930353847498145E-8;
162                f2 = f2*w+6.150865933977338354138E-10;
163                f = f1/f2;
164                result = Math.Exp(x)*w*(1+w*f);
165                return result;
166            }
167            if( (double)(x)<(double)(32) )
168            {
169                w = 1/x;
170                f1 = -2.458119367674020323359E-1;
171                f1 = f1*w-1.483382253322077687183E-1;
172                f1 = f1*w+7.248291795735551591813E-2;
173                f1 = f1*w-1.348315687380940523823E-2;
174                f1 = f1*w+1.342775069788636972294E-3;
175                f1 = f1*w-7.942465637159712264564E-5;
176                f1 = f1*w+2.644179518984235952241E-6;
177                f1 = f1*w-4.239473659313765177195E-8;
178                f2 = 1.000000000000000000000;
179                f2 = f2*w-1.044225908443871106315E-1;
180                f2 = f2*w-2.676453128101402655055E-1;
181                f2 = f2*w+9.695000254621984627876E-2;
182                f2 = f2*w-1.601745692712991078208E-2;
183                f2 = f2*w+1.496414899205908021882E-3;
184                f2 = f2*w-8.462452563778485013756E-5;
185                f2 = f2*w+2.728938403476726394024E-6;
186                f2 = f2*w-4.239462431819542051337E-8;
187                f = f1/f2;
188                result = Math.Exp(x)*w*(1+w*f);
189                return result;
190            }
191            if( (double)(x)<(double)(64) )
192            {
193                w = 1/x;
194                f1 = 1.212561118105456670844E-1;
195                f1 = f1*w-5.823133179043894485122E-1;
196                f1 = f1*w+2.348887314557016779211E-1;
197                f1 = f1*w-3.040034318113248237280E-2;
198                f1 = f1*w+1.510082146865190661777E-3;
199                f1 = f1*w-2.523137095499571377122E-5;
200                f2 = 1.000000000000000000000;
201                f2 = f2*w-1.002252150365854016662;
202                f2 = f2*w+2.928709694872224144953E-1;
203                f2 = f2*w-3.337004338674007801307E-2;
204                f2 = f2*w+1.560544881127388842819E-3;
205                f2 = f2*w-2.523137093603234562648E-5;
206                f = f1/f2;
207                result = Math.Exp(x)*w*(1+w*f);
208                return result;
209            }
210            w = 1/x;
211            f1 = -7.657847078286127362028E-1;
212            f1 = f1*w+6.886192415566705051750E-1;
213            f1 = f1*w-2.132598113545206124553E-1;
214            f1 = f1*w+3.346107552384193813594E-2;
215            f1 = f1*w-3.076541477344756050249E-3;
216            f1 = f1*w+1.747119316454907477380E-4;
217            f1 = f1*w-6.103711682274170530369E-6;
218            f1 = f1*w+1.218032765428652199087E-7;
219            f1 = f1*w-1.086076102793290233007E-9;
220            f2 = 1.000000000000000000000;
221            f2 = f2*w-1.888802868662308731041;
222            f2 = f2*w+1.066691687211408896850;
223            f2 = f2*w-2.751915982306380647738E-1;
224            f2 = f2*w+3.930852688233823569726E-2;
225            f2 = f2*w-3.414684558602365085394E-3;
226            f2 = f2*w+1.866844370703555398195E-4;
227            f2 = f2*w-6.345146083130515357861E-6;
228            f2 = f2*w+1.239754287483206878024E-7;
229            f2 = f2*w-1.086076102793126632978E-9;
230            f = f1/f2;
231            result = Math.Exp(x)*w*(1+w*f);
232            return result;
233        }
234
235
236        /*************************************************************************
237        Exponential integral En(x)
238
239        Evaluates the exponential integral
240
241                        inf.
242                          -
243                         | |   -xt
244                         |    e
245             E (x)  =    |    ----  dt.
246              n          |      n
247                       | |     t
248                        -
249                         1
250
251
252        Both n and x must be nonnegative.
253
254        The routine employs either a power series, a continued
255        fraction, or an asymptotic formula depending on the
256        relative values of n and x.
257
258        ACCURACY:
259
260                             Relative error:
261        arithmetic   domain     # trials      peak         rms
262           IEEE      0, 30       10000       1.7e-15     3.6e-16
263
264        Cephes Math Library Release 2.8:  June, 2000
265        Copyright 1985, 2000 by Stephen L. Moshier
266        *************************************************************************/
267        public static double exponentialintegralen(double x,
268            int n)
269        {
270            double result = 0;
271            double r = 0;
272            double t = 0;
273            double yk = 0;
274            double xk = 0;
275            double pk = 0;
276            double pkm1 = 0;
277            double pkm2 = 0;
278            double qk = 0;
279            double qkm1 = 0;
280            double qkm2 = 0;
281            double psi = 0;
282            double z = 0;
283            int i = 0;
284            int k = 0;
285            double big = 0;
286            double eul = 0;
287
288            eul = 0.57721566490153286060;
289            big = 1.44115188075855872*Math.Pow(10, 17);
290            if( n<0 | (double)(x)<(double)(0) | (double)(x)>(double)(170) | (double)(x)==(double)(0) & n<2 )
291            {
292                result = -1;
293                return result;
294            }
295            if( (double)(x)==(double)(0) )
296            {
297                result = (double)(1)/((double)(n-1));
298                return result;
299            }
300            if( n==0 )
301            {
302                result = Math.Exp(-x)/x;
303                return result;
304            }
305            if( n>5000 )
306            {
307                xk = x+n;
308                yk = 1/(xk*xk);
309                t = n;
310                result = yk*t*(6*x*x-8*t*x+t*t);
311                result = yk*(result+t*(t-2.0*x));
312                result = yk*(result+t);
313                result = (result+1)*Math.Exp(-x)/xk;
314                return result;
315            }
316            if( (double)(x)<=(double)(1) )
317            {
318                psi = -eul-Math.Log(x);
319                for(i=1; i<=n-1; i++)
320                {
321                    psi = psi+(double)(1)/(double)(i);
322                }
323                z = -x;
324                xk = 0;
325                yk = 1;
326                pk = 1-n;
327                if( n==1 )
328                {
329                    result = 0.0;
330                }
331                else
332                {
333                    result = 1.0/pk;
334                }
335                do
336                {
337                    xk = xk+1;
338                    yk = yk*z/xk;
339                    pk = pk+1;
340                    if( (double)(pk)!=(double)(0) )
341                    {
342                        result = result+yk/pk;
343                    }
344                    if( (double)(result)!=(double)(0) )
345                    {
346                        t = Math.Abs(yk/result);
347                    }
348                    else
349                    {
350                        t = 1;
351                    }
352                }
353                while( (double)(t)>=(double)(AP.Math.MachineEpsilon) );
354                t = 1;
355                for(i=1; i<=n-1; i++)
356                {
357                    t = t*z/i;
358                }
359                result = psi*t-result;
360                return result;
361            }
362            else
363            {
364                k = 1;
365                pkm2 = 1;
366                qkm2 = x;
367                pkm1 = 1.0;
368                qkm1 = x+n;
369                result = pkm1/qkm1;
370                do
371                {
372                    k = k+1;
373                    if( k%2==1 )
374                    {
375                        yk = 1;
376                        xk = n+((double)(k-1))/(double)(2);
377                    }
378                    else
379                    {
380                        yk = x;
381                        xk = (double)(k)/(double)(2);
382                    }
383                    pk = pkm1*yk+pkm2*xk;
384                    qk = qkm1*yk+qkm2*xk;
385                    if( (double)(qk)!=(double)(0) )
386                    {
387                        r = pk/qk;
388                        t = Math.Abs((result-r)/r);
389                        result = r;
390                    }
391                    else
392                    {
393                        t = 1;
394                    }
395                    pkm2 = pkm1;
396                    pkm1 = pk;
397                    qkm2 = qkm1;
398                    qkm1 = qk;
399                    if( (double)(Math.Abs(pk))>(double)(big) )
400                    {
401                        pkm2 = pkm2/big;
402                        pkm1 = pkm1/big;
403                        qkm2 = qkm2/big;
404                        qkm1 = qkm1/big;
405                    }
406                }
407                while( (double)(t)>=(double)(AP.Math.MachineEpsilon) );
408                result = result*Math.Exp(-x);
409            }
410            return result;
411        }
412    }
413}
Note: See TracBrowser for help on using the repository browser.