[pypy-svn] r75844 - in pypy/branch/fast-forward/pypy/module/math: . test

[benjamin at codespeak.net]

Author: benjamin
Date: Mon Jul  5 18:54:12 2010
New Revision: 75844

Modified:
   pypy/branch/fast-forward/pypy/module/math/__init__.py
   pypy/branch/fast-forward/pypy/module/math/interp_math.py
   pypy/branch/fast-forward/pypy/module/math/test/test_math.py
Log:
add math.fsum: floating point summation with partials

Modified: pypy/branch/fast-forward/pypy/module/math/__init__.py
==============================================================================
--- pypy/branch/fast-forward/pypy/module/math/__init__.py	(original)
+++ pypy/branch/fast-forward/pypy/module/math/__init__.py	Mon Jul  5 18:54:12 2010
@@ -38,5 +38,6 @@
        'isinf'          : 'interp_math.isinf',
        'isnan'          : 'interp_math.isnan',
        'trunc'          : 'interp_math.trunc',
+       'fsum'           : 'interp_math.fsum',
 }
 

Modified: pypy/branch/fast-forward/pypy/module/math/interp_math.py
==============================================================================
--- pypy/branch/fast-forward/pypy/module/math/interp_math.py	(original)
+++ pypy/branch/fast-forward/pypy/module/math/interp_math.py	Mon Jul  5 18:54:12 2010
@@ -299,3 +299,68 @@
     """
     return math1(space, math.acos, x)
 acos.unwrap_spec = [ObjSpace, float]
+
+def fsum(space, w_iterable):
+    """Sum an iterable of floats, trying to keep precision."""
+    w_iter = space.iter(w_iterable)
+    inf_sum = special_sum = 0.0
+    partials = []
+    while True:
+        try:
+            w_value = space.next(w_iter)
+        except OperationError, e:
+            if not e.match(space, space.w_StopIteration):
+                raise
+            break
+        v = space.float_w(w_value)
+        original = v
+        added = 0
+        for y in partials:
+            if abs(v) < abs(y):
+                v, y = y, v
+            hi = v + y
+            yr = hi - v
+            lo = y - yr
+            if lo != 0.0:
+                partials[added] = lo
+                added += 1
+            v = hi
+        del partials[added:]
+        if v != 0.0:
+            if rarithmetic.isinf(v) or rarithmetic.isnan(v):
+                if (not rarithmetic.isinf(original) and
+                    not rarithmetic.isnan(original)):
+                    raise OperationError(space.w_OverflowError,
+                                         space.wrap("intermediate overflow"))
+                if isinf(original):
+                    inf_sum += original
+                special_sum += original
+                del partials[:]
+            else:
+                partials.append(v)
+    if special_sum != 0.0:
+        if rarithmetic.isnan(special_sum):
+            raise OperationError(space.w_ValueError, space.wrap("-inf + inf"))
+        return space.wrap(special_sum)
+    hi = 0.0
+    if partials:
+        hi = partials[-1]
+        j = 0
+        for j in range(len(partials) - 2, -1, -1):
+            v = hi
+            y = partials[j]
+            assert abs(y) < abs(v)
+            hi = v + y
+            yr = hi - v
+            lo = y - yr
+            if lo != 0.0:
+                break
+        if j > 0 and (lo < 0.0 and partials[j - 1] < 0.0 or
+                      lo > 0.0 and partials[j - 1] > 0.0):
+            y = lo * 2.0
+            v = hi + y
+            yr = v - hi
+            if y == yr:
+                hi = v
+    return space.wrap(hi)
+fsum.unwrap_spec = [ObjSpace, W_Root]

Modified: pypy/branch/fast-forward/pypy/module/math/test/test_math.py
==============================================================================
--- pypy/branch/fast-forward/pypy/module/math/test/test_math.py	(original)
+++ pypy/branch/fast-forward/pypy/module/math/test/test_math.py	Mon Jul  5 18:54:12 2010
@@ -35,3 +35,41 @@
                 if not ok:
                     raise AssertionError("%s(%s): got %s" % (
                         fnname, ', '.join(map(str, args)), got))
+
+    def test_fsum(self):
+        # Python version of math.fsum, for comparison.  Uses a
+        # different algorithm based on frexp, ldexp and integer
+        # arithmetic.
+        import math
+
+        test_values = [
+            ([], 0.0),
+            ([0.0], 0.0),
+            ([1e100, 1.0, -1e100, 1e-100, 1e50, -1.0, -1e50], 1e-100),
+            ([2.0**53, -0.5, -2.0**-54], 2.0**53-1.0),
+            ([2.0**53, 1.0, 2.0**-100], 2.0**53+2.0),
+            ([2.0**53+10.0, 1.0, 2.0**-100], 2.0**53+12.0),
+            ([2.0**53-4.0, 0.5, 2.0**-54], 2.0**53-3.0),
+            ([1./n for n in range(1, 1001)],
+             float.fromhex('0x1.df11f45f4e61ap+2')),
+            ([(-1.)**n/n for n in range(1, 1001)],
+             float.fromhex('-0x1.62a2af1bd3624p-1')),
+            ([1.7**(i+1)-1.7**i for i in range(1000)] + [-1.7**1000], -1.0),
+            ([1e16, 1., 1e-16], 10000000000000002.0),
+            ([1e16-2., 1.-2.**-53, -(1e16-2.), -(1.-2.**-53)], 0.0),
+            # exercise code for resizing partials array
+            ([2.**n - 2.**(n+50) + 2.**(n+52) for n in range(-1074, 972, 2)] +
+             [-2.**1022],
+             float.fromhex('0x1.5555555555555p+970')),
+            ]
+
+        for i, (vals, expected) in enumerate(test_values):
+            try:
+                actual = math.fsum(vals)
+            except OverflowError:
+                py.test.fail("test %d failed: got OverflowError, expected %r "
+                          "for math.fsum(%.100r)" % (i, expected, vals))
+            except ValueError:
+                py.test.fail("test %d failed: got ValueError, expected %r "
+                          "for math.fsum(%.100r)" % (i, expected, vals))
+            assert actual == expected