diff --git a/CHANGELOG.md b/CHANGELOG.md index 97193c678..22af2493e 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 000000000..37e4a227d --- /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 4a7117429..95fd34ff2 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 1208a0f17..d5c2150ae 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 75b8c820b..73b79d776 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]))) }) }