Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/gammafunc.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: 10.8 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 gammafunc
33    {
34        /*************************************************************************
35        Gamma function
36
37        Input parameters:
38            X   -   argument
39
40        Domain:
41            0 < X < 171.6
42            -170 < X < 0, X is not an integer.
43
44        Relative error:
45         arithmetic   domain     # trials      peak         rms
46            IEEE    -170,-33      20000       2.3e-15     3.3e-16
47            IEEE     -33,  33     20000       9.4e-16     2.2e-16
48            IEEE      33, 171.6   20000       2.3e-15     3.2e-16
49
50        Cephes Math Library Release 2.8:  June, 2000
51        Original copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
52        Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
53        *************************************************************************/
54        public static double gamma(double x)
55        {
56            double result = 0;
57            double p = 0;
58            double pp = 0;
59            double q = 0;
60            double qq = 0;
61            double z = 0;
62            int i = 0;
63            double sgngam = 0;
64
65            sgngam = 1;
66            q = Math.Abs(x);
67            if( (double)(q)>(double)(33.0) )
68            {
69                if( (double)(x)<(double)(0.0) )
70                {
71                    p = (int)Math.Floor(q);
72                    i = (int)Math.Round(p);
73                    if( i%2==0 )
74                    {
75                        sgngam = -1;
76                    }
77                    z = q-p;
78                    if( (double)(z)>(double)(0.5) )
79                    {
80                        p = p+1;
81                        z = q-p;
82                    }
83                    z = q*Math.Sin(Math.PI*z);
84                    z = Math.Abs(z);
85                    z = Math.PI/(z*gammastirf(q));
86                }
87                else
88                {
89                    z = gammastirf(x);
90                }
91                result = sgngam*z;
92                return result;
93            }
94            z = 1;
95            while( (double)(x)>=(double)(3) )
96            {
97                x = x-1;
98                z = z*x;
99            }
100            while( (double)(x)<(double)(0) )
101            {
102                if( (double)(x)>(double)(-0.000000001) )
103                {
104                    result = z/((1+0.5772156649015329*x)*x);
105                    return result;
106                }
107                z = z/x;
108                x = x+1;
109            }
110            while( (double)(x)<(double)(2) )
111            {
112                if( (double)(x)<(double)(0.000000001) )
113                {
114                    result = z/((1+0.5772156649015329*x)*x);
115                    return result;
116                }
117                z = z/x;
118                x = x+1.0;
119            }
120            if( (double)(x)==(double)(2) )
121            {
122                result = z;
123                return result;
124            }
125            x = x-2.0;
126            pp = 1.60119522476751861407E-4;
127            pp = 1.19135147006586384913E-3+x*pp;
128            pp = 1.04213797561761569935E-2+x*pp;
129            pp = 4.76367800457137231464E-2+x*pp;
130            pp = 2.07448227648435975150E-1+x*pp;
131            pp = 4.94214826801497100753E-1+x*pp;
132            pp = 9.99999999999999996796E-1+x*pp;
133            qq = -2.31581873324120129819E-5;
134            qq = 5.39605580493303397842E-4+x*qq;
135            qq = -4.45641913851797240494E-3+x*qq;
136            qq = 1.18139785222060435552E-2+x*qq;
137            qq = 3.58236398605498653373E-2+x*qq;
138            qq = -2.34591795718243348568E-1+x*qq;
139            qq = 7.14304917030273074085E-2+x*qq;
140            qq = 1.00000000000000000320+x*qq;
141            result = z*pp/qq;
142            return result;
143            return result;
144        }
145
146
147        /*************************************************************************
148        Natural logarithm of gamma function
149
150        Input parameters:
151            X       -   argument
152
153        Result:
154            logarithm of the absolute value of the Gamma(X).
155
156        Output parameters:
157            SgnGam  -   sign(Gamma(X))
158
159        Domain:
160            0 < X < 2.55e305
161            -2.55e305 < X < 0, X is not an integer.
162
163        ACCURACY:
164        arithmetic      domain        # trials     peak         rms
165           IEEE    0, 3                 28000     5.4e-16     1.1e-16
166           IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
167        The error criterion was relative when the function magnitude
168        was greater than one but absolute when it was less than one.
169
170        The following test used the relative error criterion, though
171        at certain points the relative error could be much higher than
172        indicated.
173           IEEE    -200, -4             10000     4.8e-16     1.3e-16
174
175        Cephes Math Library Release 2.8:  June, 2000
176        Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
177        Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
178        *************************************************************************/
179        public static double lngamma(double x,
180            ref double sgngam)
181        {
182            double result = 0;
183            double a = 0;
184            double b = 0;
185            double c = 0;
186            double p = 0;
187            double q = 0;
188            double u = 0;
189            double w = 0;
190            double z = 0;
191            int i = 0;
192            double logpi = 0;
193            double ls2pi = 0;
194            double tmp = 0;
195
196            sgngam = 1;
197            logpi = 1.14472988584940017414;
198            ls2pi = 0.91893853320467274178;
199            if( (double)(x)<(double)(-34.0) )
200            {
201                q = -x;
202                w = lngamma(q, ref tmp);
203                p = (int)Math.Floor(q);
204                i = (int)Math.Round(p);
205                if( i%2==0 )
206                {
207                    sgngam = -1;
208                }
209                else
210                {
211                    sgngam = 1;
212                }
213                z = q-p;
214                if( (double)(z)>(double)(0.5) )
215                {
216                    p = p+1;
217                    z = p-q;
218                }
219                z = q*Math.Sin(Math.PI*z);
220                result = logpi-Math.Log(z)-w;
221                return result;
222            }
223            if( (double)(x)<(double)(13) )
224            {
225                z = 1;
226                p = 0;
227                u = x;
228                while( (double)(u)>=(double)(3) )
229                {
230                    p = p-1;
231                    u = x+p;
232                    z = z*u;
233                }
234                while( (double)(u)<(double)(2) )
235                {
236                    z = z/u;
237                    p = p+1;
238                    u = x+p;
239                }
240                if( (double)(z)<(double)(0) )
241                {
242                    sgngam = -1;
243                    z = -z;
244                }
245                else
246                {
247                    sgngam = 1;
248                }
249                if( (double)(u)==(double)(2) )
250                {
251                    result = Math.Log(z);
252                    return result;
253                }
254                p = p-2;
255                x = x+p;
256                b = -1378.25152569120859100;
257                b = -38801.6315134637840924+x*b;
258                b = -331612.992738871184744+x*b;
259                b = -1162370.97492762307383+x*b;
260                b = -1721737.00820839662146+x*b;
261                b = -853555.664245765465627+x*b;
262                c = 1;
263                c = -351.815701436523470549+x*c;
264                c = -17064.2106651881159223+x*c;
265                c = -220528.590553854454839+x*c;
266                c = -1139334.44367982507207+x*c;
267                c = -2532523.07177582951285+x*c;
268                c = -2018891.41433532773231+x*c;
269                p = x*b/c;
270                result = Math.Log(z)+p;
271                return result;
272            }
273            q = (x-0.5)*Math.Log(x)-x+ls2pi;
274            if( (double)(x)>(double)(100000000) )
275            {
276                result = q;
277                return result;
278            }
279            p = 1/(x*x);
280            if( (double)(x)>=(double)(1000.0) )
281            {
282                q = q+((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
283            }
284            else
285            {
286                a = 8.11614167470508450300*0.0001;
287                a = -(5.95061904284301438324*0.0001)+p*a;
288                a = 7.93650340457716943945*0.0001+p*a;
289                a = -(2.77777777730099687205*0.001)+p*a;
290                a = 8.33333333333331927722*0.01+p*a;
291                q = q+a/x;
292            }
293            result = q;
294            return result;
295        }
296
297
298        private static double gammastirf(double x)
299        {
300            double result = 0;
301            double y = 0;
302            double w = 0;
303            double v = 0;
304            double stir = 0;
305
306            w = 1/x;
307            stir = 7.87311395793093628397E-4;
308            stir = -2.29549961613378126380E-4+w*stir;
309            stir = -2.68132617805781232825E-3+w*stir;
310            stir = 3.47222221605458667310E-3+w*stir;
311            stir = 8.33333333333482257126E-2+w*stir;
312            w = 1+w*stir;
313            y = Math.Exp(x);
314            if( (double)(x)>(double)(143.01608) )
315            {
316                v = Math.Pow(x, 0.5*x-0.25);
317                y = v*(v/y);
318            }
319            else
320            {
321                y = Math.Pow(x, x-0.5)/y;
322            }
323            result = 2.50662827463100050242*y*w;
324            return result;
325        }
326    }
327}
Note: See TracBrowser for help on using the repository browser.