Merge pull request #4415 from juliantaylor/log1p-expm1-no-c99

BUG: fix non-c99 fallback for np.inf input to log1p/expm1

2 files changed, 39 insertions, 13 deletions

diff --git a/numpy/core/src/npymath/npy_math.c.src b/numpy/core/src/npymath/npy_math.c.src
index 61f1d79ab..1ca7033af 100644
--- a/numpy/core/src/npymath/npy_math.c.src
+++ b/numpy/core/src/npymath/npy_math.c.src
@@ -65,14 +65,19 @@
 #ifndef HAVE_EXPM1
 double npy_expm1(double x)
 {
-    const double u = npy_exp(x);
-
-    if (u == 1.0) {
+    if (npy_isinf(x) && x > 0) {
         return x;
-    } else if (u - 1.0 == -1.0) {
-        return -1;
-    } else {
-        return (u - 1.0) * x/npy_log(u);
+    }
+    else {
+        const double u = npy_exp(x);
+
+        if (u == 1.0) {
+            return x;
+        } else if (u - 1.0 == -1.0) {
+            return -1;
+        } else {
+            return (u - 1.0) * x/npy_log(u);
+        }
     }
 }
 #endif
@@ -80,13 +85,18 @@ double npy_expm1(double x)
 #ifndef HAVE_LOG1P
 double npy_log1p(double x)
 {
-    const double u = 1. + x;
-    const double d = u - 1.;
-
-    if (d == 0) {
+    if (npy_isinf(x) && x > 0) {
         return x;
-    } else {
-        return npy_log(u) * x / d;
+    }
+    else {
+        const double u = 1. + x;
+        const double d = u - 1.;
+
+        if (d == 0) {
+            return x;
+        } else {
+            return npy_log(u) * x / d;
+        }
     }
 }
 #endif
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
index 8f816169e..3646fd2a9 100644
--- a/numpy/core/tests/test_umath.py
+++ b/numpy/core/tests/test_umath.py
@@ -310,12 +310,28 @@ class TestLog1p(TestCase):
         assert_almost_equal(ncu.log1p(0.2), ncu.log(1.2))
         assert_almost_equal(ncu.log1p(1e-6), ncu.log(1+1e-6))
 
+    def test_special(self):
+        assert_equal(ncu.log1p(np.nan), np.nan)
+        assert_equal(ncu.log1p(np.inf), np.inf)
+        with np.errstate(divide="ignore"):
+            assert_equal(ncu.log1p(-1.), -np.inf)
+        with np.errstate(invalid="ignore"):
+            assert_equal(ncu.log1p(-2.), np.nan)
+            assert_equal(ncu.log1p(-np.inf), np.nan)
+
 
 class TestExpm1(TestCase):
     def test_expm1(self):
         assert_almost_equal(ncu.expm1(0.2), ncu.exp(0.2)-1)
         assert_almost_equal(ncu.expm1(1e-6), ncu.exp(1e-6)-1)
 
+    def test_special(self):
+        assert_equal(ncu.expm1(np.inf), np.inf)
+        assert_equal(ncu.expm1(0.), 0.)
+        assert_equal(ncu.expm1(-0.), -0.)
+        assert_equal(ncu.expm1(np.inf), np.inf)
+        assert_equal(ncu.expm1(-np.inf), -1.)
+
 
 class TestHypot(TestCase, object):
     def test_simple(self):