Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/igammaf.cs @ 9280

Last change on this file since 9280 was 2806, checked in by gkronber, 15 years ago

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)
File size: 13.3 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 igammaf
33    {
34        /*************************************************************************
35        Incomplete gamma integral
36
37        The function is defined by
38
39                                  x
40                                   -
41                          1       | |  -t  a-1
42         igam(a,x)  =   -----     |   e   t   dt.
43                         -      | |
44                        | (a)    -
45                                  0
46
47
48        In this implementation both arguments must be positive.
49        The integral is evaluated by either a power series or
50        continued fraction expansion, depending on the relative
51        values of a and x.
52
53        ACCURACY:
54
55                             Relative error:
56        arithmetic   domain     # trials      peak         rms
57           IEEE      0,30       200000       3.6e-14     2.9e-15
58           IEEE      0,100      300000       9.9e-14     1.5e-14
59
60        Cephes Math Library Release 2.8:  June, 2000
61        Copyright 1985, 1987, 2000 by Stephen L. Moshier
62        *************************************************************************/
63        public static double incompletegamma(double a,
64            double x)
65        {
66            double result = 0;
67            double igammaepsilon = 0;
68            double ans = 0;
69            double ax = 0;
70            double c = 0;
71            double r = 0;
72            double tmp = 0;
73
74            igammaepsilon = 0.000000000000001;
75            if( (double)(x)<=(double)(0) | (double)(a)<=(double)(0) )
76            {
77                result = 0;
78                return result;
79            }
80            if( (double)(x)>(double)(1) & (double)(x)>(double)(a) )
81            {
82                result = 1-incompletegammac(a, x);
83                return result;
84            }
85            ax = a*Math.Log(x)-x-gammafunc.lngamma(a, ref tmp);
86            if( (double)(ax)<(double)(-709.78271289338399) )
87            {
88                result = 0;
89                return result;
90            }
91            ax = Math.Exp(ax);
92            r = a;
93            c = 1;
94            ans = 1;
95            do
96            {
97                r = r+1;
98                c = c*x/r;
99                ans = ans+c;
100            }
101            while( (double)(c/ans)>(double)(igammaepsilon) );
102            result = ans*ax/a;
103            return result;
104        }
105
106
107        /*************************************************************************
108        Complemented incomplete gamma integral
109
110        The function is defined by
111
112
113         igamc(a,x)   =   1 - igam(a,x)
114
115                                   inf.
116                                     -
117                            1       | |  -t  a-1
118                      =   -----     |   e   t   dt.
119                           -      | |
120                          | (a)    -
121                                    x
122
123
124        In this implementation both arguments must be positive.
125        The integral is evaluated by either a power series or
126        continued fraction expansion, depending on the relative
127        values of a and x.
128
129        ACCURACY:
130
131        Tested at random a, x.
132                       a         x                      Relative error:
133        arithmetic   domain   domain     # trials      peak         rms
134           IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
135           IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15
136
137        Cephes Math Library Release 2.8:  June, 2000
138        Copyright 1985, 1987, 2000 by Stephen L. Moshier
139        *************************************************************************/
140        public static double incompletegammac(double a,
141            double x)
142        {
143            double result = 0;
144            double igammaepsilon = 0;
145            double igammabignumber = 0;
146            double igammabignumberinv = 0;
147            double ans = 0;
148            double ax = 0;
149            double c = 0;
150            double yc = 0;
151            double r = 0;
152            double t = 0;
153            double y = 0;
154            double z = 0;
155            double pk = 0;
156            double pkm1 = 0;
157            double pkm2 = 0;
158            double qk = 0;
159            double qkm1 = 0;
160            double qkm2 = 0;
161            double tmp = 0;
162
163            igammaepsilon = 0.000000000000001;
164            igammabignumber = 4503599627370496.0;
165            igammabignumberinv = 2.22044604925031308085*0.0000000000000001;
166            if( (double)(x)<=(double)(0) | (double)(a)<=(double)(0) )
167            {
168                result = 1;
169                return result;
170            }
171            if( (double)(x)<(double)(1) | (double)(x)<(double)(a) )
172            {
173                result = 1-incompletegamma(a, x);
174                return result;
175            }
176            ax = a*Math.Log(x)-x-gammafunc.lngamma(a, ref tmp);
177            if( (double)(ax)<(double)(-709.78271289338399) )
178            {
179                result = 0;
180                return result;
181            }
182            ax = Math.Exp(ax);
183            y = 1-a;
184            z = x+y+1;
185            c = 0;
186            pkm2 = 1;
187            qkm2 = x;
188            pkm1 = x+1;
189            qkm1 = z*x;
190            ans = pkm1/qkm1;
191            do
192            {
193                c = c+1;
194                y = y+1;
195                z = z+2;
196                yc = y*c;
197                pk = pkm1*z-pkm2*yc;
198                qk = qkm1*z-qkm2*yc;
199                if( (double)(qk)!=(double)(0) )
200                {
201                    r = pk/qk;
202                    t = Math.Abs((ans-r)/r);
203                    ans = r;
204                }
205                else
206                {
207                    t = 1;
208                }
209                pkm2 = pkm1;
210                pkm1 = pk;
211                qkm2 = qkm1;
212                qkm1 = qk;
213                if( (double)(Math.Abs(pk))>(double)(igammabignumber) )
214                {
215                    pkm2 = pkm2*igammabignumberinv;
216                    pkm1 = pkm1*igammabignumberinv;
217                    qkm2 = qkm2*igammabignumberinv;
218                    qkm1 = qkm1*igammabignumberinv;
219                }
220            }
221            while( (double)(t)>(double)(igammaepsilon) );
222            result = ans*ax;
223            return result;
224        }
225
226
227        /*************************************************************************
228        Inverse of complemented imcomplete gamma integral
229
230        Given p, the function finds x such that
231
232         igamc( a, x ) = p.
233
234        Starting with the approximate value
235
236                3
237         x = a t
238
239         where
240
241         t = 1 - d - ndtri(p) sqrt(d)
242
243        and
244
245         d = 1/9a,
246
247        the routine performs up to 10 Newton iterations to find the
248        root of igamc(a,x) - p = 0.
249
250        ACCURACY:
251
252        Tested at random a, p in the intervals indicated.
253
254                       a        p                      Relative error:
255        arithmetic   domain   domain     # trials      peak         rms
256           IEEE     0.5,100   0,0.5       100000       1.0e-14     1.7e-15
257           IEEE     0.01,0.5  0,0.5       100000       9.0e-14     3.4e-15
258           IEEE    0.5,10000  0,0.5        20000       2.3e-13     3.8e-14
259
260        Cephes Math Library Release 2.8:  June, 2000
261        Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
262        *************************************************************************/
263        public static double invincompletegammac(double a,
264            double y0)
265        {
266            double result = 0;
267            double igammaepsilon = 0;
268            double iinvgammabignumber = 0;
269            double x0 = 0;
270            double x1 = 0;
271            double x = 0;
272            double yl = 0;
273            double yh = 0;
274            double y = 0;
275            double d = 0;
276            double lgm = 0;
277            double dithresh = 0;
278            int i = 0;
279            int dir = 0;
280            double tmp = 0;
281
282            igammaepsilon = 0.000000000000001;
283            iinvgammabignumber = 4503599627370496.0;
284            x0 = iinvgammabignumber;
285            yl = 0;
286            x1 = 0;
287            yh = 1;
288            dithresh = 5*igammaepsilon;
289            d = 1/(9*a);
290            y = 1-d-normaldistr.invnormaldistribution(y0)*Math.Sqrt(d);
291            x = a*y*y*y;
292            lgm = gammafunc.lngamma(a, ref tmp);
293            i = 0;
294            while( i<10 )
295            {
296                if( (double)(x)>(double)(x0) | (double)(x)<(double)(x1) )
297                {
298                    d = 0.0625;
299                    break;
300                }
301                y = incompletegammac(a, x);
302                if( (double)(y)<(double)(yl) | (double)(y)>(double)(yh) )
303                {
304                    d = 0.0625;
305                    break;
306                }
307                if( (double)(y)<(double)(y0) )
308                {
309                    x0 = x;
310                    yl = y;
311                }
312                else
313                {
314                    x1 = x;
315                    yh = y;
316                }
317                d = (a-1)*Math.Log(x)-x-lgm;
318                if( (double)(d)<(double)(-709.78271289338399) )
319                {
320                    d = 0.0625;
321                    break;
322                }
323                d = -Math.Exp(d);
324                d = (y-y0)/d;
325                if( (double)(Math.Abs(d/x))<(double)(igammaepsilon) )
326                {
327                    result = x;
328                    return result;
329                }
330                x = x-d;
331                i = i+1;
332            }
333            if( (double)(x0)==(double)(iinvgammabignumber) )
334            {
335                if( (double)(x)<=(double)(0) )
336                {
337                    x = 1;
338                }
339                while( (double)(x0)==(double)(iinvgammabignumber) )
340                {
341                    x = (1+d)*x;
342                    y = incompletegammac(a, x);
343                    if( (double)(y)<(double)(y0) )
344                    {
345                        x0 = x;
346                        yl = y;
347                        break;
348                    }
349                    d = d+d;
350                }
351            }
352            d = 0.5;
353            dir = 0;
354            i = 0;
355            while( i<400 )
356            {
357                x = x1+d*(x0-x1);
358                y = incompletegammac(a, x);
359                lgm = (x0-x1)/(x1+x0);
360                if( (double)(Math.Abs(lgm))<(double)(dithresh) )
361                {
362                    break;
363                }
364                lgm = (y-y0)/y0;
365                if( (double)(Math.Abs(lgm))<(double)(dithresh) )
366                {
367                    break;
368                }
369                if( (double)(x)<=(double)(0.0) )
370                {
371                    break;
372                }
373                if( (double)(y)>=(double)(y0) )
374                {
375                    x1 = x;
376                    yh = y;
377                    if( dir<0 )
378                    {
379                        dir = 0;
380                        d = 0.5;
381                    }
382                    else
383                    {
384                        if( dir>1 )
385                        {
386                            d = 0.5*d+0.5;
387                        }
388                        else
389                        {
390                            d = (y0-yl)/(yh-yl);
391                        }
392                    }
393                    dir = dir+1;
394                }
395                else
396                {
397                    x0 = x;
398                    yl = y;
399                    if( dir>0 )
400                    {
401                        dir = 0;
402                        d = 0.5;
403                    }
404                    else
405                    {
406                        if( dir<-1 )
407                        {
408                            d = 0.5*d;
409                        }
410                        else
411                        {
412                            d = (y0-yl)/(yh-yl);
413                        }
414                    }
415                    dir = dir-1;
416                }
417                i = i+1;
418            }
419            result = x;
420            return result;
421        }
422    }
423}
Note: See TracBrowser for help on using the repository browser.