Merge pull request #4790 from juliantaylor/log2-windows

BUG: improve log2 windows compiler fallback of log2

2 files changed, 51 insertions, 1 deletions

diff --git a/numpy/core/src/npymath/npy_math.c.src b/numpy/core/src/npymath/npy_math.c.src
index 1ca7033af..05af0b132 100644
--- a/numpy/core/src/npymath/npy_math.c.src
+++ b/numpy/core/src/npymath/npy_math.c.src
@@ -299,7 +299,35 @@ double npy_exp2(double x)
 #ifndef HAVE_LOG2
 double npy_log2(double x)
 {
-    return NPY_LOG2E*npy_log(x);
+#ifdef HAVE_FREXP
+    if (!npy_isfinite(x) || x <= 0.) {
+        /* special value result */
+        return npy_log(x);
+    }
+    else {
+        /*
+         * fallback implementation copied from python3.4 math.log2
+         * provides int(log(2**i)) == i for i 1-64 in default rounding mode.
+         *
+         * We want log2(m * 2**e) == log(m) / log(2) + e.  Care is needed when
+         * x is just greater than 1.0: in that case e is 1, log(m) is negative,
+         * and we get significant cancellation error from the addition of
+         * log(m) / log(2) to e.  The slight rewrite of the expression below
+         * avoids this problem.
+         */
+        int e;
+        double m = frexp(x, &e);
+        if (x >= 1.0) {
+            return log(2.0 * m) / log(2.0) + (e - 1);
+        }
+        else {
+            return log(m) / log(2.0) + e;
+        }
+    }
+#else
+    /* does not provide int(log(2**i)) == i */
+    return NPY_LOG2E * npy_log(x);
+#endif
 }
 #endif
 
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
index 5f490d838..b3ddc2398 100644
--- a/numpy/core/tests/test_umath.py
+++ b/numpy/core/tests/test_umath.py
@@ -2,6 +2,7 @@ from __future__ import division, absolute_import, print_function
 
 import sys
 import platform
+import warnings
 
 from numpy.testing import *
 from numpy.testing.utils import _gen_alignment_data
@@ -189,6 +190,27 @@ class TestLog2(TestCase):
             yf = np.array(y, dtype=dt)
             assert_almost_equal(np.log2(xf), yf)
 
+    def test_log2_ints(self):
+        # a good log2 implementation should provide this,
+        # might fail on OS with bad libm
+        for i in range(1, 65):
+            v = np.log2(2.**i)
+            assert_equal(v, float(i), err_msg='at exponent %d' % i)
+
+    def test_log2_special(self):
+        assert_equal(np.log2(1.), 0.)
+        assert_equal(np.log2(np.inf), np.inf)
+        assert_(np.isnan(np.log2(np.nan)))
+
+        with warnings.catch_warnings(record=True) as w:
+            warnings.filterwarnings('always', '', RuntimeWarning)
+            assert_(np.isnan(np.log2(-1.)))
+            assert_(np.isnan(np.log2(-np.inf)))
+            assert_equal(np.log2(0.), -np.inf)
+            assert_(w[0].category is RuntimeWarning)
+            assert_(w[1].category is RuntimeWarning)
+            assert_(w[2].category is RuntimeWarning)
+
 
 class TestExp2(TestCase):
     def test_exp2_values(self) :