fix expint for negative order (#260)

Author: stevengj

src/expint.jl

@@ -286,13 +286,9 @@ end
 
 # series about origin, general ν
 # https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/06/01/04/01/01/0003/
-function En_expand_origin(ν::Number, z::Number)
-    if isinteger(ν)
-        # go to special case for integer ν
-        return En_expand_origin(Int(ν), z)
-    end
+function En_expand_origin_general(ν::Number, z::Number)
     gammaterm = En_safe_gamma_term(ν, z)
-    frac = 1
+    frac = one(real(z))
     sumterm = frac / (1 - ν)
     k, maxiter = 1, 100
     ϵ = 10*eps(real(sumterm))
@@ -305,21 +301,20 @@ function En_expand_origin(ν::Number, z::Number)
         end
         k += 1
     end
-
     return gammaterm - sumterm
 end
 
-# series about the origin, special case for integer n
+# series about the origin, special case for integer n > 0
 # https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/06/01/04/01/02/0005/
-function En_expand_origin(n::Integer, z::Number)
+function En_expand_origin_posint(n::Integer, z::Number)
     gammaterm = 1 # (-z)^(n-1) / (n-1)!
     for i = 1:n-1
         gammaterm *= -z / i
     end
 
-    gammaterm *= digamma(n) - log(z)
-    sumterm = float(n == 1 ? 0 : 1 / (1 - n))
-    frac = 1
+    frac = one(real(z))
+    gammaterm *= digamma(oftype(frac,n)) - log(z)
+    sumterm = n == 1 ? zero(frac) : frac / (1 - n)
     k, maxiter = 1, 100
     ϵ = 10*eps(real(sumterm))
     while k < maxiter
@@ -337,6 +332,14 @@ function En_expand_origin(n::Integer, z::Number)
     return gammaterm - sumterm
 end
 
+function En_expand_origin(ν::Number, z::Number)
+    if isinteger(ν) && real(ν) > 0
+        return En_expand_origin_posint(Integer(ν), z)
+    else
+        return En_expand_origin_general(ν, z)
+    end
+end
+
 # can find imaginary part of E_ν(x) for x on negative real axis analytically
 # https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/04/05/01/0003/
 function En_imagbranchcut(ν::Number, z::Number)

test/expint.jl

@@ -98,6 +98,18 @@ using Base.MathConstants
             abs(ν) ≤ 50 ? (@test expint(ν, z) ⩭ E) : (@test_throws ArgumentError expint(ν, z) ⩭ E)
         end
     end
+
+    # issue #259: positive vs. negative integer order for small/larger args
+    @test expint(-4, 0.7) ≈ 142.685471227614252544 rtol=1e-14
+    @test expint(+4, 0.7) ≈ 0.1267808300929215658  rtol=1e-14
+    @test expint(-4, 5.7) ≈ 0.0013051784816041558 rtol=1e-14
+    @test expint(+4, 5.7) ≈ 0.0003585161998373381 rtol=1e-13
+    setprecision(BigFloat, 256) do
+        @test expint(-4, big"0.7") ≈ big"142.68547122761425254403329916322251663998829971602071820534499366398779205282" rtol=1e-76
+        @test expint(+4, big"0.7") ≈ big"0.12678083009292156580635697928870772316601035828975150623334903100558297079388567" rtol=1e-76
+        @test expint(-4, big"5.7") ≈ big"0.0013051784816041557674219259305064008704427556690521012788554962045721241678128274" rtol=1e-76
+        @test expint(+4, big"5.7") ≈ big"0.00035851619983733810063957822090566023258844652925986297947736355731994111784745604432" rtol=1e-74
+    end
 end
 
 expinti_real(x) = invoke(expinti, Tuple{Real}, x)