Skip to content

Commit a9d8ae6

Browse files
devmotionandreasnoack
authored andcommitted
Avoid computing log(p) twice
(cherry picked from commit ca104bf)
1 parent 1f026ed commit a9d8ae6

File tree

1 file changed

+4
-5
lines changed

1 file changed

+4
-5
lines changed

src/gamma_inc.jl

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -725,14 +725,13 @@ function gamma_inc_inv_qsmall(a::Float64, q::Float64)
725725
end
726726

727727
"""
728-
gamma_inc_inv_asmall(a, minpq, pcase)
728+
gamma_inc_inv_asmall(a, logp)
729729
730730
Compute `x0` - initial approximation when `a` is small.
731-
Here the solution `x` of ``P(a,x)=p`` satisfies ``x_{l} < x < x_{u}``
731+
Here the solution `x` of ``P(a,x)=p=\\exp(logp)`` satisfies ``x_{l} < x < x_{u}``
732732
where ``x_{l} = (p\\Gamma(a+1))^{1/a}`` and ``x_{u} = -\\log{(1 - p\\Gamma(a+1))}``, and is used as starting value for Newton iteration.
733733
"""
734-
function gamma_inc_inv_asmall(a::Float64, minpq::Float64, pcase::Bool)
735-
logp = pcase ? log(minpq) : log1p(-minpq)
734+
function gamma_inc_inv_asmall(a::Float64, logp::Float64)
736735
return exp((logp + loggamma1p(a)) / a)
737736
end
738737

@@ -956,7 +955,7 @@ function __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
956955
elseif abs(a - 1.0) < 1.0e-4
957956
x0 = pcase ? -log1p(-minpq) : -log(minpq)
958957
elseif a < 1.0 # small value of a
959-
x0 = gamma_inc_inv_asmall(a, minpq, pcase)
958+
x0 = gamma_inc_inv_asmall(a, logp)
960959
else #large a
961960
haseta = true
962961
x0, fp = gamma_inc_inv_alarge(a, minpq, pcase)

0 commit comments

Comments
 (0)