1 /*----------------------------------------------------------------------------
2 gamma.cc (gamma function)
3 This file is a part of topaz systems
4 This module is based on the codes freely distributed by Mr. Okumura.
5 This module is free from copyright protection.
6
7
8 This program is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2 of the License, or
11 (at your option) any later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
22 ----------------------------------------------------------------------------*/
23
24 #include <stdio.h>
25 #include <math.h>
26 #include "naninf.h"
27
lngamma(double x)28 double lngamma(double x)
29 {
30 double fa, fb, fc;
31
32 fa = 1.0;
33 while (x < 8)
34 {
35 fa *= x;
36 x++;
37 }
38 fb = 1.0 / x / x;
39
40 fc = - 3617.0 / 510.0 / 16.0 / 15.0 * fb;
41 fc = (fc + 7.0 / 6.0 / 14.0 / 13.0) * fb;
42 fc = (fc - 691.0 / 2730.0 / 12.0 / 11.0) * fb;
43 fc = (fc + 5.0 / 66.0 / 10.0 / 9.0) * fb;
44 fc = (fc - 1.0 / 30.0 / 8.0 / 7.0) * fb;
45 fc = (fc + 1.0 / 42.0 / 6.0 / 5.0) * fb;
46 fc = (fc - 1.0 / 30.0 / 4.0 / 3.0) * fb;
47 fc = (fc + 1.0 / 6.0 / 2.0 / 1.0) / x;
48 fc += 0.5 * log(2.0 * M_PI) - log(fa) - x + (x - 0.5) * log(x);
49
50 return fc;
51
52 }
53
54
55 double q_gamma(double a, double x, double b);
56
p_gamma(double a,double x,double b)57 double p_gamma(double a, double x, double b)
58 {
59 double sum, value, osum;
60
61 // trivial
62 if (x == 0.0)
63 return 0.0;
64
65 if (x >= 1.0 + a)
66 return 1.0 - q_gamma(a, x, b);
67
68 value = exp(a * log(x) - x - b) / a;
69 sum = value;
70
71 for (int i = 1; i < 10000; i++)
72 {
73 osum = sum;
74 value *= x / (a + (double)i);
75 sum += value;
76 // converged !!
77 if (sum == osum)
78 return sum;
79 }
80 // not converged
81 return infvalue();
82 }
83
q_gamma(double a,double x,double b)84 double q_gamma(double a, double x, double b)
85 {
86 double sum, osum, fa, fb, fc, fd;
87
88 fc = 1.0;
89 fd = 1.0 + x - a;
90
91 if (x < 1 + a)
92 return 1 - p_gamma(a, x, b);
93
94 fa = exp(a * log(x) - x - b);
95 sum = fa / fd;
96
97 for (int i = 2; i < 10000; i++)
98 {
99 osum = sum;
100 fa *= ((double)i - 1.0 - a) / (double)i;
101 fb = (((double)i + x) * fd + ((double)i - 1.0 - a) * (fd - fc)) / (double)i;
102 fc = fd;
103 fd = fb;
104 fb = fa / (fc * fd);
105 sum += fb;
106
107 // converged !!
108 if (sum == osum)
109 return sum;
110 }
111 // not converged
112 return infvalue();
113 }
114
erf(double x)115 double erf(double x)
116 {
117 if (x >= 0.0)
118 return p_gamma(0.5, x * x, log(M_PI) / 2.0);
119 else
120 return - p_gamma(0.5, x * x, log(M_PI) / 2.0);
121 }
122
erfc(double x)123 double erfc(double x)
124 {
125 if (x >= 0.0)
126 return q_gamma(0.5, x * x, log(M_PI) / 2.0);
127 else
128 return 1.0 + p_gamma(0.5, x * x, log(M_PI) / 2.0);
129 }
130
p_normal(double x)131 double p_normal(double x)
132 {
133 if (x >= 0.0)
134 return 0.5 * (1.0 + p_gamma(0.5, 0.5 * x * x, log(M_PI) / 2.0));
135 else
136 return 0.5 * q_gamma(0.5, 0.5 * x * x, log(M_PI) / 2.0);
137 }
138
q_normal(double x)139 double q_normal(double x)
140 {
141 if (x >= 0.0)
142 return 0.5 * q_gamma(0.5, 0.5 * x * x, log(M_PI) / 2.0);
143 else
144 return 0.5 * (1.0 + p_gamma(0.5, 0.5 * x * x, log(M_PI) / 2.0));
145 }
146
gamma(double x)147 double gamma(double x)
148 {
149 if (x >= 0.0)
150 return exp(lngamma(x));
151 else
152 return M_PI / (sin(M_PI * x) * exp(lngamma(1.0 - x)));
153
154 }
155
beta(double x,double y)156 double beta(double x, double y)
157 {
158 return exp(lngamma(x) + lngamma(y) - lngamma(x + y));
159 }
160