Skip to content

Commit f94dd83

Browse files
refactor 2
1 parent 2a452e7 commit f94dd83

1 file changed

Lines changed: 9 additions & 8 deletions

File tree

src/PoissonRandom.jl

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -84,29 +84,30 @@ end
8484

8585
# Procedure F
8686
function procf(λ, K::Int, s::Float64)
87-
INV_SQRT_2PI = 0.3989422804014327 # 1/sqrt(2π)
87+
# can be pre-computed, but does not seem to affect performance
88+
INV_SQRT_2PI = inv(sqrt(2pi))
8889
ω = INV_SQRT_2PI / s
89-
b1 = 1 / (24 * λ)
90-
b2 = 0.3 * b1^2
91-
c3 = b1 * b2 / 7
90+
b1 = inv(24) / λ
91+
b2 = 0.3 * b1 * b1
92+
c3 = inv(7) * b1 * b2
9293
c2 = b2 - 15 * c3
9394
c1 = b1 - 6 * b2 + 45 * c3
9495
c0 = 1 - b1 + 3 * b2 - 15 * c3
9596

9697
if K < 10
9798
px = -λ
98-
log_py = K * log(λ) - loggamma(K + 1) # log(K!) via loggamma
99+
log_py = K * log(λ) - loggamma(K + 1) # log(K!) via loggamma
99100
py = exp(log_py)
100101
else
101-
δ = 1 / (12 * K)
102+
δ = inv(12) / K
102103
δ -= 4.8 * δ^3
103104
V =- K) / K
104-
px = K * log1pmx(V) - δ # avoids need for table
105+
px = K * log1pmx(V) - δ # avoids need for table
105106
py = INV_SQRT_2PI / sqrt(K)
106107
end
107108
X = (K - λ + 0.5) / s
108109
X2 = X^2
109-
fx = -X2 / 2 # missing negation in pseudo-algorithm, but appears in fortran code.
110+
fx = X2 / -2 # missing negation in pseudo-algorithm, but appears in fortran code.
110111
fy = ω * (((c3 * X2 + c2) * X2 + c1) * X2 + c0)
111112
return px, py, fx, fy
112113
end

0 commit comments

Comments
 (0)