* log-gamma function to support some of these distributions. The * algorithm comes from SPECFUN by Shanjie Zhang and Jianming Jin and their * book "Computation of Special Functions", 1996, John Wiley & Sons, Inc. * * If random_loggam(k+1) is being used to compute log(k!) for an integer k, consider * using logfactorial(k) instead. */
| 342 | * using logfactorial(k) instead. |
| 343 | */ |
| 344 | double random_loggam(double x) { |
| 345 | double x0, x2, lg2pi, gl, gl0; |
| 346 | RAND_INT_TYPE k, n; |
| 347 | |
| 348 | static double a[10] = {8.333333333333333e-02, -2.777777777777778e-03, |
| 349 | 7.936507936507937e-04, -5.952380952380952e-04, |
| 350 | 8.417508417508418e-04, -1.917526917526918e-03, |
| 351 | 6.410256410256410e-03, -2.955065359477124e-02, |
| 352 | 1.796443723688307e-01, -1.39243221690590e+00}; |
| 353 | |
| 354 | if ((x == 1.0) || (x == 2.0)) { |
| 355 | return 0.0; |
| 356 | } else if (x < 7.0) { |
| 357 | n = (RAND_INT_TYPE)(7 - x); |
| 358 | } else { |
| 359 | n = 0; |
| 360 | } |
| 361 | x0 = x + n; |
| 362 | x2 = (1.0 / x0) * (1.0 / x0); |
| 363 | /* log(2 * M_PI) */ |
| 364 | lg2pi = 1.8378770664093453e+00; |
| 365 | gl0 = a[9]; |
| 366 | for (k = 8; k >= 0; k--) { |
| 367 | gl0 *= x2; |
| 368 | gl0 += a[k]; |
| 369 | } |
| 370 | gl = gl0 / x0 + 0.5 * lg2pi + (x0 - 0.5) * log(x0) - x0; |
| 371 | if (x < 7.0) { |
| 372 | for (k = 1; k <= n; k++) { |
| 373 | gl -= log(x0 - 1.0); |
| 374 | x0 -= 1.0; |
| 375 | } |
| 376 | } |
| 377 | return gl; |
| 378 | } |
| 379 | |
| 380 | /* |
| 381 | double random_normal(bitgen_t *bitgen_state, double loc, double scale) { |
no test coverage detected