Fix precision loss in variance, standard deviation and related statistics (HF-491)

[marcin-kordas-hoc]

Context

The variance functions compute sum(x^2) - sum(x)^2 / n in one pass. When the mean is large relative to the spread (typical of measurement data), the two terms are nearly equal and cancel, so every significant digit is lost:

Formula	Before	Exact (stored doubles)	After
=STDEV.S({10000000.000,10000000.001,10000000.002})	#NUM!	0.0010000001639127731	0.0010000001639127731
=VAR.S( same )	-0.03125 (negative, no error)	1.0000003278255731e-06	1.0000003278255731e-06
=STDEV.S({100000000000.003,100000000000.004,100000000000.005})	0	0.0009994603901554338	0.0009994603901554338

smartRounding: false and precisionEpsilon: 0 do not change the old results.

DEVSQ, COVARIANCE.P, COVARIANCE.S, SLOPE, STEYX and the database functions DSTDEV, DSTDEVP, DVAR, DVARP compute the mean first and lose digits in the same situation when the values use all 15 significant digits. For example, for three values near 100000000000, STEYX returned 0.0187 instead of 0.00623 and SLOPE returned 1000.37 instead of 1000.53.

Affected:

• VAR.S, VAR.P, STDEV.S, STDEV.P, VARA, VARPA, STDEVA, STDEVPA, their legacy names (VAR, VARP, VARS, STDEV, STDEVP, STDEVS), and SUBTOTAL modes 7, 8, 10 and 11 (and 107, 108, 110, 111), all backed by MomentsAggregate
• DEVSQ, COVARIANCE.P, COVARIANCE.S (and COVAR, COVARIANCEP, COVARIANCES), SLOPE, STEYX, DSTDEV, DSTDEVP, DVAR, DVARP

The change

Double-double helpers (src/interpreter/doubleDouble.ts): TwoSum (in the Fast2Sum form, taking the error from the larger addend, so it cannot overflow while the sum is finite) and Dekker's TwoProduct with Veltkamp's split (no FMA needed), plus addition, subtraction, multiplication and division of double-doubles, and rounding to a double. Near the top of the double range they rescale by exact powers of two: twoProduct scales the larger factor by 2^-53 when a factor exceeds 2^996 or the product exceeds 2^1023, so its error term is exact for every finite product above 2^-969, and the divisions scale dividends above 2^1000 by 2^-64.

MomentsAggregate (src/interpreter/plugin/MomentsAggregate.ts) stays composable, so the per-range cache keeps working. Besides the count, it keeps the sums of deviations S1 and of squared deviations S2 from a shift (the first value added), stored as double-double, and rebases them when two aggregates compose. The deviations are stored multiplied by 2^-exponent, an exact power of two kept in the aggregate, so the shifted sums cannot overflow: the exponent is 0 until the scaled sum of squares exceeds 2^960 (deviations of about 1e144 and above), and two aggregates compose at the larger exponent, raised further when the difference of their shifts needs it. Adding one value, the common case when a range is folded cell by cell, takes a fast path. The sum of squared deviations S2 - S1^2 / n is evaluated in double-double as S2 - S1 * (S1 / n), so the S1^2 term cannot overflow, the variance is rounded once, and it is then multiplied by 2^(2 * exponent). The variance is #NUM! only when it exceeds the largest double, and it does not depend on the order of the arguments: =VAR.S(0, 1e154, 1e154) and =VAR.S(1e154, 1e154, 0) both return 3.33e307 (develop returns #NUM! for both), and =VAR.S(-1.2e154, 0, 1.2e154) returns 1.44e308 although the sum of squared deviations exceeds the largest double. The exponent can also go negative (down to -600) when tiny differences follow equal values, and the standard deviation is taken at the scale of the sums before scaling back, so STDEV.* is finite when only the variance overflows and non-zero for spreads below about 1e-154 (=STDEV.S(1e-170, 2e-170) returns 7.07e-171).

