diff --git a/CHANGELOG.md b/CHANGELOG.md index 787c642ad9..4bdba941dd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,10 @@ 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 `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!` 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) 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 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/deviationSums.ts b/src/interpreter/deviationSums.ts new file mode 100644 index 0000000000..f5c16c10ed --- /dev/null +++ b/src/interpreter/deviationSums.ts @@ -0,0 +1,306 @@ +/** + * @license + * Copyright (c) 2025 Handsoncode. All rights reserved. + */ + +import { + addDoubleDouble, + divideDoubleDouble, + divideDoubleDoubles, + DOUBLE_DOUBLE_ZERO, + DoubleDouble, + multiplyByPowerOfTwo, + multiplyByTwoToThe, + multiplyDoubleDouble, + roundDoubleDouble, + subtractDoubleDouble, +} from './doubleDouble' + +/** + * Sums of squared deviations and of products of deviations from the mean, for the functions built + * 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 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) - 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 */ + readonly yExponent: number, +} + +/** + * The sum of the squared deviations of the values from their mean. + * + * 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 {number} `sum((x - mean)^2)`, rounded once to a double + */ +export function sumOfSquaredDeviations(values: number[]): number { + const {deviations, exponent} = scaledDeviationsFromMean(values, false) + return multiplyByTwoToThe(roundDoubleDouble(sumOfProducts(deviations, deviations)), 2 * exponent) +} + +/** + * The covariance of two arrays: the sum of products of deviations divided by `n - deltaDegreesOfFreedom`. + * + * 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 + */ +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 sums of squared deviations and of products of deviations of a simple linear regression, each + * 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 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} 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[], 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, 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) +} + +/** + * 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 + * @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 + */ +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 largest absolute value of the values. + * + * 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 + */ +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 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 = 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) +} + +/** + * 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 new file mode 100644 index 0000000000..c2497a1fb9 --- /dev/null +++ b/src/interpreter/doubleDouble.ts @@ -0,0 +1,235 @@ +/** + * @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 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 + 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). + * + * 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 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 && 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 ((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)} +} + +/** + * 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. + * + * @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 { + 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) +} + +/** + * 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} + } + 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} +} + +/** + * 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. + * + * @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/AverageResult.ts b/src/interpreter/plugin/AverageResult.ts new file mode 100644 index 0000000000..fefceae88d --- /dev/null +++ b/src/interpreter/plugin/AverageResult.ts @@ -0,0 +1,56 @@ +/** + * @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) + + /** + * @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 + } 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/DatabasePlugin.ts b/src/interpreter/plugin/DatabasePlugin.ts index 1208a0f17f..b91c802b13 100644 --- a/src/interpreter/plugin/DatabasePlugin.ts +++ b/src/interpreter/plugin/DatabasePlugin.ts @@ -11,6 +11,7 @@ import {EmptyValue, getRawValue, InternalScalarValue, InterpreterValue, isExtend 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. @@ -337,13 +338,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return values } - if (values.length <= 1) { - 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 MomentsAggregate.of(values).stdevSValue() ?? new CellError(ErrorType.DIV_BY_ZERO) }) } @@ -362,13 +357,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return values } - if (values.length === 0) { - 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 MomentsAggregate.of(values).stdevPValue() ?? new CellError(ErrorType.DIV_BY_ZERO) }) } @@ -386,12 +375,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return values } - if (values.length <= 1) { - 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 MomentsAggregate.of(values).varSValue() ?? new CellError(ErrorType.DIV_BY_ZERO) }) } @@ -410,12 +394,7 @@ export class DatabasePlugin extends FunctionPlugin implements FunctionPluginType return values } - if (values.length === 0) { - 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 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..b962c3b485 --- /dev/null +++ b/src/interpreter/plugin/MomentsAggregate.ts @@ -0,0 +1,354 @@ +/** + * @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, 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 + */ + 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 5a675e9360..6f426c6080 100644 --- a/src/interpreter/plugin/NumericAggregationPlugin.ts +++ b/src/interpreter/plugin/NumericAggregationPlugin.ts @@ -14,6 +14,8 @@ import {coerceBooleanToNumber} from '../ArithmeticHelper' 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' @@ -31,50 +33,6 @@ function zeroForInfinite(value: InternalScalarValue) { } } -class MomentsAggregate { - - public static empty = new MomentsAggregate(0, 0, 0) - - constructor( - public readonly sumsq: number, - public readonly sum: number, - public readonly count: number, - ) { - } - - public static single(arg: number): MomentsAggregate { - return new MomentsAggregate(arg * arg, arg, 1) - } - - public compose(other: MomentsAggregate) { - return new MomentsAggregate(this.sumsq + other.sumsq, this.sum + other.sum, this.count + other.count) - } - - public averageValue(): Maybe { - if (this.count > 0) { - return this.sum / this.count - } else { - return undefined - } - } - - public varSValue(): Maybe { - if (this.count > 1) { - return (this.sumsq - (this.sum * this.sum) / this.count) / (this.count - 1) - } else { - return undefined - } - } - - public varPValue(): Maybe { - if (this.count > 0) { - return (this.sumsq - (this.sum * this.sum) / this.count) / this.count - } else { - return undefined - } - } -} - export class NumericAggregationPlugin extends FunctionPlugin implements FunctionPluginTypecheck { public static implementedFunctions: ImplementedFunctions = { 'SUM': { @@ -298,17 +256,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.doAverageA(ast.args, state) } public vars(ast: ProcedureAst, state: InterpreterState): InternalScalarValue { @@ -353,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) } } @@ -364,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) } } @@ -439,13 +385,33 @@ 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) + } + + 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. + * + * @param {Ast[]} args - the function arguments + * @param {InterpreterState} state - interpreter state + * @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 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 - } else { - return result.averageValue() ?? new CellError(ErrorType.DIV_BY_ZERO) } + return result.averageValue() ?? new CellError(ErrorType.DIV_BY_ZERO) } private doVarS(args: Ast[], state: InterpreterState): InternalScalarValue { @@ -474,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) } } @@ -485,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 75b8c820b7..c363bc1351 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,12 @@ import { sumsqerr, variance } from './3rdparty/jstat/jstat' +import {covariance, regressionSums, RegressionSums, sumOfSquaredDeviations} from '../deviationSums' +import { + divideDoubleDoubles, + multiplyByTwoToThe, + roundDoubleDouble, +} from '../doubleDouble' import {FunctionArgumentType, FunctionPlugin, FunctionPluginTypecheck, ImplementedFunctions} from './FunctionPlugin' export class StatisticalAggregationPlugin extends FunctionPlugin implements FunctionPluginTypecheck { @@ -185,7 +190,7 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (coerced.length === 0) { return 0 } - return sumsqerr(coerced) + return sumOfSquaredDeviations(coerced) }) } @@ -282,7 +287,7 @@ export class StatisticalAggregationPlugin extends FunctionPlugin implements Func if (n === 1) { return 0 } - return covariance(ret[0], ret[1]) * (n - 1) / n + return covariance(ret[0], ret[1], 0) }) } @@ -301,7 +306,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 covariance(ret[0], ret[1], 1) }) } @@ -370,7 +375,11 @@ 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 sums = nonDegenerateRegressionSums(ret[0], ret[1], true) + if (sums instanceof CellError) { + return sums + } + return multiplyByTwoToThe(Math.sqrt(roundDoubleDouble(sums.residualSumOfSquares) / (n - 2)), sums.yExponent) }) } @@ -389,7 +398,12 @@ 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]) + 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) }) } @@ -520,6 +534,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} withResidualSumOfSquares - whether to compute `residualSumOfSquares` + * @returns {RegressionSums | CellError} the scaled sums, or the error + */ +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) + } + return sums +} + function parseTwoArrays(dataX: SimpleRangeValue, dataY: SimpleRangeValue): CellError | [number[], number[]] { const xit = dataX.iterateValuesFromTopLeftCorner() const yit = dataY.iterateValuesFromTopLeftCorner()