[pypy-commit] pypy py3.3: Add rfloat.log2(), which tries to be exact for powers of two.

[amauryfa]

Author: Amaury Forgeot d'Arc <amauryfa at gmail.com>
Branch: py3.3
Changeset: r75026:d3c7a156cc70
Date: 2014-12-18 22:35 +0100
http://bitbucket.org/pypy/pypy/changeset/d3c7a156cc70/

Log:	Add rfloat.log2(), which tries to be exact for powers of two. Use it
	in math.log2.

diff --git a/pypy/module/math/interp_math.py b/pypy/module/math/interp_math.py
--- a/pypy/module/math/interp_math.py
+++ b/pypy/module/math/interp_math.py
@@ -202,6 +202,8 @@
             x = _get_double(space, w_x)
             if base == 10.0:
                 result = math.log10(x)
+            elif base == 2.0:
+                result = rfloat.log2(x)
             else:
                 result = math.log(x)
                 if base != 0.0:
diff --git a/rpython/rlib/rfloat.py b/rpython/rlib/rfloat.py
--- a/rpython/rlib/rfloat.py
+++ b/rpython/rlib/rfloat.py
@@ -281,6 +281,35 @@
             return (u - 1.) * x / math.log(u)
         return math.exp(x) - 1.
 
+def log2(x):
+    # Uses an algorithm that should:
+    #   (a) produce exact results for powers of 2, and
+    #   (b) be monotonic, assuming that the system log is monotonic.
+    if not isfinite(x):
+        if isnan(x):
+            return x  # log2(nan) = nan
+        elif x > 0.0:
+            return x  # log2(+inf) = +inf
+        else:
+            # log2(-inf) = nan, invalid-operation
+            raise ValueError("math domain error")
+
+    if x > 0.0:
+        if 0:  # HAVE_LOG2
+            return math.log2(x)
+        m, e = math.frexp(x)
+        # 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.
+        if x >= 1.0:
+            return math.log(2.0 * m) / math.log(2.0) + (e - 1)
+        else:
+            return math.log(m) / math.log(2.0) + e
+    else:
+        raise ValueError("math domain error")
+
 def round_away(x):
     # round() from libm, which is not available on all platforms!
     absx = abs(x)
diff --git a/rpython/rlib/test/test_rfloat.py b/rpython/rlib/test/test_rfloat.py
--- a/rpython/rlib/test/test_rfloat.py
+++ b/rpython/rlib/test/test_rfloat.py
@@ -265,3 +265,12 @@
                 if s.strip(): # empty s raises OperationError directly
                     py.test.raises(ParseStringError, string_to_float, s)
     py.test.raises(ParseStringError, string_to_float, "")
+
+def test_log2():
+    from rpython.rlib import rfloat
+    assert rfloat.log2(1.0) == 0.0
+    assert rfloat.log2(2.0) == 1.0
+    assert rfloat.log2(2.0**1023) == 1023.0
+    assert 1.584 < rfloat.log2(3.0) < 1.585
+    py.test.raises(ValueError, rfloat.log2, 0)
+    py.test.raises(ValueError, rfloat.log2, -1)