How to work around this error in expint?

Specific Domains Numerics
1

I ran into this error while evaluating the Gamma function from SpecialFunctions:

julia> using SpecialFunctions

julia> gamma(exp(19.39)+1, 0.0125*0.6)
ERROR: ArgumentError: Unsupported order |ν| > 50 off the positive real axis
Stacktrace:
 [1] _expint(ν::Float64, z::Float64, niter::Int64, ::Val{false})
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/9pXme/src/expint.jl:400
 [2] expint (repeats 2 times)
   @ ~/.julia/packages/SpecialFunctions/9pXme/src/expint.jl:516 [inlined]
 [3] _gamma(a::Float64, x::Float64)
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/9pXme/src/gamma_inc.jl:1059
 [4] gamma(a::Float64, x::Float64)
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/9pXme/src/gamma_inc.jl:1036
 [5] top-level scope
   @ REPL[43]:1

I traced it to this evaluation of expint

julia> expint(1-(exp(19.39)+1), 0.0125*0.6)
ERROR: ArgumentError: Unsupported order |ν| > 50 off the positive real axis
Stacktrace:
 [1] _expint(ν::Float64, z::Float64, niter::Int64, ::Val{false})
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/9pXme/src/expint.jl:400
 [2] expint (repeats 2 times)
   @ ~/.julia/packages/SpecialFunctions/9pXme/src/expint.jl:516 [inlined]
 [3] top-level scope
   @ REPL[41]:1

which in turn arises from this definition of _expint:

function _expint(ν::Number, z::Number, niter::Int=1000, ::Val{expscaled}=Val{false}()) where {expscaled}
    if abs(ν) > 50 && !(isreal(ν) && real(ν) > 0)
        throw(ArgumentError("Unsupported order |ν| > 50 off the positive real axis"))
    end
# rest of code here
end

How can I work around this?

2

According to the documentation gamma(a, x) computes the upper incomplete gamma function.
SpecialFunctions also provides gamma_inc(a, x) which computes the ratios of the lower and upper incomplete gamma(a, x) functions divided by gamma(a) and returns a tuple of both values.
Thus, you can compute your function as

mygamma(a, x) = gamma_inc(a, x)[2] * gamma(a)

For your arguments the result is just the same as gamma(exp(19.39) + 1) and way larger than floatmax(Float64). You get a meaningful result on the log scale though:

logmygamma(a, x) = log(gamma_inc(a, x)[2]) + loggamma(a)

julia> logmygamma(exp(19.39)+1, 0.0125*0.6)
4.847878583807443e9

2 Likes