add log-gamma function
@@ -179,6 +179,39 @@ __device__ double rk_beta(rk_state* state, double a, double b) {
}
}
+__device__ double loggam(double x) {
+ double x0, x2, xp, gl, gl0;
+ long k, n;
+ double a[10] = {8.333333333333333e-02,-2.777777777777778e-03,
+ 7.936507936507937e-04,-5.952380952380952e-04,
+ 8.417508417508418e-04,-1.917526917526918e-03,
+ 6.410256410256410e-03,-2.955065359477124e-02,
+ 1.796443723688307e-01,-1.39243221690590e+00};
+ x0 = x;
+ n = 0;
+ if ((x == 1.0) || (x == 2.0)) {
+ return 0.0;
+ } else if (x <= 7.0) {
+ n = (long)(7 - x);
+ x0 = x + n;
+ }
+ x2 = 1.0/(x0*x0);
+ xp = 2*M_PI;
+ gl0 = a[9];
+ for (k=8; k>=0; k--) {
+ gl0 *= x2;
+ gl0 += a[k];
+ }
+ gl = gl0/x0 + 0.5*log(xp) + (x0-0.5)*log(x0) - x0;
+ if (x <= 7.0) {
+ for (k=1; k<=n; k++) {
+ gl -= log(x0-1.0);
+ x0 -= 1.0;
+ }
+ }
+ return gl;
+}
+
__device__ int64_t rk_poisson_mult(rk_state *state, double lam) {
int64_t X;
double prod, U, enlam;
@@ -206,10 +239,10 @@ __device__ int64_t rk_poisson_ptrs(rk_state *state, double lam) {
invalpha = 1.1239 + 1.1328/(b-3.4);
vr = 0.9277 - 3.6224/(b-2);
while (1) {
- U = state->rk_double(); - 0.5;
+ U = state->rk_double() - 0.5;
V = state->rk_double();
us = 0.5 - fabs(U);
- k = (long)floor((2*a/us + b)*U + lam + 0.43);
+ k = (int64_t)floor((2*a/us + b)*U + lam + 0.43);
if ((us >= 0.07) && (V <= vr)) {
return k;
}
