From 54904ea2136f5f0580448978da83e1461d0a3856 Mon Sep 17 00:00:00 2001 From: marcin-kordas-hoc Date: Wed, 7 Oct 2026 08:22:24 +0000 Subject: [PATCH] Compute the deviation sums of DEVSQ, COVARIANCE, SLOPE, STEYX and the D-variance functions in double-double These functions subtract the mean from each value and sum the products. New src/interpreter/deviationSums.ts keeps the mean, the deviations and their sums in double-double (src/interpreter/doubleDouble.ts) and rounds once at the end, so the results stay accurate when the mean is large relative to the spread. It is used by DEVSQ, COVARIANCE.P, COVARIANCE.S, SLOPE and STEYX (StatisticalAggregationPlugin) and by DSTDEV, DSTDEVP, DVAR and DVARP (DatabasePlugin). STEYX evaluates its residual in double-double too. A division by a non-finite divisor returns the plain double quotient, so results that are finite in 3.3.0 stay finite for values up to about 1e300. SLOPE and STEYX return #NUM! when all the x values are equal, and STEYX returns 0 for points that lie exactly on a line. 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]))) }) }