Implement Math.sumPrecise (#1966)

Jint.Tests.Test262/Test262Harness.settings.json

@@ -11,7 +11,6 @@
     "decorators",
     "explicit-resource-management",
     "iterator-helpers",
-    "Math.sumPrecise",
     "regexp-lookbehind",
     "regexp-modifiers",
     "regexp-unicode-property-escapes",

Jint/Native/Math/MathInstance.cs

@@ -1154,7 +1154,7 @@ private JsValue SumPrecise(JsValue thisObject, JsValue[] arguments)
             return state;
         }
 
-        return sum.FSum();
+        return Math.SumPrecise.Sum(sum);
     }
 
     private static double[] Coerced(JsValue[] arguments)

Jint/Native/Math/SumPrecise.cs

@@ -0,0 +1,194 @@
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Jint.Native.Math;
+
+/// <summary>
+/// https://raw.githubusercontent.com/es-shims/Math.sumPrecise/main/sum.js
+/// adapted from https://github.com/tc39/proposal-math-sum/blob/f4286d0a9d8525bda61be486df964bf2527c8789/polyfill/polyfill.mjs
+/// https://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
+/// Shewchuk's algorithm for exactly floating point addition
+/// as implemented in Python's fsum: https://github.com/python/cpython/blob/48dfd74a9db9d4aa9c6f23b4a67b461e5d977173/Modules/mathmodule.c#L1359-L1474
+/// adapted to handle overflow via an additional "biased" partial, representing 2**1024 times its actual value
+/// </summary>
+internal static class SumPrecise
+{
+    // exponent 11111111110, significand all 1s
+    private const double MaxDouble = 1.79769313486231570815e+308; // i.e. (2**1024 - 2**(1023 - 52))
+
+    // exponent 11111111110, significand all 1s except final 0
+    private const double PenultimateDouble = 1.79769313486231550856e+308; // i.e. (2**1024 - 2 * 2**(1023 - 52))
+
+    private const double Two1023 = 8.98846567431158e+307; // 2 ** 1023
+
+    // exponent 11111001010, significand all 0s
+    private const double MaxUlp = MaxDouble - PenultimateDouble; // 1.99584030953471981166e+292, i.e. 2**(1023 - 52)
+
+    [StructLayout(LayoutKind.Auto)]
+    private readonly record struct TwoSumResult(double Hi, double Lo);
+
+    // prerequisite: $abs(x) >= $abs(y)
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static TwoSumResult TwoSum(double x, double y)
+    {
+        var hi = x + y;
+        var lo = y - (hi - x);
+        return new TwoSumResult(hi, lo);
+    }
+
+    // preconditions:
+    //  - array only contains numbers
+    //  - none of them are -0, NaN, or ±Infinity
+    //  - all of them are finite
+    [MethodImpl(512)]
+    internal static double Sum(List<double> array)
+    {
+        double hi, lo;
+
+        List<double> partials = [];
+
+        var overflow = 0; // conceptually 2**1024 times this value; the final partial is biased by this amount
+
+        var index = -1;
+
+        // main loop
+        while (index + 1 < array.Count)
+        {
+            var x = +array[++index];
+
+            var actuallyUsedPartials = 0;
+            for (var j = 0; j < partials.Count; j += 1)
+            {
+                var y = partials[j];
+
+                if (System.Math.Abs(x) < System.Math.Abs(y))
+                {
+                    var tmp = x;
+                    x = y;
+                    y = tmp;
+                }
+
+                (hi, lo) = TwoSum(x, y);
+
+                if (double.IsPositiveInfinity(System.Math.Abs(hi)))
+                {
+                    var sign = double.IsPositiveInfinity(hi) ? 1 : -1;
+                    overflow += sign;
+
+                    x = x - sign * Two1023 - sign * Two1023;
+                    if (System.Math.Abs(x) < System.Math.Abs(y))
+                    {
+                        var tmp2 = x;
+                        x = y;
+                        y = tmp2;
+                    }
+
+                    var s2 = TwoSum(x, y);
+                    hi = s2.Hi;
+                    lo = s2.Lo;
+                }
+
+                if (lo != 0)
+                {
+                    partials[actuallyUsedPartials] = lo;
+                    actuallyUsedPartials += 1;
+                }
+
+                x = hi;
+            }
+
+            while (partials.Count > actuallyUsedPartials)
+            {
+                partials.RemoveAt(partials.Count - 1);
+            }
+
+            if (x != 0)
+            {
+                partials.Add(x);
+            }
+        }
+
+        // compute the exact sum of partials, stopping once we lose precision
+        var n = partials.Count - 1;
+        hi = 0;
+        lo = 0;
+
+        if (overflow != 0)
+        {
+            var next = n >= 0 ? partials[n] : 0;
+            n -= 1;
+            if (System.Math.Abs(overflow) > 1 || (overflow > 0 && next > 0) || (overflow < 0 && next < 0))
+            {
+                return overflow > 0 ? double.PositiveInfinity : double.NegativeInfinity;
+            }
+
+            // here we actually have to do the arithmetic
+            // drop a factor of 2 so we can do it without overflow
+            // assert($abs(overflow) === 1)
+            var s = TwoSum(overflow * Two1023, next / 2);
+            hi = s.Hi;
+            lo = s.Lo;
+            lo *= 2;
+            if (double.IsPositiveInfinity(System.Math.Abs(2 * hi)))
+            {
+                // stupid edge case: rounding to the maximum value
+                // MAX_DOUBLE has a 1 in the last place of its significand, so if we subtract exactly half a ULP from 2**1024, the result rounds away from it (i.e. to infinity) under ties-to-even
+                // but if the next partial has the opposite sign of the current value, we need to round towards MAX_DOUBLE instead
+                // this is the same as the "handle rounding" case below, but there's only one potentially-finite case we need to worry about, so we just hardcode that one
+                if (hi > 0)
+                {
+                    if (hi == Two1023 && lo == -(MaxUlp / 2) && n >= 0 && partials[n] < 0)
+                    {
+                        return MaxDouble;
+                    }
+
+                    return double.PositiveInfinity;
+                }
+
+                if (hi == -Two1023 && lo == (MaxUlp / 2) && n >= 0 && partials[n] > 0)
+                {
+                    return -MaxDouble;
+                }
+
+                return double.NegativeInfinity;
+            }
+
+            if (lo != 0)
+            {
+                partials[n + 1] = lo;
+                n += 1;
+                lo = 0;
+            }
+
+            hi *= 2;
+        }
+
+        while (n >= 0)
+        {
+            var x1 = hi;
+            var y1 = partials[n];
+            n -= 1;
+            // assert: $abs(x1) > $abs(y1)
+            (hi, lo) = TwoSum(x1, y1);
+            if (lo != 0)
+            {
+                break; // eslint-disable-line no-restricted-syntax
+            }
+        }
+
+        // handle rounding
+        // when the roundoff error is exactly half of the ULP for the result, we need to check one more partial to know which way to round
+        if (n >= 0 && ((lo < 0.0 && partials[n] < 0.0) || (lo > 0.0 && partials[n] > 0.0)))
+        {
+            var y2 = lo * 2.0;
+            var x2 = hi + y2;
+            var yr = x2 - hi;
+            if (y2 == yr)
+            {
+                hi = x2;
+            }
+        }
+
+        return hi;
+    }
+}

README.md

@@ -133,6 +133,7 @@ and many more.
 - ✔ Float16Array (Requires NET 6 or higher)
 - ✔ Import attributes
 - ✔ JSON modules
+- ✔ Math.sumPrecise
 - ✔ `Promise.try`
 - ✔ Set methods (`intersection`, `union`, `difference`, `symmetricDifference`, `isSubsetOf`, `isSupersetOf`, `isDisjointFrom`)
 - ✔ ShadowRealm