Compute DEVSQ, COVARIANCE, SLOPE, STEYX and the D-variance functions in double-double (HF-491)

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

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<DoubleDouble>((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<DoubleDouble>((sum, deviation, index) => addDoubleDouble(sum, multiplyDoubleDouble(deviation, second[index])), DOUBLE_DOUBLE_ZERO)
+}

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.
  *

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))
     })
   }
 

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<StatisticalAggregationPlugin> {
@@ -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])))
       })
   }