bug fixes in expintx (#285)

* bug fixes in expintx

* fix spurious underflow in expintx for real(z) >> 1

Author: stevengj

src/expint.jl

@@ -169,26 +169,19 @@ function En_cf_nogamma(ν::Number, z::Number, n::Int=1000)
     scale = sqrt(floatmax(real(A)))
 
     # two recurrence steps / loop
-    iters = 0
+    iters = 1
     for i = 2:n
         iters += 1
 
-        A′ = A
-        A = z*A + (i-1) * Aprev
-        Aprev = A′
-        B′ = B
-        B = z*B + (i-1) * Bprev
-        Bprev = B′
+        A, Aprev = z*A + (i-1) * Aprev, A
+        B, Bprev = z*B + (i-1) * Bprev, B
 
-        A′ = A
-        A = A + (ν+i-1) * Aprev
-        Aprev = A′
-        B′ = B
-        B = B + (ν+i-1) * Bprev
-        Bprev = B′
+        (isinf(A) | isinf(B)) && return Aprev/Bprev, iters-1
 
-        conv = abs(Aprev*B - A*Bprev) < ϵ*abs(B*Bprev)
-        conv && i > 4 && break
+        A, Aprev = A + (ν+i-1) * Aprev, A
+        B, Bprev = B + (ν+i-1) * Bprev, B
+
+        i > 4 && abs(Aprev*B - A*Bprev) < ϵ*abs(B*Bprev) && break
 
         # rescale
         if fastabs(A) > scale
@@ -221,18 +214,14 @@ function En_cf_gamma(ν::Number, z::Number, n::Int=1000)
     ϵ = 10*eps(real(B))
     scale = sqrt(floatmax(real(A)))
 
-    iters = 0
+    iters = 1
     j = 1
     for i = 2:n
         iters += 1
 
-        A′ = A
         term = iseven(i) ? -(i÷2 - 1 - ν)*z : z*(i÷2)
-        A = (i - ν)*A + term * Aprev
-        Aprev = A′
-        B′ = B
-        B = (i - ν)*B + term * Bprev
-        Bprev = B′
+        A, Aprev = (i - ν)*A + term * Aprev, A
+        B, Bprev = (i - ν)*B + term * Bprev, B
 
         conv = abs(Aprev*B - A*Bprev) < ϵ*abs(B*Bprev)
         conv && break
@@ -458,21 +447,25 @@ function _expint(ν::Number, z::Number, niter::Int=1000, ::Val{expscaled}=Val{fa
         #   16(2), 169–177. https://doi.org/10.1145/78928.78933
         doconj = imag(z) < 0
         rez, imz = real(z), abs(imag(z))
+        zorig = z
         z = doconj ? conj(z) : z
         ν = doconj ? conj(ν) : ν
 
-        quick_niter, nmax = 50, 45
+        quick_niter = niter >> 4
         # start with small imaginary part if exactly on negative real axis
         imstart = (imz == 0) ? abs(z)*sqrt(eps(typeof(real(z)))) : imz
         z₀ = rez + imstart*im
         g_start, cf_start, i, _ = En_cf(ν, z₀, quick_niter)
         E_start = g_start + En_safeexpmult(-z₀, cf_start)
         isnan(E_start) && return E_start
-        if imz > 0 && i < nmax
+        if imz > 0 && i < quick_niter
             # didn't need to take any steps
+            if expscaled
+                E_start = En_safeexpmult(z₀, g_start) + cf_start
+            end
             return doconj ? conj(E_start) : E_start
         end
-        while i > nmax
+        while i == quick_niter
             # double imaginary part until in region with fast convergence
             imstart *= 2
             z₀ = rez + imstart*im
@@ -504,7 +497,7 @@ function _expint(ν::Number, z::Number, niter::Int=1000, ::Val{expscaled}=Val{fa
                 En = bit ? En : En + bc
             end
         end
-        return expscaled ? exp(z) * En : En
+        return expscaled ? En_safeexpmult(zorig, En) : En
     end
 end
 

test/expint.jl

@@ -93,7 +93,7 @@ using Base.MathConstants
         for (i, ans) in zip(1:15, near_pole5_ans)
             @test expint(5 - 10.0^(-i), 2.5) ⩭ ans
         end
-        
+
         @test expint(2.000001, 4.3) ⩭ 0.0022470372637239675640390482141893761501016099509821867746
         @test expint(1.999999, 4.3) ⩭ 0.0022470379319936678143503120160911823783417969321242814713
 
@@ -141,11 +141,17 @@ using Base.MathConstants
     @testset "expintx" begin
         @test expintx(1e-10) ≅ 22.448635267383787506197038031047821505309344119883144316155006107
         @test expintx(1e-100) ≅ 229.68129363450303554119263337835401832906798952693737400452702286
-        @test expintx(1e10) ≅ 9.999999999000000000199999999940000000023999999988000000007e-11 
+        @test expintx(1e10) ≅ 9.999999999000000000199999999940000000023999999988000000007e-11
         @test_broken expintx(1e10, 1e10im) ≅ 4.99999999999999999995000000001499999999700000000000000000e-11 - 4.99999999950000000004999999999999999999700000000149999999e-11im
         @test expintx(2, 1e10) ≅ 9.999999998000000000599999999760000000119999999928000000050e-11
         @test expintx(5+5im, 1e10+1e10im) ≅ 4.99999999750000000124999999940000000023624999998587499986e-11 - 4.99999999750000000149999999897500000075624999942337500043e-11im
         @test expintx(2, 1e5+1e5im) ≅ 4.999999998500059998500000003149748011339999937637483913514e-6 - 4.999900001499999998500089996850000011338866062369999513582e-6im
+
+        for ν in (1,2,-2,2+3im,2-3im,-2+3im,-2-3im), z in (0.1+0.2im, -0.1+0.2im, 4+5im, -4+5im, -4-5im, -10+1e-4im, -10-1e-4im)
+            @test expintx(ν, z) ≅ exp(z) * expint(ν, z)
+        end
+
+        @test expintx(4, 1e300) == 1e-300
     end
 
     # issue #272: more cases with real results