Handle the case when Float64(p) + Float64(q) != 1 in gamma_inc_inv

for lower precision p and q.

Author: andreasnoack

src/gamma_inc.jl

@@ -1014,8 +1014,22 @@ function _gamma_inc_inv(a::Float64, p::Float64, q::Float64)
     return x
 end
 
-_gamma_inc_inv(a::Float32, p::Float32, q::Float32) = Float32(_gamma_inc_inv(Float64(a), Float64(p), Float64(q)))
-_gamma_inc_inv(a::Float16, p::Float16, q::Float16) = Float16(_gamma_inc_inv(Float64(a), Float64(p), Float64(q)))
+function _gamma_inc_inv(a::T, p::T, q::T) where {T <: Union{Float16, Float32}}
+    p64, q64 = if p < q
+        _q = 1 - Float64(p)
+        if q != T(_q)
+            throw(ArgumentError("p + q must equal one but was $(p + q)"))
+        end
+        Float64(p), _q
+    else
+        _p = 1 - Float64(q)
+        if p != T(_p)
+            throw(ArgumentError("p + q must equal one but was $(p + q)"))
+        end
+        _p, Float64(q)
+    end
+    return T(_gamma_inc_inv(Float64(a), p64, q64))
+end
 
 # like promote(x,y), but don't complexify real values
 promotereal(x::Real, y::Real) = promote(x, y)

test/gamma_inc.jl

@@ -164,7 +164,15 @@ end
     for x = .5:5.0:100.0
         @test SpecialFunctions.stirling_error(x) ≈ log(gamma(x)) - (x-.5)*log(x)+x- log(2*pi)/2.0
     end
+
     @test_throws ArgumentError("p + q must equal one but is 0.5") gamma_inc_inv(0.4, 0.2, 0.3)
+
+    @testset "Low precision with Float64(p) + Float64(q) != 1" for T in (Float16, Float32)
+        @test gamma_inc(T(1.0), gamma_inc_inv(T(1.0), T(0.1), T(0.9)))[1]::T ≈ T(0.1)
+        @test gamma_inc(T(1.0), gamma_inc_inv(T(1.0), T(0.9), T(0.1)))[2]::T ≈ T(0.1)
+        @test_throws ArgumentError("p + q must equal one but was 1.02") gamma_inc_inv(T(1.0), T(0.1), T(0.92))
+        @test_throws ArgumentError("p + q must equal one but was 1.02") gamma_inc_inv(T(1.0), T(0.92), T(0.1))
+    end
 end
 
 double(x::Real) = Float64(x)