Implementation of Math.fsum.

Author: ppisl

graalpython/com.oracle.graal.python.test/src/tests/test_math.py

@@ -1220,6 +1220,7 @@ class II(int):
         self.assertEqual(math.ldexp(FF(10), II(12)), 40960.0)
         self.assertRaises(TypeError, math.ldexp, 'Hello', 1000000)
         self.assertRaises(TypeError, math.ldexp, 1, 'Hello')
+        self.assertEqual(math.ldexp(7589167167882033, -48), 26.962138008038156)
 
     def test_trunc(self):
         self.assertEqual(math.trunc(1), 1)
@@ -1293,3 +1294,95 @@ def testmodf(name, result, expected):
         testmodf('modf(MyFloat())', math.modf(MyFloat()), (0.6, 0.0))
         self.assertRaises(TypeError, math.modf, 'ahoj')
         testmodf('modf(BIG_INT)', math.modf(BIG_INT), (0.0, 9.999992432902008e+30))
+
+    def testFsum(self):
+        # math.fsum relies on exact rounding for correct operation.
+        # There's a known problem with IA32 floating-point that causes
+        # inexact rounding in some situations, and will cause the
+        # math.fsum tests below to fail; see issue #2937.  On non IEEE
+        # 754 platforms, and on IEEE 754 platforms that exhibit the
+        # problem described in issue #2937, we simply skip the whole
+        # test.
+
+        # Python version of math.fsum, for comparison.  Uses a
+        # different algorithm based on frexp, ldexp and integer
+        # arithmetic.
+
+        from sys import float_info
+        mant_dig = float_info.mant_dig
+        etiny = float_info.min_exp - mant_dig
+
+        def msum(iterable):
+            """Full precision summation.  Compute sum(iterable) without any
+            intermediate accumulation of error.  Based on the 'lsum' function
+            at http://code.activestate.com/recipes/393090/
+
+            """
+            tmant, texp = 0, 0
+            for x in iterable:
+                mant, exp = math.frexp(x)
+                mant, exp = int(math.ldexp(mant, mant_dig)), exp - mant_dig
+                if texp > exp:
+                    tmant <<= texp-exp
+                    texp = exp
+                else:
+                    mant <<= exp-texp
+                tmant += mant
+            # Round tmant * 2**texp to a float.  The original recipe
+            # used float(str(tmant)) * 2.0**texp for this, but that's
+            # a little unsafe because str -> float conversion can't be
+            # relied upon to do correct rounding on all platforms.
+            tail = max(len(bin(abs(tmant)))-2 - mant_dig, etiny - texp)
+            if tail > 0:
+                h = 1 << (tail-1)
+                tmant = tmant // (2*h) + bool(tmant & h and tmant & 3*h-1)
+                texp += tail
+            return math.ldexp(tmant, texp)
+
+        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:
+                self.fail("test %d failed: got OverflowError, expected %r "
+                          "for math.fsum(%.100r)" % (i, expected, vals))
+            except ValueError:
+                self.fail("test %d failed: got ValueError, expected %r "
+                          "for math.fsum(%.100r)" % (i, expected, vals))
+            self.assertEqual(actual, expected)
+
+        from random import random, gauss, shuffle
+        for j in range(1000):
+            vals = [7, 1e100, -7, -1e100, -9e-20, 8e-20] * 10
+            s = 0
+            for i in range(200):
+                v = gauss(0, random()) ** 7 - s
+                s += v
+                vals.append(v)
+            shuffle(vals)
+
+            s = msum(vals)
+            self.assertEqual(msum(vals), math.fsum(vals))
+
+        self.assertRaises(ValueError, math.fsum, [1., 2, INF, NINF])
+        self.assertEqual(math.fsum([1., 2, INF, INF]), INF)

graalpython/com.oracle.graal.python/src/com/oracle/graal/python/builtins/modules/BuiltinFunctions.java

@@ -198,50 +198,50 @@ public Object absObject(Object object,
             return result;
         }
     }