Deviation sums (src/interpreter/deviationSums.ts): DEVSQ, COVARIANCE.*, SLOPE and STEYX compute the mean in double-double, then the deviations and the sums of their products in double-double, and round once. Each array is measured at its own power-of-two scale and the rounded result is scaled back: an array whose largest magnitude is below 2^-400 is scaled up to it (exact), so squared deviations cannot underflow and the sum of squared x deviations is 0 only when all the x values are equal; an array above 2^400 is scaled down to it only when the unscaled sums overflow, so on ordinary data the results are bit-identical to the unscaled computation. DEVSQ is only scaled up, because when its unscaled sum overflows, its result does too. SLOPE and STEYX share regressionSums and the check for equal x values. STEYX does not use Syy - Sxy^2 / Sxx, which cancels when the points lie almost on a line: it computes each residual dy - slope * dx in double-double, centers the residuals on their mean (which removes the shift that the rounding errors of the two means add to every residual), and sums their squares. SLOPE skips the residuals.

D-variance functions: DVAR, DVARP, DSTDEV and DSTDEVP fold their values through MomentsAggregate, so they return the same results as VAR.S, VAR.P, STDEV.S and STDEV.P, including when the sum of squared deviations overflows.

Behaviour changes besides precision: STEYX returns the standard error instead of #NUM! or 0 for points that lie almost on a line (=STEYX({0.25,0.4,0.5499999999999999},{0.1,0.2,0.3}) returns 2.83e-17, the exact value for the stored doubles, which are not evenly spaced, where develop returns #NUM!; =STEYX({100000000,200000000,300000000.00000024},{1,2,3}) returns 9.733e-8 where develop returns 0); SLOPE and STEYX return #DIV/0! with the message "All x values are equal." (the new ErrorMessage.EqualXValues) when all the x values are equal, as Excel and Google Sheets return #DIV/0! (develop returns #NUM! or an arbitrary number), and finite results where the sum of squared deviations of the x values overflows or underflows (=SLOPE({1,2},{1e-170,2e-170}) returns 1e170; develop returns #NUM!). DVAR, DSTDEV and COVARIANCE return finite results where only an intermediate sum overflows (=DSTDEV of -1.2e154, 0, 1.2e154 returns 1.2e154, like STDEV.S).

AVERAGE, AVERAGEA and SUBTOTAL 1/101 no longer use MomentsAggregate: they reduce an AverageResult (a plain sum and a count, moved from ConditionalAggregationPlugin to its own module) under their own cache keys (_AVERAGE and _AVERAGE_A), so they do not compute the double-double sums. Each argument is evaluated once, and the results are the same as before.

Accuracy

The reference is the exact result for the stored doubles, reported as LRE (digits of agreement, McCullough 1998). Exact references: BigInt rational arithmetic, Python's statistics (exact fractions) and mpmath at 60 digits. The data is 24 datasets, including all 9 NIST StRD univariate sets. Measured for VAR/STDEV at 7a88902 and not re-measured since; the variance algorithm is the same apart from the S1 * (S1 / n) evaluation and the power-of-two exponent, which is exact scaling and is 0 for all 24 datasets.

	min LRE over 24 datasets
this PR	15
Gnumeric 1.12.56	15
R 4.3.3 sd()	4.7
NumPy 2.5.3 / SciPy 1.18.1	0
plain two-pass	0
Excel (two-pass; 50 identical values → 7.5e-8, exact 0)	0
Welford / Chan pairwise	2.4
before (one-pass)	−1 (negative variance)

Double-double arithmetic does not guarantee correct rounding, so the JSDoc states the accuracy as about one unit in the last place. Split arguments, =STDEV.S(A1:A50, A51:A100) in either order, hold the same accuracy. When a range reuses the cached aggregate of a smaller range or the values come in several arguments, partial aggregates are composed in a different order than the left-to-right fold of the D-functions, so the results can differ from them in the last bits.

Compared with Excel: near 1e7 the results are equal, but near 1e11 Excel's two-pass result differs from the exact value in the 4th significant digit (0.000999538 vs 0.000999460), and for 50 identical values whose mean is not exactly representable Excel returns 7.5e-8 instead of 0. This PR returns the exact values.

