Skip to content

Commit ca104bf

Browse files
authored
Avoid computing log(p) twice
1 parent a02c080 commit ca104bf

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

@@ -952,7 +951,7 @@ function __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
952951
elseif abs(a - 1.0) < 1.0e-4
953952
x0 = pcase ? -log1p(-minpq) : -log(minpq)
954953
elseif a < 1.0 # small value of a
955-
x0 = gamma_inc_inv_asmall(a, minpq, pcase)
954+
x0 = gamma_inc_inv_asmall(a, logp)
956955
else #large a
957956
haseta = true
958957
x0, fp = gamma_inc_inv_alarge(a, minpq, pcase)

0 commit comments

Comments
 (0)