OLD | NEW |
1 /* | 1 /* |
2 "A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964) | 2 "A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964) |
3 "Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001) | 3 "Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001) |
4 "An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004) | 4 "An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004) |
5 | 5 |
6 approximation method: | 6 approximation method: |
7 | 7 |
8 (x - 0.5) S(x) | 8 (x - 0.5) S(x) |
9 Gamma(x) = (x + g - 0.5) * ---------------- | 9 Gamma(x) = (x + g - 0.5) * ---------------- |
10 exp(x + g - 0.5) | 10 exp(x + g - 0.5) |
11 | 11 |
12 with | 12 with |
13 a1 a2 a3 aN | 13 a1 a2 a3 aN |
14 S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ] | 14 S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ] |
15 x + 1 x + 2 x + 3 x + N | 15 x + 1 x + 2 x + 3 x + N |
16 | 16 |
17 with a0, a1, a2, a3,.. aN constants which depend on g. | 17 with a0, a1, a2, a3,.. aN constants which depend on g. |
18 | 18 |
19 for x < 0 the following reflection formula is used: | 19 for x < 0 the following reflection formula is used: |
20 | 20 |
21 Gamma(x)*Gamma(-x) = -pi/(x sin(pi x)) | 21 Gamma(x)*Gamma(-x) = -pi/(x sin(pi x)) |
22 | 22 |
23 most ideas and constants are from boost and python | 23 most ideas and constants are from boost and python |
24 */ | 24 */ |
25 #include "libm.h" | 25 #include "libm.h" |
26 | 26 |
27 static const double pi = 3.141592653589793238462643383279502884; | 27 static const double pi = 3.141592653589793238462643383279502884; |
28 | 28 |
29 /* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */ | 29 /* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */ |
30 static double sinpi(double x) | 30 static double sinpi(double x) { |
31 { | 31 int n; |
32 » int n; | |
33 | 32 |
34 » /* argument reduction: x = |x| mod 2 */ | 33 /* argument reduction: x = |x| mod 2 */ |
35 » /* spurious inexact when x is odd int */ | 34 /* spurious inexact when x is odd int */ |
36 » x = x * 0.5; | 35 x = x * 0.5; |
37 » x = 2 * (x - floor(x)); | 36 x = 2 * (x - floor(x)); |
38 | 37 |
39 » /* reduce x into [-.25,.25] */ | 38 /* reduce x into [-.25,.25] */ |
40 » n = 4 * x; | 39 n = 4 * x; |
41 » n = (n+1)/2; | 40 n = (n + 1) / 2; |
42 » x -= n * 0.5; | 41 x -= n * 0.5; |
43 | 42 |
44 » x *= pi; | 43 x *= pi; |
45 » switch (n) { | 44 switch (n) { |
46 » default: /* case 4 */ | 45 default: /* case 4 */ |
47 » case 0: | 46 case 0: |
48 » » return __sin(x, 0, 0); | 47 return __sin(x, 0, 0); |
49 » case 1: | 48 case 1: |
50 » » return __cos(x, 0); | 49 return __cos(x, 0); |
51 » case 2: | 50 case 2: |
52 » » return __sin(-x, 0, 0); | 51 return __sin(-x, 0, 0); |
53 » case 3: | 52 case 3: |
54 » » return -__cos(x, 0); | 53 return -__cos(x, 0); |
55 » } | 54 } |
56 } | 55 } |
57 | 56 |
58 #define N 12 | 57 #define N 12 |
59 //static const double g = 6.024680040776729583740234375; | 58 // static const double g = 6.024680040776729583740234375; |
60 static const double gmhalf = 5.524680040776729583740234375; | 59 static const double gmhalf = 5.524680040776729583740234375; |
61 static const double Snum[N+1] = { | 60 static const double Snum[N + 1] = { |
62 » 23531376880.410759688572007674451636754734846804940, | 61 23531376880.410759688572007674451636754734846804940, |
63 » 42919803642.649098768957899047001988850926355848959, | 62 42919803642.649098768957899047001988850926355848959, |
64 » 35711959237.355668049440185451547166705960488635843, | 63 35711959237.355668049440185451547166705960488635843, |
65 » 17921034426.037209699919755754458931112671403265390, | 64 17921034426.037209699919755754458931112671403265390, |
66 » 6039542586.3520280050642916443072979210699388420708, | 65 6039542586.3520280050642916443072979210699388420708, |
67 » 1439720407.3117216736632230727949123939715485786772, | 66 1439720407.3117216736632230727949123939715485786772, |
68 » 248874557.86205415651146038641322942321632125127801, | 67 248874557.86205415651146038641322942321632125127801, |
69 » 31426415.585400194380614231628318205362874684987640, | 68 31426415.585400194380614231628318205362874684987640, |
70 » 2876370.6289353724412254090516208496135991145378768, | 69 2876370.6289353724412254090516208496135991145378768, |
71 » 186056.26539522349504029498971604569928220784236328, | 70 186056.26539522349504029498971604569928220784236328, |
72 » 8071.6720023658162106380029022722506138218516325024, | 71 8071.6720023658162106380029022722506138218516325024, |
73 » 210.82427775157934587250973392071336271166969580291, | 72 210.82427775157934587250973392071336271166969580291, |
74 » 2.5066282746310002701649081771338373386264310793408, | 73 2.5066282746310002701649081771338373386264310793408, |
75 }; | 74 }; |
76 static const double Sden[N+1] = { | 75 static const double Sden[N + 1] = { |
77 » 0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535, | 76 0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535, |
78 » 2637558, 357423, 32670, 1925, 66, 1, | 77 2637558, 357423, 32670, 1925, 66, 1, |
79 }; | 78 }; |
80 /* n! for small integer n */ | 79 /* n! for small integer n */ |
81 static const double fact[] = { | 80 static const double fact[] = { |
82 » 1, 1, 2, 6, 24, 120, 720, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800
.0, | 81 1, |
83 » 479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0, 2092278988800
0.0, | 82 1, |
84 » 355687428096000.0, 6402373705728000.0, 121645100408832000.0, | 83 2, |
85 » 2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0, | 84 6, |
| 85 24, |
| 86 120, |
| 87 720, |
| 88 5040.0, |
| 89 40320.0, |
| 90 362880.0, |
| 91 3628800.0, |
| 92 39916800.0, |
| 93 479001600.0, |
| 94 6227020800.0, |
| 95 87178291200.0, |
| 96 1307674368000.0, |
| 97 20922789888000.0, |
| 98 355687428096000.0, |
| 99 6402373705728000.0, |
| 100 121645100408832000.0, |
| 101 2432902008176640000.0, |
| 102 51090942171709440000.0, |
| 103 1124000727777607680000.0, |
86 }; | 104 }; |
87 | 105 |
88 /* S(x) rational function for positive x */ | 106 /* S(x) rational function for positive x */ |
89 static double S(double x) | 107 static double S(double x) { |
90 { | 108 double_t num = 0, den = 0; |
91 » double_t num = 0, den = 0; | 109 int i; |
92 » int i; | |
93 | 110 |
94 » /* to avoid overflow handle large x differently */ | 111 /* to avoid overflow handle large x differently */ |
95 » if (x < 8) | 112 if (x < 8) |
96 » » for (i = N; i >= 0; i--) { | 113 for (i = N; i >= 0; i--) { |
97 » » » num = num * x + Snum[i]; | 114 num = num * x + Snum[i]; |
98 » » » den = den * x + Sden[i]; | 115 den = den * x + Sden[i]; |
99 » » } | 116 } |
100 » else | 117 else |
101 » » for (i = 0; i <= N; i++) { | 118 for (i = 0; i <= N; i++) { |
102 » » » num = num / x + Snum[i]; | 119 num = num / x + Snum[i]; |
103 » » » den = den / x + Sden[i]; | 120 den = den / x + Sden[i]; |
104 » » } | 121 } |
105 » return num/den; | 122 return num / den; |
106 } | 123 } |
107 | 124 |
108 double tgamma(double x) | 125 double tgamma(double x) { |
109 { | 126 union { |
110 » union {double f; uint64_t i;} u = {x}; | 127 double f; |
111 » double absx, y; | 128 uint64_t i; |
112 » double_t dy, z, r; | 129 } u = {x}; |
113 » uint32_t ix = u.i>>32 & 0x7fffffff; | 130 double absx, y; |
114 » int sign = u.i>>63; | 131 double_t dy, z, r; |
| 132 uint32_t ix = u.i >> 32 & 0x7fffffff; |
| 133 int sign = u.i >> 63; |
115 | 134 |
116 » /* special cases */ | 135 /* special cases */ |
117 » if (ix >= 0x7ff00000) | 136 if (ix >= 0x7ff00000) |
118 » » /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with inval
id */ | 137 /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */ |
119 » » return x + INFINITY; | 138 return x + INFINITY; |
120 » if (ix < (0x3ff-54)<<20) | 139 if (ix < (0x3ff - 54) << 20) |
121 » » /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */ | 140 /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */ |
122 » » return 1/x; | 141 return 1 / x; |
123 | 142 |
124 » /* integer arguments */ | 143 /* integer arguments */ |
125 » /* raise inexact when non-integer */ | 144 /* raise inexact when non-integer */ |
126 » if (x == floor(x)) { | 145 if (x == floor(x)) { |
127 » » if (sign) | 146 if (sign) |
128 » » » return 0/0.0; | 147 return 0 / 0.0; |
129 » » if (x <= sizeof fact/sizeof *fact) | 148 if (x <= sizeof fact / sizeof *fact) |
130 » » » return fact[(int)x - 1]; | 149 return fact[(int)x - 1]; |
131 » } | 150 } |
132 | 151 |
133 » /* x >= 172: tgamma(x)=inf with overflow */ | 152 /* x >= 172: tgamma(x)=inf with overflow */ |
134 » /* x =< -184: tgamma(x)=+-0 with underflow */ | 153 /* x =< -184: tgamma(x)=+-0 with underflow */ |
135 » if (ix >= 0x40670000) { /* |x| >= 184 */ | 154 if (ix >= 0x40670000) { /* |x| >= 184 */ |
136 » » if (sign) { | 155 if (sign) { |
137 » » » FORCE_EVAL((float)(0x1p-126/x)); | 156 FORCE_EVAL((float)(0x1p-126 / x)); |
138 » » » if (floor(x) * 0.5 == floor(x * 0.5)) | 157 if (floor(x) * 0.5 == floor(x * 0.5)) |
139 » » » » return 0; | 158 return 0; |
140 » » » return -0.0; | 159 return -0.0; |
141 » » } | 160 } |
142 » » x *= 0x1p1023; | 161 x *= 0x1p1023; |
143 » » return x; | 162 return x; |
144 » } | 163 } |
145 | 164 |
146 » absx = sign ? -x : x; | 165 absx = sign ? -x : x; |
147 | 166 |
148 » /* handle the error of x + g - 0.5 */ | 167 /* handle the error of x + g - 0.5 */ |
149 » y = absx + gmhalf; | 168 y = absx + gmhalf; |
150 » if (absx > gmhalf) { | 169 if (absx > gmhalf) { |
151 » » dy = y - absx; | 170 dy = y - absx; |
152 » » dy -= gmhalf; | 171 dy -= gmhalf; |
153 » } else { | 172 } else { |
154 » » dy = y - gmhalf; | 173 dy = y - gmhalf; |
155 » » dy -= absx; | 174 dy -= absx; |
156 » } | 175 } |
157 | 176 |
158 » z = absx - 0.5; | 177 z = absx - 0.5; |
159 » r = S(absx) * exp(-y); | 178 r = S(absx) * exp(-y); |
160 » if (x < 0) { | 179 if (x < 0) { |
161 » » /* reflection formula for negative x */ | 180 /* reflection formula for negative x */ |
162 » » /* sinpi(absx) is not 0, integers are already handled */ | 181 /* sinpi(absx) is not 0, integers are already handled */ |
163 » » r = -pi / (sinpi(absx) * absx * r); | 182 r = -pi / (sinpi(absx) * absx * r); |
164 » » dy = -dy; | 183 dy = -dy; |
165 » » z = -z; | 184 z = -z; |
166 » } | 185 } |
167 » r += dy * (gmhalf+0.5) * r / y; | 186 r += dy * (gmhalf + 0.5) * r / y; |
168 » z = pow(y, 0.5*z); | 187 z = pow(y, 0.5 * z); |
169 » y = r * z * z; | 188 y = r * z * z; |
170 » return y; | 189 return y; |
171 } | 190 } |
172 | 191 |
173 #if 0 | 192 #if 0 |
174 double __lgamma_r(double x, int *sign) | 193 double __lgamma_r(double x, int *sign) |
175 { | 194 { |
176 double r, absx; | 195 double r, absx; |
177 | 196 |
178 *sign = 1; | 197 *sign = 1; |
179 | 198 |
180 /* special cases */ | 199 /* special cases */ |
(...skipping 32 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
213 /* reflection formula for negative x */ | 232 /* reflection formula for negative x */ |
214 x = sinpi(absx); | 233 x = sinpi(absx); |
215 *sign = 2*!!signbit(x) - 1; | 234 *sign = 2*!!signbit(x) - 1; |
216 r = log(pi/(fabs(x)*absx)) - r; | 235 r = log(pi/(fabs(x)*absx)) - r; |
217 } | 236 } |
218 return r; | 237 return r; |
219 } | 238 } |
220 | 239 |
221 weak_alias(__lgamma_r, lgamma_r); | 240 weak_alias(__lgamma_r, lgamma_r); |
222 #endif | 241 #endif |
OLD | NEW |