Alternatives measured

• Plain two-pass for the MomentsAggregate functions: mean first, then deviations, which needs every value again. It cannot use the range cache, and a running =STDEV.S(A$1:An) over 5,000 rows became 10x slower to build and 38x slower to recalculate after one edit. The functions that already take whole arrays (DEVSQ, COVARIANCE.*, SLOPE, STEYX, D-functions) use two passes, in double-double.
• Shifted sums plus a term that reproduced Excel's rounding of the mean: composable, but it inherited Excel's error (min LRE 0) and depended on argument order.
• Compensated (Neumaier) shifted sums: min LRE 13.1 at about the same cost as double-double.

How did you test your changes?

• The tests PR (handsontable/hyperformula-tests#70) adds 118 tests to the spec of each affected function; 63 of them fail on develop. Expected values are the exact results for the stored doubles, cross-checked with Python's statistics, compared as a ratio within 1e-14 with smartRounding off. They cover every affected function, the legacy STDEV name, SUBTOTAL 7 and 11, two ranges as separate arguments, range + scalar, a range extending a cached smaller range, two cached ranges kept at different scales (tiny spreads of different sizes, and tiny values with huge ones) combined in either order, recalculation after an edit, the NIST NumAcc4 dataset, identical values (also near 1e300), variances near 1e300, STEYX on points exactly and almost on a line, SLOPE/STEYX with equal x values (including the error message), and overflow: a squared deviation close to the largest double in either order, shifted sums that overflow, two cached ranges whose combination overflows, a running total that overflows, a value equal to the largest double, products of deviations near the largest double that cancel, the sum of squared x deviations overflowing or underflowing in SLOPE/STEYX, the D-functions and COVARIANCE matching VAR/STDEV when the sum of squared deviations overflows, tiny spreads in STDEV.S and DSTDEV, products of large and tiny values that cancel in COVARIANCE and SLOPE, an overflowing slope in STEYX, the same variance for the values in any order and for a range split at any point when the squared deviations from the first value overflow, a finite variance whose sum of squared deviations overflows, #NUM! when the variance exceeds the largest double, AVERAGE and VAR.S and AVERAGEA and VARA over the same range, and AVERAGE, AVERAGEA and SUBTOTAL 1 evaluating each argument once.
• Full suite (after merging develop): 6733 passed, 0 failed, 3 skipped.
• STEYX on points almost on a line was checked against exact rational results: 4,000 datasets (n from 3 to 12, intercepts, slopes and x values over 30 to 40 orders of magnitude, residual sums of squares down to 1e-37 of the sum of squared y deviations), largest relative error 1.8e-15, plus 60k datasets of the general COVARIANCE/SLOPE/STEYX fuzz with no STEYX result off by more than 1e-12.
• The double-double helpers were checked against exact BigInt arithmetic: about 31M cases each for twoSum and twoProduct (boundary grids around 2^996, 2^1023, the largest double, subnormals, signed zeros, infinities and NaN, plus random bit patterns) with no inexact result where an exact one is representable, and about 5M cases each for the two divisions.
• An engine differential test against develop and exact results (about 2.9M formula evaluations over ordinary data, large means with small spreads, extreme magnitudes, and the same data as one range, two ranges, scalars and a cached range extended) found no case where develop is finite and this PR returns an error, apart from the ones listed under "Not in this PR".
• A differential test of the power-of-two exponent against exact BigInt rational results: 24,000 datasets (magnitudes from 1e140 to the largest double mixed with zeros and subnormals, opposite-sign values near the largest double, large means with small spreads at every exponent, identical values, variances near the overflow threshold; n from 2 to 200), each evaluated in about 20 forms (argument orders, split ranges, range + scalar, cached prefix ranges), 10.8M checks in all, for VAR.*, STDEV.*, VARA, VARPA, STDEVA, STDEVPA and SUBTOTAL 7, 8, 10 and 11. Where the exact variance is a normal double, the largest error is 1 ulp (3 results at 2 ulp just above 2^-1022; largest relative error 3.75e-16), all forms agree within 2 ulp, and #NUM! is returned exactly when the exact variance rounds to infinity.
• Benchmark against develop at 92ebfbd (Apple M4 Pro, Node 24; fresh process per run, 7 alternating runs, warm median): one STDEV.S over 100k cells 14.5 → 18.5 ms (1.28x), AVERAGE plus VAR.S over the same 20k-cell range 5.6 → 7.2 ms (1.29x), recalculation of a running STDEV.S(A$1:An) over 5k rows after one edit 12.9 → 13.9 ms (1.08x), building a running STDEV.S(A$1:An) or AVERAGE(A$1:An) over 5k rows 1.00x and 0.99x. Data that needs the exponent (values near 1e150) adds about 6% to the 100k STDEV.S.
• The power-of-two scaling of the D-functions, COVARIANCE, SLOPE and STEYX was checked against exact BigInt results and against the previous head (92ebfbd): about 18k random datasets with 23 formulas each, 18k more from an independent generator, an 80k-dataset MomentsAggregate composition fuzz and 19,664 equal-x cases. No result differs from 92ebfbd on ordinary data, the D-functions equal VAR/STDEV in every regime, and equal x values always give #DIV/0!. On 100k ordinary values, every affected function is within about ±3% of 92ebfbd, and DVAR/DSTDEV are about 5% faster.
• Heap after building a running formula over 5k rows: STDEV.S 10.49 → 11.82 MB, AVERAGE 10.48 → 10.96 MB. A cached MomentsAggregate is larger than on develop (two double-doubles, a shift and an exponent); a cached AverageResult is smaller than develop's aggregate.

Not in this PR

• VAR.* results below 2^-1022 (spreads below about 1e-154) are subnormal numbers with few significant bits; they are rounded at the scale of the sums and then scaled back, so they can be 1 subnormal ulp off.
• CORREL, PEARSON, RSQ, Z.TEST, T.TEST, F.TEST, SKEW and SKEW.P still compute their sums of deviations in plain doubles (jStat), so they keep the precision loss on data with a large mean.
• When an intermediate sum overflows, COVARIANCE and SLOPE scale large arrays down, which rounds away values about 2^1400 times smaller than the largest one. On extremely ill-conditioned data (values near the largest double that cancel, or data spanning about 600 orders of magnitude) they can then return a finite but inaccurate result where develop returned #NUM!; no double-double method can resolve these.
• STEYX rounds twice (division, then square root), so it can be about 1 ulp off.
• Double-double edge cases that predate this change: divideDoubleDoubles with a dividend below about 2^-950 can be a few ulps less accurate than plain division, and sums or quotients within half an ulp of the largest double return infinity instead of the largest double.
• The precision of AVERAGE: =AVERAGE(1e16, 1, -1e16) returns 0, as on develop and in Excel.
• INTERCEPT is not implemented in HyperFormula.

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)
• 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. Tests: handsontable/hyperformula-tests#70
2. Merged into this branch: Compute DEVSQ, COVARIANCE, SLOPE, STEYX and the D-variance functions in double-double (HF-491) #1798 (DEVSQ, COVARIANCE, SLOPE, STEYX and the D-variance functions)

