Gamma PDF throwing NaN for very small scale parameter

[wzhorton]

Not sure if this is better suited for Distributions.jl, but I think it fits here better. Here is a reproducible example:

using Distributions
pdf(Gamma(3,1e-310),1)

returns NaN when it should either return 0.0 or throw the condition θ > zero(θ) is not satisfied error.

To add a little more detail, pdf(Gamma(3,1e-308),1) correctly returns 0.0, while pdf(Gamma(3,1e-309),1) returns NaN as do even smaller scale parameters until you get to pdf(Gamma(3,1e-324),1) when the zero-argument error is thrown. I think the source of the problem is in the gammalogpdf function from StatsFuns.jl, which is included below for reference:

gammalogpdf(k::Real, θ::Real, x::Real) = gammalogpdf(promote(k, θ, x)...)
function gammalogpdf(k::T, θ::T, x::T) where {T<:Real}
    # we ensure that `log(x)` does not error if `x < 0`
    xθ = max(x, 0) / θ
    val = -loggamma(k) + xlogy(k - 1, xθ) - log(θ) - xθ
    return x < 0 ? oftype(val, -Inf) : val
end

I think the source of the problem is in computing xθ = max(x, 0) / θ. For values of θ <= 1e-309, this quantity is Inf, but Julia doesn't acknowledge θ approx 0 until θ <= 1e-324, so we have this window of admittedly absurdly small values where the the parameter is valid (not zero) but the calculation effectively overflows to Inf through division. The result is that xlogy(k - 1, xθ) and xθ both evaluate to infinity and taking their difference leads to returning NaN.

Actually, I think I found the exact values: eps(0.0) returns 5.0e-324 while floatmin() returns around 2.2e-308.
As a possible fix, xθ = max(x, floatmin()) / θ seems to work instead of xθ = max(x, 0) / θ. Otherwise it would be nice if the interface error would show up at floatmin() instead of eps(0.0).

Specs: Julia v1.9.3, Distributions v0.25.102, StatsFuns v1.3.0, Apple Silicon M2