You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
{{ message }}
Repository navigation
fix(tdigest): stabilize numerical interpolation - #287
Finite T-Digest inputs can overflow intermediate arithmetic, and inconsistent rounding of cumulative weights can make query results decrease across interval boundaries. Weighted images also retain extrema that can differ from their outer centroid means; later updates or merges must preserve a valid distribution in that state.
This PR makes those queries finite, bounded, and monotonic, reduces rounding loss when merging centroids with very unequal weights, and avoids unnecessary compression work for small buffered inputs.
Behavior and design
Queries share overflow-safe interpolation and integer cumulative-weight centers. Rank tail masses cannot pass the adjacent centroid center, so inserting a unit-weight centroid inside stored extrema no longer produces decreasing ranks or negative PMF buckets. Rank at the maximum uses the same rounded mass coordinates as interior queries.
Weighted single-centroid queries preserve stored extrema. For min=0, max=100, centroid=(50,10), the quartiles are now 18.75, 50, 81.25 instead of all 50. Mutable and frozen queries agree. Deserialization rejects total-weight-one images whose extrema differ, in both native and reference formats.
Merged means start at the heavier centroid. Same-sign means use an FMA; opposite-sign means use a scaled correction when 1 - fraction rounds to 1, avoiding an overflowing subtraction while retaining the representable contribution.
At total weight at most k / 2, no pair of positive integer weights fits the existing K_2 merge limit. Compression skips that scan while preserving sorting, tie order, and scan direction.
The interpolation sign split follows WG21 P0811R3. Exact endpoint handling remains necessary: division can round a bracketed position to an endpoint, and FMA cannot undo rounding already introduced by end - start.
Compatibility
Public signatures, the scale function, and the serialized layout are unchanged. The seven pre-change construction fixtures still round-trip byte-for-byte. Rebuilding five construction fixtures changes one centroid mean each because the mean calculation has a different rounding order; headers, weights, counts, and image lengths remain unchanged.
The final merge operands for each changed mean were checked independently with exact rational arithmetic:
Construction
Old mean → new mean
Change
Distance from correctly rounded exact merge
Full buffer, k=200, 1,641 samples
-0.866666666666667 → -0.8666666666666671
1 ULP
0 ULP
Merge 201 + 199 samples, k=10
5.568181818181818 → 5.568181818181817
2 ULP
0 ULP
Merge 8 + 201 samples
6.0588235294117645 → 6.058823529411764
1 ULP
0 ULP
Merge 201 + 8 samples
4.794117647058823 → 4.7941176470588225
1 ULP
1 ULP
Merge 201 + 401 samples
15.073333333333334 → 15.073333333333332
1 ULP
0 ULP
Validation
cargo x check, cargo x test (745 tests including doctests), and cargo x lint pass.
Pre-PR assertions and tolerances are unchanged except for the five construction fingerprints above. The 34 pre-PR T-Digest behavior-test bodies are unchanged.
Regressions cover extreme finite values, weighted single-centroid queries, both compression directions and paths, signed zeros, subnormals, large-count boundaries, and the rounded-last-center fallback. Shared fixture code handles construction and centroid extraction; numerical expectations stay in each test.
Temporary fault injection reintroduces 12 numerical failure modes. Both the original and consolidated regression suites reject all 12, including replacing the final quantile fallback with the last centroid mean and restoring the decreasing tail formulas.
A deterministic public-API driver exercises 125,248 accepted and subsequently mutated states and 22,008,711 query points with debug assertions and overflow checks enabled. It covers finite extremes, duplicate means, counts through u64::MAX, update/borrowed merge/owned collection, mutable/frozen agreement, scalar/batch ordering, and nonnegative normalized PMF. No checked invariant failed.
An exact-rational oracle checks 11,232 mean merges across signs, magnitudes, weight ratios, and directions. Sampled error stays within epsilon * weighted_absolute_mean + smallest_subnormal; all 320 sampled symmetric opposite-sign cases whose complement rounds to 1 match the correctly rounded exact mean.
Performance and limits
Small-buffer serialization improvements apply to total weights at most k / 2: earlier measurements against v0.5.0 found 53.2% less time serializing 512 buffered 64-row partials. Merge performance remains workload dependent; earlier comparisons found borrowed-merge regressions against main for some input shapes. Owned collection is not a universal remedy.
For the tail-rank fix relative to ecc4f7f, five paired release-timing blocks on an Apple M4 Max, Rust 1.99.0-nightly, LLVM 23, and k=200 gave:
Workload
Before → after
Change
Scalar rank
36.34 → 34.65 ns
-4.7%
CDF, 100 split points
3.536 → 3.515 µs
-0.6%
PMF, 100 split points
3.571 → 3.567 µs
-0.1%
Scalar quantile
82.04 → 85.87 ns
+4.7%
Three quantiles in one batch
109.15 → 108.64 ns
-0.5%
Update 100,000 values
2.945 → 2.939 ms
-0.2%
The scalar-quantile slowdown was reproduced despite its 206 machine instructions remaining identical after normalizing addresses. The changed rank function moves quantile by eight bytes. A control build with 64-byte function alignment measured 82.30 → 82.51 ns (+0.25%, paired changes -0.81% to +0.81%), supporting a code-layout explanation. No alignment tuning is added. These are library microbenchmarks, not a general speedup or application-level performance claim.
Weighted means remain floating-point approximations. For example, (1, 2^54 - 1) merged with (-smallest_subnormal, 1) can still round to 1 instead of next_down(1) at a halfway boundary. Rank and quantile retain their existing approximation conventions and are not exact inverse functions.
Require equal stored extrema when the decoded total weight is one. The existing value-range checks then require the sample to match both extrema, preventing decreasing tail ranks and negative PMF buckets.
Add regressions for compressed and buffered native f64/f32 images and both reference formats, including valid single-sample controls. Both new tests fail before the fix; all 18 T-Digest serialization tests pass afterward.
A weighted centroid can summarize several distinct samples. Only equal extrema justify returning one value for every quantile or one rank throughout the stored range.
Exercise both query forms with a ten-sample centroid spanning 0 through 100. The expected tail quantiles follow directly from interpolation between weights 1 and 5. The new regression fails on the previous commit; all 38 T-Digest behavior tests pass with the fix.
Use bounded interpolation for quantile and rank intervals. Compute inverse interpolation with half-scaled operands only when subtracting the endpoints overflows. The sign split follows WG21 P0811R3; the query and compression policies are unchanged.
Add independently derived regressions for two extreme centroids and both weighted tails. Before the fix they return an incorrect rank or an infinite quantile; afterward both pass. Supplement them with finite-bit-pattern partial-aggregation invariants without changing existing property tests. All 41 behavior tests and 132 serialization tests pass.
Accumulate counts as u64 and convert only the required centroid centers to f64. Use the same centers for the interior scan and both tails so boundary comparisons and interpolation cannot disagree after rounding.
Add regressions for quantiles that jump from 100 back to 90 near a large centroid center, and for a rank that drops by one rounding step when entering a weighted right tail. Both fail on the previous commit; all 43 behavior tests and 132 serialization tests pass after the fix. Existing large-weight tests remain unchanged.
Anchor interpolation at the heavier centroid. This avoids losing the small coefficient through cancellation when the opposite coefficient is near one, and keeps the result within the two finite means.
Add a regression that reads the compressed mean directly for weight ratios 2^20, 2^40, and 2^54 in both merge directions. The previous formula fails the independent weighted-mean reference; all 44 behavior tests and 132 serialization tests pass after the fix.
Update five construction fingerprints for changed mean rounding: four means move by one ULP and one moves by two. Lengths, headers, weights, and centroid counts are unchanged. Exact-rational checks of the final merge operands put all five new means within one ULP of the correctly rounded result, and all seven pre-change construction images still round-trip byte-for-byte.
Keep anchoring each mean at the heavier centroid. Its interpolation fraction
is at most one half, so same-sign means need neither endpoint guards nor a
result clamp. Use the fused multiply-add directly in that case and retain
the general overflow-safe interpolation for opposite signs.
Add direct mean-bound checks for adjacent values, subnormals, extreme finite
values and unequal weights in both scan directions and compression paths.
Existing behavior assertions and serialization fingerprints are unchanged.
When total weight is at most k / 2, the K_2 weight limit is below
2, so no pair of positive integer weights can merge. Retain the
sorted centroids and run the usual finalization without evaluating
the scale function or scanning each candidate pair.
Cover weighted inputs at the bound, tiny inputs with a negative
normalizer, stable signed-zero ties, and both scan directions.
Existing assertions and serialization fingerprints are unchanged.
Opposite-sign merge loses a representable centroid correction
datasketches/src/tdigest/sketch.rs:2030
The opposite-sign path can still erase a representable light-centroid contribution. For a centroid (f64::MAX, 2^54 - 1) merged with (-f64::MAX, 1), fraction is 2^-54; this branch returns f64::MAX, while the exact weighted mean rounds to f64::MAX.next_down(). The merged mean therefore remains identical to the heavy input, contrary to the stated fix. Please use an overflow-safe scaled/FMA formulation that combines both corrections before rounding, and add this sign/ratio case as a regression.
Pin the stored-maximum fallback at total weight 2^53 - 1 and its adjacent interpolated rank through mutable and frozen scalar and batch queries. The new regression fails when the leftover scan is replaced with unreachable!().
Clarify the left-tail precondition and interpolation documentation without changing query behavior or existing test assertions.
Stored extrema can outlive heavy outer centroids when new unit-weight centroids are inserted. Reserving an entire sample beyond their midpoint makes tail ranks decrease and PMF buckets negative.
Bound tail interpolation by ordered cumulative masses and round rank(max) with the shared centroid-center representation. Cover update, borrowed merge, owned collection, and large-weight rounding without changing quantile conventions or existing test expectations.
Share construction, compression, and centroid extraction across the four mean regressions while keeping their independent numerical expectations and input matrices explicit.
Combine duplicate weighted-single-centroid queries and pin rounded maximum ranks with exact up/down rounding examples. Shorten changelog entries to observable behavior. Production code and the 34 pre-PR behavior tests are unchanged.
Validation: cargo x check, cargo x test (745 tests), cargo x lint, and a before/after comparison that rejects all 12 injected numerical faults.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Finite T-Digest inputs can overflow intermediate arithmetic, and inconsistent rounding of cumulative weights can make query results decrease across interval boundaries. Weighted images also retain extrema that can differ from their outer centroid means; later updates or merges must preserve a valid distribution in that state.
This PR makes those queries finite, bounded, and monotonic, reduces rounding loss when merging centroids with very unequal weights, and avoids unnecessary compression work for small buffered inputs.
Behavior and design
min=0, max=100, centroid=(50,10), the quartiles are now18.75, 50, 81.25instead of all50. Mutable and frozen queries agree. Deserialization rejects total-weight-one images whose extrema differ, in both native and reference formats.1 - fractionrounds to1, avoiding an overflowing subtraction while retaining the representable contribution.k / 2, no pair of positive integer weights fits the existing K_2 merge limit. Compression skips that scan while preserving sorting, tie order, and scan direction.The interpolation sign split follows WG21 P0811R3. Exact endpoint handling remains necessary: division can round a bracketed position to an endpoint, and FMA cannot undo rounding already introduced by
end - start.Compatibility
Public signatures, the scale function, and the serialized layout are unchanged. The seven pre-change construction fixtures still round-trip byte-for-byte. Rebuilding five construction fixtures changes one centroid mean each because the mean calculation has a different rounding order; headers, weights, counts, and image lengths remain unchanged.
The final merge operands for each changed mean were checked independently with exact rational arithmetic:
k=200, 1,641 samples-0.866666666666667→-0.8666666666666671k=105.568181818181818→5.5681818181818176.0588235294117645→6.0588235294117644.794117647058823→4.794117647058822515.073333333333334→15.073333333333332Validation
cargo x check,cargo x test(745 tests including doctests), andcargo x lintpass.u64::MAX, update/borrowed merge/owned collection, mutable/frozen agreement, scalar/batch ordering, and nonnegative normalized PMF. No checked invariant failed.epsilon * weighted_absolute_mean + smallest_subnormal; all 320 sampled symmetric opposite-sign cases whose complement rounds to1match the correctly rounded exact mean.Performance and limits
Small-buffer serialization improvements apply to total weights at most
k / 2: earlier measurements againstv0.5.0found 53.2% less time serializing 512 buffered 64-row partials. Merge performance remains workload dependent; earlier comparisons found borrowed-merge regressions againstmainfor some input shapes. Owned collection is not a universal remedy.For the tail-rank fix relative to
ecc4f7f, five paired release-timing blocks on an Apple M4 Max, Rust1.99.0-nightly, LLVM 23, andk=200gave:The scalar-quantile slowdown was reproduced despite its 206 machine instructions remaining identical after normalizing addresses. The changed rank function moves quantile by eight bytes. A control build with 64-byte function alignment measured 82.30 → 82.51 ns (+0.25%, paired changes -0.81% to +0.81%), supporting a code-layout explanation. No alignment tuning is added. These are library microbenchmarks, not a general speedup or application-level performance claim.
Weighted means remain floating-point approximations. For example,
(1, 2^54 - 1)merged with(-smallest_subnormal, 1)can still round to1instead ofnext_down(1)at a halfway boundary. Rank and quantile retain their existing approximation conventions and are not exact inverse functions.