Skip to content

Commit b16f79f

Browse files
authored
Handle the case when Float64(p) + Float64(q) != 1 in gamma_inc_inv (#374)
* Handle the case when Float64(p) + Float64(q) != 1 in gamma_inc_inv for lower precision p and q. * Use Speve's version
1 parent f2ccad8 commit b16f79f

File tree

2 files changed

+19
-2
lines changed

2 files changed

+19
-2
lines changed

src/gamma_inc.jl

Lines changed: 11 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1014,8 +1014,17 @@ function _gamma_inc_inv(a::Float64, p::Float64, q::Float64)
10141014
return x
10151015
end
10161016

1017-
_gamma_inc_inv(a::Float32, p::Float32, q::Float32) = Float32(_gamma_inc_inv(Float64(a), Float64(p), Float64(q)))
1018-
_gamma_inc_inv(a::Float16, p::Float16, q::Float16) = Float16(_gamma_inc_inv(Float64(a), Float64(p), Float64(q)))
1017+
function _gamma_inc_inv(a::T, p::T, q::T) where {T <: Union{Float16, Float32}}
1018+
if p + q != one(T)
1019+
throw(ArgumentError("p + q must equal one but was $(p + q)"))
1020+
end
1021+
p64, q64 = if p < q
1022+
(Float64(p), 1 - Float64(p))
1023+
else
1024+
(1 - Float64(q), Float64(q))
1025+
end
1026+
return T(_gamma_inc_inv(Float64(a), p64, q64))
1027+
end
10191028

10201029
# like promote(x,y), but don't complexify real values
10211030
promotereal(x::Real, y::Real) = promote(x, y)

test/gamma_inc.jl

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -164,7 +164,15 @@ end
164164
for x = .5:5.0:100.0
165165
@test SpecialFunctions.stirling_error(x) log(gamma(x)) - (x-.5)*log(x)+x- log(2*pi)/2.0
166166
end
167+
167168
@test_throws ArgumentError("p + q must equal one but is 0.5") gamma_inc_inv(0.4, 0.2, 0.3)
169+
170+
@testset "Low precision with Float64(p) + Float64(q) != 1" for T in (Float16, Float32)
171+
@test gamma_inc(T(1.0), gamma_inc_inv(T(1.0), T(0.1), T(0.9)))[1]::T T(0.1)
172+
@test gamma_inc(T(1.0), gamma_inc_inv(T(1.0), T(0.9), T(0.1)))[2]::T T(0.1)
173+
@test_throws ArgumentError("p + q must equal one but was 1.02") gamma_inc_inv(T(1.0), T(0.1), T(0.92))
174+
@test_throws ArgumentError("p + q must equal one but was 1.02") gamma_inc_inv(T(1.0), T(0.92), T(0.1))
175+
end
168176
end
169177

170178
double(x::Real) = Float64(x)

0 commit comments

Comments
 (0)