Checklist:

• I have reviewed the guidelines about Contributing to HyperFormula and I confirm that my code follows the code style of this project.
• I have signed the Contributor License Agreement.
• My change is compliant with the OpenDocument standard.
• My change is compatible with Microsoft Excel.
• My change is compatible with Google Sheets.
• I described my changes in the CHANGELOG.md file.
• My changes require a documentation update.
• My changes require a migration guide.

🤖 Generated with Claude Code

---

Note

Medium Risk
Touches core numeric paths for many spreadsheet functions; behavior changes (errors, edge magnitudes) are intentional but affect formula results users may rely on.

Overview
Fixes catastrophic cancellation in variance, standard deviation, and related statistics when the mean is large relative to the spread (e.g. STDEV.S on values near 1e7 with tiny differences).

Introduces double-double arithmetic (doubleDouble.ts) and two aggregation layers: MomentsAggregate for composable VAR/STDEV/SUBTOTAL/database DVAR/DSTDEV (shifted sums with power-of-two scaling), and deviationSums.ts for DEVSQ, COVARIANCE, SLOPE, and STEYX. STEYX now sums squared residuals per point (with centered residuals) instead of a formula that cancels on nearly collinear data.

SLOPE and STEYX return #DIV/0! with "All x values are equal." when all x values match. AVERAGE / AVERAGEA use a lightweight AverageResult (sum + count) instead of the variance aggregate. Release script step 3 reads CURRENT_VERSION after the step label.

