Merge pull request #393 from JuliaMath/an/usexlogyinbetaincinv

Avoid NaN in beta_inc_inv when p=q=1.

Author: andreasnoack

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"

src/beta_inc.jl

@@ -1004,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)

test/beta_inc.jl

@@ -294,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