Merge pull request #394 from JuliaMath/an/backport

More bugfixes to backport
JuliaMath · Feb 25, 2022 · d6cabd8 · d6cabd8
2 parents 2ec1157 + f80715a
commit d6cabd8
Show file tree

Hide file tree

Showing 5 changed files with 120 additions and 21 deletions.
diff --git a/Project.toml b/Project.toml
@@ -13,7 +13,7 @@ OpenSpecFun_jll = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 ChainRulesCore = "0.9.44, 0.10, 1"
 ChainRulesTestUtils = "0.6.8, 0.7, 1"
 IrrationalConstants = "0.1"
-LogExpFunctions = "0.3"
+LogExpFunctions = "0.3.2"
 OpenLibm_jll = "0.7, 0.8"
 OpenSpecFun_jll = "0.5"
 julia = "1.3"

diff --git a/src/beta_inc.jl b/src/beta_inc.jl
@@ -371,19 +371,21 @@ function beta_inc_asymptotic_asymmetric(a::Float64, b::Float64, x::Float64, y::F
     end
     z = -nu*lnx
     if b*z == 0.0
-        return error("expansion can't be computed")
+        @debug("underflow: expansion can't be computed")
+        return w
     end
 
     # COMPUTATION OF THE EXPANSION
     #SET R = EXP(-Z)*Z**B/GAMMA(B)
     r = b*(1.0 + rgamma1pm1(b))*exp(b*log(z))
     r *= exp(a*lnx)*exp(0.5*bm1*lnx)
-    u = loggammadiv(b,a) + b*log(nu)
+    u = loggammadiv(b, a) + b*log(nu)
     u = r*exp(-u)
     if u == 0.0
-        return error("expansion can't be computed")
+        @debug("underflow: expansion can't be computed")
+        return w
     end
-    (p, q) = gamma_inc(b,z,0)
+    p, q = gamma_inc(b, z, 0)
     v = inv(nu)^2/4
     t2 = lnx^2/4
     l = w/u
@@ -412,7 +414,8 @@ function beta_inc_asymptotic_asymmetric(a::Float64, b::Float64, x::Float64, y::F
         dj = d[n] * j
         sm += dj
         if sm <= 0.0
-            return error("expansion can't be computed")
+            @debug("underflow: expansion can't be computed")
+            return w
         end
         if abs(dj) <= epps*(sm+l)
             break
@@ -837,7 +840,7 @@ function _beta_inc(a::Float64, b::Float64, x::Float64, y::Float64=1-x)
     elseif a0 < min(epps, epps*b0) && b0*x0 <= 1.0
         q = beta_inc_power_series1(a0, b0, x0, epps)
         p = 1.0 - q
-    elseif max(a0,b0) > 1.0
+    elseif max(a0, b0) > 1.0
         if b0 <= 1.0
             p = beta_inc_power_series(a0, b0, x0, epps)
             q = 1.0 - p
@@ -1001,7 +1004,10 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p)
     while true
         p_approx = beta_inc(a, b, x)[1]
         xin = x
-        p_approx = (p_approx - p)*min(floatmax(), exp(lb + r*log(xin) + t*log1p(-xin)))
+        p_approx = (p_approx - p)*min(
+            floatmax(),
+            exp(lb + LogExpFunctions.xlogy(r, xin) + LogExpFunctions.xlog1py(t, -xin))
+        )
 
         if p_approx * p_approx_prev <= 0.0
             prev = max(sq, fpu)

diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl
@@ -226,7 +226,15 @@ Refer Eqn (3.12) in the paper
 """
 function coeff1(eta::Float64)
     if abs(eta) < 1.0
-        coeff1 = @horner(eta, -3.333333333438e-1, -2.070740359969e-1, -5.041806657154e-2, -4.923635739372e-3, -4.293658292782e-5) / @horner(eta, 1.000000000000e+0, 7.045554412463e-1, 2.118190062224e-1, 3.048648397436e-2, 1.605037988091e-3)
+        coeff1 = @horner(
+            eta,
+            -3.333333333438e-1, -2.070740359969e-1, -5.041806657154e-2,
+            -4.923635739372e-3, -4.293658292782e-5
+        ) / @horner(
+            eta,
+            1.000000000000e+0, 7.045554412463e-1, 2.118190062224e-1,
+            3.048648397436e-2, 1.605037988091e-3
+        )
     else
         la = lambdaeta(eta)
         coeff1 = log(eta/(la - 1.0))/eta
@@ -243,16 +251,39 @@ Refer Eqn (3.12) in the paper
 function coeff2(eta::Float64)
 
     if eta < -5.0
-        x=eta*eta
+        x = eta*eta
         lnmeta = log(-eta)
         coeff2 = (12.0 - x - 6.0*lnmeta*lnmeta)/(12.0*x*eta)
     elseif eta < -2.0
-        coeff2 = @horner(eta, -1.72847633523e-2, -1.59372646475e-2, -4.64910887221e-3, -6.06834887760e-4, -6.14830384279e-6) / @horner(eta, 1.00000000000e+0, 7.64050615669e-1, 2.97143406325e-1, 5.79490176079e-2, 5.74558524851e-3)
+        coeff2 = @horner(
+            eta,
+            -1.72847633523e-2, -1.59372646475e-2, -4.64910887221e-3,
+            -6.06834887760e-4, -6.14830384279e-6
+        ) / @horner(
+            eta,
+            1.00000000000e+0, 7.64050615669e-1, 2.97143406325e-1,
+            5.79490176079e-2, 5.74558524851e-3)
     elseif eta < 2.0
-        coeff2 = @horner(eta, -1.72839517431e-2, -1.46362417966e-2, -3.57406772616e-3, -3.91032032692e-4, 2.49634036069e-6) / @horner(eta, 1.00000000000e+0, 6.90560400696e-1, 2.49962384741e-1, 4.43843438769e-2, 4.24073217211e-3)
+        coeff2 = @horner(
+            eta,
+            -1.72839517431e-2, -1.46362417966e-2, -3.57406772616e-3,
+            -3.91032032692e-4, 2.49634036069e-6
+        ) / @horner(
+            eta,
+            1.00000000000e+0, 6.90560400696e-1, 2.49962384741e-1,
+            4.43843438769e-2, 4.24073217211e-3
+        )
     elseif eta < 1000.0
         x = 1.0/eta
-        coeff2 = @horner(x, 9.99944669480e-1, 1.04649839762e+2, 8.57204033806e+2, 7.31901559577e+2, 4.55174411671e+1) / @horner(x, 1.00000000000e+0, 1.04526456943e+2, 8.23313447808e+2, 3.11993802124e+3, 3.97003311219e+3)
+        coeff2 = @horner(
+            x,
+            9.99944669480e-1, 1.04649839762e+2, 8.57204033806e+2,
+            7.31901559577e+2, 4.55174411671e+1
+        ) / @horner(
+            x,
+            1.00000000000e+0, 1.04526456943e+2, 8.23313447808e+2,
+            3.11993802124e+3, 3.97003311219e+3
+        )/(-12.0*eta)
     else
         coeff2 = -1.0/(12.0*eta)
     end
@@ -269,19 +300,65 @@ function coeff3(eta::Float64)
     if eta < -8.0
         x=eta*eta
         y = log(-eta)/eta
-        coeff3=(-30.0+eta*y*(6.0*x*y*y-12.0+x))/(12.0*eta*x*x)
+        coeff3=(-30.0 + eta*y*(6.0*x*y*y - 12.0 + x))/(12.0*eta*x*x)
     elseif eta < -4.0
-        coeff3 = (@horner(eta, 4.95346498136e-2, 2.99521337141e-2, 6.88296911516e-3, 5.12634846317e-4, -2.01411722031e-5) / @horner(eta, 1.00000000000e+0, 7.59803615283e-1, 2.61547111595e-1, 4.64854522477e-2, 4.03751193496e-3))/(eta*eta)
+        coeff3 = (
+            @horner(
+                eta,
+                4.95346498136e-2, 2.99521337141e-2, 6.88296911516e-3,
+                5.12634846317e-4, -2.01411722031e-5
+            ) / @horner(
+                eta,
+                1.00000000000e+0, 7.59803615283e-1, 2.61547111595e-1,
+                4.64854522477e-2, 4.03751193496e-3
+            )
+        )/(eta*eta)
     elseif eta < -2.0
-        coeff3 = @horner(eta, 4.52313583942e-3, 1.20744920113e-3, -7.89724156582e-5, -5.04476066942e-5, -5.35770949796e-6) / @horner(eta, 1.00000000000e+0, 9.12203410349e-1, 4.05368773071e-1, 9.01638932349e-2, 9.48935714996e-3)
+        coeff3 = @horner(
+            eta,
+            4.52313583942e-3, 1.20744920113e-3, -7.89724156582e-5,
+            -5.04476066942e-5, -5.35770949796e-6
+        ) / @horner(
+            eta,
+            1.00000000000e+0, 9.12203410349e-1, 4.05368773071e-1,
+            9.01638932349e-2, 9.48935714996e-3
+        )
     elseif eta < 2.0
-        coeff3 = @horner(eta, 4.39937562904e-3, 4.87225670639e-4, -1.28470657374e-4, 5.29110969589e-6, 1.57166771750e-7) / @horner(eta, 1.00000000000e+0, 7.94435257415e-1, 3.33094721709e-1, 7.03527806143e-2, 8.06110846078e-3)
+        coeff3 = @horner(
+            eta,
+            4.39937562904e-3, 4.87225670639e-4, -1.28470657374e-4,
+            5.29110969589e-6, 1.57166771750e-7
+        ) / @horner(
+            eta,
+            1.00000000000e+0, 7.94435257415e-1, 3.33094721709e-1,
+            7.03527806143e-2, 8.06110846078e-3
+        )
     elseif eta < 10.0
-        x=1.0/eta
-        coeff3 = (@horner(x, -1.14811912320e-3, -1.12850923276e-1, 1.51623048511e+0, -2.18472031183e-1, 7.30002451555e-2) / @horner(x, 1.00000000000e+0, 1.42482206905e+1, 6.97360396285e+1, 2.18938950816e+2, 2.77067027185e+2))/(eta*eta)
+        x = 1.0/eta
+        coeff3 = (
+            @horner(
+                x,
+                -1.14811912320e-3, -1.12850923276e-1, 1.51623048511e+0,
+                -2.18472031183e-1, 7.30002451555e-2
+            ) / @horner(
+                x,
+                1.00000000000e+0, 1.42482206905e+1, 6.97360396285e+1,
+                2.18938950816e+2, 2.77067027185e+2
+            )
+        )/(eta*eta)
     elseif eta < 100.0
-        x=1.0/eta
-        coeff3 = (@horner(x, -1.45727889667e-4, -2.90806748131e-1, -1.33085045450e+1, 1.99722374056e+2, -1.14311378756e+1) / @horner(x, 1.00000000000e+0, 1.39612587808e+2, 2.18901116348e+3, 7.11524019009e+3, 4.55746081453e+4))/(eta*eta)
+        x = 1.0/eta
+        coeff3 = (
+            @horner(
+                x,
+                -1.45727889667e-4, -2.90806748131e-1, -1.33085045450e+1,
+                1.99722374056e+2, -1.14311378756e+1
+            ) / @horner(
+                x,
+                1.00000000000e+0, 1.39612587808e+2, 2.18901116348e+3,
+                7.11524019009e+3, 4.55746081453e+4
+            )
+        )/(eta*eta)
     else
         eta3 = eta*eta*eta
         coeff3 = -log(eta)/(12.0*eta3)

diff --git a/test/beta_inc.jl b/test/beta_inc.jl
@@ -230,6 +230,11 @@
     @test SpecialFunctions.loggammadiv(13.89, 21.0001) ≈ log(gamma(big(21.0001))/gamma(big(21.0001)+big(13.89)))
     @test SpecialFunctions.loggammadiv(big(13.89), big(21.0001)) ≈ log(gamma(big(21.0001))/gamma(big(21.0001)+big(13.89)))
     @test SpecialFunctions.stirling_corr(11.99, 100.1) ≈ SpecialFunctions.stirling_error(11.99) + SpecialFunctions.stirling_error(100.1) - SpecialFunctions.stirling_error(11.99 + 100.1)
+
+    @testset "Issue 334. Underflow without erroring" begin
+        @test beta_inc(0.1, 4000, 0.2) == (1.0, 0.0)
+        @test beta_inc(4000, 0.1, 0.2) == (0.0, 1.0)
+    end
 end
 
 @testset "inverse of incomplete beta" begin
@@ -289,4 +294,9 @@ end
             @test beta_inc_inv(a, b, p, 1 - p) === beta_inc_inv(a, b, p)
         end
     end
+
+    @testset "Avoid NaN when p=q=1" begin
+        @test first(beta_inc_inv(1.0, 1.0, 1e-21)) ≈ 1e-21
+        @test beta_inc_inv(1.0e30, 1.0, 0.49) == (1.0, 0.0)
+    end
 end
diff --git a/test/gamma_inc.jl b/test/gamma_inc.jl
@@ -201,6 +201,12 @@ end
             end
         end
     end
+
+    @testset "Issue 390 part 1" begin
+        a = 1.0309015068677239
+        q = 0.020202020202020204
+        @test last(gamma_inc(a, gamma_inc_inv(a, 1 - q, q))) ≈ q
+    end
 end
 
 double(x::Real) = Float64(x)