Go to most recent revision | Details | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1906 | serge | 1 | /* lgaml() |
2 | * |
||
3 | * Natural logarithm of gamma function |
||
4 | * |
||
5 | * |
||
6 | * |
||
7 | * SYNOPSIS: |
||
8 | * |
||
9 | * long double x, y, __lgammal_r(); |
||
10 | * int* sgngaml; |
||
11 | * y = __lgammal_r( x, sgngaml ); |
||
12 | * |
||
13 | * long double x, y, lgammal(); |
||
14 | * y = lgammal( x); |
||
15 | * |
||
16 | * |
||
17 | * |
||
18 | * DESCRIPTION: |
||
19 | * |
||
20 | * Returns the base e (2.718...) logarithm of the absolute |
||
21 | * value of the gamma function of the argument. In the reentrant |
||
22 | * version, the sign (+1 or -1) of the gamma function is returned |
||
23 | * in the variable referenced by sgngaml. |
||
24 | * |
||
25 | * For arguments greater than 33, the logarithm of the gamma |
||
26 | * function is approximated by the logarithmic version of |
||
27 | * Stirling's formula using a polynomial approximation of |
||
28 | * degree 4. Arguments between -33 and +33 are reduced by |
||
29 | * recurrence to the interval [2,3] of a rational approximation. |
||
30 | * The cosecant reflection formula is employed for arguments |
||
31 | * less than -33. |
||
32 | * |
||
33 | * Arguments greater than MAXLGML (10^4928) return MAXNUML. |
||
34 | * |
||
35 | * |
||
36 | * |
||
37 | * ACCURACY: |
||
38 | * |
||
39 | * |
||
40 | * arithmetic domain # trials peak rms |
||
41 | * IEEE -40, 40 100000 2.2e-19 4.6e-20 |
||
42 | * IEEE 10^-2000,10^+2000 20000 1.6e-19 3.3e-20 |
||
43 | * The error criterion was relative when the function magnitude |
||
44 | * was greater than one but absolute when it was less than one. |
||
45 | * |
||
46 | */ |
||
47 | |||
48 | /* |
||
49 | * Copyright 1994 by Stephen L. Moshier |
||
50 | */ |
||
51 | |||
52 | /* |
||
53 | * 26-11-2002 Modified for mingw. |
||
54 | * Danny Smith |
||
55 | */ |
||
56 | |||
57 | #ifndef __MINGW32__ |
||
58 | #include "mconf.h" |
||
59 | #ifdef ANSIPROT |
||
60 | extern long double fabsl ( long double ); |
||
61 | extern long double lgaml ( long double ); |
||
62 | extern long double logl ( long double ); |
||
63 | extern long double expl ( long double ); |
||
64 | extern long double gammal ( long double ); |
||
65 | extern long double sinl ( long double ); |
||
66 | extern long double floorl ( long double ); |
||
67 | extern long double powl ( long double, long double ); |
||
68 | extern long double polevll ( long double, void *, int ); |
||
69 | extern long double p1evll ( long double, void *, int ); |
||
70 | extern int isnanl ( long double ); |
||
71 | extern int isfinitel ( long double ); |
||
72 | #else |
||
73 | long double fabsl(), lgaml(), logl(), expl(), gammal(), sinl(); |
||
74 | long double floorl(), powl(), polevll(), p1evll(), isnanl(), isfinitel(); |
||
75 | #endif |
||
76 | #ifdef INFINITIES |
||
77 | extern long double INFINITYL; |
||
78 | #endif |
||
79 | #ifdef NANS |
||
80 | extern long double NANL; |
||
81 | #endif |
||
82 | #else /* __MINGW32__ */ |
||
83 | #include "cephes_mconf.h" |
||
84 | #endif /* __MINGW32__ */ |
||
85 | |||
86 | #if UNK |
||
87 | static long double S[9] = { |
||
88 | -1.193945051381510095614E-3L, |
||
89 | 7.220599478036909672331E-3L, |
||
90 | -9.622023360406271645744E-3L, |
||
91 | -4.219773360705915470089E-2L, |
||
92 | 1.665386113720805206758E-1L, |
||
93 | -4.200263503403344054473E-2L, |
||
94 | -6.558780715202540684668E-1L, |
||
95 | 5.772156649015328608253E-1L, |
||
96 | 1.000000000000000000000E0L, |
||
97 | }; |
||
98 | #endif |
||
99 | #if IBMPC |
||
100 | static const unsigned short S[] = { |
||
101 | 0xbaeb,0xd6d3,0x25e5,0x9c7e,0xbff5, XPD |
||
102 | 0xfe9a,0xceb4,0xc74e,0xec9a,0x3ff7, XPD |
||
103 | 0x9225,0xdfef,0xb0e9,0x9da5,0xbff8, XPD |
||
104 | 0x10b0,0xec17,0x87dc,0xacd7,0xbffa, XPD |
||
105 | 0x6b8d,0x7515,0x1905,0xaa89,0x3ffc, XPD |
||
106 | 0xf183,0x126b,0xf47d,0xac0a,0xbffa, XPD |
||
107 | 0x7bf6,0x57d1,0xa013,0xa7e7,0xbffe, XPD |
||
108 | 0xc7a9,0x7db0,0x67e3,0x93c4,0x3ffe, XPD |
||
109 | 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD |
||
110 | }; |
||
111 | #endif |
||
112 | #if MIEEE |
||
113 | static long S[27] = { |
||
114 | 0xbff50000,0x9c7e25e5,0xd6d3baeb, |
||
115 | 0x3ff70000,0xec9ac74e,0xceb4fe9a, |
||
116 | 0xbff80000,0x9da5b0e9,0xdfef9225, |
||
117 | 0xbffa0000,0xacd787dc,0xec1710b0, |
||
118 | 0x3ffc0000,0xaa891905,0x75156b8d, |
||
119 | 0xbffa0000,0xac0af47d,0x126bf183, |
||
120 | 0xbffe0000,0xa7e7a013,0x57d17bf6, |
||
121 | 0x3ffe0000,0x93c467e3,0x7db0c7a9, |
||
122 | 0x3fff0000,0x80000000,0x00000000, |
||
123 | }; |
||
124 | #endif |
||
125 | |||
126 | #if UNK |
||
127 | static long double SN[9] = { |
||
128 | 1.133374167243894382010E-3L, |
||
129 | 7.220837261893170325704E-3L, |
||
130 | 9.621911155035976733706E-3L, |
||
131 | -4.219773343731191721664E-2L, |
||
132 | -1.665386113944413519335E-1L, |
||
133 | -4.200263503402112910504E-2L, |
||
134 | 6.558780715202536547116E-1L, |
||
135 | 5.772156649015328608727E-1L, |
||
136 | -1.000000000000000000000E0L, |
||
137 | }; |
||
138 | #endif |
||
139 | #if IBMPC |
||
140 | static const unsigned SN[] = { |
||
141 | 0x5dd1,0x02de,0xb9f7,0x948d,0x3ff5, XPD |
||
142 | 0x989b,0xdd68,0xc5f1,0xec9c,0x3ff7, XPD |
||
143 | 0x2ca1,0x18f0,0x386f,0x9da5,0x3ff8, XPD |
||
144 | 0x783f,0x41dd,0x87d1,0xacd7,0xbffa, XPD |
||
145 | 0x7a5b,0xd76d,0x1905,0xaa89,0xbffc, XPD |
||
146 | 0x7f64,0x1234,0xf47d,0xac0a,0xbffa, XPD |
||
147 | 0x5e26,0x57d1,0xa013,0xa7e7,0x3ffe, XPD |
||
148 | 0xc7aa,0x7db0,0x67e3,0x93c4,0x3ffe, XPD |
||
149 | 0x0000,0x0000,0x0000,0x8000,0xbfff, XPD |
||
150 | }; |
||
151 | #endif |
||
152 | #if MIEEE |
||
153 | static long SN[27] = { |
||
154 | 0x3ff50000,0x948db9f7,0x02de5dd1, |
||
155 | 0x3ff70000,0xec9cc5f1,0xdd68989b, |
||
156 | 0x3ff80000,0x9da5386f,0x18f02ca1, |
||
157 | 0xbffa0000,0xacd787d1,0x41dd783f, |
||
158 | 0xbffc0000,0xaa891905,0xd76d7a5b, |
||
159 | 0xbffa0000,0xac0af47d,0x12347f64, |
||
160 | 0x3ffe0000,0xa7e7a013,0x57d15e26, |
||
161 | 0x3ffe0000,0x93c467e3,0x7db0c7aa, |
||
162 | 0xbfff0000,0x80000000,0x00000000, |
||
163 | }; |
||
164 | #endif |
||
165 | |||
166 | |||
167 | /* A[]: Stirling's formula expansion of log gamma |
||
168 | * B[], C[]: log gamma function between 2 and 3 |
||
169 | */ |
||
170 | |||
171 | |||
172 | /* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x A(1/x^2) |
||
173 | * x >= 8 |
||
174 | * Peak relative error 1.51e-21 |
||
175 | * Relative spread of error peaks 5.67e-21 |
||
176 | */ |
||
177 | #if UNK |
||
178 | static long double A[7] = { |
||
179 | 4.885026142432270781165E-3L, |
||
180 | -1.880801938119376907179E-3L, |
||
181 | 8.412723297322498080632E-4L, |
||
182 | -5.952345851765688514613E-4L, |
||
183 | 7.936507795855070755671E-4L, |
||
184 | -2.777777777750349603440E-3L, |
||
185 | 8.333333333333331447505E-2L, |
||
186 | }; |
||
187 | #endif |
||
188 | #if IBMPC |
||
189 | static const unsigned short A[] = { |
||
190 | 0xd984,0xcc08,0x91c2,0xa012,0x3ff7, XPD |
||
191 | 0x3d91,0x0304,0x3da1,0xf685,0xbff5, XPD |
||
192 | 0x3bdc,0xaad1,0xd492,0xdc88,0x3ff4, XPD |
||
193 | 0x8b20,0x9fce,0x844e,0x9c09,0xbff4, XPD |
||
194 | 0xf8f2,0x30e5,0x0092,0xd00d,0x3ff4, XPD |
||
195 | 0x4d88,0x03a8,0x60b6,0xb60b,0xbff6, XPD |
||
196 | 0x9fcc,0xaaaa,0xaaaa,0xaaaa,0x3ffb, XPD |
||
197 | }; |
||
198 | #endif |
||
199 | #if MIEEE |
||
200 | static long A[21] = { |
||
201 | 0x3ff70000,0xa01291c2,0xcc08d984, |
||
202 | 0xbff50000,0xf6853da1,0x03043d91, |
||
203 | 0x3ff40000,0xdc88d492,0xaad13bdc, |
||
204 | 0xbff40000,0x9c09844e,0x9fce8b20, |
||
205 | 0x3ff40000,0xd00d0092,0x30e5f8f2, |
||
206 | 0xbff60000,0xb60b60b6,0x03a84d88, |
||
207 | 0x3ffb0000,0xaaaaaaaa,0xaaaa9fcc, |
||
208 | }; |
||
209 | #endif |
||
210 | |||
211 | /* log gamma(x+2) = x B(x)/C(x) |
||
212 | * 0 <= x <= 1 |
||
213 | * Peak relative error 7.16e-22 |
||
214 | * Relative spread of error peaks 4.78e-20 |
||
215 | */ |
||
216 | #if UNK |
||
217 | static long double B[7] = { |
||
218 | -2.163690827643812857640E3L, |
||
219 | -8.723871522843511459790E4L, |
||
220 | -1.104326814691464261197E6L, |
||
221 | -6.111225012005214299996E6L, |
||
222 | -1.625568062543700591014E7L, |
||
223 | -2.003937418103815175475E7L, |
||
224 | -8.875666783650703802159E6L, |
||
225 | }; |
||
226 | static long double C[7] = { |
||
227 | /* 1.000000000000000000000E0L,*/ |
||
228 | -5.139481484435370143617E2L, |
||
229 | -3.403570840534304670537E4L, |
||
230 | -6.227441164066219501697E5L, |
||
231 | -4.814940379411882186630E6L, |
||
232 | -1.785433287045078156959E7L, |
||
233 | -3.138646407656182662088E7L, |
||
234 | -2.099336717757895876142E7L, |
||
235 | }; |
||
236 | #endif |
||
237 | #if IBMPC |
||
238 | static const unsigned short B[] = { |
||
239 | 0x9557,0x4995,0x0da1,0x873b,0xc00a, XPD |
||
240 | 0xfe44,0x9af8,0x5b8c,0xaa63,0xc00f, XPD |
||
241 | 0x5aa8,0x7cf5,0x3684,0x86ce,0xc013, XPD |
||
242 | 0x259a,0x258c,0xf206,0xba7f,0xc015, XPD |
||
243 | 0xbe18,0x1ca3,0xc0a0,0xf80a,0xc016, XPD |
||
244 | 0x168f,0x2c42,0x6717,0x98e3,0xc017, XPD |
||
245 | 0x2051,0x9d55,0x92c8,0x876e,0xc016, XPD |
||
246 | }; |
||
247 | static const unsigned short C[] = { |
||
248 | /*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/ |
||
249 | 0xaa77,0xcf2f,0xae76,0x807c,0xc008, XPD |
||
250 | 0xb280,0x0d74,0xb55a,0x84f3,0xc00e, XPD |
||
251 | 0xa505,0xcd30,0x81dc,0x9809,0xc012, XPD |
||
252 | 0x3369,0x4246,0xb8c2,0x92f0,0xc015, XPD |
||
253 | 0x63cf,0x6aee,0xbe6f,0x8837,0xc017, XPD |
||
254 | 0x26bb,0xccc7,0xb009,0xef75,0xc017, XPD |
||
255 | 0x462b,0xbae8,0xab96,0xa02a,0xc017, XPD |
||
256 | }; |
||
257 | #endif |
||
258 | #if MIEEE |
||
259 | static long B[21] = { |
||
260 | 0xc00a0000,0x873b0da1,0x49959557, |
||
261 | 0xc00f0000,0xaa635b8c,0x9af8fe44, |
||
262 | 0xc0130000,0x86ce3684,0x7cf55aa8, |
||
263 | 0xc0150000,0xba7ff206,0x258c259a, |
||
264 | 0xc0160000,0xf80ac0a0,0x1ca3be18, |
||
265 | 0xc0170000,0x98e36717,0x2c42168f, |
||
266 | 0xc0160000,0x876e92c8,0x9d552051, |
||
267 | }; |
||
268 | static long C[21] = { |
||
269 | /*0x3fff0000,0x80000000,0x00000000,*/ |
||
270 | 0xc0080000,0x807cae76,0xcf2faa77, |
||
271 | 0xc00e0000,0x84f3b55a,0x0d74b280, |
||
272 | 0xc0120000,0x980981dc,0xcd30a505, |
||
273 | 0xc0150000,0x92f0b8c2,0x42463369, |
||
274 | 0xc0170000,0x8837be6f,0x6aee63cf, |
||
275 | 0xc0170000,0xef75b009,0xccc726bb, |
||
276 | 0xc0170000,0xa02aab96,0xbae8462b, |
||
277 | }; |
||
278 | #endif |
||
279 | |||
280 | /* log( sqrt( 2*pi ) ) */ |
||
281 | static const long double LS2PI = 0.91893853320467274178L; |
||
282 | #define MAXLGM 1.04848146839019521116e+4928L |
||
283 | |||
284 | |||
285 | /* Logarithm of gamma function */ |
||
286 | /* Reentrant version */ |
||
287 | |||
288 | long double __lgammal_r(long double x, int* sgngaml) |
||
289 | { |
||
290 | long double p, q, w, z, f, nx; |
||
291 | int i; |
||
292 | |||
293 | *sgngaml = 1; |
||
294 | #ifdef NANS |
||
295 | if( isnanl(x) ) |
||
296 | return(NANL); |
||
297 | #endif |
||
298 | #ifdef INFINITIES |
||
299 | if( !isfinitel(x) ) |
||
300 | return(INFINITYL); |
||
301 | #endif |
||
302 | if( x < -34.0L ) |
||
303 | { |
||
304 | q = -x; |
||
305 | w = __lgammal_r(q, sgngaml); /* note this modifies sgngam! */ |
||
306 | p = floorl(q); |
||
307 | if( p == q ) |
||
308 | { |
||
309 | lgsing: |
||
310 | _SET_ERRNO(EDOM); |
||
311 | mtherr( "lgammal", SING ); |
||
312 | #ifdef INFINITIES |
||
313 | return (INFINITYL); |
||
314 | #else |
||
315 | return (MAXNUML); |
||
316 | #endif |
||
317 | } |
||
318 | i = p; |
||
319 | if( (i & 1) == 0 ) |
||
320 | *sgngaml = -1; |
||
321 | else |
||
322 | *sgngaml = 1; |
||
323 | z = q - p; |
||
324 | if( z > 0.5L ) |
||
325 | { |
||
326 | p += 1.0L; |
||
327 | z = p - q; |
||
328 | } |
||
329 | z = q * sinl( PIL * z ); |
||
330 | if( z == 0.0L ) |
||
331 | goto lgsing; |
||
332 | /* z = LOGPI - logl( z ) - w; */ |
||
333 | z = logl( PIL/z ) - w; |
||
334 | return( z ); |
||
335 | } |
||
336 | |||
337 | if( x < 13.0L ) |
||
338 | { |
||
339 | z = 1.0L; |
||
340 | nx = floorl( x + 0.5L ); |
||
341 | f = x - nx; |
||
342 | while( x >= 3.0L ) |
||
343 | { |
||
344 | nx -= 1.0L; |
||
345 | x = nx + f; |
||
346 | z *= x; |
||
347 | } |
||
348 | while( x < 2.0L ) |
||
349 | { |
||
350 | if( fabsl(x) <= 0.03125 ) |
||
351 | goto lsmall; |
||
352 | z /= nx + f; |
||
353 | nx += 1.0L; |
||
354 | x = nx + f; |
||
355 | } |
||
356 | if( z < 0.0L ) |
||
357 | { |
||
358 | *sgngaml = -1; |
||
359 | z = -z; |
||
360 | } |
||
361 | else |
||
362 | *sgngaml = 1; |
||
363 | if( x == 2.0L ) |
||
364 | return( logl(z) ); |
||
365 | x = (nx - 2.0L) + f; |
||
366 | p = x * polevll( x, B, 6 ) / p1evll( x, C, 7); |
||
367 | return( logl(z) + p ); |
||
368 | } |
||
369 | |||
370 | if( x > MAXLGM ) |
||
371 | { |
||
372 | _SET_ERRNO(ERANGE); |
||
373 | mtherr( "lgammal", OVERFLOW ); |
||
374 | #ifdef INFINITIES |
||
375 | return( *sgngaml * INFINITYL ); |
||
376 | #else |
||
377 | return( *sgngaml * MAXNUML ); |
||
378 | #endif |
||
379 | } |
||
380 | |||
381 | q = ( x - 0.5L ) * logl(x) - x + LS2PI; |
||
382 | if( x > 1.0e10L ) |
||
383 | return(q); |
||
384 | p = 1.0L/(x*x); |
||
385 | q += polevll( p, A, 6 ) / x; |
||
386 | return( q ); |
||
387 | |||
388 | |||
389 | lsmall: |
||
390 | if( x == 0.0L ) |
||
391 | goto lgsing; |
||
392 | if( x < 0.0L ) |
||
393 | { |
||
394 | x = -x; |
||
395 | q = z / (x * polevll( x, SN, 8 )); |
||
396 | } |
||
397 | else |
||
398 | q = z / (x * polevll( x, S, 8 )); |
||
399 | if( q < 0.0L ) |
||
400 | { |
||
401 | *sgngaml = -1; |
||
402 | q = -q; |
||
403 | } |
||
404 | else |
||
405 | *sgngaml = 1; |
||
406 | q = logl( q ); |
||
407 | return(q); |
||
408 | } |
||
409 | |||
410 | /* This is the C99 version */ |
||
411 | |||
412 | long double lgammal(long double x) |
||
413 | { |
||
414 | int local_sgngaml=0; |
||
415 | return (__lgammal_r(x, &local_sgngaml)); |
||
416 | }>>>=>>>>=>=> |