1 |
|
|
/* $OpenBSD: e_tgammal.c,v 1.4 2013/11/12 20:35:19 martynas Exp $ */ |
2 |
|
|
|
3 |
|
|
/* |
4 |
|
|
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> |
5 |
|
|
* |
6 |
|
|
* Permission to use, copy, modify, and distribute this software for any |
7 |
|
|
* purpose with or without fee is hereby granted, provided that the above |
8 |
|
|
* copyright notice and this permission notice appear in all copies. |
9 |
|
|
* |
10 |
|
|
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
11 |
|
|
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
12 |
|
|
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR |
13 |
|
|
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES |
14 |
|
|
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN |
15 |
|
|
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF |
16 |
|
|
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. |
17 |
|
|
*/ |
18 |
|
|
|
19 |
|
|
/* tgammal.c |
20 |
|
|
* |
21 |
|
|
* Gamma function |
22 |
|
|
* |
23 |
|
|
* |
24 |
|
|
* |
25 |
|
|
* SYNOPSIS: |
26 |
|
|
* |
27 |
|
|
* long double x, y, tgammal(); |
28 |
|
|
* extern int signgam; |
29 |
|
|
* |
30 |
|
|
* y = tgammal( x ); |
31 |
|
|
* |
32 |
|
|
* |
33 |
|
|
* |
34 |
|
|
* DESCRIPTION: |
35 |
|
|
* |
36 |
|
|
* Returns gamma function of the argument. The result is |
37 |
|
|
* correctly signed, and the sign (+1 or -1) is also |
38 |
|
|
* returned in a global (extern) variable named signgam. |
39 |
|
|
* This variable is also filled in by the logarithmic gamma |
40 |
|
|
* function lgamma(). |
41 |
|
|
* |
42 |
|
|
* Arguments |x| <= 13 are reduced by recurrence and the function |
43 |
|
|
* approximated by a rational function of degree 7/8 in the |
44 |
|
|
* interval (2,3). Large arguments are handled by Stirling's |
45 |
|
|
* formula. Large negative arguments are made positive using |
46 |
|
|
* a reflection formula. |
47 |
|
|
* |
48 |
|
|
* |
49 |
|
|
* ACCURACY: |
50 |
|
|
* |
51 |
|
|
* Relative error: |
52 |
|
|
* arithmetic domain # trials peak rms |
53 |
|
|
* IEEE -40,+40 10000 3.6e-19 7.9e-20 |
54 |
|
|
* IEEE -1755,+1755 10000 4.8e-18 6.5e-19 |
55 |
|
|
* |
56 |
|
|
* Accuracy for large arguments is dominated by error in powl(). |
57 |
|
|
* |
58 |
|
|
*/ |
59 |
|
|
|
60 |
|
|
#include <float.h> |
61 |
|
|
#include <math.h> |
62 |
|
|
|
63 |
|
|
#include "math_private.h" |
64 |
|
|
|
65 |
|
|
/* |
66 |
|
|
tgamma(x+2) = tgamma(x+2) P(x)/Q(x) |
67 |
|
|
0 <= x <= 1 |
68 |
|
|
Relative error |
69 |
|
|
n=7, d=8 |
70 |
|
|
Peak error = 1.83e-20 |
71 |
|
|
Relative error spread = 8.4e-23 |
72 |
|
|
*/ |
73 |
|
|
|
74 |
|
|
static long double P[8] = { |
75 |
|
|
4.212760487471622013093E-5L, |
76 |
|
|
4.542931960608009155600E-4L, |
77 |
|
|
4.092666828394035500949E-3L, |
78 |
|
|
2.385363243461108252554E-2L, |
79 |
|
|
1.113062816019361559013E-1L, |
80 |
|
|
3.629515436640239168939E-1L, |
81 |
|
|
8.378004301573126728826E-1L, |
82 |
|
|
1.000000000000000000009E0L, |
83 |
|
|
}; |
84 |
|
|
static long double Q[9] = { |
85 |
|
|
-1.397148517476170440917E-5L, |
86 |
|
|
2.346584059160635244282E-4L, |
87 |
|
|
-1.237799246653152231188E-3L, |
88 |
|
|
-7.955933682494738320586E-4L, |
89 |
|
|
2.773706565840072979165E-2L, |
90 |
|
|
-4.633887671244534213831E-2L, |
91 |
|
|
-2.243510905670329164562E-1L, |
92 |
|
|
4.150160950588455434583E-1L, |
93 |
|
|
9.999999999999999999908E-1L, |
94 |
|
|
}; |
95 |
|
|
|
96 |
|
|
/* |
97 |
|
|
static long double P[] = { |
98 |
|
|
-3.01525602666895735709e0L, |
99 |
|
|
-3.25157411956062339893e1L, |
100 |
|
|
-2.92929976820724030353e2L, |
101 |
|
|
-1.70730828800510297666e3L, |
102 |
|
|
-7.96667499622741999770e3L, |
103 |
|
|
-2.59780216007146401957e4L, |
104 |
|
|
-5.99650230220855581642e4L, |
105 |
|
|
-7.15743521530849602425e4L |
106 |
|
|
}; |
107 |
|
|
static long double Q[] = { |
108 |
|
|
1.00000000000000000000e0L, |
109 |
|
|
-1.67955233807178858919e1L, |
110 |
|
|
8.85946791747759881659e1L, |
111 |
|
|
5.69440799097468430177e1L, |
112 |
|
|
-1.98526250512761318471e3L, |
113 |
|
|
3.31667508019495079814e3L, |
114 |
|
|
1.60577839621734713377e4L, |
115 |
|
|
-2.97045081369399940529e4L, |
116 |
|
|
-7.15743521530849602412e4L |
117 |
|
|
}; |
118 |
|
|
*/ |
119 |
|
|
#define MAXGAML 1755.455L |
120 |
|
|
/*static const long double LOGPI = 1.14472988584940017414L;*/ |
121 |
|
|
|
122 |
|
|
/* Stirling's formula for the gamma function |
123 |
|
|
tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) |
124 |
|
|
z(x) = x |
125 |
|
|
13 <= x <= 1024 |
126 |
|
|
Relative error |
127 |
|
|
n=8, d=0 |
128 |
|
|
Peak error = 9.44e-21 |
129 |
|
|
Relative error spread = 8.8e-4 |
130 |
|
|
*/ |
131 |
|
|
|
132 |
|
|
static long double STIR[9] = { |
133 |
|
|
7.147391378143610789273E-4L, |
134 |
|
|
-2.363848809501759061727E-5L, |
135 |
|
|
-5.950237554056330156018E-4L, |
136 |
|
|
6.989332260623193171870E-5L, |
137 |
|
|
7.840334842744753003862E-4L, |
138 |
|
|
-2.294719747873185405699E-4L, |
139 |
|
|
-2.681327161876304418288E-3L, |
140 |
|
|
3.472222222230075327854E-3L, |
141 |
|
|
8.333333333333331800504E-2L, |
142 |
|
|
}; |
143 |
|
|
|
144 |
|
|
#define MAXSTIR 1024.0L |
145 |
|
|
static const long double SQTPI = 2.50662827463100050242E0L; |
146 |
|
|
|
147 |
|
|
/* 1/tgamma(x) = z P(z) |
148 |
|
|
* z(x) = 1/x |
149 |
|
|
* 0 < x < 0.03125 |
150 |
|
|
* Peak relative error 4.2e-23 |
151 |
|
|
*/ |
152 |
|
|
|
153 |
|
|
static long double S[9] = { |
154 |
|
|
-1.193945051381510095614E-3L, |
155 |
|
|
7.220599478036909672331E-3L, |
156 |
|
|
-9.622023360406271645744E-3L, |
157 |
|
|
-4.219773360705915470089E-2L, |
158 |
|
|
1.665386113720805206758E-1L, |
159 |
|
|
-4.200263503403344054473E-2L, |
160 |
|
|
-6.558780715202540684668E-1L, |
161 |
|
|
5.772156649015328608253E-1L, |
162 |
|
|
1.000000000000000000000E0L, |
163 |
|
|
}; |
164 |
|
|
|
165 |
|
|
/* 1/tgamma(-x) = z P(z) |
166 |
|
|
* z(x) = 1/x |
167 |
|
|
* 0 < x < 0.03125 |
168 |
|
|
* Peak relative error 5.16e-23 |
169 |
|
|
* Relative error spread = 2.5e-24 |
170 |
|
|
*/ |
171 |
|
|
|
172 |
|
|
static long double SN[9] = { |
173 |
|
|
1.133374167243894382010E-3L, |
174 |
|
|
7.220837261893170325704E-3L, |
175 |
|
|
9.621911155035976733706E-3L, |
176 |
|
|
-4.219773343731191721664E-2L, |
177 |
|
|
-1.665386113944413519335E-1L, |
178 |
|
|
-4.200263503402112910504E-2L, |
179 |
|
|
6.558780715202536547116E-1L, |
180 |
|
|
5.772156649015328608727E-1L, |
181 |
|
|
-1.000000000000000000000E0L, |
182 |
|
|
}; |
183 |
|
|
|
184 |
|
|
static const long double PIL = 3.1415926535897932384626L; |
185 |
|
|
|
186 |
|
|
static long double stirf ( long double ); |
187 |
|
|
|
188 |
|
|
/* Gamma function computed by Stirling's formula. |
189 |
|
|
*/ |
190 |
|
|
static long double stirf(long double x) |
191 |
|
|
{ |
192 |
|
|
long double y, w, v; |
193 |
|
|
|
194 |
|
|
w = 1.0L/x; |
195 |
|
|
/* For large x, use rational coefficients from the analytical expansion. */ |
196 |
|
|
if( x > 1024.0L ) |
197 |
|
|
w = (((((6.97281375836585777429E-5L * w |
198 |
|
|
+ 7.84039221720066627474E-4L) * w |
199 |
|
|
- 2.29472093621399176955E-4L) * w |
200 |
|
|
- 2.68132716049382716049E-3L) * w |
201 |
|
|
+ 3.47222222222222222222E-3L) * w |
202 |
|
|
+ 8.33333333333333333333E-2L) * w |
203 |
|
|
+ 1.0L; |
204 |
|
|
else |
205 |
|
|
w = 1.0L + w * __polevll( w, STIR, 8 ); |
206 |
|
|
y = expl(x); |
207 |
|
|
if( x > MAXSTIR ) |
208 |
|
|
{ /* Avoid overflow in pow() */ |
209 |
|
|
v = powl( x, 0.5L * x - 0.25L ); |
210 |
|
|
y = v * (v / y); |
211 |
|
|
} |
212 |
|
|
else |
213 |
|
|
{ |
214 |
|
|
y = powl( x, x - 0.5L ) / y; |
215 |
|
|
} |
216 |
|
|
y = SQTPI * y * w; |
217 |
|
|
return( y ); |
218 |
|
|
} |
219 |
|
|
|
220 |
|
|
long double |
221 |
|
|
tgammal(long double x) |
222 |
|
|
{ |
223 |
|
|
long double p, q, z; |
224 |
|
|
int i; |
225 |
|
|
|
226 |
|
|
signgam = 1; |
227 |
|
|
if( isnan(x) ) |
228 |
|
|
return(NAN); |
229 |
|
|
if(x == INFINITY) |
230 |
|
|
return(INFINITY); |
231 |
|
|
if(x == -INFINITY) |
232 |
|
|
return(x - x); |
233 |
|
|
if( x == 0.0L ) |
234 |
|
|
return( 1.0L / x ); |
235 |
|
|
q = fabsl(x); |
236 |
|
|
|
237 |
|
|
if( q > 13.0L ) |
238 |
|
|
{ |
239 |
|
|
if( q > MAXGAML ) |
240 |
|
|
goto goverf; |
241 |
|
|
if( x < 0.0L ) |
242 |
|
|
{ |
243 |
|
|
p = floorl(q); |
244 |
|
|
if( p == q ) |
245 |
|
|
return (x - x) / (x - x); |
246 |
|
|
i = p; |
247 |
|
|
if( (i & 1) == 0 ) |
248 |
|
|
signgam = -1; |
249 |
|
|
z = q - p; |
250 |
|
|
if( z > 0.5L ) |
251 |
|
|
{ |
252 |
|
|
p += 1.0L; |
253 |
|
|
z = q - p; |
254 |
|
|
} |
255 |
|
|
z = q * sinl( PIL * z ); |
256 |
|
|
z = fabsl(z) * stirf(q); |
257 |
|
|
if( z <= PIL/LDBL_MAX ) |
258 |
|
|
{ |
259 |
|
|
goverf: |
260 |
|
|
return( signgam * INFINITY); |
261 |
|
|
} |
262 |
|
|
z = PIL/z; |
263 |
|
|
} |
264 |
|
|
else |
265 |
|
|
{ |
266 |
|
|
z = stirf(x); |
267 |
|
|
} |
268 |
|
|
return( signgam * z ); |
269 |
|
|
} |
270 |
|
|
|
271 |
|
|
z = 1.0L; |
272 |
|
|
while( x >= 3.0L ) |
273 |
|
|
{ |
274 |
|
|
x -= 1.0L; |
275 |
|
|
z *= x; |
276 |
|
|
} |
277 |
|
|
|
278 |
|
|
while( x < -0.03125L ) |
279 |
|
|
{ |
280 |
|
|
z /= x; |
281 |
|
|
x += 1.0L; |
282 |
|
|
} |
283 |
|
|
|
284 |
|
|
if( x <= 0.03125L ) |
285 |
|
|
goto small; |
286 |
|
|
|
287 |
|
|
while( x < 2.0L ) |
288 |
|
|
{ |
289 |
|
|
z /= x; |
290 |
|
|
x += 1.0L; |
291 |
|
|
} |
292 |
|
|
|
293 |
|
|
if( x == 2.0L ) |
294 |
|
|
return(z); |
295 |
|
|
|
296 |
|
|
x -= 2.0L; |
297 |
|
|
p = __polevll( x, P, 7 ); |
298 |
|
|
q = __polevll( x, Q, 8 ); |
299 |
|
|
z = z * p / q; |
300 |
|
|
if( z < 0 ) |
301 |
|
|
signgam = -1; |
302 |
|
|
return z; |
303 |
|
|
|
304 |
|
|
small: |
305 |
|
|
if( x == 0.0L ) |
306 |
|
|
return (x - x) / (x - x); |
307 |
|
|
else |
308 |
|
|
{ |
309 |
|
|
if( x < 0.0L ) |
310 |
|
|
{ |
311 |
|
|
x = -x; |
312 |
|
|
q = z / (x * __polevll( x, SN, 8 )); |
313 |
|
|
signgam = -1; |
314 |
|
|
} |
315 |
|
|
else |
316 |
|
|
q = z / (x * __polevll( x, S, 8 )); |
317 |
|
|
} |
318 |
|
|
return q; |
319 |
|
|
} |