expint: fixes for #263 (#268)

* Fix branch cut * use logabsgamma * Remove local var * complex nu branch cut test
JuliaMath · Nov 17, 2020 · dde82c2 · dde82c2
1 parent 9116108
commit dde82c2
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 11 deletions.
diff --git a/src/expint.jl b/src/expint.jl
@@ -344,11 +344,10 @@ end
 # https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/04/05/01/0003/
 function En_imagbranchcut(ν::Number, z::Number)
     a = real(z)
-    e1 = exp(π*imag(ν))
+    e1 = exp(oftype(a, π) * imag(ν))
     e2 = Complex(cospi(real(ν)), -sinpi(real(ν)))
-    impart = π * im * e1 * e2 * exp((ν-1)*log(complex(a)) - loggamma(ν))
-    impart *= signbit(imag(z)) ? -1 : 1
-    return imag(impart) * im # get rid of any real error
+    lgamma, lgammasign = ν isa Real ? logabsgamma(ν) : (loggamma(ν), 1)
+    return -2 * lgammasign * e1 * π * e2 * exp((ν-1)*log(complex(a)) - lgamma) * im
 end
 
 """
@@ -375,7 +374,8 @@ function expint(ν::Number, z::Number, niter::Int=1000)
     end
 
     if z == 0
-        return oftype(z, real(ν) > 0 ? 1/(ν-1) : Inf)
+        typ = typeof(real(z))
+        return oftype(z, real(ν) > 0 ? one(typ)/(ν-1) : convert(typ, Inf))
     end
 
     if ν == 0
@@ -394,8 +394,8 @@ function expint(ν::Number, z::Number, niter::Int=1000)
         return En_expand_origin(ν, z)
     end
 
-    if z isa Real
-        return first(z > 0 ? En_cf(ν, z, niter) : En_cf_nogamma(ν, z, niter))
+    if z isa Real || real(z) > 0
+        return first(real(z) > 0 ? En_cf(ν, z, niter) : En_cf_nogamma(ν, z, niter))
     else
         # Procedure for z near the negative real axis based on
         # (without looking at the accompanying source code):
@@ -434,12 +434,22 @@ function expint(ν::Number, z::Number, niter::Int=1000)
             z₀ += Δ
         end
 
-        # more exact imaginary part available for non-integer ν
+        En = doconj ? conj(E_start) : E_start
+
+        # handle branch cut
         if imz == 0
-            E_start = real(E_start) + En_imagbranchcut(ν, z)
+            bc = En_imagbranchcut(ν, z)
+            bit = !signbit(imag(z))
+            sign = bit ? 1 : -1
+            if isreal(ν)
+                # can separate real/im in case of real ν
+                return real(En) - sign * imag(bc)/2 * im
+            else
+                return bit ? En : En + bc
+            end
+        else
+            return En
         end
-
-        return doconj ? conj(E_start) : E_start
     end
     throw("unreachable")
 end

diff --git a/test/expint.jl b/test/expint.jl
@@ -76,6 +76,12 @@ using Base.MathConstants
         @test expint(-2, -3.3) ≅ -4.745485121488663825992363318400378114254141655964565963
         @test expint(-2, -3.3) isa Real
 
+        @test expint(-2.1+3im, 4.2 +0im) ⩭ 0.003387565626497080536855067510744646235748203589175196673362104 - 0.003355838600814825542988272289126728631565803782052141768321840im
+        @test expint(-4.3, -10+0.0im) ⩭ -1505.2167858114105916067190408869152177781925766866637213104 + 0.00015439438628349460776476190593323691084185406009267711345841im
+        @test expint(-4.3, -10-0.0im) ⩭ -1505.2167858114105916067190408869152177781925766866637213104 - 0.00015439438628349460776476190593323691084185406009267711345841im
+        @test expint(-4.3+2im, -10+0.0im) ⩭ -1473.52621376586428327992147406861600035297263354271049988247182 - 210.712752794950368400059514411922967744944073260391136148642008im
+        @test expint(-4.3+2im, -10-0.0im) ⩭ -1473.47816311715810169554659685200228583538709701122499753349598 - 210.761182017526309473296237005277228612276963738826032570499367im
+
         @test isnan(expint(NaN))
         @test isnan(expint(NaN+NaN*im))
         @test isnan(expint(NaN, 1.0))