^{Reviewed by Cursor Bugbot for commit 91456eb. Bugbot is set up for automated code reviews on this repo. Configure here.}

[cloudflare-workers-and-pages]

Deploying with    Cloudflare Workers

The latest updates on your project. Learn more about integrating Git with Workers.

Status	Name	Latest Commit	Preview URL	Updated (UTC)
✅ Deployment successful! View logs	hyperformula-docs	91456eb	Commit Preview URL Branch Preview URL	Oct 08 2026, 07:32 AM

[github-actions]

Performance comparison of head (91456eb) vs base (3a9c34b)

                                     testName |    base |    head | change
--------------------------------------------------------------------------
                                      Sheet A |   507.1 |  499.66 | -1.47%
                                      Sheet B |  160.06 |  159.16 | -0.56%
                                      Sheet T |  142.05 |  141.34 | -0.50%
                                Column ranges |  472.74 |  475.08 | +0.49%
                                Sorted lookup | 13757.7 | 14190.3 | +3.14%
Sheet A:  change value, add/remove row/column |   15.75 |   16.02 | +1.71%
 Sheet B: change value, add/remove row/column |  140.76 |  134.88 | -4.18%
                   Column ranges - add column |  149.09 |  150.49 | +0.94%
                Column ranges - without batch |  462.82 |  467.26 | +0.96%
                        Column ranges - batch |  115.37 |  117.08 | +1.48%

[cursor]

Cursor Bugbot has reviewed your changes using default effort and found 2 potential issues.

^{❌ Bugbot Autofix is OFF. To automatically fix reported issues with cloud agents, enable autofix in the Cursor dashboard.}

^{Reviewed by Cursor Bugbot for commit bb43a62. Configure here.}

[codecov]

Codecov Report

❌ Patch coverage is 99.26471% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 97.28%. Comparing base (3a9c34b) to head (91456eb).

Files with missing lines	Patch %	Lines
src/interpreter/plugin/MomentsAggregate.ts	97.72%	2 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1784      +/-   ##
===========================================
+ Coverage    97.26%   97.28%   +0.02%     
===========================================
  Files          203      207       +4     
  Lines        16022    16218     +196     
  Branches      3555     3588      +33     
===========================================
+ Hits         15583    15777     +194     
- Misses         431      433       +2     
  Partials         8        8              

Files with missing lines	Coverage Δ
src/error-message.ts	100.00% <100.00%> (ø)
src/interpreter/deviationSums.ts	100.00% <100.00%> (ø)
src/interpreter/doubleDouble.ts	100.00% <100.00%> (ø)
src/interpreter/plugin/AverageResult.ts	100.00% <100.00%> (ø)
...interpreter/plugin/ConditionalAggregationPlugin.ts	100.00% <100.00%> (ø)
src/interpreter/plugin/DatabasePlugin.ts	94.97% <100.00%> (-0.29%)	⬇️
src/interpreter/plugin/NumericAggregationPlugin.ts	97.23% <100.00%> (-0.23%)	⬇️
...interpreter/plugin/StatisticalAggregationPlugin.ts	96.41% <100.00%> (+0.18%)	⬆️
src/interpreter/plugin/MomentsAggregate.ts	97.72% <97.72%> (ø)

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.
• 📦 JS Bundle Analysis: Save yourself from yourself by tracking and limiting bundle sizes in JS merges.