Fix gamma_inc_inv for some subnormal results (#404) (#406)

Author: devmotion

Project.toml

@@ -1,6 +1,6 @@
 name = "SpecialFunctions"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "1.8.6"
+version = "1.8.7"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/gamma_inc.jl

@@ -864,16 +864,13 @@ function _gamma_inc(a::Float64, x::Float64, ind::Integer)
         throw(DomainError((a, x, ind), "`a` and `x` must be greater than 0 ---- Domain : (0, Inf)"))
     elseif a == 0.0 && x == 0.0
         throw(DomainError((a, x, ind), "`a` and `x` must be greater than 0 ---- Domain : (0, Inf)"))
-    elseif a*x == 0.0
-        if x <= a
-            return (0.0, 1.0)
-        else
-            return (1.0, 0.0)
-        end
     elseif isnan(a) || isnan(x)
-        return (a*x, a*x)
-    elseif isinf(x)
+        ax = a*x
+        return (ax, ax)
+    elseif a == 0.0 || isinf(x)
         return (1.0, 0.0)
+    elseif x == 0.0
+        return (0.0, 1.0)
     end
 
     if a >= 1.0
@@ -1042,8 +1039,6 @@ function __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
     logabsgam = logabsgamma(a)[1]
     # Newton-like higher order iteration with upper limit as 15.
     while t > 1.0e-15 && n < 15
-        x = x0
-        x² = x*x
         if !haseta
             dlnr = (1.0 - a)*log(x) + x + logabsgam
             if dlnr > log(floatmax(Float64)/1000.0)
@@ -1068,22 +1063,23 @@ function __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
             end
 
             if a > 0.1
-                ck3 = (@horner(x, @horner(a, 1, -3, 2), @horner(a, 4, -4), 2))/(6*x²)
+                ck3 = (@horner(x, @horner(a, 1, -3, 2), @horner(a, 4, -4), 2))/(6*x^2)
 
                 # This check is not in the invincgam subroutine from IncgamFI but it probably
                 # should be since the routine fails to compute a finite value for very small p
                 if !isfinite(ck3)
                     break
                 end
-                x0 = @horner(ck1, x, 1.0, ck2, ck3)
+                Δx = ck1 * @horner(ck1, 1.0, ck2, ck3)
             else
-                x0 = @horner(ck1, x, 1.0, ck2)
+                Δx = ck1 * @horner(ck1, 1.0, ck2)
             end
         else
-            x0 = x + ck1
+            Δx = ck1
         end
 
-        t = abs(x/x0 - 1.0)
+        x0 = x + Δx
+        t = abs(Δx/x0)
         n += 1
         x = x0
     end

test/gamma_inc.jl

@@ -202,10 +202,22 @@ end
         end
     end
 
-    @testset "Issue 390 part 1" begin
+    @testset "Issue 390 + 403" begin
         a = 1.0309015068677239
         q = 0.020202020202020204
         @test last(gamma_inc(a, gamma_inc_inv(a, 1 - q, q))) ≈ q
+
+        a = 0.0016546512046778552
+        q = 0.7070707070707071
+        # Mathematica: InverseGammaRegularized[0.0016546512046778552, 0, 1-0.7070707070707071] = 3.04992226601142476643093`9.786272979013901*^-323
+        @test gamma_inc_inv(a, 1 - q, q) ≈ 3e-323
+    end
+
+    @testset "Distributions.jl: Issue 1567" begin
+        a = 0.0030345129757232197
+        p = 0.106
+        # Mathematica: InverseGammaRegularized[0.0030345129757232197, 0, 0.106] = 3.5283979699566549210055643`10.308610719322063*^-322
+        @test gamma_inc_inv(a, p, 1 - p) ≈ 3.5e-322
     end
 end