MAINT: special: fix expi after changing branch cut behavior of exp1

Also add tests for the branch cut.
jpreiss · Sep 8, 2019 · 7e0b254 · 7e0b254
1 parent 77d9f96
commit 7e0b254
Show file tree

Hide file tree

Showing 3 changed files with 43 additions and 6 deletions.
diff --git a/scipy/special/specfun/specfun.f b/scipy/special/specfun/specfun.f
@@ -5276,7 +5276,7 @@ SUBROUTINE EIXZ(Z,CEI)
 C       ============================================
 C
         IMPLICIT NONE
-        DOUBLE COMPLEX Z, CEI, IMF
+        DOUBLE COMPLEX Z, CEI
         DOUBLE PRECISION PI
         PI=3.141592653589793D0
         CALL E1Z(-Z, CEI)
@@ -5287,7 +5287,7 @@ SUBROUTINE EIXZ(Z,CEI)
            CEI = CEI - (0d0,1d0)*PI
         ELSE IF (DIMAG(Z).EQ.0) THEN
            IF (DBLE(Z).GT.0) THEN
-              CEI = CEI - (0d0,1d0)*PI
+              CEI = CEI + (0d0,1d0)*DSIGN(PI,DIMAG(Z))
            ENDIF
         ENDIF
         RETURN

diff --git a/scipy/special/tests/test_basic.py b/scipy/special/tests/test_basic.py
@@ -321,9 +321,6 @@ def test_erf_symmetry(self):
     def test_erfc(self):
         assert_equal(cephes.erfc(0), 1.0)
 
-    def test_expi(self):
-        cephes.expi(1)
-
     def test_expn(self):
         cephes.expn(1,1)
 

diff --git a/scipy/special/tests/test_exponential_integrals.py b/scipy/special/tests/test_exponential_integrals.py
@@ -1,10 +1,14 @@
+import pytest
+
+import numpy as np
 from numpy.testing import assert_allclose
 import scipy.special as sc
 
 
 class TestExp1(object):
 
     def test_branch_cut(self):
+        assert np.isnan(sc.exp1(-1))
         assert sc.exp1(complex(-1, 0)).imag == (
             -sc.exp1(complex(-1, -0.0)).imag
         )
@@ -22,8 +26,44 @@ def test_branch_cut(self):
             rtol=1e-15
         )
 
-    def test_exp1_834(self):
+    def test_834(self):
         # Regression test for #834
         a = sc.exp1(-complex(19.9999990))
         b = sc.exp1(-complex(19.9999991))
         assert_allclose(a.imag, b.imag, atol=0, rtol=1e-15)
+
+
+class TestExpi(object):
+
+    @pytest.mark.parametrize('result', [
+        sc.expi(complex(-1, 0)),
+        sc.expi(complex(-1, -0.0)),
+        sc.expi(-1)
+    ])
+    def test_branch_cut(self, result):
+        desired = -0.21938393439552027368  # Computed using Mpmath
+        assert_allclose(result, desired, atol=0, rtol=1e-14)
+
+    def test_near_branch_cut(self):
+        lim_from_above = sc.expi(-1 + 1e-20j)
+        lim_from_below = sc.expi(-1 - 1e-20j)
+        assert_allclose(
+            lim_from_above.real,
+            lim_from_below.real,
+            atol=0,
+            rtol=1e-15
+        )
+        assert_allclose(
+            lim_from_above.imag,
+            -lim_from_below.imag,
+            atol=0,
+            rtol=1e-15
+        )
+
+    def test_continuity_on_positive_real_axis(self):
+        assert_allclose(
+            sc.expi(complex(1, 0)),
+            sc.expi(complex(1, -0.0)),
+            atol=0,
+            rtol=1e-15
+        )