[fba9ca0] | 1 | #if FLOAT_SIZE>4 // double precision |
---|
| 2 | // based on cephes/double/igam.c |
---|
| 3 | double gammaln(double x); |
---|
| 4 | double gammainc(double a, double x); |
---|
| 5 | double gammaincc(double a, double x); |
---|
| 6 | |
---|
| 7 | double cephes_igamc(double a, double x); |
---|
| 8 | double cephes_igam(double a, double x); |
---|
| 9 | double cephes_lgam(double x); |
---|
| 10 | double cephes_lgam2(double x); |
---|
| 11 | |
---|
| 12 | double gammainc(double a, double x) |
---|
| 13 | { |
---|
| 14 | if ((x <= 0) || ( a <= 0)) return 0.0; |
---|
| 15 | if ((x > 1.0) && (x > a)) return 1.0 - cephes_igamc(a, x); |
---|
| 16 | return cephes_igam(a, x); |
---|
| 17 | } |
---|
| 18 | |
---|
| 19 | double gammaincc(double a, double x) |
---|
| 20 | { |
---|
| 21 | if ((x <= 0) || (a <= 0)) return 1.0; |
---|
| 22 | if ((x < 1.0) || (x < a)) return 1.0 - cephes_igam(a, x); |
---|
| 23 | return cephes_igamc(a, x); |
---|
| 24 | } |
---|
| 25 | |
---|
| 26 | double gammaln(double x) |
---|
| 27 | { |
---|
| 28 | if (isnan(x)) return(x); |
---|
| 29 | if (!isfinite(x)) return(INFINITY); |
---|
| 30 | return cephes_lgam(x); |
---|
| 31 | } |
---|
| 32 | |
---|
| 33 | double cephes_igamc(double a, double x) |
---|
| 34 | { |
---|
| 35 | const double MACHEP = 1.11022302462515654042E-16; // IEEE 2**-53 |
---|
| 36 | const double MAXLOG = 7.09782712893383996843E2; // IEEE log(2**1024) denormalized |
---|
| 37 | const double BIG = 4.503599627370496e15; |
---|
| 38 | const double BIGINV = 2.22044604925031308085e-16; |
---|
| 39 | double ans, ax, c, yc, r, t, y, z; |
---|
| 40 | double pk, pkm1, pkm2, qk, qkm1, qkm2; |
---|
| 41 | |
---|
| 42 | /* Compute x**a * exp(-x) / gamma(a) */ |
---|
| 43 | ax = a * log(x) - x - cephes_lgam(a); |
---|
| 44 | if (ax < -MAXLOG) return 0.0; // underflow |
---|
| 45 | ax = exp(ax); |
---|
| 46 | |
---|
| 47 | /* continued fraction */ |
---|
| 48 | y = 1.0 - a; |
---|
| 49 | z = x + y + 1.0; |
---|
| 50 | c = 0.0; |
---|
| 51 | pkm2 = 1.0; |
---|
| 52 | qkm2 = x; |
---|
| 53 | pkm1 = x + 1.0; |
---|
| 54 | qkm1 = z * x; |
---|
| 55 | ans = pkm1/qkm1; |
---|
| 56 | |
---|
| 57 | do { |
---|
| 58 | c += 1.0; |
---|
| 59 | y += 1.0; |
---|
| 60 | z += 2.0; |
---|
| 61 | yc = y * c; |
---|
| 62 | pk = pkm1 * z - pkm2 * yc; |
---|
| 63 | qk = qkm1 * z - qkm2 * yc; |
---|
| 64 | if (qk != 0) { |
---|
| 65 | r = pk/qk; |
---|
| 66 | t = fabs( (ans - r)/r ); |
---|
| 67 | ans = r; |
---|
| 68 | } else { |
---|
| 69 | t = 1.0; |
---|
| 70 | } |
---|
| 71 | pkm2 = pkm1; |
---|
| 72 | pkm1 = pk; |
---|
| 73 | qkm2 = qkm1; |
---|
| 74 | qkm1 = qk; |
---|
| 75 | if (fabs(pk) > BIG) { |
---|
| 76 | pkm2 *= BIGINV; |
---|
| 77 | pkm1 *= BIGINV; |
---|
| 78 | qkm2 *= BIGINV; |
---|
| 79 | qkm1 *= BIGINV; |
---|
| 80 | } |
---|
| 81 | } while( t > MACHEP ); |
---|
| 82 | |
---|
| 83 | return( ans * ax ); |
---|
| 84 | } |
---|
| 85 | |
---|
| 86 | |
---|
| 87 | double cephes_igam(double a, double x) |
---|
| 88 | { |
---|
| 89 | const double MACHEP = 1.11022302462515654042E-16; // IEEE 2**-53 |
---|
| 90 | const double MAXLOG = 7.09782712893383996843E2; // IEEE log(2**1024) denormalized |
---|
| 91 | double ans, ax, c, r; |
---|
| 92 | |
---|
| 93 | /* Compute x**a * exp(-x) / gamma(a) */ |
---|
| 94 | ax = a * log(x) - x - cephes_lgam(a); |
---|
| 95 | if (ax < -MAXLOG) return 0.0; // underflow |
---|
| 96 | ax = exp(ax); |
---|
| 97 | |
---|
| 98 | /* power series */ |
---|
| 99 | r = a; |
---|
| 100 | c = 1.0; |
---|
| 101 | ans = 1.0; |
---|
| 102 | |
---|
| 103 | do { |
---|
| 104 | r += 1.0; |
---|
| 105 | c *= x/r; |
---|
| 106 | ans += c; |
---|
| 107 | } while (c/ans > MACHEP); |
---|
| 108 | |
---|
| 109 | return ans * ax/a; |
---|
| 110 | } |
---|
| 111 | |
---|
| 112 | /* Logarithm of gamma function */ |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | double cephes_lgam(double x) |
---|
| 116 | { |
---|
| 117 | const double LOGPI = 1.14472988584940017414; // log(pi) |
---|
| 118 | int sgngam = 1; |
---|
| 119 | double p, q, w, z; |
---|
| 120 | int i; |
---|
| 121 | |
---|
| 122 | if (x < -34.0) { |
---|
| 123 | q = -x; |
---|
| 124 | w = cephes_lgam2(q); |
---|
| 125 | p = floor(q); |
---|
| 126 | if (p == q) return INFINITY; |
---|
| 127 | i = p; |
---|
| 128 | sgngam = ((i&1) == 0) ? -1 : 1; |
---|
| 129 | z = q - p; |
---|
| 130 | if (z > 0.5) { |
---|
| 131 | p += 1.0; |
---|
| 132 | z = p - q; |
---|
| 133 | } |
---|
| 134 | z = q * sin(M_PI * z); |
---|
| 135 | if (z == 0.0) return INFINITY; |
---|
| 136 | z = LOGPI - log(z) - w; |
---|
| 137 | return z; |
---|
| 138 | } else { |
---|
| 139 | return cephes_lgam2(x); |
---|
| 140 | } |
---|
| 141 | } |
---|
| 142 | |
---|
| 143 | double cephes_lgam2(double x) |
---|
| 144 | { |
---|
| 145 | const double LS2PI = 0.91893853320467274178; // log(sqrt(2*pi)) |
---|
| 146 | const double MAXLGM = 2.556348e305; |
---|
| 147 | int sgngam = 1; |
---|
| 148 | double p, q, u, z; |
---|
| 149 | double A, B, C; |
---|
| 150 | |
---|
| 151 | if (x < 13.0) { |
---|
| 152 | z = 1.0; |
---|
| 153 | p = 0.0; |
---|
| 154 | u = x; |
---|
| 155 | while (u >= 3.0) { |
---|
| 156 | p -= 1.0; |
---|
| 157 | u = x + p; |
---|
| 158 | z *= u; |
---|
| 159 | } |
---|
| 160 | while (u < 2.0) { |
---|
| 161 | if (u == 0.0) return INFINITY; |
---|
| 162 | z /= u; |
---|
| 163 | p += 1.0; |
---|
| 164 | u = x + p; |
---|
| 165 | } |
---|
| 166 | if (z < 0.0) { |
---|
| 167 | sgngam = -1; |
---|
| 168 | z = -z; |
---|
| 169 | } else { |
---|
| 170 | sgngam = 1; |
---|
| 171 | } |
---|
| 172 | if (u == 2.0) return log(z); |
---|
| 173 | p -= 2.0; |
---|
| 174 | x = x + p; |
---|
| 175 | B = (((((( |
---|
| 176 | -1.37825152569120859100E3)*x |
---|
| 177 | -3.88016315134637840924E4)*x |
---|
| 178 | -3.31612992738871184744E5)*x |
---|
| 179 | -1.16237097492762307383E6)*x |
---|
| 180 | -1.72173700820839662146E6)*x |
---|
| 181 | -8.53555664245765465627E5); |
---|
| 182 | C = (((((( |
---|
| 183 | /* 1.00000000000000000000E0)* */x |
---|
| 184 | -3.51815701436523470549E2)*x |
---|
| 185 | -1.70642106651881159223E4)*x |
---|
| 186 | -2.20528590553854454839E5)*x |
---|
| 187 | -1.13933444367982507207E6)*x |
---|
| 188 | -2.53252307177582951285E6)*x |
---|
| 189 | -2.01889141433532773231E6); |
---|
| 190 | p = x * B / C; |
---|
| 191 | return log(z) + p; |
---|
| 192 | } |
---|
| 193 | |
---|
| 194 | if (x > MAXLGM) return sgngam * INFINITY; |
---|
| 195 | |
---|
| 196 | q = (x - 0.5) * log(x) - x + LS2PI; |
---|
| 197 | if (x > 1.0e8) return q; |
---|
| 198 | |
---|
| 199 | p = 1.0/(x*x); |
---|
| 200 | if (x >= 1000.0) { |
---|
| 201 | q += ((7.9365079365079365079365e-4 * p |
---|
| 202 | - 2.7777777777777777777778e-3) * p |
---|
| 203 | + 0.0833333333333333333333) / x; |
---|
| 204 | } else { |
---|
| 205 | A = ((((( |
---|
| 206 | +8.11614167470508450300E-4)*p |
---|
| 207 | -5.95061904284301438324E-4)*p |
---|
| 208 | +7.93650340457716943945E-4)*p |
---|
| 209 | -2.77777777730099687205E-3)*p |
---|
| 210 | +8.33333333333331927722E-2); |
---|
| 211 | q += A / x; |
---|
| 212 | } |
---|
| 213 | return q; |
---|
| 214 | } |
---|
| 215 | |
---|
| 216 | #else // single precision |
---|
| 217 | |
---|
| 218 | // based on cephes/float/igam.c |
---|
| 219 | float gammalnf(float x); |
---|
| 220 | float gammaincf(float a, float x); |
---|
| 221 | float gammainccf(float a, float x); |
---|
| 222 | |
---|
| 223 | float cephes_igamcf(float a, float x); |
---|
| 224 | float cephes_igamf(float a, float x); |
---|
| 225 | float cephes_lgamf(float x); |
---|
| 226 | float cephes_lgam2f(float x); |
---|
| 227 | |
---|
| 228 | // Note: original uses logf, fabsf, floorf and sinf |
---|
| 229 | |
---|
| 230 | float gammaincf(float a, float x) |
---|
| 231 | { |
---|
| 232 | if ((x <= 0) || ( a <= 0)) return 0.0; |
---|
| 233 | if ((x > 1.0) && (x > a)) return 1.0 - cephes_igamcf(a, x); |
---|
| 234 | return cephes_igamf(a, x); |
---|
| 235 | } |
---|
| 236 | |
---|
| 237 | float gammainccf(float a, float x) |
---|
| 238 | { |
---|
| 239 | if ((x <= 0) || (a <= 0)) return 1.0; |
---|
| 240 | if ((x < 1.0) || (x < a)) return 1.0 - cephes_igamf(a, x); |
---|
| 241 | return cephes_igamcf(a, x); |
---|
| 242 | } |
---|
| 243 | |
---|
| 244 | float gammalnf(float x) |
---|
| 245 | { |
---|
| 246 | if (isnan(x)) return(x); |
---|
| 247 | if (!isfinite(x)) return(INFINITY); |
---|
| 248 | return cephes_lgamf(x); |
---|
| 249 | } |
---|
| 250 | |
---|
| 251 | float cephes_igamcf(float a, float x) |
---|
| 252 | { |
---|
| 253 | const float MAXLOGF = 88.72283905206835; |
---|
| 254 | const float MACHEPF = 5.9604644775390625E-8; |
---|
| 255 | const float BIG = 16777216.; |
---|
| 256 | float ans, c, yc, ax, y, z; |
---|
| 257 | float pk, pkm1, pkm2, qk, qkm1, qkm2; |
---|
| 258 | float r, t; |
---|
| 259 | |
---|
| 260 | ax = a * log(x) - x - cephes_lgamf(a); |
---|
| 261 | if (ax < -MAXLOGF) return 0.0; // underflow |
---|
| 262 | ax = expf(ax); |
---|
| 263 | |
---|
| 264 | /* continued fraction */ |
---|
| 265 | y = 1.0 - a; |
---|
| 266 | z = x + y + 1.0; |
---|
| 267 | c = 0.0; |
---|
| 268 | pkm2 = 1.0; |
---|
| 269 | qkm2 = x; |
---|
| 270 | pkm1 = x + 1.0; |
---|
| 271 | qkm1 = z * x; |
---|
| 272 | ans = pkm1/qkm1; |
---|
| 273 | |
---|
| 274 | do { |
---|
| 275 | c += 1.0; |
---|
| 276 | y += 1.0; |
---|
| 277 | z += 2.0; |
---|
| 278 | yc = y * c; |
---|
| 279 | pk = pkm1 * z - pkm2 * yc; |
---|
| 280 | qk = qkm1 * z - qkm2 * yc; |
---|
| 281 | if (qk != 0) { |
---|
| 282 | r = pk/qk; |
---|
| 283 | t = fabs((ans - r)/r); |
---|
| 284 | ans = r; |
---|
| 285 | } else { |
---|
| 286 | t = 1.0; |
---|
| 287 | } |
---|
| 288 | pkm2 = pkm1; |
---|
| 289 | pkm1 = pk; |
---|
| 290 | qkm2 = qkm1; |
---|
| 291 | qkm1 = qk; |
---|
| 292 | if (fabs(pk) > BIG) { |
---|
| 293 | pkm2 *= MACHEPF; |
---|
| 294 | pkm1 *= MACHEPF; |
---|
| 295 | qkm2 *= MACHEPF; |
---|
| 296 | qkm1 *= MACHEPF; |
---|
| 297 | } |
---|
| 298 | } while (t > MACHEPF); |
---|
| 299 | |
---|
| 300 | return ans * ax; |
---|
| 301 | } |
---|
| 302 | |
---|
| 303 | float cephes_igamf(float a, float x) |
---|
| 304 | { |
---|
| 305 | const float MAXLOGF = 88.72283905206835; |
---|
| 306 | const float MACHEPF = 5.9604644775390625E-8; |
---|
| 307 | float ans, ax, c, r; |
---|
| 308 | |
---|
| 309 | /* Compute x**a * exp(-x) / gamma(a) */ |
---|
| 310 | ax = a * log(x) - x - cephes_lgamf(a); |
---|
| 311 | if (ax < -MAXLOGF) return 0.0; // underflow |
---|
| 312 | ax = expf(ax); |
---|
| 313 | |
---|
| 314 | /* power series */ |
---|
| 315 | r = a; |
---|
| 316 | c = 1.0; |
---|
| 317 | ans = 1.0; |
---|
| 318 | |
---|
| 319 | do { |
---|
| 320 | r += 1.0; |
---|
| 321 | c *= x/r; |
---|
| 322 | ans += c; |
---|
| 323 | } while (c/ans > MACHEPF); |
---|
| 324 | |
---|
| 325 | return ans * ax/a; |
---|
| 326 | } |
---|
| 327 | |
---|
| 328 | float cephes_lgamf(float x) |
---|
| 329 | { |
---|
| 330 | const float PIINV = 0.318309886183790671538; |
---|
| 331 | float p, q, w, z; |
---|
| 332 | int i; |
---|
| 333 | int sgngamf; |
---|
| 334 | |
---|
| 335 | if (x < 0.0) { |
---|
| 336 | q = -x; |
---|
| 337 | w = cephes_lgam2f(q); |
---|
| 338 | p = floor(q); |
---|
| 339 | if (p == q) return INFINITY; |
---|
| 340 | i = p; |
---|
| 341 | sgngamf = ((i&1) == 0) ? -1 : 1; |
---|
| 342 | z = q - p; |
---|
| 343 | if (z > 0.5) { |
---|
| 344 | p += 1.0; |
---|
| 345 | z = p - q; |
---|
| 346 | } |
---|
| 347 | z = q * sin(M_PI * z); |
---|
| 348 | if (z == 0.0) return sgngamf * INFINITY; |
---|
| 349 | z = -log(PIINV*z) - w; |
---|
| 350 | return z; |
---|
| 351 | } else { |
---|
| 352 | return cephes_lgam2f(x); |
---|
| 353 | } |
---|
| 354 | } |
---|
| 355 | |
---|
| 356 | float cephes_lgam2f(float x) |
---|
| 357 | { |
---|
| 358 | const float LS2PI = 0.91893853320467274178; // log(sqrt(2*pi)) |
---|
| 359 | const float MAXLGM = 2.035093e36; |
---|
| 360 | float p, q, z; |
---|
| 361 | float nx, tx; |
---|
| 362 | int direction; |
---|
| 363 | int sgngamf = 1; |
---|
| 364 | |
---|
| 365 | if (x < 6.5) { |
---|
| 366 | direction = 0; |
---|
| 367 | z = 1.0; |
---|
| 368 | tx = x; |
---|
| 369 | nx = 0.0; |
---|
| 370 | if (x >= 1.5) { |
---|
| 371 | while (tx > 2.5) { |
---|
| 372 | nx -= 1.0; |
---|
| 373 | tx = x + nx; |
---|
| 374 | z *=tx; |
---|
| 375 | } |
---|
| 376 | x += nx - 2.0; |
---|
| 377 | } else if (x >= 1.25) { |
---|
| 378 | z *= x; |
---|
| 379 | x -= 1.0; /* x + 1 - 2 */ |
---|
| 380 | direction = 1; |
---|
| 381 | } else if (x >= 0.75) { |
---|
| 382 | x -= 1.0; |
---|
| 383 | /* log gamma(x+1), -.25 < x < .25 */ |
---|
| 384 | // p = x * polevlf( x, C, 7 ); |
---|
| 385 | p = (((((((( |
---|
| 386 | +1.369488127325832E-001)*x |
---|
| 387 | -1.590086327657347E-001)*x |
---|
| 388 | +1.692415923504637E-001)*x |
---|
| 389 | -2.067882815621965E-001)*x |
---|
| 390 | +2.705806208275915E-001)*x |
---|
| 391 | -4.006931650563372E-001)*x |
---|
| 392 | +8.224670749082976E-001)*x |
---|
| 393 | -5.772156501719101E-001)*x; |
---|
| 394 | q = 0.0; |
---|
| 395 | return p + q; |
---|
| 396 | } else { |
---|
| 397 | while (tx < 1.5) { |
---|
| 398 | if (tx == 0.0) return sgngamf * INFINITY; |
---|
| 399 | z *=tx; |
---|
| 400 | nx += 1.0; |
---|
| 401 | tx = x + nx; |
---|
| 402 | } |
---|
| 403 | direction = 1; |
---|
| 404 | x += nx - 2.0; |
---|
| 405 | } |
---|
| 406 | |
---|
| 407 | /* log gamma(x+2), -.5 < x < .5 */ |
---|
| 408 | // p = x * polevlf( x, B, 7 ); |
---|
| 409 | p = (((((((( |
---|
| 410 | +6.055172732649237E-004)*x |
---|
| 411 | -1.311620815545743E-003)*x |
---|
| 412 | +2.863437556468661E-003)*x |
---|
| 413 | -7.366775108654962E-003)*x |
---|
| 414 | +2.058355474821512E-002)*x |
---|
| 415 | -6.735323259371034E-002)*x |
---|
| 416 | +3.224669577325661E-001)*x |
---|
| 417 | +4.227843421859038E-001)*x; |
---|
| 418 | |
---|
| 419 | if (z < 0.0) { |
---|
| 420 | sgngamf = -1; |
---|
| 421 | z = -z; |
---|
| 422 | } else { |
---|
| 423 | sgngamf = 1; |
---|
| 424 | } |
---|
| 425 | q = log(z); |
---|
| 426 | if (direction) q = -q; |
---|
| 427 | return p + q; |
---|
| 428 | } else if (x > MAXLGM) { |
---|
| 429 | return sgngamf * INFINITY; |
---|
| 430 | } else { |
---|
| 431 | /* Note, though an asymptotic formula could be used for x >= 3, |
---|
| 432 | * there is cancellation error in the following if x < 6.5. |
---|
| 433 | */ |
---|
| 434 | q = LS2PI - x; |
---|
| 435 | q += (x - 0.5) * log(x); |
---|
| 436 | |
---|
| 437 | if (x <= 1.0e4) { |
---|
| 438 | z = 1.0/x; |
---|
| 439 | p = z * z; |
---|
| 440 | q += ((6.789774945028216E-004 * p |
---|
| 441 | - 2.769887652139868E-003 ) * p |
---|
| 442 | + 8.333316229807355E-002 ) * z; |
---|
| 443 | } |
---|
| 444 | return q; |
---|
| 445 | } |
---|
| 446 | } |
---|
| 447 | |
---|
| 448 | #endif // !single precision |
---|
| 449 | |
---|
| 450 | #if FLOAT_SIZE>4 |
---|
| 451 | #define sas_gammaln gammaln |
---|
| 452 | #define sas_gammainc gammainc |
---|
| 453 | #define sas_gammaincc gammaincc |
---|
| 454 | #else |
---|
| 455 | #define sas_gammaln gammalnf |
---|
| 456 | #define sas_gammainc gammaincf |
---|
| 457 | #define sas_gammaincc gammainccf |
---|
| 458 | #endif |
---|