fix spurious domainerror for expint(3.2, 1.3) (#280)

Author: stevengj

src/expint.jl

@@ -196,7 +196,11 @@ function En_cf_nogamma(ν::Number, z::Number, n::Int=1000)
 end
 
 # Calculate Γ(1 - ν) * z^(ν-1) safely
-En_safe_gamma_term(ν::Number, z::Number) = exp((ν - 1)*log(z) + loggamma(1 - oftype(z, ν)))
+function En_safe_gamma_term(ν::Number, z::Number)
+    ν1 = 1 - oftype(z, ν)
+    lgamma, lgammasign = ν1 isa Real ? logabsgamma(ν1) : (loggamma(ν1), 1)
+    return lgammasign * exp((ν - 1)*log(z) + lgamma)
+end
 En_safe_gamma_term(ν::Integer, z::Real) = (z ≥ 0 || isodd(ν) ? 1 : -1) * exp((ν - 1)*log(abs(z)) + loggamma(1 - oftype(z, ν)))
 
 # continued fraction for En(ν, z) that uses the gamma function:

test/expint.jl

@@ -128,7 +128,8 @@ using Base.MathConstants
     @test expint(-2, 2.2) ≈ 0.11696351427428933590118175681379776 rtol=1e-14
     @test expint(-2, -2.2) ≈ -2.0680909972407264331884989 rtol=1e-14
     @test expint(-2.2, 3.2) ≈ 0.024950173497409329191241353395358 rtol=1e-14
-    #@test expint(+2.2, 3.2) ≈ 0.008044603700773423319087602010 rtol=1e-14
+    @test expint(+2.2, 3.2) ≈ 0.008044603700773423319087602010 rtol=1e-12
+    @test expint(3.2, 1.3) ≈ 0.070147692224611216675759479422283315452559216337905 rtol=1e-14
 end
 
 expinti_real(x) = invoke(expinti, Tuple{Real}, x)