| /** |
| * This file has no copyright assigned and is placed in the Public Domain. |
| * This file is part of the mingw-w64 runtime package. |
| * No warranty is given; refer to the file DISCLAIMER.PD within this package. |
| */ |
| #include "cephes_mconf.h" |
| |
| /* define MAXGAM 34.84425627277176174 */ |
| |
| /* Stirling's formula for the gamma function |
| * gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) ( 1 + 1/x P(1/x) ) |
| * .028 < 1/x < .1 |
| * relative error < 1.9e-11 |
| */ |
| static const float STIR[] = { |
| -2.705194986674176E-003, |
| 3.473255786154910E-003, |
| 8.333331788340907E-002, |
| }; |
| static const float MAXSTIR = 26.77; |
| static const float SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */ |
| |
| static float stirf(float); |
| |
| /* Gamma function computed by Stirling's formula, |
| * sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) |
| * The polynomial STIR is valid for 33 <= x <= 172. |
| */ |
| static float stirf( float x ) |
| { |
| float y, w, v; |
| |
| w = 1.0/x; |
| w = 1.0 + w * polevlf(w, STIR, 2); |
| y = expf(-x); |
| if (x > MAXSTIR) |
| { /* Avoid overflow in pow() */ |
| v = powf(x, 0.5 * x - 0.25); |
| y *= v; |
| y *= v; |
| } |
| else |
| { |
| y = powf(x, x - 0.5) * y; |
| } |
| y = SQTPIF * y * w; |
| return (y); |
| } |
| |
| |
| /* gamma(x+2), 0 < x < 1 */ |
| static const float P[] = { |
| 1.536830450601906E-003, |
| 5.397581592950993E-003, |
| 4.130370201859976E-003, |
| 7.232307985516519E-002, |
| 8.203960091619193E-002, |
| 4.117857447645796E-001, |
| 4.227867745131584E-001, |
| 9.999999822945073E-001, |
| }; |
| |
| float __tgammaf_r( float x, int* sgngamf); |
| |
| float __tgammaf_r( float x, int* sgngamf) |
| { |
| float p, q, z, nz; |
| int i, direction, negative; |
| |
| #ifdef NANS |
| if (isnan(x)) |
| return (x); |
| #endif |
| #ifdef INFINITIES |
| #ifdef NANS |
| if (x == INFINITYF) |
| return (x); |
| if (x == -INFINITYF) |
| return (NANF); |
| #else |
| if (!isfinite(x)) |
| return (x); |
| #endif |
| #endif |
| |
| *sgngamf = 1; |
| negative = 0; |
| nz = 0.0; |
| if (x < 0.0) |
| { |
| negative = 1; |
| q = -x; |
| p = floorf(q); |
| if (p == q) |
| { |
| gsing: |
| _SET_ERRNO(EDOM); |
| mtherr("tgammaf", SING); |
| #ifdef INFINITIES |
| return (INFINITYF); |
| #else |
| return (MAXNUMF); |
| #endif |
| } |
| i = p; |
| if ((i & 1) == 0) |
| *sgngamf = -1; |
| nz = q - p; |
| if (nz > 0.5) |
| { |
| p += 1.0; |
| nz = q - p; |
| } |
| nz = q * sinf(PIF * nz); |
| if (nz == 0.0) |
| { |
| _SET_ERRNO(ERANGE); |
| mtherr("tgamma", OVERFLOW); |
| #ifdef INFINITIES |
| return(*sgngamf * INFINITYF); |
| #else |
| return(*sgngamf * MAXNUMF); |
| #endif |
| } |
| if (nz < 0) |
| nz = -nz; |
| x = q; |
| } |
| if (x >= 10.0) |
| { |
| z = stirf(x); |
| } |
| if (x < 2.0) |
| direction = 1; |
| else |
| direction = 0; |
| z = 1.0; |
| while (x >= 3.0) |
| { |
| x -= 1.0; |
| z *= x; |
| } |
| /* |
| while (x < 0.0) |
| { |
| if (x > -1.E-4) |
| goto Small; |
| z *=x; |
| x += 1.0; |
| } |
| */ |
| while (x < 2.0) |
| { |
| if (x < 1.e-4) |
| goto Small; |
| z *=x; |
| x += 1.0; |
| } |
| |
| if (direction) |
| z = 1.0/z; |
| |
| if (x == 2.0) |
| return (z); |
| |
| x -= 2.0; |
| p = z * polevlf(x, P, 7); |
| |
| gdone: |
| if (negative) |
| { |
| p = *sgngamf * PIF/(nz * p ); |
| } |
| return (p); |
| |
| Small: |
| if (x == 0.0) |
| { |
| goto gsing; |
| } |
| else |
| { |
| p = z / ((1.0 + 0.5772156649015329 * x) * x); |
| goto gdone; |
| } |
| } |
| |
| /* This is the C99 version */ |
| float tgammaf(float x) |
| { |
| int local_sgngamf = 0; |
| return (__tgammaf_r(x, &local_sgngamf)); |
| } |
| |