diff --git a/src/fsum.js b/src/fsum.js index 4a8e1cf..5b20c96 100644 --- a/src/fsum.js +++ b/src/fsum.js @@ -5,15 +5,20 @@ export class Adder { this._n = 0; } add(x) { - const p = this._partials; + let p = this._partials; let i = 0; - for (let j = 0; j < this._n && j < 32; j++) { + for (let j = 0; j < this._n; j++) { const y = p[j], hi = x + y, lo = Math.abs(x) < Math.abs(y) ? x - (hi - y) : y - (hi - x); if (lo) p[i++] = lo; x = hi; } + if (i === p.length) { + this._partials = new Float64Array(p.length * 2); + this._partials.set(p); + p = this._partials; + } p[i] = x; this._n = i + 1; return this; diff --git a/test/fcumsum-test.js b/test/fcumsum-test.js index 21cd3fd..d68198f 100644 --- a/test/fcumsum-test.js +++ b/test/fcumsum-test.js @@ -9,6 +9,12 @@ it("fcumsum(array) returns a Float64Array of the expected length", () => { assert.strictEqual(R.length, A.length); }); +it("fcumsum(array) preserves widely separated terms through cancellation", () => { + const values = Array.from({length: 33}, (_, i) => 2 ** (-1074 + i * 54)); + const negatives = values.slice(1).reverse().map(x => -x); + assert.deepStrictEqual(Array.from(fcumsum([...values, ...negatives])), [...values, ...values.slice(0, -1).reverse()]); +}); + it("fcumsum(array) is an exact cumsum", () => { assert.strictEqual(lastc([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]), 1); assert.strictEqual(lastc([0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, -0.3, -0.3, -0.3, -0.3, -0.3, -0.3, -0.3, -0.3, -0.3, -0.3]), 0); diff --git a/test/fsum-test.js b/test/fsum-test.js index ed9d575..a8fdbbc 100644 --- a/test/fsum-test.js +++ b/test/fsum-test.js @@ -14,6 +14,36 @@ it("+adder can be applied several times", () => { assert.strictEqual(+adder, 1); }); +it("Adder preserves widely separated terms through cancellation", () => { + const values = Array.from({length: 33}, (_, i) => 2 ** (-1074 + i * 54)); + const adder = new Adder(); + for (const value of values) adder.add(value); + // All smaller terms together are less than half an ulp of the largest. + assert.strictEqual(+adder, values[32]); + for (let i = 32; i > 0; --i) { + adder.add(-values[i]); + assert.strictEqual(+adder, values[i - 1]); + } +}); + +it("Adder preserves partials through repeated storage growth", () => { + const adder = new Adder(); + for (let e = 0; e <= 128; e += 2) { + const x = 2 ** (e + 53); + // Cancellation leaves 2 ** e in a separate partial. + adder.add(x).add(-x + 2 ** e); + } + assert(Number.isFinite(+adder)); + for (let e = 128; e > 0; e -= 2) adder.add(-(2 ** e)); + assert.strictEqual(+adder, 1); +}); + +it("fsum(array) sums widely separated finite terms", () => { + const values = Array.from({length: 33}, (_, i) => 2 ** (-1074 + i * 54)); + assert.strictEqual(fsum(values), values[32]); + assert.strictEqual(fsum([...values, ...values.slice(1).reverse().map(x => -x)]), values[0]); +}); + it("fsum(array) is an exact sum", () => { assert.strictEqual(fsum([.1, .1, .1, .1, .1, .1, .1, .1, .1, .1]), 1); assert.strictEqual(fsum([.3, .3, .3, .3, .3, .3, .3, .3, .3, .3, -.3, -.3, -.3, -.3, -.3, -.3, -.3, -.3, -.3, -.3]), 0);