Advance to stage 1 and make summation precise (#2)

.github/workflows/build.yml

@@ -11,10 +11,11 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v2
-    - uses: actions/setup-node@v3
+    - uses: actions/checkout@v4
+    - uses: actions/setup-node@v4
       with:
-        node-version: 16
+        node-version: 20
     - run: npm ci
     - run: npm run build
     - run: npm run check-format
+    - run: npm run test

.github/workflows/deploy.yml

@@ -12,13 +12,13 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v2
-    - uses: actions/setup-node@v3
+    - uses: actions/checkout@v4
+    - uses: actions/setup-node@v4
       with:
-        node-version: 16
+        node-version: 20
     - run: npm ci
     - run: npm run build
-    - uses: JamesIves/github-pages-deploy-action@v4.4.3
+    - uses: JamesIves/github-pages-deploy-action@v4.5.0
       with:
         branch: gh-pages
         folder: dist

README.md

@@ -8,31 +8,33 @@ Authors: Kevin Gibbons
 
 Champions: Kevin Gibbons
 
-This proposal is at Stage 0 of [the TC39 process](https://tc39.es/process-document/) - it has not yet been presented to committee.
+This proposal is at Stage 1 of [the TC39 process](https://tc39.es/process-document/): the proposal is under consideration.
 
 ## Motivation
 
 Summing a list is a very common operation and is one of the few remaining use cases for `Array.prototype.reduce`. Better to let users express that operation directy.
 
-Also, floating point summation can get more precise results than the naive `.reduce((a, b) => a + b, 0)` approach using [more clever algorithms](https://en.wikipedia.org/wiki/Kahan_summation_algorithm) with very little overhead, a fact which few JavaScript programmers are aware of (and even among those who are, most wouldn't bother doing it). We can make it easy to reach for the better option.
+Also, summing a list of floating point numbers can be done more precisely than the naive `.reduce((a, b) => a + b, 0)` approach using more clever algorithms, a fact which few JavaScript programmers are aware of (and even among those who are, most wouldn't bother doing it). We can make it easy to reach for the better option.
 
 ## Proposal
 
-Add a variadic `Math.sum` method which returns the sum of its arguments using a more clever algorithm than naive summation.
+Add a variadic `Math.sum` method which returns the sum of its arguments using a more precise algorithm than naive summation.
 
 ```js
-let values = Array.from({ length: 10 }).fill(0.1);
+let values = [1e20, 0.1, -1e20];
 
-values.reduce((a, b) => a + b, 0); // 0.9999999999999999
+values.reduce((a, b) => a + b, 0); // 0
 
-Math.sum(...values); // 1
+Math.sum(...values); // 0.1
 ````
 
 ## Questions
 
 ### Which algorithm?
 
-I'm hoping to get away with leaving it up to implementations. If we have to pick one, I would go with Neumaier's variant of Kahan summation. Python [adopted it](https://github.com/python/cpython/issues/100425) for their built-in `sum` in 3.12, with about 25-30% slowdown relative to naive summation. I think that's easily worth it.
+Instead of specifying any particular algorithm, this proposal requires the maximally correct answer - that is, the answer you'd get if you did arbitrary-precision arithmetic, then converted the result back to floats. This can be done without actually needing arbitrary-precision arithmetic. One way of doing this is given in [Shewchuk '96](./Shewchuk-robust-arithmetic.pdf), which I've [implemented in JS](./polyfill/polyfill.mjs) (plus some details to handle intermediate overflow). Other strategies are possible.
+
+Python's [`math.fsum`](https://docs.python.org/3/library/math.html#math.fsum) is currently implemented using the same algorithm (though without handling intermediate overflow).
 
 ### What if I have a long list, such that I can't reasonably put it on the stack with `...args`?
 

package.json

@@ -4,10 +4,12 @@
     "build": "mkdir -p dist && ecmarkup --lint-spec --strict --load-biblio @tc39/ecma262-biblio --verbose --assets-dir dist spec.html dist/index.html",
     "build-biblio": "ecmarkup --write-biblio proposal-iterator-helpers-biblio.json --lint-spec --strict --load-biblio @tc39/ecma262-biblio --verbose spec.html /dev/null",
     "format": "emu-format --write spec.html",
-    "check-format": "emu-format --check spec.html"
+    "check-format": "emu-format --check spec.html",
+    "test": "node --test"
   },
   "dependencies": {
     "@tc39/ecma262-biblio": "2.1.2657",
-    "ecmarkup": "18.1.0"
+    "ecmarkup": "18.1.0",
+    "xoshiro128": "^0.1.0"
   }
 }

polyfill/float-tools.mjs

@@ -0,0 +1,26 @@
+let bits = new ArrayBuffer(8);
+let view = new DataView(bits);
+export function floatFromParts({ sign, exponent, significand }) {
+  view.setBigUint64(0, BigInt(significand));
+  let top16Bits = (sign << 15) | (exponent << 4) | (view.getUint8(1) & 0xf);
+  view.setUint16(0, top16Bits);
+  return view.getFloat64(0);
+}
+
+export function partsFromFloat(double) {
+  view.setFloat64(0, double);
+  let sign = view.getUint8(0) >>> 7;
+  let exponent = (view.getUint16(0) >>> 4) & 0x7ff; // 11 bits. don't forget the bias!
+  let significand = Number(view.getBigUint64(0) & 0xfffffffffffffn); // 52 bits
+  return { sign, exponent, significand };
+}
+
+let interestingExponents = [2046, 2045, 1994, 1995, 1993, 0, 1, 2, 1021, 1022, 1023, 1024, 1025, 1026];
+let interestingSignificands = [0b1111111111111111111111111111111111111111111111111111, 0b1000000000000000000000000000000000000000000000000000, 0b1000000000000000000000000000000000000000000000000001, 0b1111111111111111111111111111111111111111111111111110, 0b111111111111111111111111111111111111111111111111111, 0, 1, 2];
+
+export function randomFloat(random) {
+  let sign = random.int(0, 1);
+  let exponent = random.boolean() ? random.pick(interestingExponents) : random.int(0, 2046);
+  let significand = random.boolean() ? random.pick(interestingSignificands) : Math.floor(random.float(0, 2**52));
+  return floatFromParts({ sign, exponent, significand });
+}

polyfill/full-precision-with-bigints.mjs

@@ -0,0 +1,44 @@
+import { floatFromParts, partsFromFloat } from './float-tools.mjs';
+
+export function floatToBigInt(double) {
+  if (!Number.isFinite(double)) throw new Error('todo');
+  let { sign, exponent, significand } = partsFromFloat(double);
+  let magnitude;
+  if (exponent === 0) {
+    magnitude = BigInt(significand);
+  } else {
+    significand = 2n ** 52n + BigInt(significand);
+    magnitude = 2n ** (BigInt(exponent - 1)) * significand;
+  }
+  return sign ? -magnitude : magnitude;
+}
+
+export function bigIntToFloat(bigint) {
+  let sign = bigint < 0 ? 1 : 0;
+  let magnitude = bigint < 0 ? -bigint : bigint;
+  let binary = magnitude.toString(2);
+  if (binary.length <= 52) {
+    // subnormal
+    let significand = parseInt(binary, 2);
+    let exponent = 0;
+    return floatFromParts({ sign, exponent, significand });
+  }
+
+  let significandString = binary.slice(1, 53);
+  let significand = parseInt(significandString, 2);
+  let exponent = binary.length - 52;
+  // round up if we are at least half way. if up is towards-even we always round up; otherwise we round up iff it is above the halfway point, i.e. there is a non-zero bit after the first
+  let roundUp = binary[53] === '1' && (binary[52] === '1' || binary.slice(54).includes('1'));
+  if (roundUp) {
+    significand += 1;
+    if (significand === 2 ** 52) {
+      significand = 0;
+      exponent += 1;
+    }
+  }
+  if (exponent > 2046) {
+    return sign ? -Infinity : Infinity;
+  }
+
+  return floatFromParts({ sign, exponent, significand });
+}

polyfill/polyfill.mjs

@@ -0,0 +1,183 @@
+/*
+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
+*/
+
+// exponent 11111111110, significand all 1s
+const MAX_DOUBLE = 1.79769313486231570815e+308; // i.e. (2**1024 - 2**(1023 - 52))
+
+// exponent 11111111110, significand all 1s except final 0
+const PENULTIMATE_DOUBLE = 1.79769313486231550856e+308; // i.e. (2**1024 - 2 * 2**(1023 - 52))
+
+// exponent 11111001010, significand all 0s
+const MAX_ULP = MAX_DOUBLE - PENULTIMATE_DOUBLE; // 1.99584030953471981166e+292, i.e. 2**(1023 - 52)
+
+// prerequisite: Math.abs(x) >= Math.abs(y)
+function twosum(x, y) {
+  let hi = x + y;
+  let lo = y - (hi - x);
+  return { hi, lo };
+}
+
+export function sum(iterable) {
+  let partials = [];
+
+  let overflow = 0; // conceptually 2**1024 times this value; the final partial
+
+  // for purposes of the polyfill we're going to ignore closing the iterator, sorry
+  let iterator = iterable[Symbol.iterator]();
+  let next = iterator.next.bind(iterator);
+
+  // in C this would be done using a goto
+  function drainNonFiniteValue(current) {
+    while (true) {
+      let { done, value } = next();
+      if (done) {
+        return current;
+      }
+      if (!Number.isFinite(value)) {
+        // summing any distinct two of the three non-finite values gives NaN
+        // summing any one of them with itself gives itself
+        if (!Object.is(value, current)) {
+          current = NaN;
+        }
+      }
+    }
+  }
+
+  // handle list of -0 special case
+  while (true) {
+    let { done, value } = next();
+    if (done) {
+      return -0;
+    }
+    if (!Object.is(value, -0)) {
+      if (!Number.isFinite(value)) {
+        return drainNonFiniteValue(val);
+      }
+      partials.push(value);
+      break;
+    }
+  }
+
+  // main loop
+  while (true) {
+    let { done, value } = next();
+    if (done) {
+      break;
+    }
+    let x = +value;
+    if (!Number.isFinite(x)) {
+      return drainNonFiniteValue(x);
+    }
+
+    // we're updating partials in place, but it is maybe easier to understand if you think of it as making a new copy
+    let actuallyUsedPartials = 0;
+    // let newPartials = [];
+    for (let y of partials) {
+      if (Math.abs(x) < Math.abs(y)) {
+        [x, y] = [y, x];
+      }
+      let { hi, lo } = twosum(x, y);
+      if (Math.abs(hi) === Infinity) {
+        let sign = hi === Infinity ? 1 : -1;
+        overflow += sign;
+        if (Math.abs(overflow) >= 2**53) {
+          throw new RangeError('overflow');
+        }
+
+        x = (x - sign * 2**1023) - sign * 2**1023;
+        if (Math.abs(x) < Math.abs(y)) {
+          [x, y] = [y, x];
+        }
+        0, { hi, lo } = twosum(x, y);
+      }
+      if (lo !== 0) {
+        partials[actuallyUsedPartials] = lo;
+        ++actuallyUsedPartials;
+        // newPartials.push(lo);
+      }
+      x = hi;
+    }
+    partials.length = actuallyUsedPartials;
+    // assert.deepStrictEqual(partials, newPartials)
+    // partials = newPartials
+
+    if (x !== 0) {
+      partials.push(x);
+    }
+  }
+
+  // compute the exact sum of partials, stopping once we lose precision
+  let n = partials.length - 1;
+  let hi = 0;
+  let lo = 0;
+
+  if (overflow !== 0) {
+    let next = n >= 0 ? partials[n] : 0;
+    --n;
+    if (Math.abs(overflow) > 1 || (overflow > 0 && next > 0 || overflow < 0 && next < 0)) {
+      return overflow > 0 ? Infinity : -Infinity;
+    }
+    // here we actually have to do the arithmetic
+    // drop a factor of 2 so we can do it without overflow
+    // assert(Math.abs(overflow) === 1)
+    0, { hi, lo } = twosum(overflow * 2**1023, next / 2);
+    lo *= 2;
+    if (Math.abs(2 * hi) === Infinity) {
+      // 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 === 2**1023 && lo === -(MAX_ULP / 2) && n >= 0 && partials[n] < 0) {
+          return MAX_DOUBLE;
+        }
+        return Infinity;
+      } else {
+        if (hi === -(2**1023) && lo === (MAX_ULP / 2) && n >= 0 && partials[n] > 0) {
+          return -MAX_DOUBLE;
+        }
+        return -Infinity;
+      }
+    }
+    if (lo !== 0) {
+      partials[n + 1] = lo;
+      ++n;
+      lo = 0;
+    }
+    hi *= 2;
+  }
+
+  while (n >= 0) {
+    let x = hi;
+    let y = partials[n];
+    --n;
+    // assert: Math.abs(x) > Math.abs(y)
+    0, { hi, lo } = twosum(x, y);
+    if (lo !== 0) {
+      break;
+    }
+  }
+
+  // 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))) {
+    let y = lo * 2.0;
+    let x = hi + y;
+    let yr = x - hi;
+    if (y === yr) {
+      hi = x;
+    }
+  }
+
+  return hi;
+}
+
+let naiveSum = ar => ar.reduce((a, b) => a+b);
+let v = [0, 1, 2, 1000, 0.77, 0.4123, 1274123]
+console.log(sum(v));
+console.log(naiveSum(v));