[ad9af31] | 1 | /* |
---|
| 2 | The wrapper for gamma function from OpenCL and standard libraries |
---|
[98cb4d7] | 3 | The OpenCL gamma function fails miserably on values lower than 1.0 |
---|
[ad9af31] | 4 | while works fine on larger values. |
---|
| 5 | We use gamma definition Gamma(t + 1) = t * Gamma(t) to compute |
---|
| 6 | to function for values lower than 1.0. Namely Gamma(t) = 1/t * Gamma(t + 1) |
---|
| 7 | */ |
---|
| 8 | |
---|
[6ce9048] | 9 | #if defined(NEED_TGAMMA) |
---|
| 10 | static double cephes_stirf(double x) |
---|
| 11 | { |
---|
| 12 | const double MAXSTIR=143.01608; |
---|
| 13 | const double SQTPI=2.50662827463100050242E0; |
---|
| 14 | double y, w, v; |
---|
| 15 | |
---|
| 16 | w = 1.0 / x; |
---|
| 17 | |
---|
| 18 | w = (((( |
---|
| 19 | 7.87311395793093628397E-4*w |
---|
| 20 | - 2.29549961613378126380E-4)*w |
---|
| 21 | - 2.68132617805781232825E-3)*w |
---|
| 22 | + 3.47222221605458667310E-3)*w |
---|
| 23 | + 8.33333333333482257126E-2)*w |
---|
| 24 | + 1.0; |
---|
| 25 | y = exp(x); |
---|
| 26 | if (x > MAXSTIR) |
---|
| 27 | { /* Avoid overflow in pow() */ |
---|
| 28 | v = pow(x, 0.5 * x - 0.25); |
---|
| 29 | y = v * (v / y); |
---|
| 30 | } |
---|
| 31 | else |
---|
| 32 | { |
---|
| 33 | y = pow(x, x - 0.5) / y; |
---|
| 34 | } |
---|
| 35 | y = SQTPI * y * w; |
---|
| 36 | return(y); |
---|
| 37 | } |
---|
| 38 | |
---|
| 39 | static double tgamma(double x) { |
---|
| 40 | double p, q, z; |
---|
| 41 | int sgngam; |
---|
| 42 | int i; |
---|
| 43 | |
---|
| 44 | sgngam = 1; |
---|
| 45 | if (isnan(x)) |
---|
| 46 | return(x); |
---|
| 47 | q = fabs(x); |
---|
| 48 | |
---|
| 49 | if (q > 33.0) |
---|
| 50 | { |
---|
| 51 | if (x < 0.0) |
---|
| 52 | { |
---|
| 53 | p = floor(q); |
---|
| 54 | if (p == q) |
---|
| 55 | { |
---|
| 56 | return (NAN); |
---|
| 57 | } |
---|
| 58 | i = p; |
---|
| 59 | if ((i & 1) == 0) |
---|
| 60 | sgngam = -1; |
---|
| 61 | z = q - p; |
---|
| 62 | if (z > 0.5) |
---|
| 63 | { |
---|
| 64 | p += 1.0; |
---|
| 65 | z = q - p; |
---|
| 66 | } |
---|
| 67 | z = q * sin(M_PI * z); |
---|
| 68 | if (z == 0.0) |
---|
| 69 | { |
---|
| 70 | return(NAN); |
---|
| 71 | } |
---|
| 72 | z = fabs(z); |
---|
| 73 | z = M_PI / (z * cephes_stirf(q)); |
---|
| 74 | } |
---|
| 75 | else |
---|
| 76 | { |
---|
| 77 | z = cephes_stirf(x); |
---|
| 78 | } |
---|
| 79 | return(sgngam * z); |
---|
| 80 | } |
---|
| 81 | |
---|
| 82 | z = 1.0; |
---|
| 83 | while (x >= 3.0) |
---|
| 84 | { |
---|
| 85 | x -= 1.0; |
---|
| 86 | z *= x; |
---|
| 87 | } |
---|
| 88 | |
---|
| 89 | while (x < 0.0) |
---|
| 90 | { |
---|
| 91 | if (x > -1.E-9) |
---|
| 92 | goto small; |
---|
| 93 | z /= x; |
---|
| 94 | x += 1.0; |
---|
| 95 | } |
---|
| 96 | |
---|
| 97 | while (x < 2.0) |
---|
| 98 | { |
---|
| 99 | if (x < 1.e-9) |
---|
| 100 | goto small; |
---|
| 101 | z /= x; |
---|
| 102 | x += 1.0; |
---|
| 103 | } |
---|
| 104 | |
---|
| 105 | if (x == 2.0) |
---|
| 106 | return(z); |
---|
| 107 | |
---|
| 108 | x -= 2.0; |
---|
| 109 | p = ((((( |
---|
| 110 | +1.60119522476751861407E-4*x |
---|
| 111 | + 1.19135147006586384913E-3)*x |
---|
| 112 | + 1.04213797561761569935E-2)*x |
---|
| 113 | + 4.76367800457137231464E-2)*x |
---|
| 114 | + 2.07448227648435975150E-1)*x |
---|
| 115 | + 4.94214826801497100753E-1)*x |
---|
| 116 | + 9.99999999999999996796E-1; |
---|
| 117 | q = (((((( |
---|
| 118 | -2.31581873324120129819E-5*x |
---|
| 119 | + 5.39605580493303397842E-4)*x |
---|
| 120 | - 4.45641913851797240494E-3)*x |
---|
| 121 | + 1.18139785222060435552E-2)*x |
---|
| 122 | + 3.58236398605498653373E-2)*x |
---|
| 123 | - 2.34591795718243348568E-1)*x |
---|
| 124 | + 7.14304917030273074085E-2)*x |
---|
| 125 | + 1.00000000000000000320E0; |
---|
| 126 | return(z * p / q); |
---|
| 127 | |
---|
| 128 | small: |
---|
| 129 | if (x == 0.0) |
---|
| 130 | { |
---|
| 131 | return (NAN); |
---|
| 132 | } |
---|
| 133 | else |
---|
| 134 | return(z / ((1.0 + 0.5772156649015329 * x) * x)); |
---|
| 135 | } |
---|
[1596de3] | 136 | #endif // NEED_TGAMMA |
---|
[6ce9048] | 137 | |
---|
| 138 | |
---|
| 139 | inline double sas_gamma( double x) { return tgamma(x+1)/x; } |
---|