From e2f3d7621d5b6a4a68afec7254920f9957ed358a Mon Sep 17 00:00:00 2001 From: jmh530 Date: Sat, 17 Oct 2026 20:11:01 -0400 Subject: [PATCH 1/7] Preserve summation state when merging mean accumulators Merge matching-element-type summators directly instead of converting the incoming accumulator to a scalar sum. The conversion discarded saved corrections and precise overflow information before subsequent observations could cancel the large terms. Use the summator's existing merge operation without changing public fields, scalar insertion, or batch accumulation. Retain the existing behavior for other element types. Add regressions for compensated low-order contributions, precise overflow followed by cancellation, source preservation, and empty and self merges. AI-Assisted-By: Codex (GPT-6 Astra Medium) --- .../mir/stat/descriptive/univariate/means.d | 61 ++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/source/mir/stat/descriptive/univariate/means.d b/source/mir/stat/descriptive/univariate/means.d index 9b59949..bca422c 100644 --- a/source/mir/stat/descriptive/univariate/means.d +++ b/source/mir/stat/descriptive/univariate/means.d @@ -103,7 +103,11 @@ struct MeanAccumulator(T, Summation summation) void put(F = T)(MeanAccumulator!(F, summation) m) { count += m.count; - summator.put(cast(T) m.summator); + static if (is(F == T)) + // Merge saved partials and corrections before rounding to a scalar. + summator += m.summator; + else + summator.put(cast(T) m.summator); } } @@ -150,6 +154,61 @@ unittest assert(m0.mean == 29.25 / 12); } +// Merging compensated sums must retain contributions below the rounded sum. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + import std.meta : AliasSeq; + + static foreach (mode; AliasSeq!(Summation.kbn, Summation.kb2)) + {{ + MeanAccumulator!(double, mode) source, merged; + source.put(0x1p54); + source.put(1.0); + merged.put(source); + merged.put(-0x1p54); + assert(merged.count == 3); + assert(merged.sum == 1); + assert(merged.mean == 1.0 / 3); + + // The merge must not consume or alter the source's corrections. + source.put(-0x1p54); + assert(source.sum == 1); + assert(source.count == 3); + }} +} + +// Precise summation also retains overflow internally until later cancellation. +// Copying this summator by value currently requires @system in mir-algorithm. +version(mir_stat_test) +pure nothrow @nogc +unittest +{ + alias Accumulator = MeanAccumulator!(double, Summation.precise); + Accumulator source, merged; + source.put(double.max); + source.put(double.max); + source.put(1.0); + merged.put(source); + merged.put(-double.max); + merged.put(-double.max); + assert(merged.count == 5); + assert(merged.sum == 1); + assert(merged.mean == 0.2); + + source.put(-double.max); + source.put(-double.max); + assert(source.sum == 1); + + // Empty and self merges preserve the accumulated state and count. + merged.put(Accumulator.init); + assert(merged.count == 5 && merged.sum == 1); + merged.put(merged); + assert(merged.count == 10 && merged.sum == 2); + assert(merged.mean == 0.2); +} + /++ Computes the mean of the input. From 66744cba17c31c9ffe00532fa59079155029099f Mon Sep 17 00:00:00 2001 From: jmh530 Date: Sat, 17 Oct 2026 21:04:20 -0400 Subject: [PATCH 2/7] Recover overflowing means for replayable floating-point inputs Retry nonfinite means when the input supports an independent saved traversal. Finite observations such as two maximum values can have a representable mean even when their raw sum overflows. Accumulate large observations with power-of-two scaling and retain tiny observations in a separate unscaled precise sum. Combine the subtotals before rounding when the unscaled subtotal is finite, preserving small remainders from cancellation between the groups. Otherwise divide before restoring the scale. Keep the selected summation algorithm for ordinary finite results and preserve its handling of nonfinite converted observations. Save forward ranges before consumption and support flattened slice traversal. Document that single-pass inputs, streaming accumulators, and finite cancellation errors are outside this recovery path. Add regressions for extreme finite values, subnormal cancellation, arrays and strided slices, shared-cursor forward ranges, single-pass inputs, and nonfinite observations across float, double, and real. AI-Assisted-By: Codex (GPT-6 Astra Medium) --- .../mir/stat/descriptive/univariate/means.d | 187 +++++++++++++++++- 1 file changed, 184 insertions(+), 3 deletions(-) diff --git a/source/mir/stat/descriptive/univariate/means.d b/source/mir/stat/descriptive/univariate/means.d index bca422c..5214873 100644 --- a/source/mir/stat/descriptive/univariate/means.d +++ b/source/mir/stat/descriptive/univariate/means.d @@ -209,6 +209,61 @@ unittest assert(merged.mean == 0.2); } +// Retry a nonfinite floating-point mean without overflowing the raw sum. +// The caller supplies an independently saved traversal of the same observations. +private T recoverOverflowingMean(T, Range)(scope Range values, size_t count, T original) +{ + import std.math : abs, isFinite; + + if (!count) + return original; + + // A power of two at least twice the count leaves headroom for the large + // terms' sum. Multiplication by its reciprocal is exact unless it enters + // the subnormal range; keep those small terms in a separate unscaled sum. + T scale = 2; + for (size_t n = count - 1; n; n >>= 1) + scale *= 2; + const T inverse = 1 / scale; + const T threshold = T.min_normal * scale; + Summator!(T, Summation.precise) high, low; + T lower = T.infinity, upper = -T.infinity; + foreach (value; values) + { + const T x = cast(T) value; + // Genuine NaN/Inf inputs (including overflow on conversion to T) are + // not finite-input sum overflow. Preserve the original result. + if (!isFinite(x)) + return original; + if (x < lower) lower = x; + if (x > upper) upper = x; + if (abs(x) < threshold) + low.put(x); + else + high.put(x * inverse); + } + + const T highSum = high.sum; + const T unscaled = highSum * scale; + T result; + if (isFinite(unscaled)) + { + // Cancellation can cross the two sums. Add the large subtotal to the + // still-precise small sum before rounding, so a tiny remainder is not + // lost by rounding the two opposing subtotals independently. + low.put(unscaled); + result = low.sum / cast(T) count; + } + else + { + // The unscaled subtotal itself overflows: divide before restoring its + // scale. At this magnitude the small sum cannot cancel the large one. + result = highSum * (scale / cast(T) count) + low.sum / cast(T) count; + } + // Rounding on rescaling must not leave the finite input interval. + return result < lower ? lower : result > upper ? upper : result; +} + /++ Computes the mean of the input. @@ -216,6 +271,18 @@ By default, if `F` is not floating point type or complex type, then the result will have a `double` type if `F` is implicitly convertible to a floating point type or a type for which `isComplex!F` is true. +For built-in floating-point results, a nonfinite result is retried with +scaled precise summation when the input can be traversed again. This supports +arrays, slices whose flattened view is a forward range, and other forward +ranges. For example, the mean of two maximum finite values is still that +maximum value, even though their sum overflows. A forward range's saved +traversal must yield the same observations as the first traversal. + +The retry preserves the original result if an observation is nonfinite after +conversion to the result type. Finite results use the selected summation +algorithm without a retry; this is not a general correction for cancellation +error. Single-pass inputs and $(LREF MeanAccumulator) do not retry. + Params: F = controls type of output summation = algorithm for calculating sums (default: Summation.appropriate) @@ -228,7 +295,7 @@ See_also: template mean(F, Summation summation = Summation.appropriate) { import core.lifetime: move; - import std.traits: isIterable; + import std.traits: isIterable, isFloatingPoint; /++ Params: @@ -239,8 +306,28 @@ template mean(F, Summation summation = Summation.appropriate) { alias G = typeof(return); MeanAccumulator!(G, ResolveSummationType!(summation, Range, G)) mean; + static if (isFloatingPoint!G) + { + import std.range.primitives : isForwardRange, save; + static if (isSlice!Range) + { + import mir.ndslice.topology : flattened; + auto flat = r.flattened; + static if (isForwardRange!(typeof(flat))) + auto saved = flat.save; + } + else static if (isForwardRange!Range) + auto saved = r.save; + } mean.put(r.move); - return mean.mean; + const result = mean.mean; + static if (__traits(compiles, saved)) + { + import std.math : isFinite; + if (!isFinite(result)) + return recoverOverflowingMean!G(saved.move, mean.count, result); + } + return result; } /++ @@ -252,7 +339,14 @@ template mean(F, Summation summation = Summation.appropriate) alias G = typeof(return); MeanAccumulator!(G, ResolveSummationType!(summation, const(G)[], G)) mean; mean.put(ar); - return mean.mean; + const result = mean.mean; + static if (isFloatingPoint!G) + { + import std.math : isFinite; + if (!isFinite(result)) + return recoverOverflowingMean!G(ar, mean.count, result); + } + return result; } } @@ -296,6 +390,93 @@ template mean(string summation) mixin("alias mean = .mean!(Summation." ~ summation ~ ");"); } +// Finite observations can have a representable mean but an overflowing sum. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + import std.meta : AliasSeq; + import std.math : isNaN; + import mir.ndslice.slice : sliced; + import mir.ndslice.dynamic : transposed; + + static foreach (T; AliasSeq!(float, double, real)) + static foreach (mode; AliasSeq!(Summation.naive, Summation.pairwise, + Summation.kbn, Summation.kb2, Summation.precise)) + {{ + T[4] values = [T.max, T.max, T.max, T.max]; + assert(mean!(T, mode)(values[]) == T.max); + assert(mean!(T, mode)(T.max, T.max) == T.max); + assert(mean!(T, mode)(values[].sliced) == T.max); + assert(mean!(T, mode)(values[].sliced(2, 2).transposed) == T.max); + assert(mean!(T, mode)(values[].sliced(2, 2).transposed[0]) == T.max); + values[] = -T.max; + assert(mean!(T, mode)(values[]) == -T.max); + + T[] empty; + assert(isNaN(mean!(T, mode)(empty))); + assert(isNaN(mean!(T, mode)(T.nan, T(1)))); + // Preserve the selected summator's handling of infinities. + MeanAccumulator!(T, mode) original; + original.put(T.infinity); + original.put(T(1)); + const actual = mean!(T, mode)(T.infinity, T(1)); + assert(isNaN(original.mean) ? isNaN(actual) : original.mean == actual); + }} + + static foreach (T; AliasSeq!(float, double, real)) + {{ + // Dividing every input first would erase this subnormal remainder. + const T tiny = T.min_normal * T.epsilon; + T[9] values = [T.max, T.max, -T.max, -T.max, + tiny, tiny, tiny, tiny, tiny]; + assert(mean!(T, Summation.naive)(values[]) == tiny); + + // Cancellation between the scaled and unscaled groups must happen + // before their subtotals are rounded. The exact mean is 5/9 * tiny, + // which rounds to tiny rather than zero. + const T boundary = T.min_normal * 32; + T[9] cross = [T.max, T.max, -T.max, -T.max, + boundary, -boundary / 2, -boundary / 2, 5 * tiny, 0]; + assert(mean!(T, Summation.naive)(cross[]) == tiny); + cross[] = -cross[]; + assert(mean!(T, Summation.naive)(cross[]) == -tiny); + }} + // Conversion overflow is not repaired as if it were sum overflow. + assert(mean!float(double.max, double.max) == float.infinity); +} + +// Save forward ranges before the first pass advances a shared cursor. +version(mir_stat_test) +unittest +{ + static struct Forward + { + const(double)[] data; + size_t* position; + size_t ownPosition; + bool empty() const { return (position ? *position : ownPosition) == data.length; } + double front() const { return data[position ? *position : ownPosition]; } + void popFront() { if (position) ++*position; else ++ownPosition; } + auto save() { return Forward(data, null, position ? *position : ownPosition); } + } + static struct Single + { + const(double)[] data; + size_t* position; + bool empty() const { return *position == data.length; } + double front() const { return data[*position]; } + void popFront() { ++*position; } + } + double[2] values = [double.max, double.max]; + size_t position; + assert(mean(Forward(values[], &position)) == double.max); + assert(position == 2); + position = 0; + assert(mean(Single(values[], &position)) == double.infinity); + assert(position == 2); +} + /// version(mir_stat_test) @safe pure nothrow From 29d82c4d532df42f09a62079df13f75f1e306a52 Mon Sep 17 00:00:00 2001 From: jmh530 Date: Sat, 17 Oct 2026 21:34:47 -0400 Subject: [PATCH 3/7] Recover standard deviations when variance exceeds floating-point range Computing sqrt(variance) can return zero or infinity even when the standard deviation is representable. Retry nonfinite, negative, zero, and subnormal variances for replayable floating-point inputs. Normalize observations by a power of two and use compensated two-pass accumulation. Subtract a fixed normalized observation before centering to avoid introducing spread into identical large values. Preserve the zero-centered definition of assumeZeroMean and use a cheaper equality check for zero-variance constant inputs. Keep the selected algorithm for positive normal variances, preserve existing contracts and nonfinite-input behavior, and leave single-pass inputs unchanged. Use direct loops to retain allocation-free operation. This does not repair inaccurate positive normal variances. Add regressions for extreme and subnormal magnitudes, sample and population results, constant inputs, varargs, strided slices, and saved forward-range cursors. Document recovery behavior and update the naive algorithm example to account for recoverable cancellation. AI-Assisted-By: Codex (GPT-6 Astra Medium) --- .../stat/descriptive/univariate/dispersions.d | 206 +++++++++++++++++- .../mir/stat/descriptive/univariate/means.d | 2 + 2 files changed, 201 insertions(+), 7 deletions(-) diff --git a/source/mir/stat/descriptive/univariate/dispersions.d b/source/mir/stat/descriptive/univariate/dispersions.d index ee3a2a3..ff437bd 100644 --- a/source/mir/stat/descriptive/univariate/dispersions.d +++ b/source/mir/stat/descriptive/univariate/dispersions.d @@ -1430,12 +1430,96 @@ unittest assert(x.sliced.variance!float.approxEqual(54.76562 / 11)); } +/* +Recover a representable standard deviation when squaring the observations or +residuals overflowed or underflowed. Scaling by a power of two keeps nearby +observations distinct; subtracting a fixed normalized observation avoids a +rounded mean creating a spread in constant inputs. The original variance is +computed first, so its contracts and ordinary results remain unchanged. +*/ +private T recoverStandardDeviation(T, VarianceAlgo algo, Range)( + Range values, bool isPopulation, T originalVariance) +{ + import std.math : abs, frexp, ldexp, isFinite, sqrt; + import std.range.primitives : empty, front, save; + + if (values.empty) + return sqrt(originalVariance); + const first = cast(T) values.front; + // Zero can mean either identical observations or underflowed squares. + // Equality is cheaper to check than computing the scale for every value. + if (originalVariance == 0 && isFinite(first)) + { + static if (algo == VarianceAlgo.assumeZeroMean) + const canBeConstant = first == 0; + else + enum canBeConstant = true; + if (canBeConstant) + { + bool equal = true; + foreach (x; values.save) + if (cast(T) x != first) + { + equal = false; + break; + } + if (equal) + return 0; + } + } + + T maximum = 0; + foreach (x; values.save) + { + const value = cast(T) x; + if (!isFinite(value)) + return sqrt(originalVariance); + if (abs(value) > maximum) + maximum = abs(value); + } + if (maximum == 0) + return sqrt(originalVariance); + int exponent; + frexp(maximum, exponent); + const scale = ldexp(T(0.5), exponent); + static if (algo == VarianceAlgo.assumeZeroMean) + enum anchor = T(0); + else + const anchor = first / scale; + size_t count; + Summator!(T, Summation.kbn) sum; + foreach (x; values.save) + { + sum.put(cast(T) x / scale - anchor); + ++count; + } + static if (algo == VarianceAlgo.assumeZeroMean) + enum center = T(0); + else + const center = sum.sum / count; + Summator!(T, Summation.kbn) squares; + foreach (x; values) + { + const residual = (cast(T) x / scale - anchor) - center; + squares.put(residual * residual); + } + const v = squares.sum / (count + isPopulation - 1); + return sqrt(v) * scale; +} + /++ Calculates the standard deviation of the input By default, if `F` is not floating point type, then the result will have a `double` type if `F` is implicitly convertible to a floating point type. +For replayable floating-point inputs, a nonfinite, negative, zero, or +subnormal variance triggers a scaled recalculation. This can preserve a +standard deviation even when its square is not representable. Recovery uses +compensated two-pass variance, except that `assumeZeroMean` retains its +zero-centered definition. Ordinary positive normal variances use the selected +algorithm unchanged. Single-pass inputs do not receive this recovery. + Params: F = controls type of output varianceAlgo = algorithm for calculating variance (default: VarianceAlgo.hybrid) @@ -1460,7 +1544,25 @@ template standardDeviation( { import core.lifetime: move; alias G = typeof(return); - return r.move.variance!(G, varianceAlgo, ResolveSummationType!(summation, Range, G))(isPopulation).sqrt; + import std.range.primitives : isForwardRange, ElementType, save; + import std.traits : isFloatingPoint; + static if (isSlice!Range) + { + import mir.ndslice.topology : flattened; + auto flat = r.flattened; + static if (isForwardRange!(typeof(flat)) && isFloatingPoint!(ElementType!(typeof(flat)))) + auto saved = flat.save; + } + else static if (isForwardRange!Range && isFloatingPoint!(ElementType!Range)) + auto saved = r.save; + const v = r.move.variance!(G, varianceAlgo, ResolveSummationType!(summation, Range, G))(isPopulation); + static if (__traits(compiles, saved)) + { + import std.math : isFinite; + if (!isFinite(v) || v < G.min_normal) + return recoverStandardDeviation!(G, varianceAlgo)(saved.move, isPopulation, v); + } + return v.sqrt; } /++ @@ -1470,10 +1572,102 @@ template standardDeviation( @fmamath stdevType!F standardDeviation(scope const F[] ar...) { alias G = typeof(return); - return ar.variance!(G, varianceAlgo, ResolveSummationType!(summation, const(G)[], G)).sqrt; + const v = ar.variance!(G, varianceAlgo, ResolveSummationType!(summation, const(G)[], G)); + static if (isFloatingPoint!F) + { + import std.math : isFinite; + if (!isFinite(v) || v < G.min_normal) + return recoverStandardDeviation!(G, varianceAlgo)(ar, false, v); + } + return v.sqrt; } } +// Save the forward cursor independently before computing variance. +version(mir_stat_test) +unittest +{ + static struct Forward + { + const(double)[] data; + size_t* position; + size_t ownPosition; + bool empty() const { return (position ? *position : ownPosition) == data.length; } + double front() const { return data[position ? *position : ownPosition]; } + void popFront() { if (position) ++*position; else ++ownPosition; } + auto save() { return Forward(data, null, position ? *position : ownPosition); } + } + static struct Single + { + const(double)[] data; + size_t* position; + bool empty() const { return *position == data.length; } + double front() const { return data[*position]; } + void popFront() { ++*position; } + } + size_t position; + // Standard-deviation recovery also needs an independent saved cursor. + double[3] spread = [-1e200, 0, 1e200]; + position = 0; + assert(standardDeviation!(VarianceAlgo.online)(Forward(spread[], &position)) == 1e200); + assert(position == 3); + position = 0; + assert(standardDeviation!(VarianceAlgo.online)(Single(spread[], &position)) == double.infinity); + assert(position == 3); +} + +// Recover standard deviations whose squared values are outside the type's range. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + import std.meta : AliasSeq; + import std.math : sqrt, abs, isNaN; + import mir.ndslice.slice : sliced; + import mir.ndslice.dynamic : transposed; + + static foreach (T; AliasSeq!(float, double, real)) + {{ + foreach (magnitude; [T.max / 4, T.min_normal, T.min_normal * T.epsilon]) + { + T[3] values = [-magnitude, 0, magnitude]; + static foreach (algo; AliasSeq!(VarianceAlgo.naive, VarianceAlgo.online, + VarianceAlgo.twoPass, VarianceAlgo.hybrid, VarianceAlgo.assumeZeroMean)) + {{ + assert(values[].standardDeviation!(T, algo) == magnitude); + const population = values[].standardDeviation!(T, algo)(true); + const expected = magnitude * sqrt(T(2) / 3); + assert(abs(population - expected) <= 4 * T.epsilon * magnitude + + T.min_normal * T.epsilon); + }} + } + T[3] extremes = [-T.max, 0, T.max]; + assert(extremes[].standardDeviation!(T, VarianceAlgo.naive) == T.max); + assert(standardDeviation!T(-T.max, T(0), T.max) == T.max); + T[3] largeConstant; + largeConstant[] = T.max / 4; + const rms = largeConstant[].standardDeviation!(T, VarianceAlgo.assumeZeroMean)(true); + assert(abs(rms / (T.max / 4) - 1) <= 2 * T.epsilon); + T[5] identical; + identical[] = T.max; + assert(identical[].standardDeviation == 0); + T[3] zero; + zero[] = 0; + assert(zero[].standardDeviation == 0); + T[3] constant = [T(3), T(3), T(3)]; + assert(constant[].standardDeviation!(T, VarianceAlgo.assumeZeroMean)(true) == 3); + T[2] invalid = [T.infinity, T.infinity]; + assert(invalid[].standardDeviation.isNaN); + }} + double[3] tiny = [-2e-162, 0, 2e-162]; + assert(abs(tiny[].standardDeviation / 2e-162 - 1) <= 2 * double.epsilon); + const double[6] matrix = [-1e200, 0, 1e200, -1e200, 0, 1e200]; + assert(matrix[].sliced(2, 3).transposed[0].standardDeviation == 0); + const expected = 1e200 * sqrt(4.0 / 5); + assert(abs(matrix[].sliced(2, 3).transposed.standardDeviation / expected - 1) + <= 4 * double.epsilon); +} + /// ditto template standardDeviation( VarianceAlgo varianceAlgo = VarianceAlgo.hybrid, @@ -1605,11 +1799,9 @@ unittest auto y = x.standardDeviation; assert(y.approxEqual(sqrt(54.76562 / 11))); - // The naive algorithm is numerically unstable in this case - auto z0 = x.standardDeviation!"naive"; - assert(!z0.approxEqual(y)); - - // But the two-pass algorithm provides a consistent answer + // Prefer two-pass for large offsets. Naive variance can lose the spread; + // recovery catches zero or negative results, but not an inaccurate positive result. + // The two-pass algorithm provides a consistent answer auto z1 = x.standardDeviation!"twoPass"; assert(z1.approxEqual(y)); } diff --git a/source/mir/stat/descriptive/univariate/means.d b/source/mir/stat/descriptive/univariate/means.d index 5214873..600e709 100644 --- a/source/mir/stat/descriptive/univariate/means.d +++ b/source/mir/stat/descriptive/univariate/means.d @@ -475,6 +475,8 @@ unittest position = 0; assert(mean(Single(values[], &position)) == double.infinity); assert(position == 2); + + } /// From 28cab100c2953f080de9be32f432755c95545622 Mon Sep 17 00:00:00 2001 From: jmh530 Date: Sat, 17 Oct 2026 22:29:01 -0400 Subject: [PATCH 4/7] Document integer precision limits in variance and standard deviation Explain how converting large integer observations to floating-point accumulators can discard small differences and produce zero variance. Illustrate the limitation with five consecutive long values near 2^60, whose sample variance is 2.5. Describe subtracting a known common offset before conversion as a workaround, provided the integer subtraction does not overflow. Clarify that this translation is unsuitable for assumeZeroMean and that standard deviation's floating-point recovery does not repair lost integer precision. AI-Assisted-By: Codex (GPT-6 Astra Medium) --- .../stat/descriptive/univariate/dispersions.d | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/source/mir/stat/descriptive/univariate/dispersions.d b/source/mir/stat/descriptive/univariate/dispersions.d index ff437bd..e7419cc 100644 --- a/source/mir/stat/descriptive/univariate/dispersions.d +++ b/source/mir/stat/descriptive/univariate/dispersions.d @@ -1087,6 +1087,19 @@ By default, if `F` is not floating point type or complex type, then the result will have a `double` type if `F` is implicitly convertible to a floating point type or a type for which `isComplex!F` is true. +Large integer observations that differ by small amounts can lose those +differences when converted to the floating-point accumulator type. For example, +the five `long` values `(1L << 60) + [0, 1, 2, 3, 4]` have sample variance +`2.5`, but the default calculation returns zero because their conversions to +`double` are identical. A more accurate summation algorithm cannot recover +differences already lost during conversion. + +When the data has a known common offset, subtracting it in integer arithmetic +before calling `variance` can preserve the differences. Ensure that the +subtractions fit the integer type. This translation preserves variance about +the mean, but must not be used with `VarianceAlgo.assumeZeroMean`, which +calculates squared distances from zero. + Params: F = controls type of output varianceAlgo = algorithm for calculating variance (default: VarianceAlgo.hybrid) @@ -1520,6 +1533,11 @@ compensated two-pass variance, except that `assumeZeroMean` retains its zero-centered definition. Ordinary positive normal variances use the selected algorithm unchanged. Single-pass inputs do not receive this recovery. +Integer inputs retain the conversion limitation described for $(LREF variance): +large integers separated by small differences can produce a zero standard +deviation. The scaled floating-point recovery does not repair this loss of +integer precision. + Params: F = controls type of output varianceAlgo = algorithm for calculating variance (default: VarianceAlgo.hybrid) From 23a7f486b84607a236bfdff576147290f71d02e9 Mon Sep 17 00:00:00 2001 From: jmh530 Date: Sat, 17 Oct 2026 23:28:18 -0400 Subject: [PATCH 5/7] Centralize summator merging with a capability-based fallback Extract the summator merge policy into an internal mergeSummator helper. Merge matching element types directly when the operation compiles, preserving saved partial sums and corrections. Otherwise retain the existing scalar conversion path and its constraints. Use the helper in MeanAccumulator to restore compatibility with decimal summation, whose summator currently lacks a direct merge operation. Keep capability detection entirely at compile time, with counts and statistical correction formulas remaining the caller's responsibility. Document the distinct-source requirement for self merges. Add regressions for const sources, compensated-state preservation, decimal merging, and empty and self merges. Register the helper module with Meson. AI-Assisted-By: Codex (GPT-6 Astra Medium) --- meson.build | 1 + .../mir/stat/descriptive/univariate/means.d | 31 +++++++--- source/mir/stat/internal/summation.d | 58 +++++++++++++++++++ 3 files changed, 83 insertions(+), 7 deletions(-) create mode 100644 source/mir/stat/internal/summation.d diff --git a/meson.build b/meson.build index 8160418..babc82d 100644 --- a/meson.build +++ b/meson.build @@ -43,6 +43,7 @@ sources_list = [ 'mir/stat/internal/one_minus_exp', 'mir/stat/internal/shape_transform', 'mir/stat/internal/standardize', + 'mir/stat/internal/summation', 'mir/stat/descriptive/histogram/package', 'mir/stat/descriptive/histogram/api/rc', 'mir/stat/descriptive/histogram/api/gc', diff --git a/source/mir/stat/descriptive/univariate/means.d b/source/mir/stat/descriptive/univariate/means.d index 600e709..d214f50 100644 --- a/source/mir/stat/descriptive/univariate/means.d +++ b/source/mir/stat/descriptive/univariate/means.d @@ -102,12 +102,10 @@ struct MeanAccumulator(T, Summation summation) /// void put(F = T)(MeanAccumulator!(F, summation) m) { + import mir.stat.internal.summation : mergeSummator; + count += m.count; - static if (is(F == T)) - // Merge saved partials and corrections before rounding to a scalar. - summator += m.summator; - else - summator.put(cast(T) m.summator); + mergeSummator(summator, m.summator); } } @@ -179,6 +177,27 @@ unittest }} } +// Retain scalar merging for summators without a direct merge operation. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + alias Accumulator = MeanAccumulator!(double, Summation.decimal); + Accumulator source, merged; + source.put(1.0); + source.put(3.0); + merged.put(2.0); + merged.put(source); + assert(merged.count == 3 && merged.sum == 6); + assert(merged.mean == 2); + assert(source.count == 2 && source.sum == 4); + merged.put(Accumulator.init); + assert(merged.count == 3 && merged.sum == 6); + merged.put(merged); + assert(merged.count == 6 && merged.sum == 12); + assert(merged.mean == 2); +} + // Precise summation also retains overflow internally until later cancellation. // Copying this summator by value currently requires @system in mir-algorithm. version(mir_stat_test) @@ -475,8 +494,6 @@ unittest position = 0; assert(mean(Single(values[], &position)) == double.infinity); assert(position == 2); - - } /// diff --git a/source/mir/stat/internal/summation.d b/source/mir/stat/internal/summation.d new file mode 100644 index 0000000..0faf1a1 --- /dev/null +++ b/source/mir/stat/internal/summation.d @@ -0,0 +1,58 @@ +/++ +Internal summator merging for statistical accumulators. + +License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0) +Authors: John Michael Hall +Copyright: 2026 Mir Stat Authors. ++/ +module mir.stat.internal.summation; + +import mir.math.sum : Summator, Summation; + +/* +Preserve partial sums and corrections when matching summator types support a +direct merge. Otherwise retain scalar conversion to the destination element +type. In particular, decimal summators currently lack a direct merge, and +cross-type merges retain the existing conversion path and its constraints. + +Capability detection is compile-time only. Counts and any statistical merge +corrections remain the caller's responsibility. Source and target must be +distinct summator objects; callers supporting self merges must first copy the +source, as MeanAccumulator's by-value merge parameter does. The source is not +modified. +*/ +package(mir.stat) void mergeSummator(T, U, Summation mode)( + ref Summator!(T, mode) target, ref const Summator!(U, mode) source) +{ + static if (is(T == U) && __traits(compiles, target += source)) + target += source; + else + target.put(cast(T) source); +} + +// Both paths accept a const source without changing its accumulated state. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + Summator!(double, Summation.kbn) source; + source.put(0x1p54); + source.put(1.0); + const saved = source; + Summator!(double, Summation.kbn) target; + mergeSummator(target, saved); + target.put(-0x1p54); + assert(target.sum == 1); + // The source's saved correction remains intact. + source.put(-0x1p54); + assert(source.sum == 1); + + Summator!(double, Summation.decimal) decimalSource; + decimalSource.put(2.0); + const savedDecimal = decimalSource; + Summator!(double, Summation.decimal) decimalTarget; + decimalTarget.put(1.0); + mergeSummator(decimalTarget, savedDecimal); + assert(decimalTarget.sum == 3); + assert(savedDecimal.sum == 2); +} From 8393f0c1948d988afa0e716c31e1edac9b5eabcb Mon Sep 17 00:00:00 2001 From: jmh530 Date: Sun, 18 Oct 2026 17:42:08 -0400 Subject: [PATCH 6/7] Clarify the internal summation helper's purpose and visibility Explain how mergeSummator preserves summation state when direct merging is supported and falls back to scalar conversion otherwise. Clarify that callers manage observation counts and statistical corrections. Follow the univariate internal modules' documentation style by stating the helper's package visibility and linking to the public statistical functions. AI-Assisted-By: Codex (GPT-6 Astra Medium) --- source/mir/stat/internal/summation.d | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/source/mir/stat/internal/summation.d b/source/mir/stat/internal/summation.d index 0faf1a1..756ed6b 100644 --- a/source/mir/stat/internal/summation.d +++ b/source/mir/stat/internal/summation.d @@ -1,6 +1,16 @@ /++ Internal summator merging for statistical accumulators. +`mergeSummator` combines accumulated summation state. Matching summator types +are merged directly when supported, preserving partial sums and compensation +terms. Otherwise, the source sum is converted to the destination element type +before being added. Callers manage observation counts and any additional +statistical corrections. + +The helper has `package(mir.stat)` visibility and is intended for use inside +mir.stat; it is not a public application API. See +$(MREF mir,stat,descriptive,univariate) for the public statistical functions. + License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0) Authors: John Michael Hall Copyright: 2026 Mir Stat Authors. From a2fed35d1f4298ff562b710d3fec61b71fe9fc1d Mon Sep 17 00:00:00 2001 From: jmh530 Date: Sun, 18 Oct 2026 18:06:23 -0400 Subject: [PATCH 7/7] Cover standard-deviation recovery without a usable scale Exercise the recovery fallbacks through the public assumeZeroMean API, which has no minimum-count contract. Verify that empty population input and a zero-valued singleton sample preserve their NaN results rather than attempting recovery without observations or a nonzero scale. Also verify that the singleton's population standard deviation is zero. These regressions cover both previously missed fallback returns without changing the implementation or input contracts. AI-Assisted-By: Codex (GPT-6 Astra Medium) --- .../stat/descriptive/univariate/dispersions.d | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/source/mir/stat/descriptive/univariate/dispersions.d b/source/mir/stat/descriptive/univariate/dispersions.d index e7419cc..4ba5d90 100644 --- a/source/mir/stat/descriptive/univariate/dispersions.d +++ b/source/mir/stat/descriptive/univariate/dispersions.d @@ -1634,6 +1634,24 @@ unittest assert(position == 3); } +// Preserve assumeZeroMean's existing boundary results when recovery has no scale. +// Unlike the centered algorithms, this algorithm has no minimum-count contract. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + import std.math : isNaN; + + double[] empty; + // Population variance divides zero by zero; recovery must not read front. + assert(empty.standardDeviation!(VarianceAlgo.assumeZeroMean)(true).isNaN); + // A singleton sample also divides zero by zero. All observations are zero, + // so scaling cannot recover a standard deviation from the NaN variance. + double[1] zero = [0.0]; + assert(zero[].standardDeviation!(VarianceAlgo.assumeZeroMean).isNaN); + assert(zero[].standardDeviation!(VarianceAlgo.assumeZeroMean)(true) == 0); +} + // Recover standard deviations whose squared values are outside the type's range. version(mir_stat_test) @safe pure nothrow @nogc