Infinite loop for quantile of InverseGaussian for specific value

[mentics]

I found this with Distributions v0.25.104, Julia 1.10.0, two machines: Ubuntu and Windows.

The following code will hang forever (ctrl-c will interrupt it):

using Distributions
dig = InverseGaussian{Float64}(μ=17.84806245738152, λ=163.707062977564)
quantile(dig, 0.9999981908772995)

Breaking it shows a stack trace of:

 [1] exp(x::Float64)
   @ Base.Math .\special\exp.jl:327
 [2] cdf(d::InverseGaussian{Float64}, x::Float64)
   @ Distributions C:\Users\joel\.julia\packages\Distributions\km0AH\src\univariate\continuous\inversegaussian.jl:102       
 [3] quantile_newton(d::InverseGaussian{Float64}, p::Float64, xs::Float64, tol::Float64)
   @ Distributions C:\Users\joel\.julia\packages\Distributions\km0AH\src\quantilealgs.jl:57
 [4] quantile_newton(d::InverseGaussian{Float64}, p::Float64)
   @ Distributions C:\Users\joel\.julia\packages\Distributions\km0AH\src\quantilealgs.jl:51
 [5] quantile(d::InverseGaussian{Float64}, p::Float64)
   @ Distributions C:\Users\joel\.julia\packages\Distributions\km0AH\src\quantilealgs.jl:146
 [6] top-level scope
   @ REPL[64]:1

It seems to me that the likely culprit is in the loop around Distributions:quantilealgs.jl:57.

I tested it against values just above and below it, only changing the final digit: quantile(dig, 0.9999981908772994) quantile(dig, 0.9999981908772996) and they returned values normally, so it appears I managed to get surprisingly "lucky" in finding this particular value that has a problem (though perhaps there is a family of them that would?).

[devmotion]

Possibly a duplicate of https://github.com/JuliaStats/Distributions.jl/issues/1571?

[mentics]

Yes, it does appear to be a duplicate. The examples in that issue break in the same loop. I apologize for not searching sufficiently in advance.

[devmotion]

No worries 🙂