-    
+
     // bin(object)
     @Builtin(name = BIN, fixedNumOfArguments = 1)
     @TypeSystemReference(PythonArithmeticTypes.class)
     @GenerateNodeFactory
     public abstract static class BinNode extends PythonUnaryBuiltinNode {
-        
+
         public abstract String executeObject(Object x);
-        
+
         private String buildString(boolean isNegative, String number) {
             StringBuilder sb = new StringBuilder();
-            if(isNegative) {
+            if (isNegative) {
                 sb.append('-');
             }
             sb.append("0b");
             sb.append(number);
             return sb.toString();
         }
-        
+
         @Specialization
         public String doL(long x) {
             return buildString(x < 0, Long.toBinaryString(Math.abs(x)));
         }
-        
+
         @Specialization
         public String doD(double x) {
             throw raise(TypeError, "'%p' object cannot be interpreted as an integer", x);
         }
-        
+
         @Specialization
         @TruffleBoundary
         public String doPI(PInt x) {
             BigInteger value = x.getValue();
             return buildString(value.compareTo(BigInteger.ZERO) == -1, value.abs().toString(2));
         }
-        
+
         @Specialization
         public String doO(Object x,
-                @Cached("create()") MathModuleBuiltins.ConvertToIntNode toIntNode,
-                @Cached("create()") BinNode recursiveNode) {
+                        @Cached("create()") MathModuleBuiltins.ConvertToIntNode toIntNode,
+                        @Cached("create()") BinNode recursiveNode) {
             Object value = toIntNode.execute(x);
             return recursiveNode.executeObject(value);
         }
-        
+
         protected BinNode create() {
             return BuiltinFunctionsFactory.BinNodeFactory.create();
         }

graalpython/com.oracle.graal.python/src/com/oracle/graal/python/builtins/modules/MathModuleBuiltins.java

@@ -46,11 +46,13 @@
 import com.oracle.graal.python.nodes.PGuards;
 import com.oracle.graal.python.nodes.SpecialMethodNames;
 import com.oracle.graal.python.nodes.call.special.LookupAndCallUnaryNode;
+import com.oracle.graal.python.nodes.control.GetIteratorNode;
 import com.oracle.graal.python.nodes.function.PythonBuiltinBaseNode;
 import com.oracle.graal.python.nodes.function.PythonBuiltinNode;
 import com.oracle.graal.python.nodes.function.builtins.PythonBinaryBuiltinNode;
 import com.oracle.graal.python.nodes.function.builtins.PythonUnaryBuiltinNode;
 import com.oracle.graal.python.nodes.truffle.PythonArithmeticTypes;
+import com.oracle.graal.python.runtime.exception.PException;
 import com.oracle.graal.python.runtime.exception.PythonErrorType;
 import static com.oracle.graal.python.runtime.exception.PythonErrorType.NotImplementedError;
 import com.oracle.graal.python.runtime.object.PythonObjectFactory;
@@ -65,6 +67,7 @@
 import com.oracle.truffle.api.dsl.Specialization;
 import com.oracle.truffle.api.dsl.TypeSystemReference;
 import com.oracle.truffle.api.profiles.ConditionProfile;
+import java.util.Arrays;
 
 @CoreFunctions(defineModule = "math")
 public class MathModuleBuiltins extends PythonBuiltins {
@@ -968,7 +971,7 @@ public double ldexpLD(long mantissa, double exp) {
 
         @Specialization
         public double ldexpLL(long mantissa, long exp) {
-            return exceptInfinity(Math.scalb(mantissa, makeInt(exp)), mantissa);
+            return exceptInfinity(Math.scalb((double) mantissa, makeInt(exp)), mantissa);
         }
 
         @Specialization
@@ -980,7 +983,7 @@ public double ldexpDPI(double mantissa, PInt exp) {
         @Specialization
         @TruffleBoundary
         public double ldexpLPI(long mantissa, PInt exp) {
-            return exceptInfinity(Math.scalb(mantissa, makeInt(exp)), mantissa);
+            return exceptInfinity(Math.scalb((double) mantissa, makeInt(exp)), mantissa);
         }
 
         @Specialization
@@ -1048,6 +1051,127 @@ public PTuple frexpO(Object value,
         }
     }
 
+    @Builtin(name = "fsum", fixedNumOfArguments = 1)
+    @ImportStatic(PGuards.class)
+    @GenerateNodeFactory
+    public abstract static class FsumNode extends PythonUnaryBuiltinNode {
+
+        /*
+         * This implementation is taken from CPython. The performance is not good. Should be faster.
+         * It can be easily replace with much simpler code based on BigDecimal:
+         * 
+         * BigDecimal result = BigDecimal.ZERO;
+         * 
+         * in cycle just: result = result.add(BigDecimal.valueof(x); ... The current implementation
+         * is little bit faster. The testFSum in test_math.py takes in different implementations:
+         * CPython ~0.6s CurrentImpl: ~14.3s Using BigDecimal: ~15.1
+         */
+        @Specialization
+        @TruffleBoundary
+        public double doIt(Object iterable,
+                        @Cached("create()") GetIteratorNode getIterator,
+                        @Cached("create(__NEXT__)") LookupAndCallUnaryNode next,
+                        @Cached("create()") ConvertToFloatNode toFloat,
+                        @Cached("createBinaryProfile()") ConditionProfile stopProfile) {
+            Object iterator = getIterator.executeWith(iterable);
+            double x, y, t, hi, lo = 0, yr, inf_sum = 0, special_sum = 0, sum;
+            double xsave;
+            int i, j, n = 0, arayLength = 32;
+            double[] p = new double[arayLength];
+            while (true) {
+                try {
+                    x = toFloat.execute(next.executeObject(iterator));
+                } catch (PException e) {
+                    e.expectStopIteration(getCore(), stopProfile);
+                    break;
+                }
+                xsave = x;
+                for (i = j = 0; j < n; j++) { /* for y in partials */
+                    y = p[j];
+                    if (Math.abs(x) < Math.abs(y)) {
+                        t = x;
+                        x = y;
+                        y = t;
+                    }
+                    hi = x + y;
+                    yr = hi - x;
+                    lo = y - yr;
+                    if (lo != 0.0) {
+                        p[i++] = lo;
+                    }
+                    x = hi;
+                }
+
+                n = i;
+                if (x != 0.0) {
+                    if (!Double.isFinite(x)) {
+                        /*
+                         * a nonfinite x could arise either as a result of intermediate overflow, or
+                         * as a result of a nan or inf in the summands
+                         */
+                        if (Double.isFinite(xsave)) {
+                            throw raise(OverflowError, "intermediate overflow in fsum");
+                        }
+                        if (Double.isInfinite(xsave)) {
+                            inf_sum += xsave;
+                        }
+                        special_sum += xsave;
+                        /* reset partials */
+                        n = 0;
+                    } else if (n >= arayLength) {
+                        arayLength += arayLength;
+                        p = Arrays.copyOf(p, arayLength);
+                    } else {
+                        p[n++] = x;
+                    }
+                }
+            }
+
+            if (special_sum != 0.0) {
+                if (Double.isNaN(inf_sum)) {
+                    throw raise(ValueError, "-inf + inf in fsum");
+                } else {
+                    sum = special_sum;
+                    return sum;
+                }
+            }
+
+            hi = 0.0;
+            if (n > 0) {
+                hi = p[--n];
+                /*
+                 * sum_exact(ps, hi) from the top, stop when the sum becomes inexact.
+                 */
+                while (n > 0) {
+                    x = hi;
+                    y = p[--n];
+                    assert (Math.abs(y) < Math.abs(x));
+                    hi = x + y;
+                    yr = hi - x;
+                    lo = y - yr;
+                    if (lo != 0.0)
+                        break;
+                }
+                /*
+                 * Make half-even rounding work across multiple partials. Needed so that sum([1e-16,
+                 * 1, 1e16]) will round-up the last digit to two instead of down to zero (the 1e-16
+                 * makes the 1 slightly closer to two). With a potential 1 ULP rounding error
+                 * fixed-up, math.fsum() can guarantee commutativity.
+                 */
+                if (n > 0 && ((lo < 0.0 && p[n - 1] < 0.0) ||
+                                (lo > 0.0 && p[n - 1] > 0.0))) {
+                    y = lo * 2.0;
+                    x = hi + y;
+                    yr = x - hi;
+                    if (y == yr) {
+                        hi = x;
+                    }
+                }
+            }
+            return hi;
+        }
+    }
+
     @Builtin(name = "gcd", fixedNumOfArguments = 2)
     @TypeSystemReference(PythonArithmeticTypes.class)
     @GenerateNodeFactory