From 7a88902a674bd126b6bbc167ea783e9166c201bd Mon Sep 17 00:00:00 2001 From: marcin-kordas-hoc Date: Wed, 7 Oct 2026 09:46:01 +0000 Subject: [PATCH 01/13] Fix precision loss in VAR/STDEV on data with a large mean The variance functions computed sum(x^2) - sum(x)^2 / n in one pass. When the mean is large relative to the spread, the two terms are nearly equal and cancel, so the result was 0, #NUM!, or a negative variance. MomentsAggregate now keeps the sums of deviations and of squared deviations from a shift (the first value added) in double-double arithmetic (new src/interpreter/doubleDouble.ts) and rebases them when two aggregates compose, so the range cache keeps working. The sum of squared deviations is evaluated and divided in double-double and rounded once. Affects VAR.S, VAR.P, STDEV.S, STDEV.P, VARA, VARPA, STDEVA, STDEVPA, their legacy names, and SUBTOTAL modes 7, 8, 10, 11, 107, 108, 110 and 111. A variance that is finite in 3.3.0 stays finite for values up to about 1e300. Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 1 + src/interpreter/doubleDouble.ts | 130 ++++++++++++++++++ .../plugin/NumericAggregationPlugin.ts | 84 ++++++++++- 3 files changed, 208 insertions(+), 7 deletions(-) create mode 100644 src/interpreter/doubleDouble.ts diff --git a/CHANGELOG.md b/CHANGELOG.md index aeea8dd511..97193c678f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Fixed the `AVERAGEIF` function returning a division-by-zero error when the calculated average was `0`. [#1733](https://github.com/handsontable/hyperformula/pull/1733) - Fixed the localized names of `VSTACK` and `HSTACK` in 14 language packs to match Microsoft Excel. [#1748](https://github.com/handsontable/hyperformula/pull/1748) - Fixed the MAXPOOL and MEDIANPOOL functions throwing an uncaught `TypeError` instead of returning the `#VALUE!` error when the range dimensions are not a whole multiple of the window size and the stride. [#1718](https://github.com/handsontable/hyperformula/pull/1718) +- Fixed the variance and standard deviation functions (`VAR.S`, `VAR.P`, `STDEV.S`, `STDEV.P`, `VARA`, `VARPA`, `STDEVA`, `STDEVPA`, their legacy names, and `SUBTOTAL` modes 7, 8, 10, 11, 107, 108, 110 and 111) losing all precision on data with a large mean and a small spread, where they returned `0`, `#NUM!`, or a negative variance. They now return the variance of the stored values to full double precision. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed the `MOD` function returning a remainder with the sign of the dividend instead of the sign of the divisor, which made the results differ from Excel and Google Sheets for arguments with opposite signs (e.g. `=MOD(-3, 12)` now returns `9` instead of `-3`). [#1747](https://github.com/handsontable/hyperformula/issues/1747) ## [3.4.0] - 2026-08-10 diff --git a/src/interpreter/doubleDouble.ts b/src/interpreter/doubleDouble.ts new file mode 100644 index 0000000000..4a71174292 --- /dev/null +++ b/src/interpreter/doubleDouble.ts @@ -0,0 +1,130 @@ +/** + * @license + * Copyright (c) 2025 Handsoncode. All rights reserved. + */ + +/** + * A number represented as the unevaluated sum `hi + lo` of two doubles (double-double arithmetic). + * It carries about 106 bits of significand instead of 53, so sums and products keep the low-order + * digits that plain doubles round away. After `twoSum`, `twoProduct`, `addDoubleDouble` and + * `divideDoubleDouble`, `|lo|` is at most half an ulp of `hi`; products of double-doubles are left + * unnormalized. + * + * Used where a result is the small difference of large accumulated quantities, such as a sum of + * squared deviations. The building blocks are the error-free transformations TwoSum and TwoProduct + * (splitting a factor by 2^27 + 1), which need no fused multiply-add. + */ +export interface DoubleDouble { + readonly hi: number, + readonly lo: number, +} + +export const DOUBLE_DOUBLE_ZERO: DoubleDouble = {hi: 0, lo: 0} + +/** Splitting constant for binary64: 2^27 + 1. */ +const SPLITTER = 134217729 + +/** The largest magnitude that can be multiplied by `SPLITTER` without overflowing (2^996). */ +const SPLIT_LIMIT = 2 ** 996 + +/** + * The exact sum of two doubles, as a double-double (TwoSum). + * + * @param {number} a - first addend + * @param {number} b - second addend + * @returns {DoubleDouble} `a + b` with no rounding error + */ +export function twoSum(a: number, b: number): DoubleDouble { + const hi = a + b + const bVirtual = hi - a + return {hi, lo: (a - (hi - bVirtual)) + (b - bVirtual)} +} + +/** + * The exact product of two doubles, as a double-double (TwoProduct). + * + * Factors up to 2^996 (about 6.7e299) are split and multiplied exactly. Larger factors are multiplied + * directly, which keeps the result finite whenever the plain product is finite. + * + * @param {number} a - first factor + * @param {number} b - second factor + * @returns {DoubleDouble} `a * b`, with no rounding error unless a factor exceeds 2^996 or the product + * overflows + */ +export function twoProduct(a: number, b: number): DoubleDouble { + const hi = a * b + if (Math.abs(a) > SPLIT_LIMIT || Math.abs(b) > SPLIT_LIMIT) { + return {hi, lo: 0} + } + const aScaled = SPLITTER * a + const aHigh = aScaled - (aScaled - a) + const aLow = a - aHigh + const bScaled = SPLITTER * b + const bHigh = bScaled - (bScaled - b) + const bLow = b - bHigh + return {hi, lo: ((aHigh * bHigh - hi) + aHigh * bLow + aLow * bHigh) + aLow * bLow} +} + +/** + * Sum of two double-doubles. + * + * @param {DoubleDouble} x - first addend + * @param {DoubleDouble} y - second addend + * @returns {DoubleDouble} `x + y`; a sum that is infinite or NaN is the plain double sum + */ +export function addDoubleDouble(x: DoubleDouble, y: DoubleDouble): DoubleDouble { + const sum = twoSum(x.hi, y.hi) + if (!Number.isFinite(sum.hi)) { + return {hi: sum.hi, lo: 0} + } + const lo = sum.lo + x.lo + y.lo + const hi = sum.hi + lo + return {hi, lo: lo - (hi - sum.hi)} +} + +/** + * Product of two double-doubles. + * + * @param {DoubleDouble} x - first factor + * @param {DoubleDouble} y - second factor + * @returns {DoubleDouble} `x * y` + */ +export function multiplyDoubleDouble(x: DoubleDouble, y: DoubleDouble): DoubleDouble { + const product = twoProduct(x.hi, y.hi) + return {hi: product.hi, lo: product.lo + x.hi * y.lo + x.lo * y.hi} +} + +/** + * Product of a double-double and a double. + * + * @param {DoubleDouble} x - double-double factor + * @param {number} k - double factor + * @returns {DoubleDouble} `x * k` + */ +export function scaleDoubleDouble(x: DoubleDouble, k: number): DoubleDouble { + const product = twoProduct(x.hi, k) + return {hi: product.hi, lo: product.lo + x.lo * k} +} + +/** + * Quotient of a double-double and a double. + * + * @param {DoubleDouble} x - dividend + * @param {number} k - divisor + * @returns {DoubleDouble} `x / k` + */ +export function divideDoubleDouble(x: DoubleDouble, k: number): DoubleDouble { + const quotient = x.hi / k + const product = twoProduct(quotient, k) + return twoSum(quotient, ((x.hi - product.hi) - product.lo + x.lo) / k) +} + +/** + * Rounds a double-double to the nearest double. + * + * @param {DoubleDouble} x - the value to round + * @returns {number} `x` as a double + */ +export function roundDoubleDouble(x: DoubleDouble): number { + return x.hi + x.lo +} diff --git a/src/interpreter/plugin/NumericAggregationPlugin.ts b/src/interpreter/plugin/NumericAggregationPlugin.ts index 5a675e9360..4cc618ea46 100644 --- a/src/interpreter/plugin/NumericAggregationPlugin.ts +++ b/src/interpreter/plugin/NumericAggregationPlugin.ts @@ -11,6 +11,16 @@ import {Maybe} from '../../Maybe' import {Ast, AstNodeType, CellRangeAst, ProcedureAst} from '../../parser' import {ColumnRangeAst, RowRangeAst} from '../../parser/Ast' import {coerceBooleanToNumber} from '../ArithmeticHelper' +import { + addDoubleDouble, + divideDoubleDouble, + DOUBLE_DOUBLE_ZERO, + DoubleDouble, + multiplyDoubleDouble, + roundDoubleDouble, + scaleDoubleDouble, + twoSum, +} from '../doubleDouble' import {InterpreterState} from '../InterpreterState' import {EmptyValue, ExtendedNumber, getRawValue, InternalScalarValue, isExtendedNumber} from '../InterpreterValue' import {SimpleRangeValue} from '../../SimpleRangeValue' @@ -31,23 +41,73 @@ function zeroForInfinite(value: InternalScalarValue) { } } +/** + * Moments of a set of numbers, composable so that the value of a range can be cached and reused + * for a larger range. + * + * Besides the sum and the count, it keeps the sums of deviations and of squared deviations from a + * `shift`: the first value added. Both sums are double-double. Measured from a nearby value, the + * deviations stay small, which avoids the cancellation of the textbook one-pass form + * `sum(x^2) - sum(x)^2 / n` for data with a large mean and a small spread (for example 10000000.001, + * 10000000.002, ...). The variance is the correctly rounded value for the stored values, and the + * standard deviation, its square root, is within one unit in the last place of the exact value. + */ class MomentsAggregate { - public static empty = new MomentsAggregate(0, 0, 0) + public static empty = new MomentsAggregate(0, 0, 0, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) constructor( - public readonly sumsq: number, public readonly sum: number, public readonly count: number, + public readonly shift: number, + public readonly shiftedSum: DoubleDouble, + public readonly shiftedSumOfSquares: DoubleDouble, ) { } public static single(arg: number): MomentsAggregate { - return new MomentsAggregate(arg * arg, arg, 1) + return new MomentsAggregate(arg, 1, arg, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) } - public compose(other: MomentsAggregate) { - return new MomentsAggregate(this.sumsq + other.sumsq, this.sum + other.sum, this.count + other.count) + /** + * Combines two aggregates. The result keeps this aggregate's `shift`; the other one's shifted sums + * are re-expressed relative to it. An empty aggregate is the identity: composing with it returns the + * other aggregate unchanged, with its own shift. + */ + public compose(other: MomentsAggregate): MomentsAggregate { + if (this.count === 0) { + return other + } + if (other.count === 0) { + return this + } + + const shiftDifference = twoSum(other.shift, -this.shift) + + if (other.count === 1) { + // the common case of adding one value: its deviation is the shift difference itself + return new MomentsAggregate( + this.sum + other.sum, + this.count + 1, + this.shift, + addDoubleDouble(this.shiftedSum, shiftDifference), + addDoubleDouble(this.shiftedSumOfSquares, multiplyDoubleDouble(shiftDifference, shiftDifference)), + ) + } + + // other's sums rebased: S1 + n*d and S2 + 2*d*S1 + n*d^2 + const rebasedSum = addDoubleDouble(other.shiftedSum, scaleDoubleDouble(shiftDifference, other.count)) + const rebasedSumOfSquares = addDoubleDouble( + addDoubleDouble(other.shiftedSumOfSquares, scaleDoubleDouble(multiplyDoubleDouble(shiftDifference, other.shiftedSum), 2)), + scaleDoubleDouble(multiplyDoubleDouble(shiftDifference, shiftDifference), other.count), + ) + return new MomentsAggregate( + this.sum + other.sum, + this.count + other.count, + this.shift, + addDoubleDouble(this.shiftedSum, rebasedSum), + addDoubleDouble(this.shiftedSumOfSquares, rebasedSumOfSquares), + ) } public averageValue(): Maybe { @@ -60,7 +120,7 @@ class MomentsAggregate { public varSValue(): Maybe { if (this.count > 1) { - return (this.sumsq - (this.sum * this.sum) / this.count) / (this.count - 1) + return roundDoubleDouble(divideDoubleDouble(this.sumOfSquaredDeviations(), this.count - 1)) } else { return undefined } @@ -68,11 +128,21 @@ class MomentsAggregate { public varPValue(): Maybe { if (this.count > 0) { - return (this.sumsq - (this.sum * this.sum) / this.count) / this.count + return roundDoubleDouble(divideDoubleDouble(this.sumOfSquaredDeviations(), this.count)) } else { return undefined } } + + /** + * The sum of squared deviations from the mean, `S2 - S1^2 / n` with `S1`, `S2` the shifted sums + * and `n` the count, as a double-double: the variance divides it before rounding once. + */ + private sumOfSquaredDeviations(): DoubleDouble { + const squaredSumOverCount = divideDoubleDouble(multiplyDoubleDouble(this.shiftedSum, this.shiftedSum), this.count) + const result = addDoubleDouble(this.shiftedSumOfSquares, {hi: -squaredSumOverCount.hi, lo: -squaredSumOverCount.lo}) + return result + } } export class NumericAggregationPlugin extends FunctionPlugin implements FunctionPluginTypecheck { From bb43a62e5e1a0d2fb93d1cc992593b982e0a89ea Mon Sep 17 00:00:00 2001 From: marcin-kordas-hoc Date: Thu, 8 Oct 2026 01:22:52 +0800 Subject: [PATCH 02/13] Compute DEVSQ, COVARIANCE, SLOPE, STEYX and the D-variance functions in double-double (HF-491) (#1798) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Stacked on #1784 (base `fix/variance-precision`): it uses `doubleDouble.ts` from there. Retarget to `develop` once #1784 is merged. ### Context `DEVSQ`, `COVARIANCE.P`, `COVARIANCE.S`, `SLOPE`, `STEYX` and the database functions `DSTDEV`, `DSTDEVP`, `DVAR`, `DVARP` already subtract the mean before squaring, so they are not affected by the cancellation fixed in #1784. They still round the mean to a double before subtracting it, and `STEYX` subtracts two nearly equal sums. When the values use all 15 significant digits and their mean is large relative to their spread, this costs most of the digits: | Formula (three values near 1e11, `y = {1,2,3}`) | Before | Exact (stored doubles) | After | |---|---|---|---| | `=SLOPE({1,2,3}, {100000000000.003,100000000000.004,100000000000.005})` | `1000.37` | 1000.5301841348769 | 1000.5301841348769 | | `=STEYX(` same `)` | `0.0187` | 0.006232737795987631 | 0.006232737795987631 | | `=DEVSQ({100000000000.003,100000000000.004,100000000000.005})` | `1.99815e-06` | 1.997842142979304e-06 | 1.997842142979304e-06 | ### The change New `src/interpreter/deviationSums.ts` computes the mean, the deviations and their sums of squares and products in double-double (`doubleDouble.ts` from #1784, plus a double-double division and a rounding helper) and rounds once at the end. `DEVSQ`, `COVARIANCE.P/S`, `SLOPE` and `STEYX` call it instead of the jStat helpers, and `DSTDEV`, `DSTDEVP`, `DVAR` and `DVARP` use it instead of their own plain-double loops. A division by zero or any other non-finite result is returned as before. ### Accuracy The reference is the exact result for the stored doubles (BigInt rational arithmetic, cross-checked with Python's `statistics` and mpmath). Minimum digits of agreement over 23 datasets, including the 9 NIST StRD univariate sets: | Function | Before | After | |---|---|---| | `DEVSQ` | 3.4 | **15** | | `COVARIANCE.S` | 14.3 | **15** | | `SLOPE` | 3.4 | **15** | | `STEYX` | 0 | **15** | | `DSTDEV` | 4.1 | **15** | On ordinary data the old code was already at 11.5–15 digits; the loss shows on 15-digit inputs and on near-perfect fits. A plain-double version of the same code reproduces the "before" numbers, so the comparison is meaningful. ### Cost Medians of 5 alternating process pairs, each the median of 7 runs, building an engine and evaluating the formula (so the figures include parsing and loading the data): `DEVSQ` over 100k cells 1.08x, `SLOPE` over 100k pairs 1.12x, `STEYX` 1.14x, `COVARIANCE.S` 1.05x, 2000 `DEVSQ` formulas over 100 cells each 1.25x, `DSTDEV` over 20k rows 0.97x. On running ranges, `=DEVSQ(A$1:An)` over 4,000 rows takes 1.4x as long as before (1.8 s to 2.5 s, medians over 6 processes on a quiet machine), because every formula recomputes the whole range and the work per value is larger. `STDEV.S` on running ranges is unaffected by the first pull request (+4%, inside the noise). ### Behavior changes - `SLOPE` and `STEYX` now return `#NUM!` when all the x values are equal. Before they returned an arbitrary number (`SLOPE` 0, `STEYX` 1.414 for y = 1, 2, 3). - `STEYX` returns 0 instead of `#NUM!` for points that lie exactly on a line, where the residual rounds to a tiny negative number. - Results that were finite in 3.3.0 stay finite for values up to about 1e300. Both behaviors are covered by tests and mentioned in the changelog entry. ### How did you test your changes? - 30 new tests in the tests repo, added to the spec of each function. Expected values are exact for the stored doubles and compared as a ratio within 1e-14 with `smartRounding` off. They cover a case near 1e7 and one near 1e11 per function, covariance with both series near 1e11 (where the old mean rounding showed in the 4th digit), values near 1e300 and 1e150 that stay finite, `STEYX` on points that lie exactly on a line, and `SLOPE`/`STEYX` with all x values equal. Against develop 12 fail, against the earlier head of this PR 7 fail (the values that overflowed), and with this change all pass. - Full suite: 6291 passed, 0 failed (3 skipped, as on develop). - Lint: no new warnings (104 in the changed files against 106 on develop); the new files have none. ### Not in this PR `CORREL`, `PEARSON`, `RSQ`, `F.TEST`, `T.TEST`, `Z.TEST`, `SKEW` and `SKEW.P` still use the jStat helpers. `INTERCEPT` does not exist yet (HF-413). ### Types of changes - [ ] Breaking change (a fix or a feature because of which an existing functionality doesn't work as expected anymore) - [ ] New feature or improvement (a non-breaking change that adds functionality) - [x] Bug fix (a non-breaking change that fixes an issue) - [ ] Additional language file, or a change to an existing language file (translations) - [ ] Change to the documentation ### Related issues: 1. Stacked on #1784 2. Tests: handsontable/hyperformula-tests#74 ### Checklist: - [x] I have reviewed the guidelines about [Contributing to HyperFormula](https://hyperformula.handsontable.com/docs/guide/contributing.html) and I confirm that my code follows the code style of this project. - [x] I have signed the [Contributor License Agreement](https://goo.gl/forms/yuutGuN0RjsikVpM2). - [ ] My change is compliant with the [OpenDocument](https://docs.oasis-open.org/office/OpenDocument/v1.3/os/part4-formula/OpenDocument-v1.3-os-part4-formula.html) standard. - [ ] My change is compatible with Microsoft Excel. - [ ] My change is compatible with Google Sheets. - [x] I described my changes in the [CHANGELOG.md](https://github.com/handsontable/hyperformula/blob/master/CHANGELOG.md) file. - [ ] My changes require a documentation update. - [ ] My changes require a migration guide. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --- > [!NOTE] > **Medium Risk** > Changes numeric results for affected formulas (intentional bug fix) and touches shared double-double math used by multiple statistical and database functions; edge-case error behavior for SLOPE/STEYX also changes. > > **Overview** > Fixes **numerical precision** for deviation-based statistics when inputs have a **large mean and tiny spread** (e.g. values near `1e11` differing only in the last digits). > > Adds **`deviationSums.ts`** to compute means, squared deviations, and covariance-style products in **double-double** (rounding once at the end), plus **`divideDoubleDoubles`** in `doubleDouble.ts`. **`DEVSQ`**, **`COVARIANCE.P` / `.S`** (and aliases), **`SLOPE`**, **`STEYX`**, and database **`DSTDEV`**, **`DSTDEVP`**, **`DVAR`**, **`DVARP`** now use this path instead of jStat / plain-double loops. > > **`STEYX`** no longer yields **`#NUM!`** for a perfect linear fit; **`SLOPE`** and **`STEYX`** return **`#NUM!`** when all **x** values are equal. **CHANGELOG** updated. > > Reviewed by [Cursor Bugbot](https://cursor.com/bugbot) for commit 54904ea2136f5f0580448978da83e1461d0a3856. Bugbot is set up for automated code reviews on this repo. Configure [here](https://www.cursor.com/dashboard/bugbot). Co-authored-by: Claude Sonnet 5.5 --- CHANGELOG.md | 1 + src/interpreter/deviationSums.ts | 76 +++++++++++++++++++ src/interpreter/doubleDouble.ts | 19 +++++ src/interpreter/plugin/DatabasePlugin.ts | 16 ++-- .../plugin/StatisticalAggregationPlugin.ts | 24 ++++-- 5 files changed, 120 insertions(+), 16 deletions(-) create mode 100644 src/interpreter/deviationSums.ts diff --git a/CHANGELOG.md b/CHANGELOG.md index 97193c678f..22af2493e4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Fixed the localized names of `VSTACK` and `HSTACK` in 14 language packs to match Microsoft Excel. [#1748](https://github.com/handsontable/hyperformula/pull/1748) - Fixed the MAXPOOL and MEDIANPOOL functions throwing an uncaught `TypeError` instead of returning the `#VALUE!` error when the range dimensions are not a whole multiple of the window size and the stride. [#1718](https://github.com/handsontable/hyperformula/pull/1718) - Fixed the variance and standard deviation functions (`VAR.S`, `VAR.P`, `STDEV.S`, `STDEV.P`, `VARA`, `VARPA`, `STDEVA`, `STDEVPA`, their legacy names, and `SUBTOTAL` modes 7, 8, 10, 11, 107, 108, 110 and 111) losing all precision on data with a large mean and a small spread, where they returned `0`, `#NUM!`, or a negative variance. They now return the variance of the stored values to full double precision. [#1784](https://github.com/handsontable/hyperformula/pull/1784) +- Fixed `DEVSQ`, `COVARIANCE.P`, `COVARIANCE.S` (and their aliases `COVAR`, `COVARIANCEP` and `COVARIANCES`), `SLOPE`, `STEYX`, `DSTDEV`, `DSTDEVP`, `DVAR` and `DVARP` losing digits when the values use all 15 significant digits and their mean is large relative to their spread (for example, `STEYX` returned `0.0187` instead of `0.00623` and `SLOPE` returned `1000.37` instead of `1000.53` for three values near `100000000000`). They now return results accurate to the last digits of the stored values. `STEYX` no longer returns `#NUM!` for points that lie exactly on a line, and `SLOPE` and `STEYX` return `#NUM!` instead of an arbitrary number when all the x values are equal. [#1798](https://github.com/handsontable/hyperformula/pull/1798) - Fixed the `MOD` function returning a remainder with the sign of the dividend instead of the sign of the divisor, which made the results differ from Excel and Google Sheets for arguments with opposite signs (e.g. `=MOD(-3, 12)` now returns `9` instead of `-3`). [#1747](https://github.com/handsontable/hyperformula/issues/1747) ## [3.4.0] - 2026-08-10 diff --git a/src/interpreter/deviationSums.ts b/src/interpreter/deviationSums.ts new file mode 100644 index 0000000000..37e4a227d2 --- /dev/null +++ b/src/interpreter/deviationSums.ts @@ -0,0 +1,76 @@ +/** + * @license + * Copyright (c) 2025 Handsoncode. All rights reserved. + */ + +import { + addDoubleDouble, + divideDoubleDouble, + DOUBLE_DOUBLE_ZERO, + DoubleDouble, + multiplyDoubleDouble, +} from './doubleDouble' + +/** + * Sums of squared deviations and of products of deviations from the mean, for the functions built + * on them (DEVSQ, COVARIANCE, SLOPE, STEYX, DSTDEV, DVAR and their variants). + * + * The mean, the deviations and their sums are all kept in double-double and rounded once at the end, + * so the result is accurate to the last digits of the stored values even when the mean is large + * relative to the spread. + */ + +/** + * The mean of the values. + * + * @param {number[]} values - a non-empty array of numbers + * @returns {DoubleDouble} the mean, without rounding to a double + */ +function mean(values: number[]): DoubleDouble { + const total = values.reduce((sum, value) => addDoubleDouble(sum, {hi: value, lo: 0}), DOUBLE_DOUBLE_ZERO) + return divideDoubleDouble(total, values.length) +} + +/** + * The deviations of the values from their mean. + * + * @param {number[]} values - a non-empty array of numbers + * @returns {DoubleDouble[]} `value - mean` for each value, without rounding + */ +function deviationsFromMean(values: number[]): DoubleDouble[] { + const center = mean(values) + return values.map((value) => addDoubleDouble({hi: value, lo: 0}, {hi: -center.hi, lo: -center.lo})) +} + +/** + * The sum of the squared deviations of the values from their mean. + * + * @param {number[]} values - a non-empty array of numbers + * @returns {DoubleDouble} `sum((x - mean)^2)` + */ +export function sumOfSquaredDeviations(values: number[]): DoubleDouble { + const deviations = deviationsFromMean(values) + return sumOfProducts(deviations, deviations) +} + +/** + * The sum of the products of the paired deviations of two arrays from their means. + * + * @param {number[]} first - a non-empty array of numbers + * @param {number[]} second - an array of the same length as `first` + * @returns {DoubleDouble} `sum((x - mean(x)) * (y - mean(y)))` + */ +export function sumOfProductsOfDeviations(first: number[], second: number[]): DoubleDouble { + return sumOfProducts(deviationsFromMean(first), deviationsFromMean(second)) +} + +/** + * The sum of the products of paired double-doubles. + * + * @param {DoubleDouble[]} first - the first factors + * @param {DoubleDouble[]} second - the second factors, of the same length + * @returns {DoubleDouble} `sum(first[i] * second[i])` + */ +function sumOfProducts(first: DoubleDouble[], second: DoubleDouble[]): DoubleDouble { + return first.reduce((sum, deviation, index) => addDoubleDouble(sum, multiplyDoubleDouble(deviation, second[index])), DOUBLE_DOUBLE_ZERO) +} diff --git a/src/interpreter/doubleDouble.ts b/src/interpreter/doubleDouble.ts index 4a71174292..95fd34ff23 100644 --- a/src/interpreter/doubleDouble.ts +++ b/src/interpreter/doubleDouble.ts @@ -119,6 +119,25 @@ export function divideDoubleDouble(x: DoubleDouble, k: number): DoubleDouble { return twoSum(quotient, ((x.hi - product.hi) - product.lo + x.lo) / k) } +/** + * Quotient of two double-doubles. + * + * A zero or non-finite divisor, or a non-finite quotient, gives the plain double quotient, so infinities, + * zeros and NaNs are the same as in ordinary arithmetic. + * + * @param {DoubleDouble} x - dividend + * @param {DoubleDouble} y - divisor + * @returns {DoubleDouble} `x / y` + */ +export function divideDoubleDoubles(x: DoubleDouble, y: DoubleDouble): DoubleDouble { + const quotient = x.hi / y.hi + if (!Number.isFinite(quotient) || !Number.isFinite(y.hi)) { + return {hi: quotient, lo: 0} + } + const remainder = addDoubleDouble(x, scaleDoubleDouble({hi: -y.hi, lo: -y.lo}, quotient)) + return twoSum(quotient, remainder.hi / y.hi) +} + /** * Rounds a double-double to the nearest double. * diff --git a/src/interpreter/plugin/DatabasePlugin.ts b/src/interpreter/plugin/DatabasePlugin.ts index 1208a0f17f..d5c2150ae3 100644 --- a/src/interpreter/plugin/DatabasePlugin.ts +++ b/src/interpreter/plugin/DatabasePlugin.ts @@ -6,6 +6,8 @@ import {CellError, ErrorType} from '../../Cell' import {ErrorMessage} from '../../error-message' import {ProcedureAst} from '../../parser' +import {sumOfSquaredDeviations} from '../deviationSums' +import {divideDoubleDouble, roundDoubleDouble} from '../doubleDouble' import {InterpreterState} from '../InterpreterState' import {EmptyValue, getRawValue, InternalScalarValue, InterpreterValue, isExtendedNumber, RawScalarValue} from '../InterpreterValue' import {SimpleRangeValue} from '../../SimpleRangeValue' @@ -341,9 +343,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return new CellError(ErrorType.DIV_BY_ZERO) } - const mean = values.reduce((a, b) => a + b, 0) / values.length - const variance = values.reduce((sum, v) => sum + (v - mean) ** 2, 0) / (values.length - 1) - return Math.sqrt(variance) + return Math.sqrt(roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations(values), values.length - 1))) }) } @@ -366,9 +366,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return new CellError(ErrorType.DIV_BY_ZERO) } - const mean = values.reduce((a, b) => a + b, 0) / values.length - const variance = values.reduce((sum, v) => sum + (v - mean) ** 2, 0) / values.length - return Math.sqrt(variance) + return Math.sqrt(roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations(values), values.length))) }) } @@ -390,8 +388,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return new CellError(ErrorType.DIV_BY_ZERO) } - const mean = values.reduce((a, b) => a + b, 0) / values.length - return values.reduce((sum, v) => sum + (v - mean) ** 2, 0) / (values.length - 1) + return roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations(values), values.length - 1)) }) } @@ -414,8 +411,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return new CellError(ErrorType.DIV_BY_ZERO) } - const mean = values.reduce((a, b) => a + b, 0) / values.length - return values.reduce((sum, v) => sum + (v - mean) ** 2, 0) / values.length + return roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations(values), values.length)) }) } diff --git a/src/interpreter/plugin/StatisticalAggregationPlugin.ts b/src/interpreter/plugin/StatisticalAggregationPlugin.ts index 75b8c820b7..73b79d776d 100644 --- a/src/interpreter/plugin/StatisticalAggregationPlugin.ts +++ b/src/interpreter/plugin/StatisticalAggregationPlugin.ts @@ -19,7 +19,6 @@ import { centralF, chisquare, corrcoeff, - covariance, geomean, mean, normal, @@ -28,6 +27,14 @@ import { sumsqerr, variance } from './3rdparty/jstat/jstat' +import {sumOfProductsOfDeviations, sumOfSquaredDeviations} from '../deviationSums' +import { + addDoubleDouble, + divideDoubleDouble, + divideDoubleDoubles, + multiplyDoubleDouble, + roundDoubleDouble, +} from '../doubleDouble' import {FunctionArgumentType, FunctionPlugin, FunctionPluginTypecheck, ImplementedFunctions} from './FunctionPlugin' export class StatisticalAggregationPlugin extends FunctionPlugin implements FunctionPluginTypecheck { @@ -185,7 +192,7 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (coerced.length === 0) { return 0 } - return sumsqerr(coerced) + return roundDoubleDouble(sumOfSquaredDeviations(coerced)) }) } @@ -282,7 +289,7 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n === 1) { return 0 } - return covariance(ret[0], ret[1]) * (n - 1) / n + return roundDoubleDouble(divideDoubleDouble(sumOfProductsOfDeviations(ret[0], ret[1]), n)) }) } @@ -301,7 +308,7 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 1) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.TwoValues) } - return covariance(ret[0], ret[1]) + return roundDoubleDouble(divideDoubleDouble(sumOfProductsOfDeviations(ret[0], ret[1]), n - 1)) }) } @@ -370,7 +377,12 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 2) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.ThreeValues) } - return Math.sqrt((sumsqerr(ret[0]) - Math.pow(covariance(ret[0], ret[1]) * (n - 1), 2) / sumsqerr(ret[1])) / (n - 2)) + const sumOfProducts = sumOfProductsOfDeviations(ret[0], ret[1]) + const slope = divideDoubleDoubles(sumOfProducts, sumOfSquaredDeviations(ret[1])) + const explained = multiplyDoubleDouble(slope, sumOfProducts) + const residual = addDoubleDouble(sumOfSquaredDeviations(ret[0]), {hi: -explained.hi, lo: -explained.lo}) + // the residual is non-negative; a tiny negative rounding remainder counts as 0 + return Math.sqrt(Math.max(0, roundDoubleDouble(residual)) / (n - 2)) }) } @@ -389,7 +401,7 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 1) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.TwoValues) } - return covariance(ret[0], ret[1]) * (n - 1) / sumsqerr(ret[1]) + return roundDoubleDouble(divideDoubleDoubles(sumOfProductsOfDeviations(ret[0], ret[1]), sumOfSquaredDeviations(ret[1]))) }) } From f42263e7610a0bfcadeb2968bd61378631207911 Mon Sep 17 00:00:00 2001 From: Kuba Sekowski Date: Wed, 7 Oct 2026 20:29:55 +0200 Subject: [PATCH 03/13] Fix overflow regressions in the double-double variance code (HF-491) Values of about 1e153 and above made VAR, STDEV, COVARIANCE, SLOPE, STEYX and the D-variance functions return #NUM! where develop returned a finite result. - twoProduct: scale the larger factor by 2^-53 when a factor exceeds 2^996 or the product exceeds 2^1023, so the Veltkamp split cannot overflow; the error term is now exact for every finite product above 2^-969. - twoSum: take the error from the larger-magnitude addend (Fast2Sum), which cannot overflow while the sum is finite. - divideDoubleDouble, divideDoubleDoubles: scale dividends above 2^1000 down by 2^64 so the quotient times the divisor cannot round to infinity. - MomentsAggregate: evaluate S1^2 / n as S1 * (S1 / n), and keep the plain sum of squares so the variance falls back to the previous one-pass form when the shifted sums overflow. - deviationSums mean(): sum scaled values when the running total overflows. - SLOPE, STEYX: return #NUM! when the sum of squared x deviations overflows, instead of a slope of 0. Also: - AVERAGE, AVERAGEA and SUBTOTAL 1/101 reduce a plain sum and count under their own cache keys instead of MomentsAggregate. - Shared variance/covariance and subtractDoubleDouble helpers; SLOPE and STEYX compute each array's deviations once. - Correct the accuracy claims in the JSDoc and the CHANGELOG. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 4 +- src/interpreter/deviationSums.ts | 54 ++++++++++-- src/interpreter/doubleDouble.ts | 88 +++++++++++++++++-- src/interpreter/plugin/DatabasePlugin.ts | 11 ++- .../plugin/NumericAggregationPlugin.ts | 86 +++++++++++------- .../plugin/StatisticalAggregationPlugin.ts | 30 ++++--- 6 files changed, 210 insertions(+), 63 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 22af2493e4..2a674b9ac3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,8 +12,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Fixed the `AVERAGEIF` function returning a division-by-zero error when the calculated average was `0`. [#1733](https://github.com/handsontable/hyperformula/pull/1733) - Fixed the localized names of `VSTACK` and `HSTACK` in 14 language packs to match Microsoft Excel. [#1748](https://github.com/handsontable/hyperformula/pull/1748) - Fixed the MAXPOOL and MEDIANPOOL functions throwing an uncaught `TypeError` instead of returning the `#VALUE!` error when the range dimensions are not a whole multiple of the window size and the stride. [#1718](https://github.com/handsontable/hyperformula/pull/1718) -- Fixed the variance and standard deviation functions (`VAR.S`, `VAR.P`, `STDEV.S`, `STDEV.P`, `VARA`, `VARPA`, `STDEVA`, `STDEVPA`, their legacy names, and `SUBTOTAL` modes 7, 8, 10, 11, 107, 108, 110 and 111) losing all precision on data with a large mean and a small spread, where they returned `0`, `#NUM!`, or a negative variance. They now return the variance of the stored values to full double precision. [#1784](https://github.com/handsontable/hyperformula/pull/1784) -- Fixed `DEVSQ`, `COVARIANCE.P`, `COVARIANCE.S` (and their aliases `COVAR`, `COVARIANCEP` and `COVARIANCES`), `SLOPE`, `STEYX`, `DSTDEV`, `DSTDEVP`, `DVAR` and `DVARP` losing digits when the values use all 15 significant digits and their mean is large relative to their spread (for example, `STEYX` returned `0.0187` instead of `0.00623` and `SLOPE` returned `1000.37` instead of `1000.53` for three values near `100000000000`). They now return results accurate to the last digits of the stored values. `STEYX` no longer returns `#NUM!` for points that lie exactly on a line, and `SLOPE` and `STEYX` return `#NUM!` instead of an arbitrary number when all the x values are equal. [#1798](https://github.com/handsontable/hyperformula/pull/1798) +- Fixed the variance and standard deviation functions (`VAR.S`, `VAR.P`, `STDEV.S`, `STDEV.P`, `VARA`, `VARPA`, `STDEVA`, `STDEVPA`, their legacy names, and `SUBTOTAL` modes 7, 8, 10, 11, 107, 108, 110 and 111) losing all precision on data with a large mean and a small spread, where they returned `0`, `#NUM!`, or a negative variance. They now return the variance of the stored values accurate to about the last digit. [#1784](https://github.com/handsontable/hyperformula/pull/1784) +- Fixed `DEVSQ`, `COVARIANCE.P`, `COVARIANCE.S` (and their aliases `COVAR`, `COVARIANCEP` and `COVARIANCES`), `SLOPE`, `STEYX`, `DSTDEV`, `DSTDEVP`, `DVAR` and `DVARP` losing digits when the values use all 15 significant digits and their mean is large relative to their spread (for example, `STEYX` returned `0.0187` instead of `0.00623` and `SLOPE` returned `1000.37` instead of `1000.53` for three values near `100000000000`). They now return results accurate to the last digits of the stored values. `STEYX` no longer returns `#NUM!` for points that lie exactly on a line, and `SLOPE` and `STEYX` return `#NUM!` instead of an arbitrary number when all the x values are equal or the sum of their squared deviations overflows. [#1798](https://github.com/handsontable/hyperformula/pull/1798) - Fixed the `MOD` function returning a remainder with the sign of the dividend instead of the sign of the divisor, which made the results differ from Excel and Google Sheets for arguments with opposite signs (e.g. `=MOD(-3, 12)` now returns `9` instead of `-3`). [#1747](https://github.com/handsontable/hyperformula/issues/1747) ## [3.4.0] - 2026-08-10 diff --git a/src/interpreter/deviationSums.ts b/src/interpreter/deviationSums.ts index 37e4a227d2..6bf9a0c8a5 100644 --- a/src/interpreter/deviationSums.ts +++ b/src/interpreter/deviationSums.ts @@ -8,7 +8,10 @@ import { divideDoubleDouble, DOUBLE_DOUBLE_ZERO, DoubleDouble, + multiplyByPowerOfTwo, multiplyDoubleDouble, + roundDoubleDouble, + subtractDoubleDouble, } from './doubleDouble' /** @@ -23,12 +26,30 @@ import { /** * The mean of the values. * + * When the running total overflows, the values are summed scaled down by a power of two no smaller + * than their count, so that no partial sum can exceed the largest double, and the mean is scaled back. + * * @param {number[]} values - a non-empty array of numbers * @returns {DoubleDouble} the mean, without rounding to a double */ function mean(values: number[]): DoubleDouble { - const total = values.reduce((sum, value) => addDoubleDouble(sum, {hi: value, lo: 0}), DOUBLE_DOUBLE_ZERO) - return divideDoubleDouble(total, values.length) + const total = sum(values) + if (Number.isFinite(total.hi)) { + return divideDoubleDouble(total, values.length) + } + const scale = 2 ** Math.ceil(Math.log2(values.length)) + const scaledTotal = sum(values.map((value) => value / scale)) + return multiplyByPowerOfTwo(divideDoubleDouble(scaledTotal, values.length), scale) +} + +/** + * The sum of the values. + * + * @param {number[]} values - an array of numbers + * @returns {DoubleDouble} the sum, without rounding to a double + */ +function sum(values: number[]): DoubleDouble { + return values.reduce((total, value) => addDoubleDouble(total, {hi: value, lo: 0}), DOUBLE_DOUBLE_ZERO) } /** @@ -37,9 +58,9 @@ function mean(values: number[]): DoubleDouble { * @param {number[]} values - a non-empty array of numbers * @returns {DoubleDouble[]} `value - mean` for each value, without rounding */ -function deviationsFromMean(values: number[]): DoubleDouble[] { +export function deviationsFromMean(values: number[]): DoubleDouble[] { const center = mean(values) - return values.map((value) => addDoubleDouble({hi: value, lo: 0}, {hi: -center.hi, lo: -center.lo})) + return values.map((value) => subtractDoubleDouble({hi: value, lo: 0}, center)) } /** @@ -64,6 +85,29 @@ export function sumOfProductsOfDeviations(first: number[], second: number[]): Do return sumOfProducts(deviationsFromMean(first), deviationsFromMean(second)) } +/** + * The variance of the values: the sum of squared deviations divided by `n - deltaDegreesOfFreedom`. + * + * @param {number[]} values - a non-empty array of numbers + * @param {number} deltaDegreesOfFreedom - 0 for the population variance, 1 for the sample variance + * @returns {number} the variance, rounded once to a double + */ +export function variance(values: number[], deltaDegreesOfFreedom: number): number { + return roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations(values), values.length - deltaDegreesOfFreedom)) +} + +/** + * The covariance of two arrays: the sum of products of deviations divided by `n - deltaDegreesOfFreedom`. + * + * @param {number[]} first - a non-empty array of numbers + * @param {number[]} second - an array of the same length as `first` + * @param {number} deltaDegreesOfFreedom - 0 for the population covariance, 1 for the sample covariance + * @returns {number} the covariance, rounded once to a double + */ +export function covariance(first: number[], second: number[], deltaDegreesOfFreedom: number): number { + return roundDoubleDouble(divideDoubleDouble(sumOfProductsOfDeviations(first, second), first.length - deltaDegreesOfFreedom)) +} + /** * The sum of the products of paired double-doubles. * @@ -71,6 +115,6 @@ export function sumOfProductsOfDeviations(first: number[], second: number[]): Do * @param {DoubleDouble[]} second - the second factors, of the same length * @returns {DoubleDouble} `sum(first[i] * second[i])` */ -function sumOfProducts(first: DoubleDouble[], second: DoubleDouble[]): DoubleDouble { +export function sumOfProducts(first: DoubleDouble[], second: DoubleDouble[]): DoubleDouble { return first.reduce((sum, deviation, index) => addDoubleDouble(sum, multiplyDoubleDouble(deviation, second[index])), DOUBLE_DOUBLE_ZERO) } diff --git a/src/interpreter/doubleDouble.ts b/src/interpreter/doubleDouble.ts index 95fd34ff23..e7e5cf0555 100644 --- a/src/interpreter/doubleDouble.ts +++ b/src/interpreter/doubleDouble.ts @@ -27,42 +27,84 @@ const SPLITTER = 134217729 /** The largest magnitude that can be multiplied by `SPLITTER` without overflowing (2^996). */ const SPLIT_LIMIT = 2 ** 996 +/** + * The largest product magnitude for which the product of the split high halves (at most `|a * b|` + * times (1 + 2^-26)^2) cannot overflow (2^1023). + */ +const PRODUCT_LIMIT = 2 ** 1023 + +/** + * An exact power-of-two scale (2^53) for the larger factor of a product above `PRODUCT_LIMIT` or with + * a factor above `SPLIT_LIMIT`, so that splitting it cannot overflow. + */ +const FACTOR_SCALE = 2 ** 53 + +/** + * The largest dividend magnitude that `divideDoubleDouble` and `divideDoubleDoubles` divide directly + * (2^1000). Above it, the quotient times the divisor can round to infinity, so the dividend is scaled + * down by `DIVIDEND_SCALE` first. + */ +const DIVIDEND_LIMIT = 2 ** 1000 + +/** An exact power-of-two scale for dividends above `DIVIDEND_LIMIT` (2^64). */ +const DIVIDEND_SCALE = 2 ** 64 + /** * The exact sum of two doubles, as a double-double (TwoSum). * + * The rounding error is computed from the addend of larger magnitude (Fast2Sum), which is exact and, + * unlike the branch-free form, cannot overflow while the sum is finite. + * * @param {number} a - first addend * @param {number} b - second addend * @returns {DoubleDouble} `a + b` with no rounding error */ export function twoSum(a: number, b: number): DoubleDouble { const hi = a + b - const bVirtual = hi - a - return {hi, lo: (a - (hi - bVirtual)) + (b - bVirtual)} + return {hi, lo: Math.abs(a) >= Math.abs(b) ? b - (hi - a) : a - (hi - b)} } /** * The exact product of two doubles, as a double-double (TwoProduct). * - * Factors up to 2^996 (about 6.7e299) are split and multiplied exactly. Larger factors are multiplied - * directly, which keeps the result finite whenever the plain product is finite. + * When a factor exceeds 2^996 or the product exceeds 2^1023, the larger factor is scaled by 2^-53 + * (exactly) so the splitting cannot overflow, and the error term is scaled back. * * @param {number} a - first factor * @param {number} b - second factor - * @returns {DoubleDouble} `a * b`, with no rounding error unless a factor exceeds 2^996 or the product - * overflows + * @returns {DoubleDouble} `a * b`, with no rounding error whenever the product is finite and at least + * 2^-969 in magnitude; an infinite or NaN product has `lo` 0 */ export function twoProduct(a: number, b: number): DoubleDouble { const hi = a * b - if (Math.abs(a) > SPLIT_LIMIT || Math.abs(b) > SPLIT_LIMIT) { + if (Math.abs(a) <= SPLIT_LIMIT && Math.abs(b) <= SPLIT_LIMIT && Math.abs(hi) <= PRODUCT_LIMIT) { + return {hi, lo: productError(a, b, hi)} + } + if (!Number.isFinite(hi)) { return {hi, lo: 0} } + const [larger, smaller] = Math.abs(a) >= Math.abs(b) ? [a, b] : [b, a] + const largerScaled = larger / FACTOR_SCALE + return {hi, lo: productError(largerScaled, smaller, largerScaled * smaller) * FACTOR_SCALE} +} + +/** + * The rounding error of `a * b`, where `hi` is `a * b` rounded (Dekker's product with Veltkamp + * splitting). Exact when both factors are at most `SPLIT_LIMIT` and `|hi|` is at most `PRODUCT_LIMIT`. + * + * @param {number} a - first factor + * @param {number} b - second factor + * @param {number} hi - `a * b` rounded to a double + * @returns {number} `a * b - hi` + */ +function productError(a: number, b: number, hi: number): number { const aScaled = SPLITTER * a const aHigh = aScaled - (aScaled - a) const aLow = a - aHigh const bScaled = SPLITTER * b const bHigh = bScaled - (bScaled - b) const bLow = b - bHigh - return {hi, lo: ((aHigh * bHigh - hi) + aHigh * bLow + aLow * bHigh) + aLow * bLow} + return ((aHigh * bHigh - hi) + aHigh * bLow + aLow * bHigh) + aLow * bLow } /** @@ -82,6 +124,17 @@ export function addDoubleDouble(x: DoubleDouble, y: DoubleDouble): DoubleDouble return {hi, lo: lo - (hi - sum.hi)} } +/** + * Difference of two double-doubles. + * + * @param {DoubleDouble} x - minuend + * @param {DoubleDouble} y - subtrahend + * @returns {DoubleDouble} `x - y`; a difference that is infinite or NaN is the plain double difference + */ +export function subtractDoubleDouble(x: DoubleDouble, y: DoubleDouble): DoubleDouble { + return addDoubleDouble(x, {hi: -y.hi, lo: -y.lo}) +} + /** * Product of two double-doubles. * @@ -114,6 +167,9 @@ export function scaleDoubleDouble(x: DoubleDouble, k: number): DoubleDouble { * @returns {DoubleDouble} `x / k` */ export function divideDoubleDouble(x: DoubleDouble, k: number): DoubleDouble { + if (Math.abs(x.hi) > DIVIDEND_LIMIT && Number.isFinite(x.hi)) { + return multiplyByPowerOfTwo(divideDoubleDouble(multiplyByPowerOfTwo(x, 1 / DIVIDEND_SCALE), k), DIVIDEND_SCALE) + } const quotient = x.hi / k const product = twoProduct(quotient, k) return twoSum(quotient, ((x.hi - product.hi) - product.lo + x.lo) / k) @@ -134,10 +190,24 @@ export function divideDoubleDoubles(x: DoubleDouble, y: DoubleDouble): DoubleDou if (!Number.isFinite(quotient) || !Number.isFinite(y.hi)) { return {hi: quotient, lo: 0} } - const remainder = addDoubleDouble(x, scaleDoubleDouble({hi: -y.hi, lo: -y.lo}, quotient)) + if (Math.abs(x.hi) > DIVIDEND_LIMIT) { + return multiplyByPowerOfTwo(divideDoubleDoubles(multiplyByPowerOfTwo(x, 1 / DIVIDEND_SCALE), y), DIVIDEND_SCALE) + } + const remainder = subtractDoubleDouble(x, scaleDoubleDouble(y, quotient)) return twoSum(quotient, remainder.hi / y.hi) } +/** + * Product of a double-double and a power of two, which is exact unless it overflows or underflows. + * + * @param {DoubleDouble} x - double-double factor + * @param {number} powerOfTwo - a power of two + * @returns {DoubleDouble} `x * powerOfTwo` + */ +export function multiplyByPowerOfTwo(x: DoubleDouble, powerOfTwo: number): DoubleDouble { + return {hi: x.hi * powerOfTwo, lo: x.lo * powerOfTwo} +} + /** * Rounds a double-double to the nearest double. * diff --git a/src/interpreter/plugin/DatabasePlugin.ts b/src/interpreter/plugin/DatabasePlugin.ts index d5c2150ae3..462266d6fc 100644 --- a/src/interpreter/plugin/DatabasePlugin.ts +++ b/src/interpreter/plugin/DatabasePlugin.ts @@ -6,8 +6,7 @@ import {CellError, ErrorType} from '../../Cell' import {ErrorMessage} from '../../error-message' import {ProcedureAst} from '../../parser' -import {sumOfSquaredDeviations} from '../deviationSums' -import {divideDoubleDouble, roundDoubleDouble} from '../doubleDouble' +import {variance} from '../deviationSums' import {InterpreterState} from '../InterpreterState' import {EmptyValue, getRawValue, InternalScalarValue, InterpreterValue, isExtendedNumber, RawScalarValue} from '../InterpreterValue' import {SimpleRangeValue} from '../../SimpleRangeValue' @@ -343,7 +342,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return new CellError(ErrorType.DIV_BY_ZERO) } - return Math.sqrt(roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations(values), values.length - 1))) + return Math.sqrt(variance(values, 1)) }) } @@ -366,7 +365,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return new CellError(ErrorType.DIV_BY_ZERO) } - return Math.sqrt(roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations(values), values.length))) + return Math.sqrt(variance(values, 0)) }) } @@ -388,7 +387,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return new CellError(ErrorType.DIV_BY_ZERO) } - return roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations(values), values.length - 1)) + return variance(values, 1) }) } @@ -411,7 +410,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return new CellError(ErrorType.DIV_BY_ZERO) } - return roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations(values), values.length)) + return variance(values, 0) }) } diff --git a/src/interpreter/plugin/NumericAggregationPlugin.ts b/src/interpreter/plugin/NumericAggregationPlugin.ts index 4cc618ea46..8e76a7ace6 100644 --- a/src/interpreter/plugin/NumericAggregationPlugin.ts +++ b/src/interpreter/plugin/NumericAggregationPlugin.ts @@ -19,6 +19,7 @@ import { multiplyDoubleDouble, roundDoubleDouble, scaleDoubleDouble, + subtractDoubleDouble, twoSum, } from '../doubleDouble' import {InterpreterState} from '../InterpreterState' @@ -49,15 +50,23 @@ function zeroForInfinite(value: InternalScalarValue) { * `shift`: the first value added. Both sums are double-double. Measured from a nearby value, the * deviations stay small, which avoids the cancellation of the textbook one-pass form * `sum(x^2) - sum(x)^2 / n` for data with a large mean and a small spread (for example 10000000.001, - * 10000000.002, ...). The variance is the correctly rounded value for the stored values, and the - * standard deviation, its square root, is within one unit in the last place of the exact value. + * 10000000.002, ...). The variance and the standard deviation are accurate to about one unit in the + * last place of the exact values for the stored numbers. Double-double arithmetic does not guarantee + * correct rounding. + * + * It also keeps the plain sum and sum of squares. When the deviations from the shift are so large + * (about 1e154 and above) that the shifted sums overflow, the variance falls back to the one-pass form + * used before, with the same accuracy as before: it loses about `log10(2n)` digits to cancellation. + * It is `#NUM!` when the plain sum of squares overflows too, as before, although the shifted sums of + * the same values in another order can be finite. */ class MomentsAggregate { - public static empty = new MomentsAggregate(0, 0, 0, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) + public static empty = new MomentsAggregate(0, 0, 0, 0, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) constructor( public readonly sum: number, + public readonly sumOfSquares: number, public readonly count: number, public readonly shift: number, public readonly shiftedSum: DoubleDouble, @@ -65,8 +74,14 @@ class MomentsAggregate { ) { } + /** + * The moments of one value, which is also the shift. + * + * @param {number} arg - the value + * @returns {MomentsAggregate} an aggregate of the single value + */ public static single(arg: number): MomentsAggregate { - return new MomentsAggregate(arg, 1, arg, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) + return new MomentsAggregate(arg, arg * arg, 1, arg, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) } /** @@ -88,6 +103,7 @@ class MomentsAggregate { // the common case of adding one value: its deviation is the shift difference itself return new MomentsAggregate( this.sum + other.sum, + this.sumOfSquares + other.sumOfSquares, this.count + 1, this.shift, addDoubleDouble(this.shiftedSum, shiftDifference), @@ -103,6 +119,7 @@ class MomentsAggregate { ) return new MomentsAggregate( this.sum + other.sum, + this.sumOfSquares + other.sumOfSquares, this.count + other.count, this.shift, addDoubleDouble(this.shiftedSum, rebasedSum), @@ -110,14 +127,6 @@ class MomentsAggregate { ) } - public averageValue(): Maybe { - if (this.count > 0) { - return this.sum / this.count - } else { - return undefined - } - } - public varSValue(): Maybe { if (this.count > 1) { return roundDoubleDouble(divideDoubleDouble(this.sumOfSquaredDeviations(), this.count - 1)) @@ -137,11 +146,20 @@ class MomentsAggregate { /** * The sum of squared deviations from the mean, `S2 - S1^2 / n` with `S1`, `S2` the shifted sums * and `n` the count, as a double-double: the variance divides it before rounding once. + * + * `S1^2 / n` is evaluated as `S1 * (S1 / n)`, so it cannot overflow while `S2` is finite + * (`S1^2 / n <= S2`), even when `S1^2` alone would exceed the largest double. When the shifted sums + * overflow, it falls back to `sum(x^2) - sum(x)^2 / n` in plain doubles, evaluated as before unless + * `sum(x)^2` overflows. */ private sumOfSquaredDeviations(): DoubleDouble { - const squaredSumOverCount = divideDoubleDouble(multiplyDoubleDouble(this.shiftedSum, this.shiftedSum), this.count) - const result = addDoubleDouble(this.shiftedSumOfSquares, {hi: -squaredSumOverCount.hi, lo: -squaredSumOverCount.lo}) - return result + const squaredSumOverCount = multiplyDoubleDouble(this.shiftedSum, divideDoubleDouble(this.shiftedSum, this.count)) + const result = subtractDoubleDouble(this.shiftedSumOfSquares, squaredSumOverCount) + if (Number.isFinite(roundDoubleDouble(result))) { + return result + } + const onePass = this.sumOfSquares - this.sum * this.sum / this.count + return {hi: Number.isFinite(onePass) ? onePass : this.sumOfSquares - this.sum * (this.sum / this.count), lo: 0} } } @@ -368,17 +386,7 @@ export class NumericAggregationPlugin extends FunctionPlugin implements Function } public averagea(ast: ProcedureAst, state: InterpreterState): InternalScalarValue { - const result = this.reduce(ast.args, state, MomentsAggregate.empty, '_AGGREGATE_A', - (left, right) => left.compose(right), - (arg): MomentsAggregate => MomentsAggregate.single(getRawValue(arg)), - numbersBooleans - ) - - if (result instanceof CellError) { - return result - } else { - return result.averageValue() ?? new CellError(ErrorType.DIV_BY_ZERO) - } + return this.averageOf(ast.args, state, '_AVERAGE_A', numbersBooleans) } public vars(ast: ProcedureAst, state: InterpreterState): InternalScalarValue { @@ -509,13 +517,29 @@ export class NumericAggregationPlugin extends FunctionPlugin implements Function } private doAverage(args: Ast[], state: InterpreterState): InternalScalarValue { - const result = this.reduceAggregate(args, state) + return this.averageOf(args, state, '_AVERAGE', strictlyNumbers) + } - if (result instanceof CellError) { - return result - } else { - return result.averageValue() ?? new CellError(ErrorType.DIV_BY_ZERO) + /** + * The plain sum of the values divided by their count, as in AVERAGE and AVERAGEA. The sum and the + * count are cached per range under their own keys, so AVERAGE does not compute the variance sums. + * + * @param {Ast[]} args - the function arguments + * @param {InterpreterState} state - interpreter state + * @param {string} cacheKey - prefix of the range cache keys for the sum and the count + * @param {coercionOperation} coercion - which values count + * @returns {InternalScalarValue} the average, `#DIV/0!` when there are no values, or the first error + */ + private averageOf(args: Ast[], state: InterpreterState, cacheKey: string, coercion: coercionOperation): InternalScalarValue { + const sum = this.reduce(args, state, 0, `${cacheKey}_SUM`, (left, right) => left + right, getRawValue, coercion) + if (sum instanceof CellError) { + return sum + } + const count = this.reduce(args, state, 0, `${cacheKey}_COUNT`, (left, right) => left + right, () => 1, coercion) + if (count instanceof CellError) { + return count } + return count > 0 ? sum / count : new CellError(ErrorType.DIV_BY_ZERO) } private doVarS(args: Ast[], state: InterpreterState): InternalScalarValue { diff --git a/src/interpreter/plugin/StatisticalAggregationPlugin.ts b/src/interpreter/plugin/StatisticalAggregationPlugin.ts index 73b79d776d..5fee2d208d 100644 --- a/src/interpreter/plugin/StatisticalAggregationPlugin.ts +++ b/src/interpreter/plugin/StatisticalAggregationPlugin.ts @@ -27,13 +27,12 @@ import { sumsqerr, variance } from './3rdparty/jstat/jstat' -import {sumOfProductsOfDeviations, sumOfSquaredDeviations} from '../deviationSums' +import {covariance, deviationsFromMean, sumOfProducts, sumOfSquaredDeviations} from '../deviationSums' import { - addDoubleDouble, - divideDoubleDouble, divideDoubleDoubles, multiplyDoubleDouble, roundDoubleDouble, + subtractDoubleDouble, } from '../doubleDouble' import {FunctionArgumentType, FunctionPlugin, FunctionPluginTypecheck, ImplementedFunctions} from './FunctionPlugin' @@ -289,7 +288,7 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n === 1) { return 0 } - return roundDoubleDouble(divideDoubleDouble(sumOfProductsOfDeviations(ret[0], ret[1]), n)) + return covariance(ret[0], ret[1], 0) }) } @@ -308,7 +307,7 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 1) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.TwoValues) } - return roundDoubleDouble(divideDoubleDouble(sumOfProductsOfDeviations(ret[0], ret[1]), n - 1)) + return covariance(ret[0], ret[1], 1) }) } @@ -377,10 +376,16 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 2) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.ThreeValues) } - const sumOfProducts = sumOfProductsOfDeviations(ret[0], ret[1]) - const slope = divideDoubleDoubles(sumOfProducts, sumOfSquaredDeviations(ret[1])) - const explained = multiplyDoubleDouble(slope, sumOfProducts) - const residual = addDoubleDouble(sumOfSquaredDeviations(ret[0]), {hi: -explained.hi, lo: -explained.lo}) + const yDeviations = deviationsFromMean(ret[0]) + const xDeviations = deviationsFromMean(ret[1]) + const xSumOfSquares = sumOfProducts(xDeviations, xDeviations) + if (!Number.isFinite(roundDoubleDouble(xSumOfSquares))) { + return new CellError(ErrorType.NUM, ErrorMessage.NaN) + } + const productsSum = sumOfProducts(yDeviations, xDeviations) + const slope = divideDoubleDoubles(productsSum, xSumOfSquares) + const explained = multiplyDoubleDouble(slope, productsSum) + const residual = subtractDoubleDouble(sumOfProducts(yDeviations, yDeviations), explained) // the residual is non-negative; a tiny negative rounding remainder counts as 0 return Math.sqrt(Math.max(0, roundDoubleDouble(residual)) / (n - 2)) }) @@ -401,7 +406,12 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 1) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.TwoValues) } - return roundDoubleDouble(divideDoubleDoubles(sumOfProductsOfDeviations(ret[0], ret[1]), sumOfSquaredDeviations(ret[1]))) + const xDeviations = deviationsFromMean(ret[1]) + const xSumOfSquares = sumOfProducts(xDeviations, xDeviations) + if (!Number.isFinite(roundDoubleDouble(xSumOfSquares))) { + return new CellError(ErrorType.NUM, ErrorMessage.NaN) + } + return roundDoubleDouble(divideDoubleDoubles(sumOfProducts(deviationsFromMean(ret[0]), xDeviations), xSumOfSquares)) }) } From d3a1bec56068aed4f3e3e118acf119364a3c4f56 Mon Sep 17 00:00:00 2001 From: Kuba Sekowski Date: Wed, 7 Oct 2026 21:19:03 +0200 Subject: [PATCH 04/13] Fix STEYX returning 0 when the slope overflows, and return #DIV/0! from SLOPE and STEYX for equal x values (HF-491) STEYX computed the explained sum of squares as slope * Sxy. When the slope overflowed, the residual became -Infinity and was clamped to 0. The explained sum is now computed as Sxy^2 / Sxx in that case, and a non-finite residual returns #NUM! instead of 0. When all the x values are equal, SLOPE and STEYX return #DIV/0!, as in Excel and Google Sheets, instead of #NUM!. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 2 +- .../plugin/StatisticalAggregationPlugin.ts | 18 +++++++++++++++--- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a674b9ac3..c45fcb8028 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,7 +13,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Fixed the localized names of `VSTACK` and `HSTACK` in 14 language packs to match Microsoft Excel. [#1748](https://github.com/handsontable/hyperformula/pull/1748) - Fixed the MAXPOOL and MEDIANPOOL functions throwing an uncaught `TypeError` instead of returning the `#VALUE!` error when the range dimensions are not a whole multiple of the window size and the stride. [#1718](https://github.com/handsontable/hyperformula/pull/1718) - Fixed the variance and standard deviation functions (`VAR.S`, `VAR.P`, `STDEV.S`, `STDEV.P`, `VARA`, `VARPA`, `STDEVA`, `STDEVPA`, their legacy names, and `SUBTOTAL` modes 7, 8, 10, 11, 107, 108, 110 and 111) losing all precision on data with a large mean and a small spread, where they returned `0`, `#NUM!`, or a negative variance. They now return the variance of the stored values accurate to about the last digit. [#1784](https://github.com/handsontable/hyperformula/pull/1784) -- Fixed `DEVSQ`, `COVARIANCE.P`, `COVARIANCE.S` (and their aliases `COVAR`, `COVARIANCEP` and `COVARIANCES`), `SLOPE`, `STEYX`, `DSTDEV`, `DSTDEVP`, `DVAR` and `DVARP` losing digits when the values use all 15 significant digits and their mean is large relative to their spread (for example, `STEYX` returned `0.0187` instead of `0.00623` and `SLOPE` returned `1000.37` instead of `1000.53` for three values near `100000000000`). They now return results accurate to the last digits of the stored values. `STEYX` no longer returns `#NUM!` for points that lie exactly on a line, and `SLOPE` and `STEYX` return `#NUM!` instead of an arbitrary number when all the x values are equal or the sum of their squared deviations overflows. [#1798](https://github.com/handsontable/hyperformula/pull/1798) +- Fixed `DEVSQ`, `COVARIANCE.P`, `COVARIANCE.S` (and their aliases `COVAR`, `COVARIANCEP` and `COVARIANCES`), `SLOPE`, `STEYX`, `DSTDEV`, `DSTDEVP`, `DVAR` and `DVARP` losing digits when the values use all 15 significant digits and their mean is large relative to their spread (for example, `STEYX` returned `0.0187` instead of `0.00623` and `SLOPE` returned `1000.37` instead of `1000.53` for three values near `100000000000`). They now return results accurate to the last digits of the stored values. `STEYX` no longer returns `#NUM!` for points that lie exactly on a line, `SLOPE` and `STEYX` return `#DIV/0!` when all the x values are equal, and `#NUM!` instead of an arbitrary number when the sum of squared deviations of the x values overflows. [#1798](https://github.com/handsontable/hyperformula/pull/1798) - Fixed the `MOD` function returning a remainder with the sign of the dividend instead of the sign of the divisor, which made the results differ from Excel and Google Sheets for arguments with opposite signs (e.g. `=MOD(-3, 12)` now returns `9` instead of `-3`). [#1747](https://github.com/handsontable/hyperformula/issues/1747) ## [3.4.0] - 2026-08-10 diff --git a/src/interpreter/plugin/StatisticalAggregationPlugin.ts b/src/interpreter/plugin/StatisticalAggregationPlugin.ts index 5fee2d208d..4b234d21d3 100644 --- a/src/interpreter/plugin/StatisticalAggregationPlugin.ts +++ b/src/interpreter/plugin/StatisticalAggregationPlugin.ts @@ -379,15 +379,24 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func const yDeviations = deviationsFromMean(ret[0]) const xDeviations = deviationsFromMean(ret[1]) const xSumOfSquares = sumOfProducts(xDeviations, xDeviations) + if (xSumOfSquares.hi === 0) { + return new CellError(ErrorType.DIV_BY_ZERO) + } if (!Number.isFinite(roundDoubleDouble(xSumOfSquares))) { return new CellError(ErrorType.NUM, ErrorMessage.NaN) } const productsSum = sumOfProducts(yDeviations, xDeviations) const slope = divideDoubleDoubles(productsSum, xSumOfSquares) - const explained = multiplyDoubleDouble(slope, productsSum) - const residual = subtractDoubleDouble(sumOfProducts(yDeviations, yDeviations), explained) + // when the slope overflows, the explained sum of squares Sxy^2 / Sxx can still be finite + const explained = Number.isFinite(slope.hi) + ? multiplyDoubleDouble(slope, productsSum) + : divideDoubleDoubles(multiplyDoubleDouble(productsSum, productsSum), xSumOfSquares) + const residual = roundDoubleDouble(subtractDoubleDouble(sumOfProducts(yDeviations, yDeviations), explained)) + if (!Number.isFinite(residual)) { + return new CellError(ErrorType.NUM, ErrorMessage.NaN) + } // the residual is non-negative; a tiny negative rounding remainder counts as 0 - return Math.sqrt(Math.max(0, roundDoubleDouble(residual)) / (n - 2)) + return Math.sqrt(Math.max(0, residual) / (n - 2)) }) } @@ -408,6 +417,9 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func } const xDeviations = deviationsFromMean(ret[1]) const xSumOfSquares = sumOfProducts(xDeviations, xDeviations) + if (xSumOfSquares.hi === 0) { + return new CellError(ErrorType.DIV_BY_ZERO) + } if (!Number.isFinite(roundDoubleDouble(xSumOfSquares))) { return new CellError(ErrorType.NUM, ErrorMessage.NaN) } From 92ebfbd495ae633261795603eea493e6c7901903 Mon Sep 17 00:00:00 2001 From: Kuba Sekowski Date: Wed, 7 Oct 2026 21:19:03 +0200 Subject: [PATCH 05/13] Keep the variance sums with a power-of-two exponent and compute AVERAGE in one pass (HF-491) MomentsAggregate stores the deviations from the shift multiplied by 2^-exponent, so the shifted sums cannot overflow. The exponent is raised when the scaled sum of squares exceeds 2^960, and two aggregates are composed at the larger exponent. The variance is rounded once and scaled back. It is #NUM! only when it exceeds the largest double, and it no longer depends on the order of the arguments. The plain sum and sum of squares and the one-pass fallback are removed. AVERAGE, AVERAGEA and SUBTOTAL 1/101 reduce a single AverageResult (sum and count), moved out of ConditionalAggregationPlugin into its own module, so each argument is evaluated once. Before, the sum and the count were reduced in two passes, which evaluated non-range arguments twice: nested AVERAGE calls took exponential time, and a volatile argument could give a sum and a count from different evaluations. Co-Authored-By: Claude Opus 5.5 --- src/interpreter/plugin/AverageResult.ts | 35 +++ .../plugin/ConditionalAggregationPlugin.ts | 26 +- .../plugin/NumericAggregationPlugin.ts | 222 +++++++++++++----- 3 files changed, 196 insertions(+), 87 deletions(-) create mode 100644 src/interpreter/plugin/AverageResult.ts diff --git a/src/interpreter/plugin/AverageResult.ts b/src/interpreter/plugin/AverageResult.ts new file mode 100644 index 0000000000..86351f9625 --- /dev/null +++ b/src/interpreter/plugin/AverageResult.ts @@ -0,0 +1,35 @@ +/** + * @license + * Copyright (c) 2025 Handsoncode. All rights reserved. + */ + +import {Maybe} from '../../Maybe' + +/** + * The sum and the count of a set of numbers, composable so that the value of a range can be cached + * and reused for a larger range. + */ +export class AverageResult { + public static empty = new AverageResult(0, 0) + + constructor( + public readonly sum: number, + public readonly count: number, + ) {} + + public static single(arg: number): AverageResult { + return new AverageResult(arg, 1) + } + + public compose(other: AverageResult) { + return new AverageResult(this.sum + other.sum, this.count + other.count) + } + + public averageValue(): Maybe { + if (this.count > 0) { + return this.sum / this.count + } else { + return undefined + } + } +} diff --git a/src/interpreter/plugin/ConditionalAggregationPlugin.ts b/src/interpreter/plugin/ConditionalAggregationPlugin.ts index 34e0fc7db6..774f418f6f 100644 --- a/src/interpreter/plugin/ConditionalAggregationPlugin.ts +++ b/src/interpreter/plugin/ConditionalAggregationPlugin.ts @@ -18,33 +18,9 @@ import { RawScalarValue } from '../InterpreterValue' import {SimpleRangeValue} from '../../SimpleRangeValue' +import {AverageResult} from './AverageResult' import {FunctionArgumentType, FunctionPlugin, FunctionPluginTypecheck, ImplementedFunctions} from './FunctionPlugin' -class AverageResult { - public static empty = new AverageResult(0, 0) - - constructor( - public readonly sum: number, - public readonly count: number, - ) {} - - public static single(arg: number): AverageResult { - return new AverageResult(arg, 1) - } - - public compose(other: AverageResult) { - return new AverageResult(this.sum + other.sum, this.count + other.count) - } - - public averageValue(): Maybe { - if (this.count > 0) { - return this.sum / this.count - } else { - return undefined - } - } -} - /** Computes key for criterion function cache */ function conditionalAggregationFunctionCacheKey(functionName: string): (conditions: Condition[]) => string { return (conditions: Condition[]): string => { diff --git a/src/interpreter/plugin/NumericAggregationPlugin.ts b/src/interpreter/plugin/NumericAggregationPlugin.ts index 8e76a7ace6..77f61cb0e4 100644 --- a/src/interpreter/plugin/NumericAggregationPlugin.ts +++ b/src/interpreter/plugin/NumericAggregationPlugin.ts @@ -16,6 +16,7 @@ import { divideDoubleDouble, DOUBLE_DOUBLE_ZERO, DoubleDouble, + multiplyByPowerOfTwo, multiplyDoubleDouble, roundDoubleDouble, scaleDoubleDouble, @@ -25,6 +26,7 @@ import { import {InterpreterState} from '../InterpreterState' import {EmptyValue, ExtendedNumber, getRawValue, InternalScalarValue, isExtendedNumber} from '../InterpreterValue' import {SimpleRangeValue} from '../../SimpleRangeValue' +import {AverageResult} from './AverageResult' import {FunctionArgumentType, FunctionPlugin, FunctionPluginTypecheck, ImplementedFunctions} from './FunctionPlugin' import {RangeVertex} from '../../DependencyGraph' @@ -42,11 +44,35 @@ function zeroForInfinite(value: InternalScalarValue) { } } +/** + * The largest scaled sum of squared deviations a `MomentsAggregate` keeps (2^960). Above it, the + * aggregate raises its exponent. Every scaled deviation is then at most 2^480, so rebasing the sums of + * two aggregates cannot overflow. + */ +const MAX_SCALED_SUM_OF_SQUARES = 2 ** 960 + +/** + * The largest scaled difference of shifts that two aggregates are rebased with (2^470). Above it, the + * exponent is raised first. + */ +const MAX_SCALED_SHIFT_DIFFERENCE = 2 ** 470 + +/** + * The smallest exponent `e >= 0` (up to one) for which `magnitude * 2^-e` is at most `limit`. + * + * @param {number} magnitude - a non-negative finite number + * @param {number} limit - a positive power of two + * @returns {number} the exponent + */ +function exponentToFit(magnitude: number, limit: number): number { + return magnitude <= limit ? 0 : Math.ceil(Math.log2(magnitude / limit)) +} + /** * Moments of a set of numbers, composable so that the value of a range can be cached and reused * for a larger range. * - * Besides the sum and the count, it keeps the sums of deviations and of squared deviations from a + * Besides the count, it keeps the sums of deviations `S1` and of squared deviations `S2` from a * `shift`: the first value added. Both sums are double-double. Measured from a nearby value, the * deviations stay small, which avoids the cancellation of the textbook one-pass form * `sum(x^2) - sum(x)^2 / n` for data with a large mean and a small spread (for example 10000000.001, @@ -54,21 +80,26 @@ function zeroForInfinite(value: InternalScalarValue) { * last place of the exact values for the stored numbers. Double-double arithmetic does not guarantee * correct rounding. * - * It also keeps the plain sum and sum of squares. When the deviations from the shift are so large - * (about 1e154 and above) that the shifted sums overflow, the variance falls back to the one-pass form - * used before, with the same accuracy as before: it loses about `log10(2n)` digits to cancellation. - * It is `#NUM!` when the plain sum of squares overflows too, as before, although the shifted sums of - * the same values in another order can be finite. + * The deviations are stored multiplied by `2^-exponent`, which is exact, so that the sums cannot + * overflow. The exponent is 0 until the scaled sum of squares exceeds `MAX_SCALED_SUM_OF_SQUARES` + * (deviations of about 1e144 and above), and grows as needed. The variance is `#NUM!` only when it + * exceeds the largest double, whatever the order of the values. */ class MomentsAggregate { - public static empty = new MomentsAggregate(0, 0, 0, 0, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) + public static empty = new MomentsAggregate(0, 0, 0, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) + /** + * @param {number} count - the number of values + * @param {number} shift - the value the deviations are measured from + * @param {number} exponent - the deviations are stored multiplied by `2^-exponent` + * @param {DoubleDouble} shiftedSum - `S1`, the sum of `(x - shift) * 2^-exponent` + * @param {DoubleDouble} shiftedSumOfSquares - `S2`, the sum of `((x - shift) * 2^-exponent)^2` + */ constructor( - public readonly sum: number, - public readonly sumOfSquares: number, public readonly count: number, public readonly shift: number, + public readonly exponent: number, public readonly shiftedSum: DoubleDouble, public readonly shiftedSumOfSquares: DoubleDouble, ) { @@ -81,13 +112,40 @@ class MomentsAggregate { * @returns {MomentsAggregate} an aggregate of the single value */ public static single(arg: number): MomentsAggregate { - return new MomentsAggregate(arg, arg * arg, 1, arg, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) + return new MomentsAggregate(1, arg, 0, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) + } + + /** + * An aggregate whose exponent is raised, if needed, so that its scaled sum of squares is at most + * `MAX_SCALED_SUM_OF_SQUARES`. + * + * @param {number} count - the number of values + * @param {number} shift - the value the deviations are measured from + * @param {number} exponent - the exponent of the given sums + * @param {DoubleDouble} shiftedSum - the scaled sum of deviations + * @param {DoubleDouble} shiftedSumOfSquares - the scaled sum of squared deviations + * @returns {MomentsAggregate} the aggregate + */ + private static normalized(count: number, shift: number, exponent: number, shiftedSum: DoubleDouble, shiftedSumOfSquares: DoubleDouble): MomentsAggregate { + const excess = exponentToFit(Math.abs(shiftedSumOfSquares.hi), MAX_SCALED_SUM_OF_SQUARES) + if (excess === 0) { + return new MomentsAggregate(count, shift, exponent, shiftedSum, shiftedSumOfSquares) + } + const raise = Math.ceil(excess / 2) + return new MomentsAggregate(count, shift, exponent + raise, + multiplyByPowerOfTwo(shiftedSum, 2 ** -raise), + multiplyByPowerOfTwo(shiftedSumOfSquares, 2 ** (-2 * raise)), + ) } /** * Combines two aggregates. The result keeps this aggregate's `shift`; the other one's shifted sums - * are re-expressed relative to it. An empty aggregate is the identity: composing with it returns the - * other aggregate unchanged, with its own shift. + * are re-expressed relative to it, at the larger exponent of the two (raised further when the + * difference of the shifts needs it). An empty aggregate is the identity: composing with it returns + * the other aggregate unchanged, with its own shift. + * + * @param {MomentsAggregate} other - the aggregate to add + * @returns {MomentsAggregate} the aggregate of the values of both */ public compose(other: MomentsAggregate): MomentsAggregate { if (this.count === 0) { @@ -97,39 +155,36 @@ class MomentsAggregate { return this } - const shiftDifference = twoSum(other.shift, -this.shift) - - if (other.count === 1) { - // the common case of adding one value: its deviation is the shift difference itself - return new MomentsAggregate( - this.sum + other.sum, - this.sumOfSquares + other.sumOfSquares, - this.count + 1, - this.shift, - addDoubleDouble(this.shiftedSum, shiftDifference), - addDoubleDouble(this.shiftedSumOfSquares, multiplyDoubleDouble(shiftDifference, shiftDifference)), - ) + // the common case: no rescaling, because both aggregates have the same exponent or the other one + // is a single value, whose shifted sums are 0 at any exponent + if (other.exponent === this.exponent || (other.count === 1 && this.exponent > 0)) { + const exponent = this.exponent + const shiftDifference = exponent === 0 + ? twoSum(other.shift, -this.shift) + : twoSum(other.shift * 2 ** -exponent, -this.shift * 2 ** -exponent) + if (Math.abs(shiftDifference.hi) <= MAX_SCALED_SHIFT_DIFFERENCE) { + return this.composeScaled(other, exponent, this.shiftedSum, this.shiftedSumOfSquares, other.shiftedSum, other.shiftedSumOfSquares, shiftDifference) + } } - // other's sums rebased: S1 + n*d and S2 + 2*d*S1 + n*d^2 - const rebasedSum = addDoubleDouble(other.shiftedSum, scaleDoubleDouble(shiftDifference, other.count)) - const rebasedSumOfSquares = addDoubleDouble( - addDoubleDouble(other.shiftedSumOfSquares, scaleDoubleDouble(multiplyDoubleDouble(shiftDifference, other.shiftedSum), 2)), - scaleDoubleDouble(multiplyDoubleDouble(shiftDifference, shiftDifference), other.count), - ) - return new MomentsAggregate( - this.sum + other.sum, - this.sumOfSquares + other.sumOfSquares, - this.count + other.count, - this.shift, - addDoubleDouble(this.shiftedSum, rebasedSum), - addDoubleDouble(this.shiftedSumOfSquares, rebasedSumOfSquares), + // halving the shifts first keeps their difference finite + const shiftDifferenceExponent = exponentToFit(Math.abs(other.shift / 2 - this.shift / 2), MAX_SCALED_SHIFT_DIFFERENCE / 2) + const exponent = Math.max(this.exponent, other.exponent, shiftDifferenceExponent) + const thisRaise = exponent - this.exponent + const otherRaise = exponent - other.exponent + const scale = 2 ** -exponent + return this.composeScaled(other, exponent, + multiplyByPowerOfTwo(this.shiftedSum, 2 ** -thisRaise), + multiplyByPowerOfTwo(this.shiftedSumOfSquares, 2 ** (-2 * thisRaise)), + multiplyByPowerOfTwo(other.shiftedSum, 2 ** -otherRaise), + multiplyByPowerOfTwo(other.shiftedSumOfSquares, 2 ** (-2 * otherRaise)), + twoSum(other.shift * scale, -this.shift * scale), ) } public varSValue(): Maybe { if (this.count > 1) { - return roundDoubleDouble(divideDoubleDouble(this.sumOfSquaredDeviations(), this.count - 1)) + return this.variance(this.count - 1) } else { return undefined } @@ -137,29 +192,72 @@ class MomentsAggregate { public varPValue(): Maybe { if (this.count > 0) { - return roundDoubleDouble(divideDoubleDouble(this.sumOfSquaredDeviations(), this.count)) + return this.variance(this.count) } else { return undefined } } /** - * The sum of squared deviations from the mean, `S2 - S1^2 / n` with `S1`, `S2` the shifted sums - * and `n` the count, as a double-double: the variance divides it before rounding once. + * Adds the other aggregate's sums, already scaled to `exponent`, to this aggregate's sums, also + * scaled to `exponent`. * - * `S1^2 / n` is evaluated as `S1 * (S1 / n)`, so it cannot overflow while `S2` is finite - * (`S1^2 / n <= S2`), even when `S1^2` alone would exceed the largest double. When the shifted sums - * overflow, it falls back to `sum(x^2) - sum(x)^2 / n` in plain doubles, evaluated as before unless - * `sum(x)^2` overflows. + * @param {MomentsAggregate} other - the aggregate to add + * @param {number} exponent - the common exponent + * @param {DoubleDouble} shiftedSum - this aggregate's `S1` at `exponent` + * @param {DoubleDouble} shiftedSumOfSquares - this aggregate's `S2` at `exponent` + * @param {DoubleDouble} otherShiftedSum - the other aggregate's `S1` at `exponent` + * @param {DoubleDouble} otherShiftedSumOfSquares - the other aggregate's `S2` at `exponent` + * @param {DoubleDouble} shiftDifference - `(other.shift - this.shift) * 2^-exponent` + * @returns {MomentsAggregate} the aggregate of the values of both */ - private sumOfSquaredDeviations(): DoubleDouble { - const squaredSumOverCount = multiplyDoubleDouble(this.shiftedSum, divideDoubleDouble(this.shiftedSum, this.count)) - const result = subtractDoubleDouble(this.shiftedSumOfSquares, squaredSumOverCount) - if (Number.isFinite(roundDoubleDouble(result))) { - return result + private composeScaled( + other: MomentsAggregate, + exponent: number, + shiftedSum: DoubleDouble, + shiftedSumOfSquares: DoubleDouble, + otherShiftedSum: DoubleDouble, + otherShiftedSumOfSquares: DoubleDouble, + shiftDifference: DoubleDouble, + ): MomentsAggregate { + const count = this.count + other.count + if (other.count === 1) { + // the common case of adding one value: its deviation is the shift difference itself + return MomentsAggregate.normalized(count, this.shift, exponent, + addDoubleDouble(shiftedSum, shiftDifference), + addDoubleDouble(shiftedSumOfSquares, multiplyDoubleDouble(shiftDifference, shiftDifference)), + ) } - const onePass = this.sumOfSquares - this.sum * this.sum / this.count - return {hi: Number.isFinite(onePass) ? onePass : this.sumOfSquares - this.sum * (this.sum / this.count), lo: 0} + + // other's sums rebased: S1 + n*d and S2 + 2*d*S1 + n*d^2 + const rebasedSum = addDoubleDouble(otherShiftedSum, scaleDoubleDouble(shiftDifference, other.count)) + const rebasedSumOfSquares = addDoubleDouble( + addDoubleDouble(otherShiftedSumOfSquares, scaleDoubleDouble(multiplyDoubleDouble(shiftDifference, otherShiftedSum), 2)), + scaleDoubleDouble(multiplyDoubleDouble(shiftDifference, shiftDifference), other.count), + ) + return MomentsAggregate.normalized(count, this.shift, exponent, + addDoubleDouble(shiftedSum, rebasedSum), + addDoubleDouble(shiftedSumOfSquares, rebasedSumOfSquares), + ) + } + + /** + * The sum of squared deviations from the mean, `S2 - S1^2 / n` with `S1`, `S2` the shifted sums and + * `n` the count, divided by `divisor` and rounded once. + * + * `S1^2 / n` is evaluated as `S1 * (S1 / n)`, which cannot overflow because `S1^2 / n <= S2`. The + * quotient is computed at the scale of the sums and then multiplied by `2^(2 * exponent)` in two + * steps, so that the power of two itself cannot overflow; the result overflows only when the + * variance exceeds the largest double. + * + * @param {number} divisor - `n - 1` for the sample variance, `n` for the population variance + * @returns {number} the variance + */ + private variance(divisor: number): number { + const squaredSumOverCount = multiplyDoubleDouble(this.shiftedSum, divideDoubleDouble(this.shiftedSum, this.count)) + const sumOfSquaredDeviations = subtractDoubleDouble(this.shiftedSumOfSquares, squaredSumOverCount) + const scale = 2 ** this.exponent + return roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations, divisor)) * scale * scale } } @@ -522,24 +620,24 @@ export class NumericAggregationPlugin extends FunctionPlugin implements Function /** * The plain sum of the values divided by their count, as in AVERAGE and AVERAGEA. The sum and the - * count are cached per range under their own keys, so AVERAGE does not compute the variance sums. + * count are folded in one pass and cached per range, so AVERAGE does not compute the variance sums. * * @param {Ast[]} args - the function arguments * @param {InterpreterState} state - interpreter state - * @param {string} cacheKey - prefix of the range cache keys for the sum and the count + * @param {string} cacheKey - the range cache key * @param {coercionOperation} coercion - which values count * @returns {InternalScalarValue} the average, `#DIV/0!` when there are no values, or the first error */ private averageOf(args: Ast[], state: InterpreterState, cacheKey: string, coercion: coercionOperation): InternalScalarValue { - const sum = this.reduce(args, state, 0, `${cacheKey}_SUM`, (left, right) => left + right, getRawValue, coercion) - if (sum instanceof CellError) { - return sum - } - const count = this.reduce(args, state, 0, `${cacheKey}_COUNT`, (left, right) => left + right, () => 1, coercion) - if (count instanceof CellError) { - return count + const result = this.reduce(args, state, AverageResult.empty, cacheKey, + (left, right) => left.compose(right), + (arg) => AverageResult.single(getRawValue(arg)), + coercion, + ) + if (result instanceof CellError) { + return result } - return count > 0 ? sum / count : new CellError(ErrorType.DIV_BY_ZERO) + return result.averageValue() ?? new CellError(ErrorType.DIV_BY_ZERO) } private doVarS(args: Ast[], state: InterpreterState): InternalScalarValue { From 1d0c8229e55fa36a1e125de77f282a78bee10f31 Mon Sep 17 00:00:00 2001 From: marcin-kordas-hoc Date: Thu, 8 Oct 2026 04:23:44 +0000 Subject: [PATCH 06/13] Shorten the variance changelog entries (HF-491) Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c45fcb8028..478d62658c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,8 +12,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Fixed the `AVERAGEIF` function returning a division-by-zero error when the calculated average was `0`. [#1733](https://github.com/handsontable/hyperformula/pull/1733) - Fixed the localized names of `VSTACK` and `HSTACK` in 14 language packs to match Microsoft Excel. [#1748](https://github.com/handsontable/hyperformula/pull/1748) - Fixed the MAXPOOL and MEDIANPOOL functions throwing an uncaught `TypeError` instead of returning the `#VALUE!` error when the range dimensions are not a whole multiple of the window size and the stride. [#1718](https://github.com/handsontable/hyperformula/pull/1718) -- Fixed the variance and standard deviation functions (`VAR.S`, `VAR.P`, `STDEV.S`, `STDEV.P`, `VARA`, `VARPA`, `STDEVA`, `STDEVPA`, their legacy names, and `SUBTOTAL` modes 7, 8, 10, 11, 107, 108, 110 and 111) losing all precision on data with a large mean and a small spread, where they returned `0`, `#NUM!`, or a negative variance. They now return the variance of the stored values accurate to about the last digit. [#1784](https://github.com/handsontable/hyperformula/pull/1784) -- Fixed `DEVSQ`, `COVARIANCE.P`, `COVARIANCE.S` (and their aliases `COVAR`, `COVARIANCEP` and `COVARIANCES`), `SLOPE`, `STEYX`, `DSTDEV`, `DSTDEVP`, `DVAR` and `DVARP` losing digits when the values use all 15 significant digits and their mean is large relative to their spread (for example, `STEYX` returned `0.0187` instead of `0.00623` and `SLOPE` returned `1000.37` instead of `1000.53` for three values near `100000000000`). They now return results accurate to the last digits of the stored values. `STEYX` no longer returns `#NUM!` for points that lie exactly on a line, `SLOPE` and `STEYX` return `#DIV/0!` when all the x values are equal, and `#NUM!` instead of an arbitrary number when the sum of squared deviations of the x values overflows. [#1798](https://github.com/handsontable/hyperformula/pull/1798) +- Fixed the variance and standard deviation functions (`VAR`, `STDEV` and their variants, and the matching `SUBTOTAL` modes) returning `0`, `#NUM!`, or a negative variance on data with a large mean and a small spread. [#1784](https://github.com/handsontable/hyperformula/pull/1784) +- Fixed `DEVSQ`, `COVARIANCE`, `SLOPE`, `STEYX`, and the `DSTDEV` and `DVAR` functions losing digits on data with a large mean and a small spread. [#1798](https://github.com/handsontable/hyperformula/pull/1798) - Fixed the `MOD` function returning a remainder with the sign of the dividend instead of the sign of the divisor, which made the results differ from Excel and Google Sheets for arguments with opposite signs (e.g. `=MOD(-3, 12)` now returns `9` instead of `-3`). [#1747](https://github.com/handsontable/hyperformula/issues/1747) ## [3.4.0] - 2026-08-10 From 4c97baefad925f742b3ed7871e4a81c09644efd0 Mon Sep 17 00:00:00 2001 From: Kuba Sekowski Date: Thu, 8 Oct 2026 08:15:54 +0200 Subject: [PATCH 07/13] Share the x-deviation checks of SLOPE and STEYX and route AVERAGEA through doAverageA (HF-491) Co-Authored-By: Claude Opus 5.5 --- .../plugin/NumericAggregationPlugin.ts | 6 ++- .../plugin/StatisticalAggregationPlugin.ts | 44 ++++++++++++------- 2 files changed, 34 insertions(+), 16 deletions(-) diff --git a/src/interpreter/plugin/NumericAggregationPlugin.ts b/src/interpreter/plugin/NumericAggregationPlugin.ts index 77f61cb0e4..956c456a51 100644 --- a/src/interpreter/plugin/NumericAggregationPlugin.ts +++ b/src/interpreter/plugin/NumericAggregationPlugin.ts @@ -484,7 +484,7 @@ export class NumericAggregationPlugin extends FunctionPlugin implements Function } public averagea(ast: ProcedureAst, state: InterpreterState): InternalScalarValue { - return this.averageOf(ast.args, state, '_AVERAGE_A', numbersBooleans) + return this.doAverageA(ast.args, state) } public vars(ast: ProcedureAst, state: InterpreterState): InternalScalarValue { @@ -618,6 +618,10 @@ export class NumericAggregationPlugin extends FunctionPlugin implements Function return this.averageOf(args, state, '_AVERAGE', strictlyNumbers) } + private doAverageA(args: Ast[], state: InterpreterState): InternalScalarValue { + return this.averageOf(args, state, '_AVERAGE_A', numbersBooleans) + } + /** * The plain sum of the values divided by their count, as in AVERAGE and AVERAGEA. The sum and the * count are folded in one pass and cached per range, so AVERAGE does not compute the variance sums. diff --git a/src/interpreter/plugin/StatisticalAggregationPlugin.ts b/src/interpreter/plugin/StatisticalAggregationPlugin.ts index 4b234d21d3..f81f12ab71 100644 --- a/src/interpreter/plugin/StatisticalAggregationPlugin.ts +++ b/src/interpreter/plugin/StatisticalAggregationPlugin.ts @@ -30,6 +30,7 @@ import { import {covariance, deviationsFromMean, sumOfProducts, sumOfSquaredDeviations} from '../deviationSums' import { divideDoubleDoubles, + DoubleDouble, multiplyDoubleDouble, roundDoubleDouble, subtractDoubleDouble, @@ -376,15 +377,12 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 2) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.ThreeValues) } - const yDeviations = deviationsFromMean(ret[0]) - const xDeviations = deviationsFromMean(ret[1]) - const xSumOfSquares = sumOfProducts(xDeviations, xDeviations) - if (xSumOfSquares.hi === 0) { - return new CellError(ErrorType.DIV_BY_ZERO) - } - if (!Number.isFinite(roundDoubleDouble(xSumOfSquares))) { - return new CellError(ErrorType.NUM, ErrorMessage.NaN) + const xMoments = xDeviationsAndSumOfSquares(ret[1]) + if (xMoments instanceof CellError) { + return xMoments } + const {xDeviations, xSumOfSquares} = xMoments + const yDeviations = deviationsFromMean(ret[0]) const productsSum = sumOfProducts(yDeviations, xDeviations) const slope = divideDoubleDoubles(productsSum, xSumOfSquares) // when the slope overflows, the explained sum of squares Sxy^2 / Sxx can still be finite @@ -415,14 +413,11 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 1) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.TwoValues) } - const xDeviations = deviationsFromMean(ret[1]) - const xSumOfSquares = sumOfProducts(xDeviations, xDeviations) - if (xSumOfSquares.hi === 0) { - return new CellError(ErrorType.DIV_BY_ZERO) - } - if (!Number.isFinite(roundDoubleDouble(xSumOfSquares))) { - return new CellError(ErrorType.NUM, ErrorMessage.NaN) + const xMoments = xDeviationsAndSumOfSquares(ret[1]) + if (xMoments instanceof CellError) { + return xMoments } + const {xDeviations, xSumOfSquares} = xMoments return roundDoubleDouble(divideDoubleDoubles(sumOfProducts(deviationsFromMean(ret[0]), xDeviations), xSumOfSquares)) }) } @@ -554,6 +549,25 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func } } +/** + * The deviations of the x values of SLOPE or STEYX from their mean and the sum of their squares. + * + * @param {number[]} values - the x values + * @returns {CellError | {xDeviations: DoubleDouble[], xSumOfSquares: DoubleDouble}} `#DIV/0!` when all the + * x values are equal, `#NUM!` when the sum of squares overflows + */ +function xDeviationsAndSumOfSquares(values: number[]): CellError | {xDeviations: DoubleDouble[], xSumOfSquares: DoubleDouble} { + const xDeviations = deviationsFromMean(values) + const xSumOfSquares = sumOfProducts(xDeviations, xDeviations) + if (xSumOfSquares.hi === 0) { + return new CellError(ErrorType.DIV_BY_ZERO) + } + if (!Number.isFinite(roundDoubleDouble(xSumOfSquares))) { + return new CellError(ErrorType.NUM, ErrorMessage.NaN) + } + return {xDeviations, xSumOfSquares} +} + function parseTwoArrays(dataX: SimpleRangeValue, dataY: SimpleRangeValue): CellError | [number[], number[]] { const xit = dataX.iterateValuesFromTopLeftCorner() const yit = dataY.iterateValuesFromTopLeftCorner() From 9d88a3887173183fa28a1ee1502eb1cad128c618 Mon Sep 17 00:00:00 2001 From: Kuba Sekowski Date: Thu, 8 Oct 2026 08:15:54 +0200 Subject: [PATCH 08/13] Simplify the changelog entries for the variance fixes (HF-491) Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 478d62658c..89dbed3e4d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,8 +12,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Fixed the `AVERAGEIF` function returning a division-by-zero error when the calculated average was `0`. [#1733](https://github.com/handsontable/hyperformula/pull/1733) - Fixed the localized names of `VSTACK` and `HSTACK` in 14 language packs to match Microsoft Excel. [#1748](https://github.com/handsontable/hyperformula/pull/1748) - Fixed the MAXPOOL and MEDIANPOOL functions throwing an uncaught `TypeError` instead of returning the `#VALUE!` error when the range dimensions are not a whole multiple of the window size and the stride. [#1718](https://github.com/handsontable/hyperformula/pull/1718) -- Fixed the variance and standard deviation functions (`VAR`, `STDEV` and their variants, and the matching `SUBTOTAL` modes) returning `0`, `#NUM!`, or a negative variance on data with a large mean and a small spread. [#1784](https://github.com/handsontable/hyperformula/pull/1784) -- Fixed `DEVSQ`, `COVARIANCE`, `SLOPE`, `STEYX`, and the `DSTDEV` and `DVAR` functions losing digits on data with a large mean and a small spread. [#1798](https://github.com/handsontable/hyperformula/pull/1798) +- Fixed the `VAR`, `STDEV`, `DEVSQ`, `COVARIANCE`, `SLOPE`, `STEYX`, `DVAR` and `DSTDEV` functions, their variants, and the matching `SUBTOTAL` modes losing precision on data with a large mean and a small spread. [#1784](https://github.com/handsontable/hyperformula/pull/1784) +- Fixed a bug where `SLOPE` and `STEYX` returned `#NUM!` or an arbitrary number instead of `#DIV/0!` when all the x values are equal. [#1784](https://github.com/handsontable/hyperformula/pull/1784) +- Fixed a bug where `STEYX` returned `#NUM!` instead of `0` for points that lie exactly on a line. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed the `MOD` function returning a remainder with the sign of the dividend instead of the sign of the divisor, which made the results differ from Excel and Google Sheets for arguments with opposite signs (e.g. `=MOD(-3, 12)` now returns `9` instead of `-3`). [#1747](https://github.com/handsontable/hyperformula/issues/1747) ## [3.4.0] - 2026-08-10 From 6d4925ed9b4c66c04269607268419700b14f9152 Mon Sep 17 00:00:00 2001 From: Kuba Sekowski Date: Thu, 8 Oct 2026 08:16:55 +0200 Subject: [PATCH 09/13] Fold the D-variance functions through MomentsAggregate and scale the deviation sums by powers of two (HF-491) DVAR, DVARP, DSTDEV and DSTDEVP now use the same MomentsAggregate as VAR and STDEV, so they return the same results, including when the sum of squared deviations exceeds the largest double. MomentsAggregate moves to its own module, takes the standard deviation at the scale of its sums, and lowers its exponent for tiny differences, so STDEV no longer returns 0 for spreads below about 1e-154. COVARIANCE, SLOPE, STEYX and DEVSQ keep two passes over their arrays. An array whose largest magnitude is below 2^-400 is scaled up to it, so squared deviations cannot underflow and SLOPE and STEYX return #DIV/0! only when all the x values are equal. An array above 2^400 is scaled down only when the unscaled sums overflow. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 1 + src/interpreter/deviationSums.ts | 245 +++++++++--- src/interpreter/doubleDouble.ts | 16 + src/interpreter/plugin/DatabasePlugin.ts | 26 +- src/interpreter/plugin/MomentsAggregate.ts | 351 ++++++++++++++++++ .../plugin/NumericAggregationPlugin.ts | 242 +----------- .../plugin/StatisticalAggregationPlugin.ts | 58 +-- 7 files changed, 590 insertions(+), 349 deletions(-) create mode 100644 src/interpreter/plugin/MomentsAggregate.ts diff --git a/CHANGELOG.md b/CHANGELOG.md index 89dbed3e4d..38b68140a8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Fixed the `VAR`, `STDEV`, `DEVSQ`, `COVARIANCE`, `SLOPE`, `STEYX`, `DVAR` and `DSTDEV` functions, their variants, and the matching `SUBTOTAL` modes losing precision on data with a large mean and a small spread. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed a bug where `SLOPE` and `STEYX` returned `#NUM!` or an arbitrary number instead of `#DIV/0!` when all the x values are equal. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed a bug where `STEYX` returned `#NUM!` instead of `0` for points that lie exactly on a line. [#1784](https://github.com/handsontable/hyperformula/pull/1784) +- Fixed a bug where `STDEV`, `DVAR`, `DSTDEV`, `COVARIANCE`, `SLOPE`, `STEYX` and their variants returned an error or `0` for very large or very small values although the result was within the range of numbers. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed the `MOD` function returning a remainder with the sign of the dividend instead of the sign of the divisor, which made the results differ from Excel and Google Sheets for arguments with opposite signs (e.g. `=MOD(-3, 12)` now returns `9` instead of `-3`). [#1747](https://github.com/handsontable/hyperformula/issues/1747) ## [3.4.0] - 2026-08-10 diff --git a/src/interpreter/deviationSums.ts b/src/interpreter/deviationSums.ts index 6bf9a0c8a5..e29eccd87f 100644 --- a/src/interpreter/deviationSums.ts +++ b/src/interpreter/deviationSums.ts @@ -9,6 +9,7 @@ import { DOUBLE_DOUBLE_ZERO, DoubleDouble, multiplyByPowerOfTwo, + multiplyByTwoToThe, multiplyDoubleDouble, roundDoubleDouble, subtractDoubleDouble, @@ -16,96 +17,244 @@ import { /** * Sums of squared deviations and of products of deviations from the mean, for the functions built - * on them (DEVSQ, COVARIANCE, SLOPE, STEYX, DSTDEV, DVAR and their variants). + * on them (DEVSQ, COVARIANCE, SLOPE and STEYX). * * The mean, the deviations and their sums are all kept in double-double and rounded once at the end, * so the result is accurate to the last digits of the stored values even when the mean is large * relative to the spread. + * + * Each array is measured at its own power-of-two scale, and the rounded results are scaled back, so + * that a result is not lost to an intermediate sum that overflows or underflows: + * - an array whose largest magnitude is below `MIN_UNSCALED_MAGNITUDE` is scaled up to it, which is + * exact, so that its squared deviations cannot underflow; + * - an array whose largest magnitude is above `MAX_UNSCALED_MAGNITUDE` is scaled down to it only when + * the sums overflow without that. Scaling down rounds away the low-order bits of values about 2^1400 + * times smaller than the largest one, which matter when they pair with large deviations of the other + * array, so it is used only where the sums could not be computed at all otherwise. */ /** - * The mean of the values. + * The binary exponent of the magnitudes that the arrays are scaled to (400). + */ +const SCALED_MAGNITUDE_EXPONENT = 400 + +/** + * The smallest largest magnitude of an array that is used without scaling up (2^-400). Distinct + * doubles near the largest magnitude differ by at least 2^-53 of it, so unless all the values are + * equal, the largest deviation from their mean is at least about 2^-454 and its square, at least + * 2^-908, is computed exactly (`twoProduct` is exact above 2^-969). So the sum of squared deviations + * is 0 only when all the values are equal. + */ +const MIN_UNSCALED_MAGNITUDE = 2 ** -SCALED_MAGNITUDE_EXPONENT + +/** + * The largest magnitude of an array that is used without scaling down when the sums overflow (2^400). + * Every deviation is then at most 2^402, and every product of two deviations at most 2^804, so no + * sum can overflow. + */ +const MAX_UNSCALED_MAGNITUDE = 2 ** SCALED_MAGNITUDE_EXPONENT + +/** + * Deviations from the mean of an array of values, multiplied by `2^-exponent`. + */ +interface ScaledDeviations { + /** `(value - mean) * 2^-exponent` for each value, without rounding */ + readonly deviations: DoubleDouble[], + /** the exponent of the power of two the deviations are scaled by */ + readonly exponent: number, +} + +/** + * Sums of products of the deviations of two paired arrays, each array at its own scale. + */ +interface PairedSums { + /** the sums, as returned by the function that computed them */ + readonly sums: DoubleDouble[], + /** the deviations of the first array are multiplied by `2^-firstExponent` */ + readonly firstExponent: number, + /** the deviations of the second array are multiplied by `2^-secondExponent` */ + readonly secondExponent: number, +} + +/** + * The sums a simple linear regression of `y` on `x` is computed from. + */ +export interface RegressionSums { + /** `sum((x - mean(x))^2) * 2^(-2 * xExponent)` */ + readonly xSumOfSquares: DoubleDouble, + /** `sum((x - mean(x)) * (y - mean(y))) * 2^(-xExponent - yExponent)` */ + readonly productsSum: DoubleDouble, + /** `sum((y - mean(y))^2) * 2^(-2 * yExponent)`, when requested */ + readonly ySumOfSquares: DoubleDouble, + /** the exponent of the power of two the x deviations are scaled by */ + readonly xExponent: number, + /** the exponent of the power of two the y deviations are scaled by */ + readonly yExponent: number, +} + +/** + * The sum of the squared deviations of the values from their mean. * - * When the running total overflows, the values are summed scaled down by a power of two no smaller - * than their count, so that no partial sum can exceed the largest double, and the mean is scaled back. + * Only a tiny array is scaled: when the unscaled sum overflows, the sum itself exceeds the largest + * double. * * @param {number[]} values - a non-empty array of numbers - * @returns {DoubleDouble} the mean, without rounding to a double + * @returns {number} `sum((x - mean)^2)`, rounded once to a double */ -function mean(values: number[]): DoubleDouble { - const total = sum(values) - if (Number.isFinite(total.hi)) { - return divideDoubleDouble(total, values.length) - } - const scale = 2 ** Math.ceil(Math.log2(values.length)) - const scaledTotal = sum(values.map((value) => value / scale)) - return multiplyByPowerOfTwo(divideDoubleDouble(scaledTotal, values.length), scale) +export function sumOfSquaredDeviations(values: number[]): number { + const {deviations, exponent} = scaledDeviationsFromMean(values, false) + return multiplyByTwoToThe(roundDoubleDouble(sumOfProducts(deviations, deviations)), 2 * exponent) } /** - * The sum of the values. + * The covariance of two arrays: the sum of products of deviations divided by `n - deltaDegreesOfFreedom`. * - * @param {number[]} values - an array of numbers - * @returns {DoubleDouble} the sum, without rounding to a double + * The quotient is rounded at the scale of the deviations and scaled back, so it is infinite only when + * the covariance exceeds the largest double. + * + * @param {number[]} first - a non-empty array of numbers + * @param {number[]} second - an array of the same length as `first` + * @param {number} deltaDegreesOfFreedom - 0 for the population covariance, 1 for the sample covariance + * @returns {number} the covariance, rounded once to a double */ -function sum(values: number[]): DoubleDouble { - return values.reduce((total, value) => addDoubleDouble(total, {hi: value, lo: 0}), DOUBLE_DOUBLE_ZERO) +export function covariance(first: number[], second: number[], deltaDegreesOfFreedom: number): number { + const {sums: [productsSum], firstExponent, secondExponent} = pairedSums(first, second, (x, y) => [sumOfProducts(x, y)]) + const scaledCovariance = roundDoubleDouble(divideDoubleDouble(productsSum, first.length - deltaDegreesOfFreedom)) + return multiplyByTwoToThe(scaledCovariance, firstExponent + secondExponent) } /** - * The deviations of the values from their mean. + * The sums of squared deviations and of products of deviations of a simple linear regression, each + * array at its own scale. * - * @param {number[]} values - a non-empty array of numbers - * @returns {DoubleDouble[]} `value - mean` for each value, without rounding + * Unless all the x values are equal, the scaled sum of squares of x is at least about 2^-908, so the + * scaled slope `productsSum / xSumOfSquares` is finite whenever `ySumOfSquares` is. + * + * @param {number[]} knownYs - a non-empty array of the dependent values + * @param {number[]} knownXs - an array of the independent values, of the same length as `knownYs` + * @param {boolean} withYSumOfSquares - whether to compute `ySumOfSquares`; otherwise it is 0 + * @returns {RegressionSums} the scaled sums and the exponents of the scales */ -export function deviationsFromMean(values: number[]): DoubleDouble[] { - const center = mean(values) - return values.map((value) => subtractDoubleDouble({hi: value, lo: 0}, center)) +export function regressionSums(knownYs: number[], knownXs: number[], withYSumOfSquares: boolean): RegressionSums { + const {sums: [xSumOfSquares, productsSum, ySumOfSquares], firstExponent, secondExponent} = pairedSums(knownXs, knownYs, + (x, y) => withYSumOfSquares + ? [sumOfProducts(x, x), sumOfProducts(y, x), sumOfProducts(y, y)] + : [sumOfProducts(x, x), sumOfProducts(y, x), DOUBLE_DOUBLE_ZERO], + ) + return {xSumOfSquares, productsSum, ySumOfSquares, xExponent: firstExponent, yExponent: secondExponent} } /** - * The sum of the squared deviations of the values from their mean. + * Sums of products of the deviations of two paired arrays from their means. + * + * The sums are computed first with only tiny arrays scaled, which gives the same results as no scaling + * at all for arrays of ordinary magnitude. Only when a sum is not finite, they are computed again with + * the arrays whose largest magnitude exceeds `MAX_UNSCALED_MAGNITUDE` scaled down. + * + * @param {number[]} first - a non-empty array of numbers + * @param {number[]} second - an array of the same length as `first` + * @param {Function} computeSums - computes the sums from the scaled deviations of `first` and `second` + * @returns {PairedSums} the scaled sums and the exponents of the scales + */ +function pairedSums(first: number[], second: number[], computeSums: (first: DoubleDouble[], second: DoubleDouble[]) => DoubleDouble[]): PairedSums { + let firstDeviations = scaledDeviationsFromMean(first, false) + let secondDeviations = scaledDeviationsFromMean(second, false) + let sums = computeSums(firstDeviations.deviations, secondDeviations.deviations) + if (sums.some((sum) => !Number.isFinite(sum.hi))) { + firstDeviations = scaledDeviationsFromMean(first, true) + secondDeviations = scaledDeviationsFromMean(second, true) + sums = computeSums(firstDeviations.deviations, secondDeviations.deviations) + } + return {sums, firstExponent: firstDeviations.exponent, secondExponent: secondDeviations.exponent} +} + +/** + * The deviations of the values from their mean, measured at a power-of-two scale. + * + * The values are first multiplied by `2^-exponent`, which brings their largest magnitude up to + * `MIN_UNSCALED_MAGNITUDE` when it is below it, and, when `scaleDown` is set, down to + * `MAX_UNSCALED_MAGNITUDE` when it is above it. Scaling up is exact. Scaling down is exact except for + * values that become subnormal, at least 2^1400 times smaller than the largest one. * * @param {number[]} values - a non-empty array of numbers - * @returns {DoubleDouble} `sum((x - mean)^2)` + * @param {boolean} scaleDown - whether to scale down values above `MAX_UNSCALED_MAGNITUDE` + * @returns {ScaledDeviations} the scaled deviations, without rounding, and the exponent of the scale + */ +function scaledDeviationsFromMean(values: number[], scaleDown: boolean): ScaledDeviations { + const exponent = scalingExponent(largestMagnitude(values), scaleDown) + const scale = 2 ** -exponent + const scaledValues = exponent === 0 ? values : values.map((value) => value * scale) + const center = mean(scaledValues) + return { + deviations: scaledValues.map((value) => subtractDoubleDouble({hi: value, lo: 0}, center)), + exponent, + } +} + +/** + * The exponent `scaledDeviationsFromMean` scales the values by. + * + * @param {number} magnitude - the largest magnitude of the values + * @param {boolean} scaleDown - whether to scale down magnitudes above `MAX_UNSCALED_MAGNITUDE` + * @returns {number} an exponent between -674 and 624: 0 when the values are not scaled, otherwise the + * binary exponent of `magnitude` minus that of the magnitude it is scaled to */ -export function sumOfSquaredDeviations(values: number[]): DoubleDouble { - const deviations = deviationsFromMean(values) - return sumOfProducts(deviations, deviations) +function scalingExponent(magnitude: number, scaleDown: boolean): number { + if (magnitude > 0 && magnitude < MIN_UNSCALED_MAGNITUDE) { + return Math.floor(Math.log2(magnitude)) + SCALED_MAGNITUDE_EXPONENT + } + if (scaleDown && magnitude > MAX_UNSCALED_MAGNITUDE && Number.isFinite(magnitude)) { + return Math.floor(Math.log2(magnitude)) - SCALED_MAGNITUDE_EXPONENT + } + return 0 } /** - * The sum of the products of the paired deviations of two arrays from their means. + * The largest absolute value of the values. * - * @param {number[]} first - a non-empty array of numbers - * @param {number[]} second - an array of the same length as `first` - * @returns {DoubleDouble} `sum((x - mean(x)) * (y - mean(y)))` + * A plain loop, which is several times faster than `reduce` with `Math.max` on long arrays. + * + * @param {number[]} values - an array of numbers + * @returns {number} the largest `|value|`, or 0 for an empty array */ -export function sumOfProductsOfDeviations(first: number[], second: number[]): DoubleDouble { - return sumOfProducts(deviationsFromMean(first), deviationsFromMean(second)) +function largestMagnitude(values: number[]): number { + let largest = 0 + for (const value of values) { + const magnitude = Math.abs(value) + if (magnitude > largest) { + largest = magnitude + } + } + return largest } /** - * The variance of the values: the sum of squared deviations divided by `n - deltaDegreesOfFreedom`. + * The mean of the values. + * + * When the running total overflows, the values are summed scaled down by a power of two no smaller + * than their count, so that no partial sum can exceed the largest double, and the mean is scaled back. * * @param {number[]} values - a non-empty array of numbers - * @param {number} deltaDegreesOfFreedom - 0 for the population variance, 1 for the sample variance - * @returns {number} the variance, rounded once to a double + * @returns {DoubleDouble} the mean, without rounding to a double */ -export function variance(values: number[], deltaDegreesOfFreedom: number): number { - return roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations(values), values.length - deltaDegreesOfFreedom)) +function mean(values: number[]): DoubleDouble { + const total = sum(values) + if (Number.isFinite(total.hi)) { + return divideDoubleDouble(total, values.length) + } + const scale = 2 ** Math.ceil(Math.log2(values.length)) + const scaledTotal = sum(values.map((value) => value / scale)) + return multiplyByPowerOfTwo(divideDoubleDouble(scaledTotal, values.length), scale) } /** - * The covariance of two arrays: the sum of products of deviations divided by `n - deltaDegreesOfFreedom`. + * The sum of the values. * - * @param {number[]} first - a non-empty array of numbers - * @param {number[]} second - an array of the same length as `first` - * @param {number} deltaDegreesOfFreedom - 0 for the population covariance, 1 for the sample covariance - * @returns {number} the covariance, rounded once to a double + * @param {number[]} values - an array of numbers + * @returns {DoubleDouble} the sum, without rounding to a double */ -export function covariance(first: number[], second: number[], deltaDegreesOfFreedom: number): number { - return roundDoubleDouble(divideDoubleDouble(sumOfProductsOfDeviations(first, second), first.length - deltaDegreesOfFreedom)) +function sum(values: number[]): DoubleDouble { + return values.reduce((total, value) => addDoubleDouble(total, {hi: value, lo: 0}), DOUBLE_DOUBLE_ZERO) } /** @@ -115,6 +264,6 @@ export function covariance(first: number[], second: number[], deltaDegreesOfFree * @param {DoubleDouble[]} second - the second factors, of the same length * @returns {DoubleDouble} `sum(first[i] * second[i])` */ -export function sumOfProducts(first: DoubleDouble[], second: DoubleDouble[]): DoubleDouble { +function sumOfProducts(first: DoubleDouble[], second: DoubleDouble[]): DoubleDouble { return first.reduce((sum, deviation, index) => addDoubleDouble(sum, multiplyDoubleDouble(deviation, second[index])), DOUBLE_DOUBLE_ZERO) } diff --git a/src/interpreter/doubleDouble.ts b/src/interpreter/doubleDouble.ts index e7e5cf0555..c2497a1fb9 100644 --- a/src/interpreter/doubleDouble.ts +++ b/src/interpreter/doubleDouble.ts @@ -208,6 +208,22 @@ export function multiplyByPowerOfTwo(x: DoubleDouble, powerOfTwo: number): Doubl return {hi: x.hi * powerOfTwo, lo: x.lo * powerOfTwo} } +/** + * Product of a double and `2^exponent`, for exponents beyond the range of a finite power of two. + * + * The power is applied in two halves of the same sign, so that neither half overflows or underflows + * and an intermediate product overflows only when the result does. The product is exact unless it + * overflows or is subnormal; a subnormal product can be rounded twice. + * + * @param {number} value - the double factor + * @param {number} exponent - an integer, at most 2046 in magnitude + * @returns {number} `value * 2^exponent` + */ +export function multiplyByTwoToThe(value: number, exponent: number): number { + const half = Math.trunc(exponent / 2) + return value * 2 ** half * 2 ** (exponent - half) +} + /** * Rounds a double-double to the nearest double. * diff --git a/src/interpreter/plugin/DatabasePlugin.ts b/src/interpreter/plugin/DatabasePlugin.ts index 462266d6fc..b91c802b13 100644 --- a/src/interpreter/plugin/DatabasePlugin.ts +++ b/src/interpreter/plugin/DatabasePlugin.ts @@ -6,12 +6,12 @@ import {CellError, ErrorType} from '../../Cell' import {ErrorMessage} from '../../error-message' import {ProcedureAst} from '../../parser' -import {variance} from '../deviationSums' import {InterpreterState} from '../InterpreterState' import {EmptyValue, getRawValue, InternalScalarValue, InterpreterValue, isExtendedNumber, RawScalarValue} from '../InterpreterValue' import {SimpleRangeValue} from '../../SimpleRangeValue' import {FunctionArgumentType, FunctionPlugin, FunctionPluginTypecheck, ImplementedFunctions} from './FunctionPlugin' import {CriterionLambda} from '../Criterion' +import {MomentsAggregate} from './MomentsAggregate' /** * Parsed criterion for a single cell in the criteria range. @@ -338,11 +338,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return values } - if (values.length <= 1) { - return new CellError(ErrorType.DIV_BY_ZERO) - } - - return Math.sqrt(variance(values, 1)) + return MomentsAggregate.of(values).stdevSValue() ?? new CellError(ErrorType.DIV_BY_ZERO) }) } @@ -361,11 +357,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return values } - if (values.length === 0) { - return new CellError(ErrorType.DIV_BY_ZERO) - } - - return Math.sqrt(variance(values, 0)) + return MomentsAggregate.of(values).stdevPValue() ?? new CellError(ErrorType.DIV_BY_ZERO) }) } @@ -383,11 +375,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return values } - if (values.length <= 1) { - return new CellError(ErrorType.DIV_BY_ZERO) - } - - return variance(values, 1) + return MomentsAggregate.of(values).varSValue() ?? new CellError(ErrorType.DIV_BY_ZERO) }) } @@ -406,11 +394,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return values } - if (values.length === 0) { - return new CellError(ErrorType.DIV_BY_ZERO) - } - - return variance(values, 0) + return MomentsAggregate.of(values).varPValue() ?? new CellError(ErrorType.DIV_BY_ZERO) }) } diff --git a/src/interpreter/plugin/MomentsAggregate.ts b/src/interpreter/plugin/MomentsAggregate.ts new file mode 100644 index 0000000000..c475658936 --- /dev/null +++ b/src/interpreter/plugin/MomentsAggregate.ts @@ -0,0 +1,351 @@ +/** + * @license + * Copyright (c) 2025 Handsoncode. All rights reserved. + */ + +import {Maybe} from '../../Maybe' +import { + addDoubleDouble, + divideDoubleDouble, + DOUBLE_DOUBLE_ZERO, + DoubleDouble, + multiplyByPowerOfTwo, + multiplyByTwoToThe, + multiplyDoubleDouble, + roundDoubleDouble, + scaleDoubleDouble, + subtractDoubleDouble, + twoSum, +} from '../doubleDouble' + +/** + * The largest scaled sum of squared deviations a `MomentsAggregate` keeps (2^960). Above it, the + * aggregate raises its exponent. Every scaled deviation is then at most 2^480, so rebasing the sums of + * two aggregates cannot overflow. + */ +const MAX_SCALED_SUM_OF_SQUARES = 2 ** 960 + +/** + * The largest scaled difference of shifts that two aggregates are rebased with (2^470). Above it, the + * exponent is raised first. + */ +const MAX_SCALED_SHIFT_DIFFERENCE = 2 ** 470 + +/** + * The smallest non-zero scaled difference of shifts that is added to an aggregate whose sum of squares + * is 0 without lowering its exponent first (2^-450). Its square, at least 2^-900, is then computed + * exactly (`twoProduct` is exact above 2^-969). + */ +const MIN_SCALED_SHIFT_DIFFERENCE = 2 ** -450 + +/** + * The lowest exponent of an aggregate (-600). Scaled by `2^600`, even the smallest difference of two + * doubles, 2^-1074, has a square above 2^-969, while `2^600` and `2^-600` are finite. + */ +const MIN_EXPONENT = -600 + +/** + * The smallest exponent `e >= 0` (up to one) for which `magnitude * 2^-e` is at most `limit`. + * + * @param {number} magnitude - a non-negative finite number + * @param {number} limit - a positive power of two + * @returns {number} the exponent + */ +function exponentToFit(magnitude: number, limit: number): number { + return magnitude <= limit ? 0 : Math.ceil(Math.log2(magnitude / limit)) +} + +/** + * The exponent at which two shifts can be subtracted and their scaled difference squared: 0 when the + * difference is within [`MIN_SCALED_SHIFT_DIFFERENCE`, `MAX_SCALED_SHIFT_DIFFERENCE`], higher when it + * is larger, lower (down to `MIN_EXPONENT`) when it is smaller, and `-Infinity` when it is 0. + * + * @param {number} shift - the first shift + * @param {number} otherShift - the second shift + * @returns {number} the exponent + */ +function shiftDifferenceExponent(shift: number, otherShift: number): number { + const difference = Math.abs(otherShift - shift) + if (difference === 0) { + return -Infinity + } + if (difference < MIN_SCALED_SHIFT_DIFFERENCE) { + return Math.max(MIN_EXPONENT, Math.floor(Math.log2(difference))) + } + // halving the shifts first keeps their difference finite + return exponentToFit(Math.abs(otherShift / 2 - shift / 2), MAX_SCALED_SHIFT_DIFFERENCE / 2) +} + +/** + * Moments of a set of numbers, composable so that the value of a range can be cached and reused + * for a larger range. + * + * Besides the count, it keeps the sums of deviations `S1` and of squared deviations `S2` from a + * `shift`: the first value added. Both sums are double-double. Measured from a nearby value, the + * deviations stay small, which avoids the cancellation of the textbook one-pass form + * `sum(x^2) - sum(x)^2 / n` for data with a large mean and a small spread (for example 10000000.001, + * 10000000.002, ...). The variance and the standard deviation are accurate to about one unit in the + * last place of the exact values for the stored numbers. Double-double arithmetic does not guarantee + * correct rounding. + * + * The deviations are stored multiplied by `2^-exponent`, which is exact, so that the sums can neither + * overflow nor underflow. The exponent is 0 until the scaled sum of squares exceeds + * `MAX_SCALED_SUM_OF_SQUARES` (deviations of about 1e144 and above), and grows as needed. When all the + * values so far are equal and the next one differs from them by less than `MIN_SCALED_SHIFT_DIFFERENCE` + * (about 1e-135), the exponent is lowered instead, so that the squared difference does not underflow. + * An aggregate whose sum of squares is 0 has exponent 0. The variance and the standard deviation are + * computed at the scale of the sums and scaled back, so each is `#NUM!` only when it exceeds the + * largest double, and 0 only when it is below the smallest one, whatever the order of the values. + */ +export class MomentsAggregate { + + public static empty = new MomentsAggregate(0, 0, 0, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) + + /** + * @param {number} count - the number of values + * @param {number} shift - the value the deviations are measured from + * @param {number} exponent - the deviations are stored multiplied by `2^-exponent` + * @param {DoubleDouble} shiftedSum - `S1`, the sum of `(x - shift) * 2^-exponent` + * @param {DoubleDouble} shiftedSumOfSquares - `S2`, the sum of `((x - shift) * 2^-exponent)^2` + */ + constructor( + public readonly count: number, + public readonly shift: number, + public readonly exponent: number, + public readonly shiftedSum: DoubleDouble, + public readonly shiftedSumOfSquares: DoubleDouble, + ) { + } + + /** + * The moments of one value, which is also the shift. + * + * @param {number} arg - the value + * @returns {MomentsAggregate} an aggregate of the single value + */ + public static single(arg: number): MomentsAggregate { + return new MomentsAggregate(1, arg, 0, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) + } + + /** + * The moments of an array of numbers, added one by one in order, the same way as the values of a + * range are folded for VAR.S, VAR.P, STDEV.S and STDEV.P, so that the results are the same. + * + * @param {number[]} values - the values + * @returns {MomentsAggregate} an aggregate of the values + */ + public static of(values: number[]): MomentsAggregate { + return values.reduce((aggregate, value) => aggregate.compose(MomentsAggregate.single(value)), MomentsAggregate.empty) + } + + /** + * An aggregate whose exponent is raised, if needed, so that its scaled sum of squares is at most + * `MAX_SCALED_SUM_OF_SQUARES`. + * + * @param {number} count - the number of values + * @param {number} shift - the value the deviations are measured from + * @param {number} exponent - the exponent of the given sums + * @param {DoubleDouble} shiftedSum - the scaled sum of deviations + * @param {DoubleDouble} shiftedSumOfSquares - the scaled sum of squared deviations + * @returns {MomentsAggregate} the aggregate + */ + private static normalized(count: number, shift: number, exponent: number, shiftedSum: DoubleDouble, shiftedSumOfSquares: DoubleDouble): MomentsAggregate { + const excess = exponentToFit(Math.abs(shiftedSumOfSquares.hi), MAX_SCALED_SUM_OF_SQUARES) + if (excess === 0) { + return new MomentsAggregate(count, shift, exponent, shiftedSum, shiftedSumOfSquares) + } + const raise = Math.ceil(excess / 2) + return new MomentsAggregate(count, shift, exponent + raise, + multiplyByPowerOfTwo(shiftedSum, 2 ** -raise), + multiplyByPowerOfTwo(shiftedSumOfSquares, 2 ** (-2 * raise)), + ) + } + + /** + * Combines two aggregates. The result keeps this aggregate's `shift`; the other one's shifted sums + * are re-expressed relative to it, at the larger exponent of the two (raised further when the + * difference of the shifts needs it). The exponent of an aggregate whose sum of squares is 0 does not + * count, so the exponent is lowered when all the values are equal but for a tiny difference of the + * shifts. An empty aggregate is the identity: composing with it returns the other aggregate + * unchanged, with its own shift. + * + * @param {MomentsAggregate} other - the aggregate to add + * @returns {MomentsAggregate} the aggregate of the values of both + */ + public compose(other: MomentsAggregate): MomentsAggregate { + if (this.count === 0) { + return other + } + if (other.count === 0) { + return this + } + + // the common case: no rescaling, because both aggregates have the same exponent or the other one + // is a single value, whose shifted sums are 0 at any exponent + if (other.exponent === this.exponent || other.count === 1) { + const exponent = this.exponent + const shiftDifference = exponent === 0 + ? twoSum(other.shift, -this.shift) + : twoSum(other.shift * 2 ** -exponent, -this.shift * 2 ** -exponent) + const magnitude = Math.abs(shiftDifference.hi) + // a difference below MIN_SCALED_SHIFT_DIFFERENCE would lose precision when squared, which matters + // only when it is not 0 and no non-zero sum of squares outweighs it; the exponent is then lowered + const squaresPrecisely = magnitude >= MIN_SCALED_SHIFT_DIFFERENCE || magnitude === 0 + || this.shiftedSumOfSquares.hi !== 0 || other.shiftedSumOfSquares.hi !== 0 + if (magnitude <= MAX_SCALED_SHIFT_DIFFERENCE && squaresPrecisely) { + return this.composeScaled(other, exponent, this.shiftedSum, this.shiftedSumOfSquares, other.shiftedSum, other.shiftedSumOfSquares, shiftDifference) + } + } + + // an aggregate whose sum of squares is 0 does not constrain the exponent + const exponent = Math.max( + this.shiftedSumOfSquares.hi === 0 ? -Infinity : this.exponent, + other.shiftedSumOfSquares.hi === 0 ? -Infinity : other.exponent, + shiftDifferenceExponent(this.shift, other.shift), + ) + const thisScale = 2 ** (this.exponent - exponent) + const otherScale = 2 ** (other.exponent - exponent) + const scale = 2 ** -exponent + return this.composeScaled(other, exponent, + multiplyByPowerOfTwo(this.shiftedSum, thisScale), + multiplyByPowerOfTwo(multiplyByPowerOfTwo(this.shiftedSumOfSquares, thisScale), thisScale), + multiplyByPowerOfTwo(other.shiftedSum, otherScale), + multiplyByPowerOfTwo(multiplyByPowerOfTwo(other.shiftedSumOfSquares, otherScale), otherScale), + twoSum(other.shift * scale, -this.shift * scale), + ) + } + + /** + * The sample variance, as in VAR.S. + * + * @returns {Maybe} the variance, or `undefined` for fewer than two values + */ + public varSValue(): Maybe { + if (this.count > 1) { + return this.variance(this.count - 1) + } else { + return undefined + } + } + + /** + * The population variance, as in VAR.P. + * + * @returns {Maybe} the variance, or `undefined` for no values + */ + public varPValue(): Maybe { + if (this.count > 0) { + return this.variance(this.count) + } else { + return undefined + } + } + + /** + * The sample standard deviation, as in STDEV.S. + * + * @returns {Maybe} the standard deviation, or `undefined` for fewer than two values + */ + public stdevSValue(): Maybe { + if (this.count > 1) { + return this.standardDeviation(this.count - 1) + } else { + return undefined + } + } + + /** + * The population standard deviation, as in STDEV.P. + * + * @returns {Maybe} the standard deviation, or `undefined` for no values + */ + public stdevPValue(): Maybe { + if (this.count > 0) { + return this.standardDeviation(this.count) + } else { + return undefined + } + } + + /** + * Adds the other aggregate's sums, already scaled to `exponent`, to this aggregate's sums, also + * scaled to `exponent`. + * + * @param {MomentsAggregate} other - the aggregate to add + * @param {number} exponent - the common exponent + * @param {DoubleDouble} shiftedSum - this aggregate's `S1` at `exponent` + * @param {DoubleDouble} shiftedSumOfSquares - this aggregate's `S2` at `exponent` + * @param {DoubleDouble} otherShiftedSum - the other aggregate's `S1` at `exponent` + * @param {DoubleDouble} otherShiftedSumOfSquares - the other aggregate's `S2` at `exponent` + * @param {DoubleDouble} shiftDifference - `(other.shift - this.shift) * 2^-exponent` + * @returns {MomentsAggregate} the aggregate of the values of both + */ + private composeScaled( + other: MomentsAggregate, + exponent: number, + shiftedSum: DoubleDouble, + shiftedSumOfSquares: DoubleDouble, + otherShiftedSum: DoubleDouble, + otherShiftedSumOfSquares: DoubleDouble, + shiftDifference: DoubleDouble, + ): MomentsAggregate { + const count = this.count + other.count + if (other.count === 1) { + // the common case of adding one value: its deviation is the shift difference itself + return MomentsAggregate.normalized(count, this.shift, exponent, + addDoubleDouble(shiftedSum, shiftDifference), + addDoubleDouble(shiftedSumOfSquares, multiplyDoubleDouble(shiftDifference, shiftDifference)), + ) + } + + // other's sums rebased: S1 + n*d and S2 + 2*d*S1 + n*d^2 + const rebasedSum = addDoubleDouble(otherShiftedSum, scaleDoubleDouble(shiftDifference, other.count)) + const rebasedSumOfSquares = addDoubleDouble( + addDoubleDouble(otherShiftedSumOfSquares, scaleDoubleDouble(multiplyDoubleDouble(shiftDifference, otherShiftedSum), 2)), + scaleDoubleDouble(multiplyDoubleDouble(shiftDifference, shiftDifference), other.count), + ) + return MomentsAggregate.normalized(count, this.shift, exponent, + addDoubleDouble(shiftedSum, rebasedSum), + addDoubleDouble(shiftedSumOfSquares, rebasedSumOfSquares), + ) + } + + /** + * The variance: the sum of squared deviations from the mean divided by `divisor`. + * + * The scaled variance is multiplied by `2^(2 * exponent)` in two steps, so that the power of two + * itself cannot overflow; the result overflows only when the variance exceeds the largest double. + * + * @param {number} divisor - `n - 1` for the sample variance, `n` for the population variance + * @returns {number} the variance + */ + private variance(divisor: number): number { + return multiplyByTwoToThe(this.scaledVariance(divisor), 2 * this.exponent) + } + + /** + * The standard deviation: the square root of the variance, taken at the scale of the sums, so that + * it is finite and non-zero whenever the exact standard deviation is, even when the variance is not. + * + * @param {number} divisor - `n - 1` for the sample standard deviation, `n` for the population one + * @returns {number} the standard deviation + */ + private standardDeviation(divisor: number): number { + return Math.sqrt(this.scaledVariance(divisor)) * 2 ** this.exponent + } + + /** + * The sum of squared deviations from the mean at the scale of the sums, `S2 - S1^2 / n` with `S1`, + * `S2` the shifted sums and `n` the count, divided by `divisor` and rounded once. + * + * `S1^2 / n` is evaluated as `S1 * (S1 / n)`, which cannot overflow because `S1^2 / n <= S2`. + * + * @param {number} divisor - `n - 1` for the sample variance, `n` for the population variance + * @returns {number} the variance multiplied by `2^(-2 * exponent)` + */ + private scaledVariance(divisor: number): number { + const squaredSumOverCount = multiplyDoubleDouble(this.shiftedSum, divideDoubleDouble(this.shiftedSum, this.count)) + const sumOfSquaredDeviations = subtractDoubleDouble(this.shiftedSumOfSquares, squaredSumOverCount) + return roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations, divisor)) + } +} diff --git a/src/interpreter/plugin/NumericAggregationPlugin.ts b/src/interpreter/plugin/NumericAggregationPlugin.ts index 956c456a51..6f426c6080 100644 --- a/src/interpreter/plugin/NumericAggregationPlugin.ts +++ b/src/interpreter/plugin/NumericAggregationPlugin.ts @@ -11,22 +11,11 @@ import {Maybe} from '../../Maybe' import {Ast, AstNodeType, CellRangeAst, ProcedureAst} from '../../parser' import {ColumnRangeAst, RowRangeAst} from '../../parser/Ast' import {coerceBooleanToNumber} from '../ArithmeticHelper' -import { - addDoubleDouble, - divideDoubleDouble, - DOUBLE_DOUBLE_ZERO, - DoubleDouble, - multiplyByPowerOfTwo, - multiplyDoubleDouble, - roundDoubleDouble, - scaleDoubleDouble, - subtractDoubleDouble, - twoSum, -} from '../doubleDouble' import {InterpreterState} from '../InterpreterState' import {EmptyValue, ExtendedNumber, getRawValue, InternalScalarValue, isExtendedNumber} from '../InterpreterValue' import {SimpleRangeValue} from '../../SimpleRangeValue' import {AverageResult} from './AverageResult' +import {MomentsAggregate} from './MomentsAggregate' import {FunctionArgumentType, FunctionPlugin, FunctionPluginTypecheck, ImplementedFunctions} from './FunctionPlugin' import {RangeVertex} from '../../DependencyGraph' @@ -44,223 +33,6 @@ function zeroForInfinite(value: InternalScalarValue) { } } -/** - * The largest scaled sum of squared deviations a `MomentsAggregate` keeps (2^960). Above it, the - * aggregate raises its exponent. Every scaled deviation is then at most 2^480, so rebasing the sums of - * two aggregates cannot overflow. - */ -const MAX_SCALED_SUM_OF_SQUARES = 2 ** 960 - -/** - * The largest scaled difference of shifts that two aggregates are rebased with (2^470). Above it, the - * exponent is raised first. - */ -const MAX_SCALED_SHIFT_DIFFERENCE = 2 ** 470 - -/** - * The smallest exponent `e >= 0` (up to one) for which `magnitude * 2^-e` is at most `limit`. - * - * @param {number} magnitude - a non-negative finite number - * @param {number} limit - a positive power of two - * @returns {number} the exponent - */ -function exponentToFit(magnitude: number, limit: number): number { - return magnitude <= limit ? 0 : Math.ceil(Math.log2(magnitude / limit)) -} - -/** - * Moments of a set of numbers, composable so that the value of a range can be cached and reused - * for a larger range. - * - * Besides the count, it keeps the sums of deviations `S1` and of squared deviations `S2` from a - * `shift`: the first value added. Both sums are double-double. Measured from a nearby value, the - * deviations stay small, which avoids the cancellation of the textbook one-pass form - * `sum(x^2) - sum(x)^2 / n` for data with a large mean and a small spread (for example 10000000.001, - * 10000000.002, ...). The variance and the standard deviation are accurate to about one unit in the - * last place of the exact values for the stored numbers. Double-double arithmetic does not guarantee - * correct rounding. - * - * The deviations are stored multiplied by `2^-exponent`, which is exact, so that the sums cannot - * overflow. The exponent is 0 until the scaled sum of squares exceeds `MAX_SCALED_SUM_OF_SQUARES` - * (deviations of about 1e144 and above), and grows as needed. The variance is `#NUM!` only when it - * exceeds the largest double, whatever the order of the values. - */ -class MomentsAggregate { - - public static empty = new MomentsAggregate(0, 0, 0, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) - - /** - * @param {number} count - the number of values - * @param {number} shift - the value the deviations are measured from - * @param {number} exponent - the deviations are stored multiplied by `2^-exponent` - * @param {DoubleDouble} shiftedSum - `S1`, the sum of `(x - shift) * 2^-exponent` - * @param {DoubleDouble} shiftedSumOfSquares - `S2`, the sum of `((x - shift) * 2^-exponent)^2` - */ - constructor( - public readonly count: number, - public readonly shift: number, - public readonly exponent: number, - public readonly shiftedSum: DoubleDouble, - public readonly shiftedSumOfSquares: DoubleDouble, - ) { - } - - /** - * The moments of one value, which is also the shift. - * - * @param {number} arg - the value - * @returns {MomentsAggregate} an aggregate of the single value - */ - public static single(arg: number): MomentsAggregate { - return new MomentsAggregate(1, arg, 0, DOUBLE_DOUBLE_ZERO, DOUBLE_DOUBLE_ZERO) - } - - /** - * An aggregate whose exponent is raised, if needed, so that its scaled sum of squares is at most - * `MAX_SCALED_SUM_OF_SQUARES`. - * - * @param {number} count - the number of values - * @param {number} shift - the value the deviations are measured from - * @param {number} exponent - the exponent of the given sums - * @param {DoubleDouble} shiftedSum - the scaled sum of deviations - * @param {DoubleDouble} shiftedSumOfSquares - the scaled sum of squared deviations - * @returns {MomentsAggregate} the aggregate - */ - private static normalized(count: number, shift: number, exponent: number, shiftedSum: DoubleDouble, shiftedSumOfSquares: DoubleDouble): MomentsAggregate { - const excess = exponentToFit(Math.abs(shiftedSumOfSquares.hi), MAX_SCALED_SUM_OF_SQUARES) - if (excess === 0) { - return new MomentsAggregate(count, shift, exponent, shiftedSum, shiftedSumOfSquares) - } - const raise = Math.ceil(excess / 2) - return new MomentsAggregate(count, shift, exponent + raise, - multiplyByPowerOfTwo(shiftedSum, 2 ** -raise), - multiplyByPowerOfTwo(shiftedSumOfSquares, 2 ** (-2 * raise)), - ) - } - - /** - * Combines two aggregates. The result keeps this aggregate's `shift`; the other one's shifted sums - * are re-expressed relative to it, at the larger exponent of the two (raised further when the - * difference of the shifts needs it). An empty aggregate is the identity: composing with it returns - * the other aggregate unchanged, with its own shift. - * - * @param {MomentsAggregate} other - the aggregate to add - * @returns {MomentsAggregate} the aggregate of the values of both - */ - public compose(other: MomentsAggregate): MomentsAggregate { - if (this.count === 0) { - return other - } - if (other.count === 0) { - return this - } - - // the common case: no rescaling, because both aggregates have the same exponent or the other one - // is a single value, whose shifted sums are 0 at any exponent - if (other.exponent === this.exponent || (other.count === 1 && this.exponent > 0)) { - const exponent = this.exponent - const shiftDifference = exponent === 0 - ? twoSum(other.shift, -this.shift) - : twoSum(other.shift * 2 ** -exponent, -this.shift * 2 ** -exponent) - if (Math.abs(shiftDifference.hi) <= MAX_SCALED_SHIFT_DIFFERENCE) { - return this.composeScaled(other, exponent, this.shiftedSum, this.shiftedSumOfSquares, other.shiftedSum, other.shiftedSumOfSquares, shiftDifference) - } - } - - // halving the shifts first keeps their difference finite - const shiftDifferenceExponent = exponentToFit(Math.abs(other.shift / 2 - this.shift / 2), MAX_SCALED_SHIFT_DIFFERENCE / 2) - const exponent = Math.max(this.exponent, other.exponent, shiftDifferenceExponent) - const thisRaise = exponent - this.exponent - const otherRaise = exponent - other.exponent - const scale = 2 ** -exponent - return this.composeScaled(other, exponent, - multiplyByPowerOfTwo(this.shiftedSum, 2 ** -thisRaise), - multiplyByPowerOfTwo(this.shiftedSumOfSquares, 2 ** (-2 * thisRaise)), - multiplyByPowerOfTwo(other.shiftedSum, 2 ** -otherRaise), - multiplyByPowerOfTwo(other.shiftedSumOfSquares, 2 ** (-2 * otherRaise)), - twoSum(other.shift * scale, -this.shift * scale), - ) - } - - public varSValue(): Maybe { - if (this.count > 1) { - return this.variance(this.count - 1) - } else { - return undefined - } - } - - public varPValue(): Maybe { - if (this.count > 0) { - return this.variance(this.count) - } else { - return undefined - } - } - - /** - * Adds the other aggregate's sums, already scaled to `exponent`, to this aggregate's sums, also - * scaled to `exponent`. - * - * @param {MomentsAggregate} other - the aggregate to add - * @param {number} exponent - the common exponent - * @param {DoubleDouble} shiftedSum - this aggregate's `S1` at `exponent` - * @param {DoubleDouble} shiftedSumOfSquares - this aggregate's `S2` at `exponent` - * @param {DoubleDouble} otherShiftedSum - the other aggregate's `S1` at `exponent` - * @param {DoubleDouble} otherShiftedSumOfSquares - the other aggregate's `S2` at `exponent` - * @param {DoubleDouble} shiftDifference - `(other.shift - this.shift) * 2^-exponent` - * @returns {MomentsAggregate} the aggregate of the values of both - */ - private composeScaled( - other: MomentsAggregate, - exponent: number, - shiftedSum: DoubleDouble, - shiftedSumOfSquares: DoubleDouble, - otherShiftedSum: DoubleDouble, - otherShiftedSumOfSquares: DoubleDouble, - shiftDifference: DoubleDouble, - ): MomentsAggregate { - const count = this.count + other.count - if (other.count === 1) { - // the common case of adding one value: its deviation is the shift difference itself - return MomentsAggregate.normalized(count, this.shift, exponent, - addDoubleDouble(shiftedSum, shiftDifference), - addDoubleDouble(shiftedSumOfSquares, multiplyDoubleDouble(shiftDifference, shiftDifference)), - ) - } - - // other's sums rebased: S1 + n*d and S2 + 2*d*S1 + n*d^2 - const rebasedSum = addDoubleDouble(otherShiftedSum, scaleDoubleDouble(shiftDifference, other.count)) - const rebasedSumOfSquares = addDoubleDouble( - addDoubleDouble(otherShiftedSumOfSquares, scaleDoubleDouble(multiplyDoubleDouble(shiftDifference, otherShiftedSum), 2)), - scaleDoubleDouble(multiplyDoubleDouble(shiftDifference, shiftDifference), other.count), - ) - return MomentsAggregate.normalized(count, this.shift, exponent, - addDoubleDouble(shiftedSum, rebasedSum), - addDoubleDouble(shiftedSumOfSquares, rebasedSumOfSquares), - ) - } - - /** - * The sum of squared deviations from the mean, `S2 - S1^2 / n` with `S1`, `S2` the shifted sums and - * `n` the count, divided by `divisor` and rounded once. - * - * `S1^2 / n` is evaluated as `S1 * (S1 / n)`, which cannot overflow because `S1^2 / n <= S2`. The - * quotient is computed at the scale of the sums and then multiplied by `2^(2 * exponent)` in two - * steps, so that the power of two itself cannot overflow; the result overflows only when the - * variance exceeds the largest double. - * - * @param {number} divisor - `n - 1` for the sample variance, `n` for the population variance - * @returns {number} the variance - */ - private variance(divisor: number): number { - const squaredSumOverCount = multiplyDoubleDouble(this.shiftedSum, divideDoubleDouble(this.shiftedSum, this.count)) - const sumOfSquaredDeviations = subtractDoubleDouble(this.shiftedSumOfSquares, squaredSumOverCount) - const scale = 2 ** this.exponent - return roundDoubleDouble(divideDoubleDouble(sumOfSquaredDeviations, divisor)) * scale * scale - } -} - export class NumericAggregationPlugin extends FunctionPlugin implements FunctionPluginTypecheck { public static implementedFunctions: ImplementedFunctions = { 'SUM': { @@ -529,8 +301,7 @@ export class NumericAggregationPlugin extends FunctionPlugin implements Function if (result instanceof CellError) { return result } else { - const val = result.varSValue() - return val === undefined ? new CellError(ErrorType.DIV_BY_ZERO) : Math.sqrt(val) + return result.stdevSValue() ?? new CellError(ErrorType.DIV_BY_ZERO) } } @@ -540,8 +311,7 @@ export class NumericAggregationPlugin extends FunctionPlugin implements Function if (result instanceof CellError) { return result } else { - const val = result.varPValue() - return val === undefined ? new CellError(ErrorType.DIV_BY_ZERO) : Math.sqrt(val) + return result.stdevPValue() ?? new CellError(ErrorType.DIV_BY_ZERO) } } @@ -670,8 +440,7 @@ export class NumericAggregationPlugin extends FunctionPlugin implements Function if (result instanceof CellError) { return result } else { - const val = result.varSValue() - return val === undefined ? new CellError(ErrorType.DIV_BY_ZERO) : Math.sqrt(val) + return result.stdevSValue() ?? new CellError(ErrorType.DIV_BY_ZERO) } } @@ -681,8 +450,7 @@ export class NumericAggregationPlugin extends FunctionPlugin implements Function if (result instanceof CellError) { return result } else { - const val = result.varPValue() - return val === undefined ? new CellError(ErrorType.DIV_BY_ZERO) : Math.sqrt(val) + return result.stdevPValue() ?? new CellError(ErrorType.DIV_BY_ZERO) } } diff --git a/src/interpreter/plugin/StatisticalAggregationPlugin.ts b/src/interpreter/plugin/StatisticalAggregationPlugin.ts index f81f12ab71..eb9223d43b 100644 --- a/src/interpreter/plugin/StatisticalAggregationPlugin.ts +++ b/src/interpreter/plugin/StatisticalAggregationPlugin.ts @@ -27,10 +27,10 @@ import { sumsqerr, variance } from './3rdparty/jstat/jstat' -import {covariance, deviationsFromMean, sumOfProducts, sumOfSquaredDeviations} from '../deviationSums' +import {covariance, regressionSums, sumOfSquaredDeviations} from '../deviationSums' import { divideDoubleDoubles, - DoubleDouble, + multiplyByTwoToThe, multiplyDoubleDouble, roundDoubleDouble, subtractDoubleDouble, @@ -192,7 +192,7 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (coerced.length === 0) { return 0 } - return roundDoubleDouble(sumOfSquaredDeviations(coerced)) + return sumOfSquaredDeviations(coerced) }) } @@ -377,24 +377,15 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 2) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.ThreeValues) } - const xMoments = xDeviationsAndSumOfSquares(ret[1]) - if (xMoments instanceof CellError) { - return xMoments - } - const {xDeviations, xSumOfSquares} = xMoments - const yDeviations = deviationsFromMean(ret[0]) - const productsSum = sumOfProducts(yDeviations, xDeviations) - const slope = divideDoubleDoubles(productsSum, xSumOfSquares) - // when the slope overflows, the explained sum of squares Sxy^2 / Sxx can still be finite - const explained = Number.isFinite(slope.hi) - ? multiplyDoubleDouble(slope, productsSum) - : divideDoubleDoubles(multiplyDoubleDouble(productsSum, productsSum), xSumOfSquares) - const residual = roundDoubleDouble(subtractDoubleDouble(sumOfProducts(yDeviations, yDeviations), explained)) - if (!Number.isFinite(residual)) { - return new CellError(ErrorType.NUM, ErrorMessage.NaN) + const sums = regressionSums(ret[0], ret[1], true) + if (sums.xSumOfSquares.hi === 0) { + return new CellError(ErrorType.DIV_BY_ZERO) } + // at the scale of the sums, the slope and the explained sum of squares Sxy^2 / Sxx <= Syy are finite + const explained = multiplyDoubleDouble(divideDoubleDoubles(sums.productsSum, sums.xSumOfSquares), sums.productsSum) + const residual = roundDoubleDouble(subtractDoubleDouble(sums.ySumOfSquares, explained)) // the residual is non-negative; a tiny negative rounding remainder counts as 0 - return Math.sqrt(Math.max(0, residual) / (n - 2)) + return multiplyByTwoToThe(Math.sqrt(Math.max(0, residual) / (n - 2)), sums.yExponent) }) } @@ -413,12 +404,12 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 1) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.TwoValues) } - const xMoments = xDeviationsAndSumOfSquares(ret[1]) - if (xMoments instanceof CellError) { - return xMoments + const sums = regressionSums(ret[0], ret[1], false) + if (sums.xSumOfSquares.hi === 0) { + return new CellError(ErrorType.DIV_BY_ZERO) } - const {xDeviations, xSumOfSquares} = xMoments - return roundDoubleDouble(divideDoubleDoubles(sumOfProducts(deviationsFromMean(ret[0]), xDeviations), xSumOfSquares)) + const scaledSlope = roundDoubleDouble(divideDoubleDoubles(sums.productsSum, sums.xSumOfSquares)) + return multiplyByTwoToThe(scaledSlope, sums.yExponent - sums.xExponent) }) } @@ -549,25 +540,6 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func } } -/** - * The deviations of the x values of SLOPE or STEYX from their mean and the sum of their squares. - * - * @param {number[]} values - the x values - * @returns {CellError | {xDeviations: DoubleDouble[], xSumOfSquares: DoubleDouble}} `#DIV/0!` when all the - * x values are equal, `#NUM!` when the sum of squares overflows - */ -function xDeviationsAndSumOfSquares(values: number[]): CellError | {xDeviations: DoubleDouble[], xSumOfSquares: DoubleDouble} { - const xDeviations = deviationsFromMean(values) - const xSumOfSquares = sumOfProducts(xDeviations, xDeviations) - if (xSumOfSquares.hi === 0) { - return new CellError(ErrorType.DIV_BY_ZERO) - } - if (!Number.isFinite(roundDoubleDouble(xSumOfSquares))) { - return new CellError(ErrorType.NUM, ErrorMessage.NaN) - } - return {xDeviations, xSumOfSquares} -} - function parseTwoArrays(dataX: SimpleRangeValue, dataY: SimpleRangeValue): CellError | [number[], number[]] { const xit = dataX.iterateValuesFromTopLeftCorner() const yit = dataY.iterateValuesFromTopLeftCorner() From 2b3c8abf7ea1054cd7610b8199dc08273097e40e Mon Sep 17 00:00:00 2001 From: Kuba Sekowski Date: Thu, 8 Oct 2026 08:49:45 +0200 Subject: [PATCH 10/13] Add a message to the equal-x #DIV/0! of SLOPE and STEYX and fix the docs of the variance changes (HF-491) SLOPE and STEYX share the equal-x check in one helper that returns #DIV/0! with ErrorMessage.EqualXValues. AverageResult gets JSDoc, the JSDoc of MomentsAggregate.of no longer claims the same results as a cached or split range, and the changelog lists VAR and the SUBTOTAL modes among the functions fixed for very large and very small values. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 2 +- src/error-message.ts | 1 + src/interpreter/plugin/AverageResult.ts | 21 +++++++++++++ src/interpreter/plugin/MomentsAggregate.ts | 7 +++-- .../plugin/StatisticalAggregationPlugin.ts | 31 ++++++++++++++----- 5 files changed, 52 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e388c366f9..9f8e8e3917 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,7 +25,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Fixed the `VAR`, `STDEV`, `DEVSQ`, `COVARIANCE`, `SLOPE`, `STEYX`, `DVAR` and `DSTDEV` functions, their variants, and the matching `SUBTOTAL` modes losing precision on data with a large mean and a small spread. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed a bug where `SLOPE` and `STEYX` returned `#NUM!` or an arbitrary number instead of `#DIV/0!` when all the x values are equal. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed a bug where `STEYX` returned `#NUM!` instead of `0` for points that lie exactly on a line. [#1784](https://github.com/handsontable/hyperformula/pull/1784) -- Fixed a bug where `STDEV`, `DVAR`, `DSTDEV`, `COVARIANCE`, `SLOPE`, `STEYX` and their variants returned an error or `0` for very large or very small values although the result was within the range of numbers. [#1784](https://github.com/handsontable/hyperformula/pull/1784) +- Fixed a bug where `VAR`, `STDEV`, `DVAR`, `DSTDEV`, `COVARIANCE`, `SLOPE`, `STEYX`, their variants, and the matching `SUBTOTAL` modes returned an error or `0` for very large or very small values although the result was within the range of numbers. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed the `MOD` function returning a remainder with the sign of the dividend instead of the sign of the divisor, which made the results differ from Excel and Google Sheets for arguments with opposite signs (e.g. `=MOD(-3, 12)` now returns `9` instead of `-3`). [#1747](https://github.com/handsontable/hyperformula/issues/1747) - Fixed a bug where moving or pasting a formula with an undefined name to another sheet incorrectly added an empty global named expression. [#1728](https://github.com/handsontable/hyperformula/pull/1728) diff --git a/src/error-message.ts b/src/error-message.ts index 36bc5231af..d42040901c 100644 --- a/src/error-message.ts +++ b/src/error-message.ts @@ -48,6 +48,7 @@ export class ErrorMessage { public static OneValue = 'Needs at least one value.' public static TwoValues = 'Range needs to contain at least two elements.' public static ThreeValues = 'Range needs to contain at least three elements.' + public static EqualXValues = 'All x values are equal.' public static IndexBounds = 'Index out of bounds.' public static IndexLarge = 'Index too large.' public static Formula = 'Expected formula.' diff --git a/src/interpreter/plugin/AverageResult.ts b/src/interpreter/plugin/AverageResult.ts index 86351f9625..fefceae88d 100644 --- a/src/interpreter/plugin/AverageResult.ts +++ b/src/interpreter/plugin/AverageResult.ts @@ -12,19 +12,40 @@ import {Maybe} from '../../Maybe' export class AverageResult { public static empty = new AverageResult(0, 0) + /** + * @param {number} sum - the sum of the values + * @param {number} count - the number of values + */ constructor( public readonly sum: number, public readonly count: number, ) {} + /** + * The sum and the count of one value. + * + * @param {number} arg - the value + * @returns {AverageResult} an aggregate of the single value + */ public static single(arg: number): AverageResult { return new AverageResult(arg, 1) } + /** + * Combines two aggregates by adding their sums and their counts. + * + * @param {AverageResult} other - the aggregate to add + * @returns {AverageResult} the aggregate of the values of both + */ public compose(other: AverageResult) { return new AverageResult(this.sum + other.sum, this.count + other.count) } + /** + * The average: the sum divided by the count. + * + * @returns {Maybe} the average, or `undefined` for no values + */ public averageValue(): Maybe { if (this.count > 0) { return this.sum / this.count diff --git a/src/interpreter/plugin/MomentsAggregate.ts b/src/interpreter/plugin/MomentsAggregate.ts index c475658936..b962c3b485 100644 --- a/src/interpreter/plugin/MomentsAggregate.ts +++ b/src/interpreter/plugin/MomentsAggregate.ts @@ -128,8 +128,11 @@ export class MomentsAggregate { } /** - * The moments of an array of numbers, added one by one in order, the same way as the values of a - * range are folded for VAR.S, VAR.P, STDEV.S and STDEV.P, so that the results are the same. + * The moments of an array of numbers, added one by one in order, with the first value as the shift. + * + * VAR.S, VAR.P, STDEV.S and STDEV.P fold a single uncached range the same way. When they reuse the + * cached aggregate of a smaller range or get several arguments, they compose partial aggregates in a + * different order, so their results can differ from this one in the last bits. * * @param {number[]} values - the values * @returns {MomentsAggregate} an aggregate of the values diff --git a/src/interpreter/plugin/StatisticalAggregationPlugin.ts b/src/interpreter/plugin/StatisticalAggregationPlugin.ts index eb9223d43b..e7ad6b921c 100644 --- a/src/interpreter/plugin/StatisticalAggregationPlugin.ts +++ b/src/interpreter/plugin/StatisticalAggregationPlugin.ts @@ -27,7 +27,7 @@ import { sumsqerr, variance } from './3rdparty/jstat/jstat' -import {covariance, regressionSums, sumOfSquaredDeviations} from '../deviationSums' +import {covariance, regressionSums, RegressionSums, sumOfSquaredDeviations} from '../deviationSums' import { divideDoubleDoubles, multiplyByTwoToThe, @@ -377,9 +377,9 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 2) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.ThreeValues) } - const sums = regressionSums(ret[0], ret[1], true) - if (sums.xSumOfSquares.hi === 0) { - return new CellError(ErrorType.DIV_BY_ZERO) + const sums = nonDegenerateRegressionSums(ret[0], ret[1], true) + if (sums instanceof CellError) { + return sums } // at the scale of the sums, the slope and the explained sum of squares Sxy^2 / Sxx <= Syy are finite const explained = multiplyDoubleDouble(divideDoubleDoubles(sums.productsSum, sums.xSumOfSquares), sums.productsSum) @@ -404,9 +404,9 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n <= 1) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.TwoValues) } - const sums = regressionSums(ret[0], ret[1], false) - if (sums.xSumOfSquares.hi === 0) { - return new CellError(ErrorType.DIV_BY_ZERO) + const sums = nonDegenerateRegressionSums(ret[0], ret[1], false) + if (sums instanceof CellError) { + return sums } const scaledSlope = roundDoubleDouble(divideDoubleDoubles(sums.productsSum, sums.xSumOfSquares)) return multiplyByTwoToThe(scaledSlope, sums.yExponent - sums.xExponent) @@ -540,6 +540,23 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func } } +/** + * The sums of a simple linear regression of `knownYs` on `knownXs`, or `#DIV/0!` when all the x values + * are equal, so that the slope is undefined. + * + * @param {number[]} knownYs - a non-empty array of the dependent values + * @param {number[]} knownXs - an array of the independent values, of the same length as `knownYs` + * @param {boolean} withYSumOfSquares - whether to compute `ySumOfSquares` + * @returns {RegressionSums | CellError} the scaled sums, or the error + */ +function nonDegenerateRegressionSums(knownYs: number[], knownXs: number[], withYSumOfSquares: boolean): RegressionSums | CellError { + const sums = regressionSums(knownYs, knownXs, withYSumOfSquares) + if (sums.xSumOfSquares.hi === 0) { + return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.EqualXValues) + } + return sums +} + function parseTwoArrays(dataX: SimpleRangeValue, dataY: SimpleRangeValue): CellError | [number[], number[]] { const xit = dataX.iterateValuesFromTopLeftCorner() const yit = dataY.iterateValuesFromTopLeftCorner() From e0d1e7d57708db4824d8c2f9a4ec557c5b95dcd9 Mon Sep 17 00:00:00 2001 From: Kuba Sekowski Date: Thu, 8 Oct 2026 09:22:38 +0200 Subject: [PATCH 11/13] Compute the STEYX residual sum of squares from the residuals themselves (HF-491) STEYX computed the residual sum of squares as Syy - Sxy^2 / Sxx, which cancels when the points lie almost on a line: for y = {1e8, 2e8, 300000000.00000024} and x = {1, 2, 3} it was 0.8% off. The residuals are now computed one by one in double-double and centered on their mean, which removes the shift that the rounding errors of the x and y means add to every residual. Co-Authored-By: Claude Opus 5.5 --- src/interpreter/deviationSums.ts | 57 +++++++++++++++---- .../plugin/StatisticalAggregationPlugin.ts | 14 ++--- 2 files changed, 51 insertions(+), 20 deletions(-) diff --git a/src/interpreter/deviationSums.ts b/src/interpreter/deviationSums.ts index e29eccd87f..f5c16c10ed 100644 --- a/src/interpreter/deviationSums.ts +++ b/src/interpreter/deviationSums.ts @@ -6,6 +6,7 @@ import { addDoubleDouble, divideDoubleDouble, + divideDoubleDoubles, DOUBLE_DOUBLE_ZERO, DoubleDouble, multiplyByPowerOfTwo, @@ -84,8 +85,11 @@ export interface RegressionSums { readonly xSumOfSquares: DoubleDouble, /** `sum((x - mean(x)) * (y - mean(y))) * 2^(-xExponent - yExponent)` */ readonly productsSum: DoubleDouble, - /** `sum((y - mean(y))^2) * 2^(-2 * yExponent)`, when requested */ - readonly ySumOfSquares: DoubleDouble, + /** + * `sum((y - mean(y) - slope * (x - mean(x)))^2) * 2^(-2 * yExponent)`, the residual sum of squares, + * when requested + */ + readonly residualSumOfSquares: DoubleDouble, /** the exponent of the power of two the x deviations are scaled by */ readonly xExponent: number, /** the exponent of the power of two the y deviations are scaled by */ @@ -128,20 +132,53 @@ export function covariance(first: number[], second: number[], deltaDegreesOfFree * array at its own scale. * * Unless all the x values are equal, the scaled sum of squares of x is at least about 2^-908, so the - * scaled slope `productsSum / xSumOfSquares` is finite whenever `ySumOfSquares` is. + * scaled slope `productsSum / xSumOfSquares` is finite whenever the sum of squared y deviations is, + * and so is every residual, whose square is at most the residual sum of squares. * * @param {number[]} knownYs - a non-empty array of the dependent values * @param {number[]} knownXs - an array of the independent values, of the same length as `knownYs` - * @param {boolean} withYSumOfSquares - whether to compute `ySumOfSquares`; otherwise it is 0 + * @param {boolean} withResidualSumOfSquares - whether to compute `residualSumOfSquares`; otherwise it is 0 * @returns {RegressionSums} the scaled sums and the exponents of the scales */ -export function regressionSums(knownYs: number[], knownXs: number[], withYSumOfSquares: boolean): RegressionSums { - const {sums: [xSumOfSquares, productsSum, ySumOfSquares], firstExponent, secondExponent} = pairedSums(knownXs, knownYs, - (x, y) => withYSumOfSquares - ? [sumOfProducts(x, x), sumOfProducts(y, x), sumOfProducts(y, y)] - : [sumOfProducts(x, x), sumOfProducts(y, x), DOUBLE_DOUBLE_ZERO], +export function regressionSums(knownYs: number[], knownXs: number[], withResidualSumOfSquares: boolean): RegressionSums { + const {sums: [xSumOfSquares, productsSum, residualSumOfSquares], firstExponent, secondExponent} = pairedSums(knownXs, knownYs, + (x, y) => { + const xSquaresSum = sumOfProducts(x, x) + const xyProductsSum = sumOfProducts(y, x) + return [xSquaresSum, xyProductsSum, withResidualSumOfSquares ? sumOfSquaredResiduals(y, x, xSquaresSum, xyProductsSum) : DOUBLE_DOUBLE_ZERO] + }, ) - return {xSumOfSquares, productsSum, ySumOfSquares, xExponent: firstExponent, yExponent: secondExponent} + return {xSumOfSquares, productsSum, residualSumOfSquares, xExponent: firstExponent, yExponent: secondExponent} +} + +/** + * The sum of squared residuals of a simple linear regression, from the deviations of y and x. + * + * The residuals are computed one by one rather than as `sum(dy^2) - productsSum^2 / xSumOfSquares`, + * which cancels when the points lie almost on a line: each residual is then the small difference of + * a deviation and its fitted value, which double-double keeps to about 2^-106 of the deviation. + * + * The rounding errors of the two means shift every residual by the same amount, the error of the y + * mean minus the slope times the error of the x mean, which can exceed the residuals themselves when + * the slope is large. The exact residuals sum to 0, so they are centered on their mean, which + * removes that shift, before they are squared. + * + * @param {DoubleDouble[]} yDeviations - the deviations of y from its mean + * @param {DoubleDouble[]} xDeviations - the deviations of x from its mean, at the same scale as the sums + * @param {DoubleDouble} xSumOfSquares - `sum(dx^2)` + * @param {DoubleDouble} productsSum - `sum(dy * dx)` + * @returns {DoubleDouble} `sum((dy - slope * dx)^2)`, or 0 when all the x values are equal + */ +function sumOfSquaredResiduals(yDeviations: DoubleDouble[], xDeviations: DoubleDouble[], xSumOfSquares: DoubleDouble, productsSum: DoubleDouble): DoubleDouble { + if (xSumOfSquares.hi === 0) { + return DOUBLE_DOUBLE_ZERO + } + const slope = divideDoubleDoubles(productsSum, xSumOfSquares) + const residuals = yDeviations.map((deviation, index) => subtractDoubleDouble(deviation, multiplyDoubleDouble(slope, xDeviations[index]))) + const residualsSum = residuals.reduce((total, residual) => addDoubleDouble(total, residual), DOUBLE_DOUBLE_ZERO) + const residualsMean = divideDoubleDouble(residualsSum, residuals.length) + const centeredResiduals = residuals.map((residual) => subtractDoubleDouble(residual, residualsMean)) + return sumOfProducts(centeredResiduals, centeredResiduals) } /** diff --git a/src/interpreter/plugin/StatisticalAggregationPlugin.ts b/src/interpreter/plugin/StatisticalAggregationPlugin.ts index e7ad6b921c..c363bc1351 100644 --- a/src/interpreter/plugin/StatisticalAggregationPlugin.ts +++ b/src/interpreter/plugin/StatisticalAggregationPlugin.ts @@ -31,9 +31,7 @@ import {covariance, regressionSums, RegressionSums, sumOfSquaredDeviations} from import { divideDoubleDoubles, multiplyByTwoToThe, - multiplyDoubleDouble, roundDoubleDouble, - subtractDoubleDouble, } from '../doubleDouble' import {FunctionArgumentType, FunctionPlugin, FunctionPluginTypecheck, ImplementedFunctions} from './FunctionPlugin' @@ -381,11 +379,7 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (sums instanceof CellError) { return sums } - // at the scale of the sums, the slope and the explained sum of squares Sxy^2 / Sxx <= Syy are finite - const explained = multiplyDoubleDouble(divideDoubleDoubles(sums.productsSum, sums.xSumOfSquares), sums.productsSum) - const residual = roundDoubleDouble(subtractDoubleDouble(sums.ySumOfSquares, explained)) - // the residual is non-negative; a tiny negative rounding remainder counts as 0 - return multiplyByTwoToThe(Math.sqrt(Math.max(0, residual) / (n - 2)), sums.yExponent) + return multiplyByTwoToThe(Math.sqrt(roundDoubleDouble(sums.residualSumOfSquares) / (n - 2)), sums.yExponent) }) } @@ -546,11 +540,11 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func * * @param {number[]} knownYs - a non-empty array of the dependent values * @param {number[]} knownXs - an array of the independent values, of the same length as `knownYs` - * @param {boolean} withYSumOfSquares - whether to compute `ySumOfSquares` + * @param {boolean} withResidualSumOfSquares - whether to compute `residualSumOfSquares` * @returns {RegressionSums | CellError} the scaled sums, or the error */ -function nonDegenerateRegressionSums(knownYs: number[], knownXs: number[], withYSumOfSquares: boolean): RegressionSums | CellError { - const sums = regressionSums(knownYs, knownXs, withYSumOfSquares) +function nonDegenerateRegressionSums(knownYs: number[], knownXs: number[], withResidualSumOfSquares: boolean): RegressionSums | CellError { + const sums = regressionSums(knownYs, knownXs, withResidualSumOfSquares) if (sums.xSumOfSquares.hi === 0) { return new CellError(ErrorType.DIV_BY_ZERO, ErrorMessage.EqualXValues) } From 6001b6e32f66bf2710cc212572795bc55738610e Mon Sep 17 00:00:00 2001 From: Kuba Sekowski Date: Thu, 8 Oct 2026 09:23:43 +0200 Subject: [PATCH 12/13] Describe the STEYX fix for points that lie almost on a line in the changelog (HF-491) Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9f8e8e3917..4bdba941dd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,7 +24,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Fixed the MAXPOOL and MEDIANPOOL functions throwing an uncaught `TypeError` instead of returning the `#VALUE!` error when the range dimensions are not a whole multiple of the window size and the stride. [#1718](https://github.com/handsontable/hyperformula/pull/1718) - Fixed the `VAR`, `STDEV`, `DEVSQ`, `COVARIANCE`, `SLOPE`, `STEYX`, `DVAR` and `DSTDEV` functions, their variants, and the matching `SUBTOTAL` modes losing precision on data with a large mean and a small spread. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed a bug where `SLOPE` and `STEYX` returned `#NUM!` or an arbitrary number instead of `#DIV/0!` when all the x values are equal. [#1784](https://github.com/handsontable/hyperformula/pull/1784) -- Fixed a bug where `STEYX` returned `#NUM!` instead of `0` for points that lie exactly on a line. [#1784](https://github.com/handsontable/hyperformula/pull/1784) +- Fixed a bug where `STEYX` returned `#NUM!` or `0` instead of the standard error for points that lie almost on a line. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed a bug where `VAR`, `STDEV`, `DVAR`, `DSTDEV`, `COVARIANCE`, `SLOPE`, `STEYX`, their variants, and the matching `SUBTOTAL` modes returned an error or `0` for very large or very small values although the result was within the range of numbers. [#1784](https://github.com/handsontable/hyperformula/pull/1784) - Fixed the `MOD` function returning a remainder with the sign of the dividend instead of the sign of the divisor, which made the results differ from Excel and Google Sheets for arguments with opposite signs (e.g. `=MOD(-3, 12)` now returns `9` instead of `-3`). [#1747](https://github.com/handsontable/hyperformula/issues/1747) - Fixed a bug where moving or pasting a formula with an undefined name to another sheet incorrectly added an empty global named expression. [#1728](https://github.com/handsontable/hyperformula/pull/1728) From 91456eb95d4efa958b8bdc55a54d64a38f2434ae Mon Sep 17 00:00:00 2001 From: Kuba Sekowski Date: Thu, 8 Oct 2026 09:27:15 +0200 Subject: [PATCH 13/13] Fix the release script --- script/release/release.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/script/release/release.sh b/script/release/release.sh index cebdc65fd9..50d7f0051d 100644 --- a/script/release/release.sh +++ b/script/release/release.sh @@ -593,9 +593,9 @@ step "2. Create release/$VERSION in hyperformula-tests, then sync the suite" run npm run test:setup-private # 3. Bump version + release date (each half skipped when already correct, so a -CURRENT_VERSION="$(node -e 'process.stdout.write(require("./package.json").version||"")' 2>/dev/null || true)" # re-run after a later failure does not rewrite files it already wrote) step "3. Bump version + HT_RELEASE_DATE" +CURRENT_VERSION="$(node -e 'process.stdout.write(require("./package.json").version||"")' 2>/dev/null || true)" if [[ "$CURRENT_VERSION" == "$VERSION" ]]; then skip "package.json already at $VERSION" elif $DRY_RUN; then