From 969aac1213bd5ad29beebe467d2ad784c14b70ea Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 14:12:35 +0800 Subject: [PATCH 01/14] fix(tdigest): reject inconsistent single-sample extrema 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. --- CHANGELOG.md | 2 +- datasketches/src/tdigest/sketch.rs | 15 +++- .../tests/serde_tests/tdigest.rs | 74 +++++++++++++++++++ 3 files changed, 89 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c5c6553b..c7ac8c29 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -46,7 +46,7 @@ All significant changes to this project will be documented in this file. * `TDigestMut` updates and merges now panic without changing the digest if the total weight would exceed `u64::MAX`, including in release builds. * `TDigestMut::merge` now preserves the true minimum and maximum from both inputs, including compressed digests. * `TDigest` and `TDigestMut` quantile queries no longer extrapolate above the stored maximum when large total weights cause rounding near rank `1.0`. -* `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. +* `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, inconsistent single-sample extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. ## v0.5.0 (2026-09-04) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index 5a7d1c65..edac70a8 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -765,7 +765,8 @@ impl TDigestMut { "malformed data: centroid means must be finite, within extrema, and nondecreasing", )); } - checked_weight_sum(compressed_weight, num_buffered as u64)?; + let total_weight = checked_weight_sum(compressed_weight, num_buffered as u64)?; + check_single_sample_extrema(min, max, total_weight)?; let mut buffered_values_valid = true; for bytes in buffered_payload.chunks_exact(buffered_value_bytes) { let value = if is_f32 { @@ -846,6 +847,7 @@ impl TDigestMut { "malformed data: centroid means in compat double format must be finite, within extrema, and nondecreasing", )); } + check_single_sample_extrema(min, max, total_weight)?; Ok(TDigestMut::make( k, false, @@ -900,6 +902,7 @@ impl TDigestMut { "malformed data: centroid means in compat float format must be finite, within extrema, and nondecreasing", )); } + check_single_sample_extrema(min, max, total_weight)?; Ok(TDigestMut::make( k, false, @@ -2064,6 +2067,16 @@ fn check_extrema(min: f64, max: f64, format: &'static str) -> Result<(), Error> Ok(()) } +fn check_single_sample_extrema(min: f64, max: f64, total_weight: u64) -> Result<(), Error> { + // Tail interpolation assumes that distinct extrema represent at least two samples. + if total_weight == 1 && min != max { + return Err(Error::deserial( + "malformed data: a single sample must have equal extrema", + )); + } + Ok(()) +} + fn check_nonzero(value: u64, tag: &'static str) -> Result { NonZeroU64::new(value) .ok_or_else(|| Error::deserial(format!("malformed data: {tag} cannot be zero"))) diff --git a/tests-integration/tests/serde_tests/tdigest.rs b/tests-integration/tests/serde_tests/tdigest.rs index 1b871ef8..22cf1b7b 100644 --- a/tests-integration/tests/serde_tests/tdigest.rs +++ b/tests-integration/tests/serde_tests/tdigest.rs @@ -362,6 +362,80 @@ fn assert_invalid_tdigest(bytes: &[u8]) { assert_eq!(error.kind(), datasketches::error::ErrorKind::InvalidData); } +#[test] +fn test_deserialize_single_sample_requires_equal_extrema() { + // The multiple-value format can store a single centroid or a single buffered value. + for buffered in [false, true] { + for is_f32 in [false, true] { + for (min, max) in [(0., 100.), (0., 50.), (50., 100.), (50., 50.)] { + let mut bytes = vec![2, 1, 20, 100, 0, 0, 0, 0]; + bytes.extend(u32::from(!buffered).to_le_bytes()); + bytes.extend(u32::from(buffered).to_le_bytes()); + if is_f32 { + for value in [min, max, 50.] { + bytes.extend((value as f32).to_le_bytes()); + } + if !buffered { + bytes.extend(1_u32.to_le_bytes()); + } + } else { + for value in [min, max, 50_f64] { + bytes.extend(value.to_le_bytes()); + } + if !buffered { + bytes.extend(1_u64.to_le_bytes()); + } + } + let result = if is_f32 { + TDigestMut::deserialize_f32(&bytes) + } else { + TDigestMut::deserialize(&bytes) + }; + if min == max { + let mut digest = result.unwrap(); + assert_eq!(digest.total_weight(), 1); + assert_eq!(digest.quantiles(&[0., 0.5, 1.]), Some(vec![50.; 3])); + } else { + assert_eq!( + result.unwrap_err().kind(), + datasketches::error::ErrorKind::InvalidData + ); + } + } + } + } +} + +#[test] +fn test_deserialize_compat_single_sample_requires_equal_extrema() { + for compact in [false, true] { + for (min, max) in [(0., 100.), (0., 50.), (50., 100.), (50., 50.)] { + let mut bytes = if compact { 2_u32 } else { 1_u32 }.to_be_bytes().to_vec(); + bytes.extend(f64::to_be_bytes(min)); + bytes.extend(f64::to_be_bytes(max)); + if compact { + bytes.extend(100_f32.to_be_bytes()); + bytes.extend(0_u32.to_be_bytes()); // unused capacities + bytes.extend(1_u16.to_be_bytes()); + bytes.extend(1_f32.to_be_bytes()); + bytes.extend(50_f32.to_be_bytes()); + } else { + bytes.extend(100_f64.to_be_bytes()); + bytes.extend(1_u32.to_be_bytes()); + bytes.extend(1_f64.to_be_bytes()); + bytes.extend(50_f64.to_be_bytes()); + } + if min == max { + let mut digest = TDigestMut::deserialize(&bytes).unwrap(); + assert_eq!(digest.total_weight(), 1); + assert_eq!(digest.quantiles(&[0., 0.5, 1.]), Some(vec![50.; 3])); + } else { + assert_invalid_tdigest(&bytes); + } + } + } +} + #[test] fn test_deserialize_rejects_unknown_or_conflicting_flags() { let mut unknown = serialized_two_value_digest(); From 08b5d557895c1b6b2fc1987e92cfe290502ed0a8 Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 14:13:37 +0800 Subject: [PATCH 02/14] fix(tdigest): retain extrema in single-centroid queries 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. --- CHANGELOG.md | 2 +- datasketches/src/tdigest/sketch.rs | 9 ++++---- .../tests/tdigest_test/sketch.rs | 22 +++++++++++++++++++ 3 files changed, 28 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c7ac8c29..bc039ddf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,7 +45,7 @@ All significant changes to this project will be documented in this file. * `ReqSketch::cdf` and its sorted view now preserve the rank at each split point and end at exactly `1.0`, avoiding extra rounding from accumulating PMF values. * `TDigestMut` updates and merges now panic without changing the digest if the total weight would exceed `u64::MAX`, including in release builds. * `TDigestMut::merge` now preserves the true minimum and maximum from both inputs, including compressed digests. -* `TDigest` and `TDigestMut` quantile queries no longer extrapolate above the stored maximum when large total weights cause rounding near rank `1.0`. +* `TDigest` and `TDigestMut` quantile queries no longer extrapolate above the stored maximum when large total weights cause rounding near rank `1.0`. Single-centroid queries also preserve stored extrema and interpolate their tails. * `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, inconsistent single-sample extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. ## v0.5.0 (2026-09-04) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index edac70a8..62626d56 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -1513,6 +1513,8 @@ impl TDigest { /// Computes the approximate quantile for the given normalized rank. /// + /// Ranks `0.0` and `1.0` return the stored minimum and maximum, respectively. + /// /// Returns `None` if this t-digest is empty. /// /// # Panics @@ -1643,8 +1645,7 @@ impl TDigestView<'_> { if value > self.max { return Some(1.0); } - // one centroid and value == min == max - if self.centroids.len() == 1 { + if self.min == self.max { return Some(0.5); } @@ -1774,8 +1775,8 @@ impl TDigestView<'_> { ranks: impl DoubleEndedIterator, quantiles: &mut [f64], ) { - if self.centroids.len() == 1 { - quantiles.fill(self.centroids[0].mean); + if self.min == self.max { + quantiles.fill(self.min); return; } diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index cae816bb..10ccf2d4 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -663,6 +663,28 @@ fn test_quantile_handles_two_sample_last_centroid() { assert_eq!(tdigest.quantile(0.75), Some(100.0)); } +#[test] +fn test_single_centroid_preserves_stored_tail_information() { + let digest = deserialize_with_centroids(100, 0., 100., &[(50., 10)]); + // Ten samples have a centroid center at weight 5. The tail spans weights 1 through 5, + // so weight 2.5 lies 3/8 of the way from the stored minimum to the mean. + assert_quantile_queries( + &digest, + &[0., 0.1, 0.25, 0.5, 0.75, 0.9, 1.], + &[0., 0., 18.75, 50., 81.25, 100., 100.], + ); + let digest = digest.freeze(); + for (value, rank) in [ + (0., 0.05), + (18.75, 0.25), + (50., 0.5), + (81.25, 0.75), + (100., 0.95), + ] { + assert_that!(digest.rank(value).unwrap(), near(rank, 1e-12)); + } +} + #[test] fn test_batch_quantiles_match_scalar_queries_in_input_order() { let mut tdigest = TDigestMut::new(100).unwrap(); From 20f893052713abf995d9ec214bdf3dfc338256df Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 14:15:11 +0800 Subject: [PATCH 03/14] fix(tdigest): avoid overflow in finite-value interpolation 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. --- CHANGELOG.md | 2 +- datasketches/src/tdigest/sketch.rs | 122 +++++++++++------- .../tests/tdigest_test/property.rs | 44 +++++++ .../tests/tdigest_test/sketch.rs | 38 ++++++ 4 files changed, 159 insertions(+), 47 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bc039ddf..e6d2f652 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,7 +45,7 @@ All significant changes to this project will be documented in this file. * `ReqSketch::cdf` and its sorted view now preserve the rank at each split point and end at exactly `1.0`, avoiding extra rounding from accumulating PMF values. * `TDigestMut` updates and merges now panic without changing the digest if the total weight would exceed `u64::MAX`, including in release builds. * `TDigestMut::merge` now preserves the true minimum and maximum from both inputs, including compressed digests. -* `TDigest` and `TDigestMut` quantile queries no longer extrapolate above the stored maximum when large total weights cause rounding near rank `1.0`. Single-centroid queries also preserve stored extrema and interpolate their tails. +* Fix `TDigest` and `TDigestMut` queries that could return non-finite or out-of-range results for extreme finite values or large total weights. Single-centroid queries also preserve stored extrema and interpolate their tails. * `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, inconsistent single-sample extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. ## v0.5.0 (2026-09-04) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index 62626d56..74022394 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -1655,32 +1655,29 @@ impl TDigestView<'_> { // left tail let first_mean = self.centroids[0].mean; if value < first_mean { - if first_mean - self.min > 0. { - return Some(if value == self.min { - 0.5 / centroids_weight - } else { - (1. + (((value - self.min) / (first_mean - self.min)) - * ((self.centroids[0].weight() / 2.) - 1.))) - / centroids_weight - }); - } - return Some(0.); // should never happen + return Some(if value == self.min { + 0.5 / centroids_weight + } else { + interpolate( + 1., + self.centroids[0].weight() / 2., + interpolation_fraction(value, self.min, first_mean), + ) / centroids_weight + }); } // right tail let last_mean = self.centroids[num_centroids - 1].mean; if value > last_mean { - if self.max - last_mean > 0. { - return Some(if value == self.max { - 1. - (0.5 / centroids_weight) - } else { - 1.0 - ((1.0 - + (((self.max - value) / (self.max - last_mean)) - * ((self.centroids[num_centroids - 1].weight() / 2.) - 1.))) - / centroids_weight) - }); - } - return Some(1.); // should never happen + return Some(if value == self.max { + 1. - (0.5 / centroids_weight) + } else { + 1. - interpolate( + 1., + self.centroids[num_centroids - 1].weight() / 2., + interpolation_fraction(value, self.max, last_mean), + ) / centroids_weight + }); } let mut lower = self @@ -1715,16 +1712,16 @@ impl TDigestView<'_> { } weight_delta -= self.centroids[lower].weight() / 2.; weight_delta += self.centroids[upper].weight() / 2.; - Some( - if self.centroids[upper].mean - self.centroids[lower].mean > 0. { - (weight_below - + (weight_delta * (value - self.centroids[lower].mean) - / (self.centroids[upper].mean - self.centroids[lower].mean))) - / centroids_weight - } else { - (weight_below + weight_delta / 2.) / centroids_weight - }, - ) + let fraction = if self.centroids[lower].mean < self.centroids[upper].mean { + interpolation_fraction( + value, + self.centroids[lower].mean, + self.centroids[upper].mean, + ) + } else { + 0.5 + }; + Some(interpolate(weight_below, weight_below + weight_delta, fraction) / centroids_weight) } fn quantile(&self, rank: f64) -> Option { @@ -1804,9 +1801,11 @@ impl TDigestView<'_> { quantiles[index] = if last_weight == 2. { self.max } else { - self.max - - (((centroids_weight - weight - 1.) / ((last_weight / 2.) - 1.)) - * (self.max - last.mean)) + interpolate( + self.max, + last.mean, + (centroids_weight - weight - 1.) / (last_weight / 2. - 1.), + ) }; } } @@ -1825,8 +1824,11 @@ impl TDigestView<'_> { while let Some((index, weight)) = queries.next_if(|&(_, weight)| weight < first_weight / 2.) { - quantiles[index] = self.min - + (((weight - 1.) / ((first_weight / 2.) - 1.)) * (first.mean - self.min)); + quantiles[index] = interpolate( + self.min, + first.mean, + (weight - 1.) / (first_weight / 2. - 1.), + ); } } if queries.peek().is_none() { @@ -1867,8 +1869,11 @@ impl TDigestView<'_> { // centroid, so the estimate approaches the nearer one. let distance_from_left = weight - weight_so_far - left_weight; let distance_to_right = weight_so_far + dw - weight - right_weight; - quantiles[index] = - weighted_average(left.mean, distance_to_right, right.mean, distance_from_left); + quantiles[index] = interpolate( + left.mean, + right.mean, + distance_from_left / (distance_from_left + distance_to_right), + ); } weight_so_far += dw; } @@ -2125,14 +2130,39 @@ mod scale_function { } } -fn weighted_average(x1: f64, w1: f64, x2: f64, w2: f64) -> f64 { - let total_weight = w1 + w2; - let ratio = w2 / total_weight; - if x1.is_sign_positive() != x2.is_sign_positive() { - // Subtracting opposite-signed finite extremes can overflow. - x1 * (1. - ratio) + x2 * ratio +/// Interpolates finite endpoints, preserving their bounds despite rounding. +/// +/// The sign split follows WG21 P0811R3: +/// . +fn interpolate(start: f64, end: f64, fraction: f64) -> f64 { + debug_assert!(start.is_finite() && end.is_finite() && fraction.is_finite()); + if fraction <= 0. { + return start; + } + if fraction >= 1. { + return end; + } + if start.is_sign_positive() != end.is_sign_positive() { + // Avoid the potentially overflowing difference of opposite-signed endpoints. + start * (1. - fraction) + end * fraction + } else { + // Same-sign subtraction is finite; rounding must not overshoot the end. + let value = (end - start).mul_add(fraction, start); + if start < end { + value.min(end) + } else { + value.max(end) + } + } +} + +/// Locates a value between distinct finite endpoints, in either direction. +fn interpolation_fraction(value: f64, start: f64, end: f64) -> f64 { + let width = end - start; + if width.is_finite() { + (value - start) / width } else { - // Same-sign subtraction is finite and avoids summing two near-maximum terms. - (x2 - x1).mul_add(ratio, x1) + // Scaling is needed only when opposite-signed endpoints overflow their difference. + (value * 0.5 - start * 0.5) / (end * 0.5 - start * 0.5) } } diff --git a/tests-integration/tests/tdigest_test/property.rs b/tests-integration/tests/tdigest_test/property.rs index 6e5fbc14..61577efd 100644 --- a/tests-integration/tests/tdigest_test/property.rs +++ b/tests-integration/tests/tdigest_test/property.rs @@ -51,6 +51,50 @@ fn assert_quantiles_are_monotonic(tdigest: &mut TDigestMut) { } } +#[test] +fn prop_finite_values_survive_partial_aggregation() { + fn property(bits: Vec) -> TestResult { + let mut values: Vec<_> = bits + .into_iter() + .map(f64::from_bits) + .filter(|value| value.is_finite()) + .collect(); + if values.is_empty() { + return TestResult::discard(); + } + let mut digest = TDigestMut::default(); + for chunk in values.chunks(64) { + let mut partial = TDigestMut::default(); + for &value in chunk { + partial.update(value); + } + digest.merge(&TDigestMut::deserialize(&partial.serialize()).unwrap()); + } + let mut digest = TDigestMut::deserialize(&digest.serialize()).unwrap(); + assert_eq!(digest.total_weight(), values.len() as u64); + assert_quantiles_are_monotonic(&mut digest); + + // Use the input samples as split points, independently of the computed quantiles. + values.sort_unstable_by(f64::total_cmp); + values.dedup(); + assert_eq!(digest.quantile(0.), values.first().copied()); + assert_eq!(digest.quantile(1.), values.last().copied()); + let cdf = digest.cdf(&values).unwrap(); + assert!(cdf.is_sorted()); + assert!(cdf.iter().all(|rank| (0.0..=1.0).contains(rank))); + let pmf = digest.pmf(&values).unwrap(); + assert!(pmf.iter().all(|p| (0.0..=1.0).contains(p))); + assert!((pmf.iter().sum::() - 1.).abs() < 1e-12); + TestResult::passed() + } + + QuickCheck::new() + .tests(128) + .min_tests_passed(128) + .rng(Gen::new(512)) + .quickcheck(property as fn(Vec) -> TestResult); +} + #[test] fn prop_quantile_is_non_decreasing_and_within_the_observed_range() { fn property(values: Vec) -> TestResult { diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index 10ccf2d4..c62c1a86 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -576,6 +576,44 @@ fn test_extreme_values_produce_finite_quantiles() { } } +#[test] +fn test_rank_interpolation_across_extreme_finite_values() { + for scale in [1., f64::MAX] { + let mut digest = TDigestMut::default(); + digest.update(-scale); + digest.update(scale); + let digest = digest.freeze(); + // The two centroid centers have ranks 1/4 and 3/4, independently of scale. + let points = [-scale, -scale / 2., 0., scale / 2., scale]; + for (point, expected) in points.into_iter().zip([0.25, 0.375, 0.5, 0.625, 0.75]) { + assert_that!( + digest.rank(point).unwrap(), + near(expected, 4. * f64::EPSILON) + ); + } + } +} + +#[test] +fn test_tail_interpolation_across_extreme_finite_values() { + let max = f64::MAX; + let left = deserialize_with_centroids(100, -max, max, &[(max / 2., 10), (max, 1)]); + let right = deserialize_with_centroids(100, -max, max, &[(-max, 1), (-max / 2., 10)]); + // Weight 3 is halfway from the left tail's endpoint (weight 1) to its center (weight 5). + // The expected value is therefore (-MAX + MAX/2)/2 = -MAX/4; the right tail mirrors it. + for (digest, rank, expected) in [(left, 3. / 11., -0.25), (right, 8. / 11., 0.25)] { + let digest = digest.freeze(); + assert_that!( + digest.quantile(rank).unwrap() / max, + near(expected, 4. * f64::EPSILON) + ); + assert_that!( + digest.rank(expected * max).unwrap(), + near(rank, 4. * f64::EPSILON) + ); + } +} + #[test] fn test_estimate_repeat_values() { let mut tdigest = TDigestMut::default(); From 548bec3821c83c0b7af33a613ba57a28f2b16f88 Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 14:18:10 +0800 Subject: [PATCH 04/14] fix(tdigest): keep query centers consistent across intervals 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. --- CHANGELOG.md | 2 +- datasketches/src/tdigest/sketch.rs | 100 +++++++++--------- .../tests/tdigest_test/sketch.rs | 43 ++++++++ 3 files changed, 93 insertions(+), 52 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e6d2f652..1dcfe198 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,7 +45,7 @@ All significant changes to this project will be documented in this file. * `ReqSketch::cdf` and its sorted view now preserve the rank at each split point and end at exactly `1.0`, avoiding extra rounding from accumulating PMF values. * `TDigestMut` updates and merges now panic without changing the digest if the total weight would exceed `u64::MAX`, including in release builds. * `TDigestMut::merge` now preserves the true minimum and maximum from both inputs, including compressed digests. -* Fix `TDigest` and `TDigestMut` queries that could return non-finite or out-of-range results for extreme finite values or large total weights. Single-centroid queries also preserve stored extrema and interpolate their tails. +* Fix `TDigest` and `TDigestMut` queries that could return non-finite, out-of-range, or non-monotonic results for extreme finite values or large total weights. Single-centroid queries also preserve stored extrema and interpolate their tails. * `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, inconsistent single-sample extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. ## v0.5.0 (2026-09-04) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index 74022394..2d03eb57 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -1660,7 +1660,7 @@ impl TDigestView<'_> { } else { interpolate( 1., - self.centroids[0].weight() / 2., + centroid_center(0, self.centroids[0].weight.get()), interpolation_fraction(value, self.min, first_mean), ) / centroids_weight }); @@ -1672,10 +1672,11 @@ impl TDigestView<'_> { return Some(if value == self.max { 1. - (0.5 / centroids_weight) } else { - 1. - interpolate( - 1., - self.centroids[num_centroids - 1].weight() / 2., - interpolation_fraction(value, self.max, last_mean), + let last_weight = self.centroids[num_centroids - 1].weight.get(); + interpolate( + centroid_center(self.centroids_weight - last_weight, last_weight), + (self.centroids_weight - 1) as f64, + interpolation_fraction(value, last_mean, self.max), ) / centroids_weight }); } @@ -1697,31 +1698,25 @@ impl TDigestView<'_> { upper -= 1; } - let mut weight_below = 0.; - let mut i = 0; - while i < lower { - weight_below += self.centroids[i].weight(); - i += 1; - } - weight_below += self.centroids[lower].weight() / 2.; - - let mut weight_delta = 0.; - while i < upper { - weight_delta += self.centroids[i].weight(); - i += 1; - } - weight_delta -= self.centroids[lower].weight() / 2.; - weight_delta += self.centroids[upper].weight() / 2.; - let fraction = if self.centroids[lower].mean < self.centroids[upper].mean { - interpolation_fraction( - value, - self.centroids[lower].mean, - self.centroids[upper].mean, - ) + let weight_below = self.centroids[..lower].iter().map(|c| c.weight.get()).sum(); + let weight_between: u64 = self.centroids[lower..upper] + .iter() + .map(|c| c.weight.get()) + .sum(); + let left = &self.centroids[lower]; + let right = &self.centroids[upper]; + let fraction = if left.mean < right.mean { + interpolation_fraction(value, left.mean, right.mean) } else { 0.5 }; - Some(interpolate(weight_below, weight_below + weight_delta, fraction) / centroids_weight) + Some( + interpolate( + centroid_center(weight_below, left.weight.get()), + centroid_center(weight_below + weight_between, right.weight.get()), + fraction, + ) / centroids_weight, + ) } fn quantile(&self, rank: f64) -> Option { @@ -1784,7 +1779,7 @@ impl TDigestView<'_> { .peekable(); // Consume the right tail from the back, so the centroid scan can omit all tail checks. - while let Some((index, _)) = queries.next_if(|&(_, weight)| weight > centroids_weight - 1.) + while let Some((index, _)) = queries.next_if(|&(_, weight)| centroids_weight - weight < 1.) { quantiles[index] = self.max; } @@ -1794,9 +1789,10 @@ impl TDigestView<'_> { let last = self.centroids.last().unwrap(); let last_weight = last.weight(); + let last_center = + centroid_center(self.centroids_weight - last.weight.get(), last.weight.get()); if last_weight > 1. { - while let Some((index, weight)) = - queries.next_if(|&(_, weight)| centroids_weight - weight <= last_weight / 2.) + while let Some((index, weight)) = queries.next_if(|&(_, weight)| weight >= last_center) { quantiles[index] = if last_weight == 2. { self.max @@ -1820,9 +1816,9 @@ impl TDigestView<'_> { let first = &self.centroids[0]; let first_weight = first.weight(); + let first_center = centroid_center(0, first.weight.get()); if first_weight > 1. { - while let Some((index, weight)) = - queries.next_if(|&(_, weight)| weight < first_weight / 2.) + while let Some((index, weight)) = queries.next_if(|&(_, weight)| weight < first_center) { quantiles[index] = interpolate( self.min, @@ -1836,46 +1832,42 @@ impl TDigestView<'_> { } // Answer the remaining queries by centroid interval; both streams advance only forward. - let mut weight_so_far = first_weight / 2.; + let mut weight_before = 0; + let mut left_center = first_center; for pair in self.centroids.windows(2) { let left = &pair[0]; let right = &pair[1]; - let dw = (left.weight() + right.weight()) / 2.; + let right_before = weight_before + left.weight.get(); + let right_center = centroid_center(right_before, right.weight.get()); loop { let Some(&(index, weight)) = queries.peek() else { return; }; - if weight_so_far + dw <= weight { + if right_center <= weight { break; } queries.next(); - let mut left_weight = 0.; + let mut start = left_center; if left.weight.get() == 1 { - if weight - weight_so_far < 0.5 { + start = right_before as f64; + if weight < start { quantiles[index] = left.mean; continue; } - left_weight = 0.5; } - let mut right_weight = 0.; + let mut end = right_center; if right.weight.get() == 1 { - if weight_so_far + dw - weight <= 0.5 { + end = right_before as f64; + if weight >= end { quantiles[index] = right.mean; continue; } - right_weight = 0.5; } - // Each centroid is weighted by the distance from the target to the *other* - // centroid, so the estimate approaches the nearer one. - let distance_from_left = weight - weight_so_far - left_weight; - let distance_to_right = weight_so_far + dw - weight - right_weight; - quantiles[index] = interpolate( - left.mean, - right.mean, - distance_from_left / (distance_from_left + distance_to_right), - ); + quantiles[index] = + interpolate(left.mean, right.mean, (weight - start) / (end - start)); } - weight_so_far += dw; + weight_before = right_before; + left_center = right_center; } // Rounding at large total weights can exhaust the scan near the maximum. @@ -2130,6 +2122,12 @@ mod scale_function { } } +/// Keeps cumulative counts exact until a centroid's midpoint is needed for interpolation. +/// All intervals and tails use this same center to avoid inconsistent boundary rounding. +fn centroid_center(weight_before: u64, weight: u64) -> f64 { + (weight_before + weight / 2) as f64 + (weight % 2) as f64 * 0.5 +} + /// Interpolates finite endpoints, preserving their bounds despite rounding. /// /// The sign split follows WG21 P0811R3: diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index c62c1a86..17eb799f 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -851,6 +851,49 @@ fn test_quantiles_stay_within_extrema_at_large_total_weights() { } } +#[test] +fn test_queries_stay_ordered_at_large_centroid_centers() { + for total_weight in [1_u64 << 53, 1_u64 << 54, u64::MAX] { + for centroids in [ + vec![(10., total_weight - 5), (50., 1), (90., 4)], + vec![(0., 1), (90., total_weight - 1)], + ] { + let digest = deserialize_with_centroids(100, 0., 100., ¢roids).freeze(); + // Probe either side of a heavy centroid's center and the rounded upper endpoint. + let ranks = [ + 0., + 0.5_f64.next_down(), + 0.5, + 0.5_f64.next_up(), + 1.0_f64.next_down(), + 1., + ]; + let quantiles = digest.quantiles(&ranks).unwrap(); + assert_eq!(quantiles[0], 0.); + assert_eq!(quantiles[5], 100.); + assert!( + quantiles.is_sorted(), + "centroids {centroids:?}: {quantiles:?}" + ); + assert!(quantiles.iter().all(|value| (0.0..=100.0).contains(value))); + + let cdf = digest.cdf(&[0., 5., 10., 50., 90., 95., 100.]).unwrap(); + assert!(cdf.is_sorted(), "centroids {centroids:?}: {cdf:?}"); + assert!(cdf.iter().all(|rank| (0.0..=1.0).contains(rank))); + } + } +} + +#[test] +fn test_rank_stays_monotonic_at_the_right_tail_boundary() { + let digest = deserialize_with_centroids(100, 0., 300., &[(0., 1), (100., 29)]).freeze(); + let points = [100., 100_f64.next_up(), 150., 300.]; + let cdf = digest.cdf(&points).unwrap(); + // The last center is at weight 1 + 29/2 = 15.5. Entering the tail must not lower its rank. + assert_eq!(cdf[0], 15.5 / 30.); + assert!(cdf.is_sorted(), "CDF at {points:?}: {cdf:?}"); +} + #[test] fn test_batch_quantiles_reject_invalid_ranks() { for values in [&[][..], &[1.0, 2.0, 3.0][..]] { From 6959b4c78f634536f089934ed8ebc20a75d23bbd Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 14:25:14 +0800 Subject: [PATCH 05/14] fix(tdigest): preserve light centroid contributions during merging 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. --- CHANGELOG.md | 1 + datasketches/src/tdigest/sketch.rs | 18 +++-------- .../tests/serde_tests/tdigest.rs | 10 +++--- .../tests/tdigest_test/sketch.rs | 32 +++++++++++++++++++ 4 files changed, 42 insertions(+), 19 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1dcfe198..712fed5d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -46,6 +46,7 @@ All significant changes to this project will be documented in this file. * `TDigestMut` updates and merges now panic without changing the digest if the total weight would exceed `u64::MAX`, including in release builds. * `TDigestMut::merge` now preserves the true minimum and maximum from both inputs, including compressed digests. * Fix `TDigest` and `TDigestMut` queries that could return non-finite, out-of-range, or non-monotonic results for extreme finite values or large total weights. Single-centroid queries also preserve stored extrema and interpolate their tails. +* T-Digest compression no longer discards the contribution of a lightly weighted centroid when merging very unequal weights. * `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, inconsistent single-sample extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. ## v0.5.0 (2026-09-04) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index 2d03eb57..a0046ef9 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -2007,22 +2007,12 @@ impl Centroid { .checked_add(other.weight.get()) .expect("weight overflow"); - let (self_mean, other_mean) = (self.mean, other.mean); - let ratio_other = other_weight / total_weight; - let delta = other_mean - self_mean; - self.mean = if delta.is_finite() { - delta.mul_add(ratio_other, self_mean) + // Start at the heavier centroid so a small contribution is not rounded away in 1 - ratio. + self.mean = if self_weight >= other_weight { + interpolate(self.mean, other.mean, other_weight / total_weight) } else { - let ratio_self = self_weight / total_weight; - self_mean.mul_add(ratio_self, other_mean * ratio_other) + interpolate(other.mean, self.mean, self_weight / total_weight) }; - - debug_assert!( - self.mean.is_finite(), - "Centroid's mean must be finite; self: {}, other: {}", - self_mean, - other_mean - ); } fn weight(&self) -> f64 { diff --git a/tests-integration/tests/serde_tests/tdigest.rs b/tests-integration/tests/serde_tests/tdigest.rs index 22cf1b7b..b90b1ca6 100644 --- a/tests-integration/tests/serde_tests/tdigest.rs +++ b/tests-integration/tests/serde_tests/tdigest.rs @@ -247,7 +247,7 @@ fn test_serialized_bytes_stable_for_full_and_merged_digests() { let mut full_buffer = patterned_digest(200, 1_641, 0); let bytes = full_buffer.serialize(); assert_eq!(bytes.len(), 2_864); - assert_eq!(fnv1a(&bytes), 0x5c01_c50d_d1c8_fdbb); + assert_eq!(fnv1a(&bytes), 0xaffe_62f2_3716_59b8); let mut left = patterned_digest(10, 201, 2); let mut right = patterned_digest(10, 199, 3); @@ -255,12 +255,12 @@ fn test_serialized_bytes_stable_for_full_and_merged_digests() { left.merge(&right); let bytes = left.serialize(); assert_eq!(bytes.len(), 272); - assert_eq!(fnv1a(&bytes), 0x7d2e_a927_9b9e_f559); + assert_eq!(fnv1a(&bytes), 0x6037_b570_d640_d31f); for &(left_len, right_len, expected_len, expected_hash) in &[ - (8, 201, 272, 0x8522_1f3f_152f_24e5), - (201, 8, 256, 0xe60d_1f6f_f4b0_73e0), - (201, 401, 288, 0x4cb8_4037_5e68_ca4b), + (8, 201, 272, 0x4ec5_2661_c556_c10e), + (201, 8, 256, 0x2b43_056d_9a57_b2a5), + (201, 401, 288, 0x1e43_34b7_ba13_4fe0), (401, 201, 288, 0x6f6d_e965_77a7_a53f), ] { let mut left = patterned_digest(10, left_len, 2); diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index 17eb799f..06b42949 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -723,6 +723,38 @@ fn test_single_centroid_preserves_stored_tail_information() { } } +#[test] +fn test_compression_preserves_a_small_centroid_weight() { + for weight in [1_u64 << 20, 1_u64 << 40, 1_u64 << 54] { + for reverse in [false, true] { + let mut digest = deserialize_with_centroids( + 10, + 0., + f64::MAX, + &[(0., weight), (1., weight), (1e300, 1), (f64::MAX, weight)], + ); + let mut bytes = digest.serialize(); + if reverse { + bytes[5] |= 1 << 2; + } + let mut digest = TDigestMut::deserialize(&bytes).unwrap(); + digest.update(0.); + + // Read the merged mean directly: quantile interpolation must not hide a bad mean. + let bytes = digest.serialize(); + let centroid = bytes[32..] + .chunks_exact(16) + .find(|c| u64::from_le_bytes(c[8..].try_into().unwrap()) == weight + 1) + .expect("the middle centroids should merge"); + let mean = f64::from_le_bytes(centroid[..8].try_into().unwrap()); + // The exact mean is (weight * 1 + 1e300) / (weight + 1). The first term is far + // below one ULP of the numerator, so this division is an independent reference. + let expected = 1e300 / (weight + 1) as f64; + assert_that!(mean / expected, near(1., 4. * f64::EPSILON)); + } + } +} + #[test] fn test_batch_quantiles_match_scalar_queries_in_input_order() { let mut tdigest = TDigestMut::new(100).unwrap(); From c802adf3287441c6adaef523403fa959226f7b43 Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 15:26:28 +0800 Subject: [PATCH 06/14] perf(tdigest): specialize centroid mean interpolation 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. --- datasketches/src/tdigest/sketch.rs | 16 ++++- .../tests/tdigest_test/sketch.rs | 70 +++++++++++++++++++ 2 files changed, 83 insertions(+), 3 deletions(-) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index a0046ef9..26e5febe 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -2008,10 +2008,20 @@ impl Centroid { .expect("weight overflow"); // Start at the heavier centroid so a small contribution is not rounded away in 1 - ratio. - self.mean = if self_weight >= other_weight { - interpolate(self.mean, other.mean, other_weight / total_weight) + let (start, end, fraction) = if self_weight >= other_weight { + (self.mean, other.mean, other_weight / total_weight) } else { - interpolate(other.mean, self.mean, self_weight / total_weight) + (other.mean, self.mean, self_weight / total_weight) + }; + debug_assert!(start.is_finite() && end.is_finite() && fraction > 0. && fraction <= 0.5); + + // Same-sign means within a factor of two in magnitude subtract exactly. For wider + // gaps, a step of at most 1/2 cannot cross either endpoint even if subtraction rounds. + // Queries can interpolate farther and still need the general bounds. + self.mean = if start.is_sign_positive() == end.is_sign_positive() { + (end - start).mul_add(fraction, start) + } else { + interpolate(start, end, fraction) }; } diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index 06b42949..724479c5 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -755,6 +755,76 @@ fn test_compression_preserves_a_small_centroid_weight() { } } +#[test] +fn test_compression_keeps_merged_means_within_their_endpoints() { + let tiny = f64::from_bits(1); + let heavy = 1_u64 << 54; + for (left, right) in [ + (0., tiny), + (tiny, 2. * tiny), + ( + f64::from_bits(f64::MIN_POSITIVE.to_bits() - 1), + f64::MIN_POSITIVE, + ), + (1., f64::from_bits(1_f64.to_bits() + 1)), + (f64::from_bits(f64::MAX.to_bits() - 1), f64::MAX), + (f64::MAX / 4., f64::MAX), + (-f64::MAX, -f64::MAX / 4.), + (-2. * tiny, -tiny), + (-tiny, tiny), + (-f64::MAX, f64::MAX), + ] { + for (left_weight, right_weight) in [(1, 1), (1, heavy), (heavy, 1)] { + for reverse in [false, true] { + for owned in [false, true] { + let mut digest = deserialize_with_centroids( + 10, + -f64::MAX, + f64::MAX, + &[ + (-f64::MAX, heavy), + (left, left_weight), + (right, right_weight), + (f64::MAX, heavy), + ], + ); + let mut bytes = digest.serialize(); + if reverse { + bytes[5] |= 1 << 2; + } + let mut digest = TDigestMut::deserialize(&bytes).unwrap(); + if owned { + let other = + deserialize_with_centroids(10, f64::MAX, f64::MAX, &[(f64::MAX, 1)]); + digest = [digest, other].into_iter().collect(); + } else { + digest.update(-f64::MAX); + } + + // Inspect compression directly; query interpolation could hide an overshoot. + let bytes = digest.serialize(); + let centroid = bytes[32..] + .chunks_exact(16) + .find(|c| { + u64::from_le_bytes(c[8..].try_into().unwrap()) + == left_weight + right_weight + }) + .unwrap_or_else(|| { + panic!( + "middle pair {left:?}:{left_weight}, {right:?}:{right_weight} did not merge (reverse={reverse}, owned={owned})" + ) + }); + let mean = f64::from_le_bytes(centroid[..8].try_into().unwrap()); + assert!( + mean.is_finite() && (left..=right).contains(&mean), + "{mean:?} is outside [{left:?}, {right:?}] (weights={left_weight}:{right_weight}, reverse={reverse}, owned={owned})" + ); + } + } + } + } +} + #[test] fn test_batch_quantiles_match_scalar_queries_in_input_order() { let mut tdigest = TDigestMut::new(100).unwrap(); From d54fac77c35824d05e2dfd1fe87a013675ad9f54 Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 15:46:33 +0800 Subject: [PATCH 07/14] perf(tdigest): skip compression scans that cannot merge 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. --- CHANGELOG.md | 1 + datasketches/src/tdigest/sketch.rs | 54 +++++++++++-------- .../tests/tdigest_test/sketch.rs | 43 +++++++++++++++ 3 files changed, 75 insertions(+), 23 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 712fed5d..7183c1ce 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,7 @@ All significant changes to this project will be documented in this file. ### Improvements +* Serializing buffered `TDigestMut` inputs is faster when their total weight is at most `k / 2`. * Improve the readability of `Debug` output for HLL and CPC sketches and unions. * `KllSketch::quantile` and `ReqSketch::quantile` are faster at ranks `0.0` and `1.0`, with query work independent of the number of retained samples. * `BloomFilter::insert` is faster for integer and raw-byte inputs. `BloomFilter::contains_and_insert` is also faster when checking already-present integer values. diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index 26e5febe..dd92788e 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -942,30 +942,38 @@ impl TDigestMut { let mut num_centroids = 1; let len = centroids.len(); - let compressed_weight = self.compressed_weight as f64; - let normalizer = scale_function::normalizer(2.0 * f64::from(self.k), compressed_weight); - let mut current = 1; - let mut weight_so_far = 0.; - while current < len { - let c = centroids[current]; - let proposed_weight = centroids[num_centroids - 1].weight() + c.weight(); - if should_merge_centroid( - current, - len, - weight_so_far, - proposed_weight, - compressed_weight, - normalizer, - ) { - // merge into existing centroid - centroids[num_centroids - 1].add(c); - } else { - // copy to a new centroid - weight_so_far += centroids[num_centroids - 1].weight(); - centroids[num_centroids] = c; - num_centroids += 1; + // No two positive integer weights fit when n <= k / 2. For a positive K_2 + // normalizer, the maximum limit is n / (4 * normalizer) = r * (ln(r) + 6), + // where r = n / (2 * k) <= 1/4: the limit stays below 1.154, hence below 2. + // A nonpositive normalizer also rejects all merges at these small counts. + if total_weight <= u64::from(self.k) / 2 { + num_centroids = len; + } else { + let compressed_weight = self.compressed_weight as f64; + let normalizer = scale_function::normalizer(2.0 * f64::from(self.k), compressed_weight); + let mut current = 1; + let mut weight_so_far = 0.; + while current < len { + let c = centroids[current]; + let proposed_weight = centroids[num_centroids - 1].weight() + c.weight(); + if should_merge_centroid( + current, + len, + weight_so_far, + proposed_weight, + compressed_weight, + normalizer, + ) { + // merge into existing centroid + centroids[num_centroids - 1].add(c); + } else { + // copy to a new centroid + weight_so_far += centroids[num_centroids - 1].weight(); + centroids[num_centroids] = c; + num_centroids += 1; + } + current += 1; } - current += 1; } centroids.truncate(num_centroids); diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index 724479c5..aee86aa1 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -723,6 +723,49 @@ fn test_single_centroid_preserves_stored_tail_information() { } } +#[test] +fn test_small_compression_preserves_weighted_centroids_and_tie_order() { + // Exercise the no-merge bound at two k values and a tiny input whose K_2 normalizer + // is negative. Existing weighted centroids must remain intact in every case. + for (k, middle_weight) in [(200, 95), (u16::MAX, 2), (u16::MAX, 32_762)] { + for reverse in [false, true] { + let mut digest = deserialize_with_centroids( + k, + -10., + 10., + &[(-10., 1), (-0., 1), (0., middle_weight), (1., 1), (10., 1)], + ); + let mut bytes = digest.serialize(); + if reverse { + bytes[5] |= 1 << 2; + } + let mut digest = TDigestMut::deserialize(&bytes).unwrap(); + digest.update(0.); + + // Buffered values precede retained centroids on ties, including signed zeros. + let mut expected = deserialize_with_centroids( + k, + -10., + 10., + &[ + (-10., 1), + (0., 1), + (-0., 1), + (0., middle_weight), + (1., 1), + (10., 1), + ], + ) + .serialize(); + if !reverse { + expected[5] |= 1 << 2; + } + assert_eq!(digest.total_weight(), middle_weight + 5); + assert_eq!(digest.serialize(), expected); + } + } +} + #[test] fn test_compression_preserves_a_small_centroid_weight() { for weight in [1_u64 << 20, 1_u64 << 40, 1_u64 << 54] { From 7257a99fb1e0424e2adf24cc8a71d009579fb16d Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 19:36:48 +0800 Subject: [PATCH 08/14] fix(tdigest): preserve weighted tails in mutable rank queries --- CHANGELOG.md | 2 +- datasketches/src/tdigest/sketch.rs | 3 +-- tests-integration/tests/tdigest_test/sketch.rs | 18 ++++++++++++++++++ 3 files changed, 20 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7183c1ce..9890b89a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -46,7 +46,7 @@ All significant changes to this project will be documented in this file. * `ReqSketch::cdf` and its sorted view now preserve the rank at each split point and end at exactly `1.0`, avoiding extra rounding from accumulating PMF values. * `TDigestMut` updates and merges now panic without changing the digest if the total weight would exceed `u64::MAX`, including in release builds. * `TDigestMut::merge` now preserves the true minimum and maximum from both inputs, including compressed digests. -* Fix `TDigest` and `TDigestMut` queries that could return non-finite, out-of-range, or non-monotonic results for extreme finite values or large total weights. Single-centroid queries also preserve stored extrema and interpolate their tails. +* Fix `TDigest` and `TDigestMut` queries that could return non-finite, out-of-range, or non-monotonic results for extreme finite values or large total weights. Mutable and frozen single-centroid queries also preserve stored extrema and interpolate their tails. * T-Digest compression no longer discards the contribution of a lightly weighted centroid when merging very unequal weights. * `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, inconsistent single-sample extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index dd92788e..b8e685a6 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -490,8 +490,7 @@ impl TDigestMut { if value > self.max { return Some(1.0); } - // one centroid and value == min == max - if self.buffer.len() == 1 { + if self.min == self.max { return Some(0.5); } diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index aee86aa1..9c5fa5fa 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -723,6 +723,24 @@ fn test_single_centroid_preserves_stored_tail_information() { } } +#[test] +fn test_mutable_rank_preserves_weighted_single_centroid_tails() { + let mut digest = deserialize_with_centroids(100, 0., 100., &[(50., 10)]); + let frozen = digest.clone().freeze(); + let points = [-1., 0., 18.75, 50., 81.25, 100., 101.]; + let expected = [0., 0.05, 0.25, 0.5, 0.75, 0.95, 1.]; + + // A quarter rank is mass 2.5: 3/8 of the tail from mass 1 at min to mass 5 at mean. + // Repeating the queries must preserve the result before and after creating a view. + for _ in 0..2 { + for (&point, &rank) in points.iter().zip(&expected) { + assert_eq!(digest.rank(point), Some(rank)); + assert_eq!(frozen.rank(point), Some(rank)); + } + assert_eq!(&digest.cdf(&points).unwrap()[..points.len()], &expected); + } +} + #[test] fn test_small_compression_preserves_weighted_centroids_and_tie_order() { // Exercise the no-merge bound at two k values and a tiny input whose K_2 normalizer From 715d1f8c5957317ded0f681c483a7eb857f35039 Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 19:41:09 +0800 Subject: [PATCH 09/14] fix(tdigest): align tail quantiles with cumulative centroid centers --- CHANGELOG.md | 2 +- datasketches/src/tdigest/sketch.rs | 29 ++++++++-------- .../tests/tdigest_test/sketch.rs | 33 +++++++++++++++++++ 3 files changed, 48 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9890b89a..ca853671 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -46,7 +46,7 @@ All significant changes to this project will be documented in this file. * `ReqSketch::cdf` and its sorted view now preserve the rank at each split point and end at exactly `1.0`, avoiding extra rounding from accumulating PMF values. * `TDigestMut` updates and merges now panic without changing the digest if the total weight would exceed `u64::MAX`, including in release builds. * `TDigestMut::merge` now preserves the true minimum and maximum from both inputs, including compressed digests. -* Fix `TDigest` and `TDigestMut` queries that could return non-finite, out-of-range, or non-monotonic results for extreme finite values or large total weights. Mutable and frozen single-centroid queries also preserve stored extrema and interpolate their tails. +* Fix `TDigest` and `TDigestMut` queries that could return non-finite, out-of-range, or non-monotonic results for extreme finite values or large total weights. Tail quantiles now use the same centroid centers as rank queries, including above `2^53` total weight. Mutable and frozen single-centroid queries also preserve stored extrema and interpolate their tails. * T-Digest compression no longer discards the contribution of a lightly weighted centroid when merging very unequal weights. * `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, inconsistent single-sample extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index b8e685a6..f0904721 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -1795,19 +1795,21 @@ impl TDigestView<'_> { } let last = self.centroids.last().unwrap(); - let last_weight = last.weight(); - let last_center = - centroid_center(self.centroids_weight - last.weight.get(), last.weight.get()); - if last_weight > 1. { + let last_weight = last.weight.get(); + let last_center = centroid_center(self.centroids_weight - last_weight, last_weight); + if last_weight > 1 { + let tail_end = (self.centroids_weight - 1) as f64; while let Some((index, weight)) = queries.next_if(|&(_, weight)| weight >= last_center) { - quantiles[index] = if last_weight == 2. { + // Integer counts can round a wider tail to zero width as well. At this + // boundary prefer the maximum, as for an exact two-sample tail. + quantiles[index] = if last_center == tail_end { self.max } else { interpolate( - self.max, last.mean, - (centroids_weight - weight - 1.) / (last_weight / 2. - 1.), + self.max, + (weight - last_center) / (tail_end - last_center), ) }; } @@ -1822,16 +1824,13 @@ impl TDigestView<'_> { } let first = &self.centroids[0]; - let first_weight = first.weight(); - let first_center = centroid_center(0, first.weight.get()); - if first_weight > 1. { + let first_weight = first.weight.get(); + let first_center = centroid_center(0, first_weight); + if first_weight > 1 { while let Some((index, weight)) = queries.next_if(|&(_, weight)| weight < first_center) { - quantiles[index] = interpolate( - self.min, - first.mean, - (weight - 1.) / (first_weight / 2. - 1.), - ); + quantiles[index] = + interpolate(self.min, first.mean, (weight - 1.) / (first_center - 1.)); } } if queries.peek().is_none() { diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index 9c5fa5fa..ba2085e6 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -1047,6 +1047,39 @@ fn test_queries_stay_ordered_at_large_centroid_centers() { } } +#[test] +fn test_quantile_right_tail_uses_the_same_center_as_rank() { + for (total_weight, last_weight, rank) in [ + (1_u64 << 53, 5, 1. - f64::EPSILON), + ((1_u64 << 53) + 1, 4, 1_f64.next_down()), + (1_u64 << 54, 9, 1. - f64::EPSILON), + (u64::MAX, 4097, 1_f64.next_down()), + ] { + let mut digest = deserialize_with_centroids( + 100, + 0., + 100., + &[(0., total_weight - last_weight), (90., last_weight)], + ); + // These ranks land exactly on the rounded mass N - last_weight/2. Quantile must + // return the centroid mean there, before interpolating toward the stored maximum. + assert_eq!(digest.rank(90.), Some(rank)); + assert_quantile_queries(&digest, &[rank, rank, 1.], &[90., 90., 100.]); + } + + // At N = 2^53 + 2 the last center and N - 1 round to the same mass. The zero-width + // tail must stay finite and ordered, just as a two-sample tail does at smaller counts. + let total_weight = (1_u64 << 53) + 2; + let digest = deserialize_with_centroids(100, 0., 100., &[(0., total_weight - 3), (90., 3)]); + let rank = 1. - f64::EPSILON; + let ranks = [rank.next_down(), rank, rank.next_up(), 1.]; + let quantiles = digest.clone().quantiles(&ranks).unwrap(); + assert!(quantiles.is_sorted()); + assert!(quantiles.iter().all(|value| (0.0..=100.0).contains(value))); + assert_eq!(quantiles[3], 100.); + assert_quantile_queries(&digest, &ranks, &quantiles); +} + #[test] fn test_rank_stays_monotonic_at_the_right_tail_boundary() { let digest = deserialize_with_centroids(100, 0., 300., &[(0., 1), (100., 29)]).freeze(); From 641abff7eb0d18f2e569ef00af521c9a75c22deb Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 20:01:07 +0800 Subject: [PATCH 10/14] test(tdigest): pin interpolation endpoints and signed-zero behavior --- datasketches/src/tdigest/sketch.rs | 23 ++++++++++- .../tests/tdigest_test/sketch.rs | 41 +++++++++++++++++++ 2 files changed, 62 insertions(+), 2 deletions(-) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index f0904721..c64dc684 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -2134,12 +2134,14 @@ fn centroid_center(weight_before: u64, weight: u64) -> f64 { (weight_before + weight / 2) as f64 + (weight % 2) as f64 * 0.5 } -/// Interpolates finite endpoints, preserving their bounds despite rounding. +/// Interpolates finite endpoints at a fraction in `[0, 1]`, preserving their bounds. /// /// The sign split follows WG21 P0811R3: /// . fn interpolate(start: f64, end: f64, fraction: f64) -> f64 { - debug_assert!(start.is_finite() && end.is_finite() && fraction.is_finite()); + debug_assert!(start.is_finite() && end.is_finite() && (0.0..=1.0).contains(&fraction)); + // Callers guarantee the fraction's range, but division may round it to an endpoint. + // Even an FMA need not reproduce `end` at 1: its subtracted difference already rounded. if fraction <= 0. { return start; } @@ -2170,3 +2172,20 @@ fn interpolation_fraction(value: f64, start: f64, end: f64) -> f64 { (value * 0.5 - start * 0.5) / (end * 0.5 - start * 0.5) } } + +#[cfg(test)] +mod tests { + use super::interpolate; + + #[test] + fn interpolation_preserves_endpoints_after_subtraction_rounds() { + let start = f64::EPSILON / 2.; + let end = 1_f64.next_up(); + // end - start rounds to 1, and even an FMA at fraction 1 then returns 1, not end. + // The endpoint branch also retains the sign of a zero endpoint. + for (start, end) in [(start, end), (end, start), (-0., 1.), (-1., 0.)] { + assert_eq!(interpolate(start, end, 0.).to_bits(), start.to_bits()); + assert_eq!(interpolate(start, end, 1.).to_bits(), end.to_bits()); + } + } +} diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index ba2085e6..d68aa82b 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -886,6 +886,47 @@ fn test_compression_keeps_merged_means_within_their_endpoints() { } } +#[test] +fn test_signed_zero_and_symmetric_centroids_survive_interpolation() { + for magnitude in [0., f64::from_bits(1), 1., f64::MAX] { + let heavy = 1_u64 << 54; + for reverse in [false, true] { + let mut digest = deserialize_with_centroids( + 10, + -f64::MAX, + f64::MAX, + &[ + (-f64::MAX, heavy), + (-magnitude, 1), + (magnitude, 1), + (f64::MAX, heavy), + ], + ); + let mut bytes = digest.serialize(); + if reverse { + bytes[5] |= 1 << 2; + } + let mut digest = TDigestMut::deserialize(&bytes).unwrap(); + digest.update(-f64::MAX); + let bytes = digest.serialize(); + let merged = bytes[32..] + .chunks_exact(16) + .find(|c| u64::from_le_bytes(c[8..].try_into().unwrap()) == 2) + .expect("the equally weighted middle centroids should merge"); + // Equal and opposite contributions cancel exactly, including subnormal values. + assert_eq!(f64::from_le_bytes(merged[..8].try_into().unwrap()), 0.); + } + } + + let mut digest = + deserialize_with_centroids(100, -1., 1., &[(-1., 1), (-0., 2), (0., 2), (1., 1)]); + assert_eq!(digest.rank(-0.), Some(0.5)); + assert_eq!(digest.rank(0.), Some(0.5)); + assert_quantile_queries(&digest, &[1. / 3., 0.5, 2. / 3.], &[0., 0., 0.]); + assert_eq!(digest.cdf(&[-0.]), Some(vec![0.5, 1.])); + assert_eq!(digest.pmf(&[0.]), Some(vec![0.5, 0.5])); +} + #[test] fn test_batch_quantiles_match_scalar_queries_in_input_order() { let mut tdigest = TDigestMut::new(100).unwrap(); From e5aced9de285db869bf6aa24f47b740b7651ddde Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 21:01:19 +0800 Subject: [PATCH 11/14] fix(tdigest): retain opposite-sign corrections at tiny weight ratios --- CHANGELOG.md | 2 +- datasketches/src/tdigest/sketch.rs | 5 ++ .../tests/tdigest_test/sketch.rs | 67 +++++++++++++++++++ 3 files changed, 73 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ca853671..72852d54 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,7 +47,7 @@ All significant changes to this project will be documented in this file. * `TDigestMut` updates and merges now panic without changing the digest if the total weight would exceed `u64::MAX`, including in release builds. * `TDigestMut::merge` now preserves the true minimum and maximum from both inputs, including compressed digests. * Fix `TDigest` and `TDigestMut` queries that could return non-finite, out-of-range, or non-monotonic results for extreme finite values or large total weights. Tail quantiles now use the same centroid centers as rank queries, including above `2^53` total weight. Mutable and frozen single-centroid queries also preserve stored extrema and interpolate their tails. -* T-Digest compression no longer discards the contribution of a lightly weighted centroid when merging very unequal weights. +* T-Digest compression reduces rounding loss in centroid means when merging very unequal weights, including opposite-sign centroids. * `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, inconsistent single-sample extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. ## v0.5.0 (2026-09-04) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index c64dc684..07adc912 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -2026,6 +2026,11 @@ impl Centroid { // Queries can interpolate farther and still need the general bounds. self.mean = if start.is_sign_positive() == end.is_sign_positive() { (end - start).mul_add(fraction, start) + } else if 1. - fraction == 1. { + // Both corrections can matter even when 1 - fraction rounds to 1. + // Halving before subtraction avoids overflow; the FMA combines the + // correction with start before rounding the result. + (end * 0.5 - start * 0.5).mul_add(2. * fraction, start) } else { interpolate(start, end, fraction) }; diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index d68aa82b..3e7b3ddd 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -816,6 +816,73 @@ fn test_compression_preserves_a_small_centroid_weight() { } } +#[test] +fn test_opposite_sign_merge_preserves_a_representable_correction() { + let outer_weight = 1_u64 << 56; + for magnitude in [1., f64::MAX] { + for exponent in [54, 55] { + let total_weight = 1_u64 << exponent; + // The exact mean is magnitude * (1 - 2 / total_weight). At 2^54 it + // rounds to the next value toward zero; at 2^55 it rounds back to magnitude. + let expected = if exponent == 54 { + magnitude.next_down() + } else { + magnitude + }; + for negative_heavy in [false, true] { + let (left_weight, right_weight, expected) = if negative_heavy { + (total_weight - 1, 1, -expected) + } else { + (1, total_weight - 1, expected) + }; + for reverse in [false, true] { + for owned in [false, true] { + let mut digest = deserialize_with_centroids( + 10, + -f64::MAX, + f64::MAX, + &[ + (-f64::MAX, outer_weight), + (-magnitude, left_weight), + (magnitude, right_weight), + (f64::MAX, outer_weight), + ], + ); + let mut bytes = digest.serialize(); + if reverse { + bytes[5] |= 1 << 2; + } + let mut digest = TDigestMut::deserialize(&bytes).unwrap(); + if owned { + let other = deserialize_with_centroids( + 10, + f64::MAX, + f64::MAX, + &[(f64::MAX, 1)], + ); + digest = [digest, other].into_iter().collect(); + } else { + digest.update(-f64::MAX); + } + let bytes = digest.serialize(); + let centroid = bytes[32..] + .chunks_exact(16) + .find(|c| u64::from_le_bytes(c[8..].try_into().unwrap()) == total_weight) + .unwrap_or_else(|| { + panic!("middle pair did not merge: magnitude={magnitude:?}, exponent={exponent}, negative_heavy={negative_heavy}, reverse={reverse}, owned={owned}") + }); + let mean = f64::from_le_bytes(centroid[..8].try_into().unwrap()); + assert_eq!( + mean, expected, + "magnitude={magnitude:?}, exponent={exponent}, negative_heavy={negative_heavy}, reverse={reverse}, owned={owned}" + ); + } + } + } + } + } +} + #[test] fn test_compression_keeps_merged_means_within_their_endpoints() { let tiny = f64::from_bits(1); From ecc4f7f32c91573dead0bb4edca2953000834ac2 Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 21:07:45 +0800 Subject: [PATCH 12/14] test(tdigest): cover quantile fallback at a rounded last center 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. --- datasketches/src/tdigest/sketch.rs | 3 ++- .../tests/tdigest_test/sketch.rs | 21 +++++++++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index 07adc912..dbfa3d2a 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -1827,6 +1827,7 @@ impl TDigestView<'_> { let first_weight = first.weight.get(); let first_center = centroid_center(0, first_weight); if first_weight > 1 { + // At first_weight == 2, the weight < 1 pre-pass already consumed this interval. while let Some((index, weight)) = queries.next_if(|&(_, weight)| weight < first_center) { quantiles[index] = @@ -2167,7 +2168,7 @@ fn interpolate(start: f64, end: f64, fraction: f64) -> f64 { } } -/// Locates a value between distinct finite endpoints, in either direction. +/// Locates a value between distinct finite endpoints. fn interpolation_fraction(value: f64, start: f64, end: f64) -> f64 { let width = end - start; if width.is_finite() { diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index 3e7b3ddd..4827764e 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -1188,6 +1188,27 @@ fn test_quantile_right_tail_uses_the_same_center_as_rank() { assert_quantile_queries(&digest, &ranks, &quantiles); } +#[test] +fn test_quantile_uses_max_when_last_center_rounds_down() { + let total_weight = (1_u64 << 53) - 1; + let first_weight = (total_weight - 1) / 2; + let middle_weight = total_weight - first_weight - 1; + let digest = deserialize_with_centroids( + 100, + 0., + 100., + &[(10., first_weight), (50., middle_weight), (90., 1)], + ); + + // The last center's +0.5 ties down to the even count 2^53 - 2. The highest rank + // below 1 lands there and exhausts the scan; the adjacent rank stays inside it. + assert_quantile_queries( + &digest, + &[1. - f64::EPSILON, 1_f64.next_down()], + &[90_f64.next_down(), 100.], + ); +} + #[test] fn test_rank_stays_monotonic_at_the_right_tail_boundary() { let digest = deserialize_with_centroids(100, 0., 300., &[(0., 1), (100., 29)]).freeze(); From 0cf1b6ee98eedb23a1c7087f0427d7973ed03470 Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 21:52:19 +0800 Subject: [PATCH 13/14] fix(tdigest): keep tail rank masses ordered after aggregation 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. --- CHANGELOG.md | 2 +- datasketches/src/tdigest/sketch.rs | 32 ++++--- .../tests/tdigest_test/sketch.rs | 92 +++++++++++++++++++ 3 files changed, 111 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 72852d54..037284dc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -46,7 +46,7 @@ All significant changes to this project will be documented in this file. * `ReqSketch::cdf` and its sorted view now preserve the rank at each split point and end at exactly `1.0`, avoiding extra rounding from accumulating PMF values. * `TDigestMut` updates and merges now panic without changing the digest if the total weight would exceed `u64::MAX`, including in release builds. * `TDigestMut::merge` now preserves the true minimum and maximum from both inputs, including compressed digests. -* Fix `TDigest` and `TDigestMut` queries that could return non-finite, out-of-range, or non-monotonic results for extreme finite values or large total weights. Tail quantiles now use the same centroid centers as rank queries, including above `2^53` total weight. Mutable and frozen single-centroid queries also preserve stored extrema and interpolate their tails. +* Fix `TDigest` and `TDigestMut` queries that could return non-finite, out-of-range, or non-monotonic results for extreme finite values or large total weights. Tail queries now use consistent cumulative weights, including above `2^53` total weight and after updating or merging deserialized digests, avoiding decreasing ranks and negative PMF buckets. Mutable and frozen single-centroid queries also preserve stored extrema and interpolate their tails. * T-Digest compression reduces rounding loss in centroid means when merging very unequal weights, including opposite-sign centroids. * `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, inconsistent single-sample extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index dbfa3d2a..ef39ffdb 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -1659,31 +1659,35 @@ impl TDigestView<'_> { let centroids_weight = self.centroids_weight as f64; let num_centroids = self.centroids.len(); - // left tail - let first_mean = self.centroids[0].mean; - if value < first_mean { + // Reserve at most one sample for each stored extremum, without passing the adjacent + // centroid center. Updates can place a unit-weight centroid inside the stored extrema, + // leaving no mass to interpolate in that tail. + let first = &self.centroids[0]; + if value < first.mean { + let center = centroid_center(0, first.weight.get()); return Some(if value == self.min { 0.5 / centroids_weight } else { interpolate( - 1., - centroid_center(0, self.centroids[0].weight.get()), - interpolation_fraction(value, self.min, first_mean), + center.min(1.), + center, + interpolation_fraction(value, self.min, first.mean), ) / centroids_weight }); } - // right tail - let last_mean = self.centroids[num_centroids - 1].mean; - if value > last_mean { + let last = &self.centroids[num_centroids - 1]; + if value > last.mean { + let center = + centroid_center(self.centroids_weight - last.weight.get(), last.weight.get()); return Some(if value == self.max { - 1. - (0.5 / centroids_weight) + // Round the maximum's mass in the same coordinates as the centroid centers. + centroid_center(self.centroids_weight - 1, 1) / centroids_weight } else { - let last_weight = self.centroids[num_centroids - 1].weight.get(); interpolate( - centroid_center(self.centroids_weight - last_weight, last_weight), - (self.centroids_weight - 1) as f64, - interpolation_fraction(value, last_mean, self.max), + center, + center.max((self.centroids_weight - 1) as f64), + interpolation_fraction(value, last.mean, self.max), ) / centroids_weight }); } diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index 4827764e..0790d415 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -741,6 +741,98 @@ fn test_mutable_rank_preserves_weighted_single_centroid_tails() { } } +#[test] +fn test_rank_tails_stay_monotonic_after_updating_weighted_centroids() { + let points = [ + 0., + 0_f64.next_up(), + 0.25, + 0.5, + 1., + 50., + 99., + 99.5, + 99.75, + 100., + ]; + for operation in 0..3 { + let mut digest = deserialize_with_centroids(100, 0., 100., &[(50., 10)]); + let mut updates = TDigestMut::new(100).unwrap(); + updates.update(1.); + updates.update(99.); + match operation { + 0 => { + digest.update(1.); + digest.update(99.); + } + 1 => digest.merge(&updates), + _ => digest = [digest, updates].into_iter().collect(), + } + let frozen = digest.clone().freeze(); + + // The new unit-weight tails lie inside the retained extrema. Their half-sample + // centers leave no mass for interpolation between min/first or last/max. + for (&point, expected) in points.iter().zip([ + 0.5 / 12., + 0.5 / 12., + 0.5 / 12., + 0.5 / 12., + 0.5 / 12., + 0.5, + 11.5 / 12., + 11.5 / 12., + 11.5 / 12., + 11.5 / 12., + ]) { + assert_eq!(digest.rank(point), Some(expected)); + assert_eq!(frozen.rank(point), Some(expected)); + } + let cdf = digest.cdf(&points).unwrap(); + assert_eq!(cdf, frozen.cdf(&points).unwrap()); + assert!(cdf.is_sorted(), "CDF at {points:?}: {cdf:?}"); + let pmf = digest.pmf(&points).unwrap(); + assert_eq!(pmf, frozen.pmf(&points).unwrap()); + assert!(pmf.iter().all(|&mass| mass >= 0.), "PMF: {pmf:?}"); + assert_that!(pmf.iter().sum::(), near(1., 1e-12)); + } +} + +#[test] +fn test_rank_tail_masses_stay_ordered_when_counts_round() { + for total_weight in [ + (1_u64 << 52) - 2, + (1_u64 << 53) - 2, + (1_u64 << 53) - 1, + 1_u64 << 53, + (1_u64 << 53) + 1, + (1_u64 << 54) - 2, + u64::MAX, + ] { + for last_weight in 1..=4 { + let mut digest = deserialize_with_centroids( + 100, + 0., + 100., + &[(0., total_weight - last_weight), (90., last_weight)], + ); + let frozen = digest.clone().freeze(); + let points = [90., 90_f64.next_up(), 95., 100_f64.next_down(), 100.]; + let cdf = digest.cdf(&points).unwrap(); + assert_eq!(cdf, frozen.cdf(&points).unwrap()); + assert!( + cdf.is_sorted(), + "N={total_weight}, last={last_weight}: {cdf:?}" + ); + assert!(cdf.iter().all(|rank| (0.0..=1.0).contains(rank))); + + // At N = 2^53 - 2, the unit-weight last center rounds up to N. Computing + // rank(max) as 1 - 0.5/N instead would round down and produce a negative bucket. + let pmf = frozen.pmf(&points).unwrap(); + assert!(pmf.iter().all(|&mass| mass >= 0.), "PMF: {pmf:?}"); + } + } +} + #[test] fn test_small_compression_preserves_weighted_centroids_and_tie_order() { // Exercise the no-merge bound at two k values and a tiny input whose K_2 normalizer From 618279830f5cda812cafd5867d0ee78828139cad Mon Sep 17 00:00:00 2001 From: tison Date: Sun, 4 Oct 2026 22:48:34 +0800 Subject: [PATCH 14/14] test(tdigest): consolidate numerical regression fixtures 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. --- CHANGELOG.md | 4 +- .../tests/tdigest_test/sketch.rs | 262 +++++------------- 2 files changed, 78 insertions(+), 188 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 037284dc..995aa648 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -46,9 +46,9 @@ All significant changes to this project will be documented in this file. * `ReqSketch::cdf` and its sorted view now preserve the rank at each split point and end at exactly `1.0`, avoiding extra rounding from accumulating PMF values. * `TDigestMut` updates and merges now panic without changing the digest if the total weight would exceed `u64::MAX`, including in release builds. * `TDigestMut::merge` now preserves the true minimum and maximum from both inputs, including compressed digests. -* Fix `TDigest` and `TDigestMut` queries that could return non-finite, out-of-range, or non-monotonic results for extreme finite values or large total weights. Tail queries now use consistent cumulative weights, including above `2^53` total weight and after updating or merging deserialized digests, avoiding decreasing ranks and negative PMF buckets. Mutable and frozen single-centroid queries also preserve stored extrema and interpolate their tails. +* Fix non-finite, out-of-range, or decreasing `TDigest` and `TDigestMut` query results for extreme finite values, large weights, and updated or merged deserialized digests. Single-centroid queries now honor stored extrema, and PMF buckets remain nonnegative. * T-Digest compression reduces rounding loss in centroid means when merging very unequal weights, including opposite-sign centroids. -* `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for invalid flags or extrema, inconsistent single-sample extrema, out-of-range values, unsorted centroids, or non-empty images with no stored values. +* `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for malformed images, including unsorted centroids, values outside stored extrema, and single-sample images with distinct extrema. ## v0.5.0 (2026-09-04) diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index 0790d415..428d96c7 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -623,11 +623,7 @@ fn test_estimate_repeat_values() { assert_eq!(tdigest.quantile(0.9), Some(1.0)); } -/// Builds a digest whose centroids carry the given weights. -/// -/// Compression never merges the extreme centroids, so digests built through `update` and `merge` -/// always keep unit-weight tails. Heavier tails arrive only through deserialization, including the -/// reference implementation format, and they select the tail interpolation branches. +/// Builds a weighted image, including tail centroids unavailable through raw updates. fn deserialize_with_centroids(k: u16, min: f64, max: f64, centroids: &[(f64, u64)]) -> TDigestMut { const PREAMBLE_LONGS: u8 = 2; const SERIAL_VERSION: u8 = 1; @@ -703,7 +699,7 @@ fn test_quantile_handles_two_sample_last_centroid() { #[test] fn test_single_centroid_preserves_stored_tail_information() { - let digest = deserialize_with_centroids(100, 0., 100., &[(50., 10)]); + let mut digest = deserialize_with_centroids(100, 0., 100., &[(50., 10)]); // Ten samples have a centroid center at weight 5. The tail spans weights 1 through 5, // so weight 2.5 lies 3/8 of the way from the stored minimum to the mean. assert_quantile_queries( @@ -711,27 +707,11 @@ fn test_single_centroid_preserves_stored_tail_information() { &[0., 0.1, 0.25, 0.5, 0.75, 0.9, 1.], &[0., 0., 18.75, 50., 81.25, 100., 100.], ); - let digest = digest.freeze(); - for (value, rank) in [ - (0., 0.05), - (18.75, 0.25), - (50., 0.5), - (81.25, 0.75), - (100., 0.95), - ] { - assert_that!(digest.rank(value).unwrap(), near(rank, 1e-12)); - } -} - -#[test] -fn test_mutable_rank_preserves_weighted_single_centroid_tails() { - let mut digest = deserialize_with_centroids(100, 0., 100., &[(50., 10)]); let frozen = digest.clone().freeze(); let points = [-1., 0., 18.75, 50., 81.25, 100., 101.]; let expected = [0., 0.05, 0.25, 0.5, 0.75, 0.95, 1.]; - // A quarter rank is mass 2.5: 3/8 of the tail from mass 1 at min to mass 5 at mean. - // Repeating the queries must preserve the result before and after creating a view. + // Repeated queries agree before and after creating a mutable view. for _ in 0..2 { for (&point, &rank) in points.iter().zip(&expected) { assert_eq!(digest.rank(point), Some(rank)); @@ -743,18 +723,8 @@ fn test_mutable_rank_preserves_weighted_single_centroid_tails() { #[test] fn test_rank_tails_stay_monotonic_after_updating_weighted_centroids() { - let points = [ - 0., - 0_f64.next_up(), - 0.25, - 0.5, - 1., - 50., - 99., - 99.5, - 99.75, - 100., - ]; + let points = [0., 0_f64.next_up(), 0.5, 1., 99., 99.5, 100.]; + let expected = [0.5, 0.5, 0.5, 0.5, 11.5, 11.5, 11.5].map(|mass| mass / 12.); for operation in 0..3 { let mut digest = deserialize_with_centroids(100, 0., 100., &[(50., 10)]); let mut updates = TDigestMut::new(100).unwrap(); @@ -772,18 +742,7 @@ fn test_rank_tails_stay_monotonic_after_updating_weighted_centroids() { // The new unit-weight tails lie inside the retained extrema. Their half-sample // centers leave no mass for interpolation between min/first or last/max. - for (&point, expected) in points.iter().zip([ - 0.5 / 12., - 0.5 / 12., - 0.5 / 12., - 0.5 / 12., - 0.5 / 12., - 0.5, - 11.5 / 12., - 11.5 / 12., - 11.5 / 12., - 11.5 / 12., - ]) { + for (&point, &expected) in points.iter().zip(&expected) { assert_eq!(digest.rank(point), Some(expected)); assert_eq!(frozen.rank(point), Some(expected)); } @@ -797,42 +756,6 @@ fn test_rank_tails_stay_monotonic_after_updating_weighted_centroids() { } } -#[test] -fn test_rank_tail_masses_stay_ordered_when_counts_round() { - for total_weight in [ - (1_u64 << 52) - 2, - (1_u64 << 53) - 2, - (1_u64 << 53) - 1, - 1_u64 << 53, - (1_u64 << 53) + 1, - (1_u64 << 54) - 2, - u64::MAX, - ] { - for last_weight in 1..=4 { - let mut digest = deserialize_with_centroids( - 100, - 0., - 100., - &[(0., total_weight - last_weight), (90., last_weight)], - ); - let frozen = digest.clone().freeze(); - let points = [90., 90_f64.next_up(), 95., 100_f64.next_down(), 100.]; - let cdf = digest.cdf(&points).unwrap(); - assert_eq!(cdf, frozen.cdf(&points).unwrap()); - assert!( - cdf.is_sorted(), - "N={total_weight}, last={last_weight}: {cdf:?}" - ); - assert!(cdf.iter().all(|rank| (0.0..=1.0).contains(rank))); - - // At N = 2^53 - 2, the unit-weight last center rounds up to N. Computing - // rank(max) as 1 - 0.5/N instead would round down and produce a negative bucket. - let pmf = frozen.pmf(&points).unwrap(); - assert!(pmf.iter().all(|&mass| mass >= 0.), "PMF: {pmf:?}"); - } - } -} - #[test] fn test_small_compression_preserves_weighted_centroids_and_tie_order() { // Exercise the no-merge bound at two k values and a tiny input whose K_2 normalizer @@ -876,30 +799,45 @@ fn test_small_compression_preserves_weighted_centroids_and_tie_order() { } } +// Exercise both compression paths and inspect the mean directly; query interpolation can hide +// rounding errors. The heavy outer centroids keep the middle pair away from the protected tails. +fn merged_centroid_mean(centroids: [(f64, u64); 4], reverse: bool, owned: bool) -> f64 { + let min = centroids[0].0; + let max = centroids[3].0; + let mut bytes = deserialize_with_centroids(10, min, max, ¢roids).serialize(); + if reverse { + bytes[5] |= 1 << 2; + } + let mut digest = TDigestMut::deserialize(&bytes).unwrap(); + if owned { + let other = deserialize_with_centroids(10, max, max, &[(max, 1)]); + digest = [digest, other].into_iter().collect(); + } else { + digest.update(min); + } + let bytes = digest.serialize(); + let merged_weight = centroids[1].1 + centroids[2].1; + let centroid = bytes[32..] + .chunks_exact(16) + .find(|c| u64::from_le_bytes(c[8..].try_into().unwrap()) == merged_weight) + .unwrap_or_else(|| { + panic!( + "middle pair {:?} did not merge (reverse={reverse}, owned={owned})", + ¢roids[1..3] + ) + }); + f64::from_le_bytes(centroid[..8].try_into().unwrap()) +} + #[test] fn test_compression_preserves_a_small_centroid_weight() { for weight in [1_u64 << 20, 1_u64 << 40, 1_u64 << 54] { for reverse in [false, true] { - let mut digest = deserialize_with_centroids( - 10, - 0., - f64::MAX, - &[(0., weight), (1., weight), (1e300, 1), (f64::MAX, weight)], + let mean = merged_centroid_mean( + [(0., weight), (1., weight), (1e300, 1), (f64::MAX, weight)], + reverse, + false, ); - let mut bytes = digest.serialize(); - if reverse { - bytes[5] |= 1 << 2; - } - let mut digest = TDigestMut::deserialize(&bytes).unwrap(); - digest.update(0.); - - // Read the merged mean directly: quantile interpolation must not hide a bad mean. - let bytes = digest.serialize(); - let centroid = bytes[32..] - .chunks_exact(16) - .find(|c| u64::from_le_bytes(c[8..].try_into().unwrap()) == weight + 1) - .expect("the middle centroids should merge"); - let mean = f64::from_le_bytes(centroid[..8].try_into().unwrap()); // The exact mean is (weight * 1 + 1e300) / (weight + 1). The first term is far // below one ULP of the numerator, so this division is an independent reference. let expected = 1e300 / (weight + 1) as f64; @@ -929,41 +867,16 @@ fn test_opposite_sign_merge_preserves_a_representable_correction() { }; for reverse in [false, true] { for owned in [false, true] { - let mut digest = deserialize_with_centroids( - 10, - -f64::MAX, - f64::MAX, - &[ + let mean = merged_centroid_mean( + [ (-f64::MAX, outer_weight), (-magnitude, left_weight), (magnitude, right_weight), (f64::MAX, outer_weight), ], + reverse, + owned, ); - let mut bytes = digest.serialize(); - if reverse { - bytes[5] |= 1 << 2; - } - let mut digest = TDigestMut::deserialize(&bytes).unwrap(); - if owned { - let other = deserialize_with_centroids( - 10, - f64::MAX, - f64::MAX, - &[(f64::MAX, 1)], - ); - digest = [digest, other].into_iter().collect(); - } else { - digest.update(-f64::MAX); - } - let bytes = digest.serialize(); - let centroid = bytes[32..] - .chunks_exact(16) - .find(|c| u64::from_le_bytes(c[8..].try_into().unwrap()) == total_weight) - .unwrap_or_else(|| { - panic!("middle pair did not merge: magnitude={magnitude:?}, exponent={exponent}, negative_heavy={negative_heavy}, reverse={reverse}, owned={owned}") - }); - let mean = f64::from_le_bytes(centroid[..8].try_into().unwrap()); assert_eq!( mean, expected, "magnitude={magnitude:?}, exponent={exponent}, negative_heavy={negative_heavy}, reverse={reverse}, owned={owned}" @@ -982,12 +895,9 @@ fn test_compression_keeps_merged_means_within_their_endpoints() { for (left, right) in [ (0., tiny), (tiny, 2. * tiny), - ( - f64::from_bits(f64::MIN_POSITIVE.to_bits() - 1), - f64::MIN_POSITIVE, - ), - (1., f64::from_bits(1_f64.to_bits() + 1)), - (f64::from_bits(f64::MAX.to_bits() - 1), f64::MAX), + (f64::MIN_POSITIVE.next_down(), f64::MIN_POSITIVE), + (1., 1_f64.next_up()), + (f64::MAX.next_down(), f64::MAX), (f64::MAX / 4., f64::MAX), (-f64::MAX, -f64::MAX / 4.), (-2. * tiny, -tiny), @@ -997,44 +907,16 @@ fn test_compression_keeps_merged_means_within_their_endpoints() { for (left_weight, right_weight) in [(1, 1), (1, heavy), (heavy, 1)] { for reverse in [false, true] { for owned in [false, true] { - let mut digest = deserialize_with_centroids( - 10, - -f64::MAX, - f64::MAX, - &[ + let mean = merged_centroid_mean( + [ (-f64::MAX, heavy), (left, left_weight), (right, right_weight), (f64::MAX, heavy), ], + reverse, + owned, ); - let mut bytes = digest.serialize(); - if reverse { - bytes[5] |= 1 << 2; - } - let mut digest = TDigestMut::deserialize(&bytes).unwrap(); - if owned { - let other = - deserialize_with_centroids(10, f64::MAX, f64::MAX, &[(f64::MAX, 1)]); - digest = [digest, other].into_iter().collect(); - } else { - digest.update(-f64::MAX); - } - - // Inspect compression directly; query interpolation could hide an overshoot. - let bytes = digest.serialize(); - let centroid = bytes[32..] - .chunks_exact(16) - .find(|c| { - u64::from_le_bytes(c[8..].try_into().unwrap()) - == left_weight + right_weight - }) - .unwrap_or_else(|| { - panic!( - "middle pair {left:?}:{left_weight}, {right:?}:{right_weight} did not merge (reverse={reverse}, owned={owned})" - ) - }); - let mean = f64::from_le_bytes(centroid[..8].try_into().unwrap()); assert!( mean.is_finite() && (left..=right).contains(&mean), "{mean:?} is outside [{left:?}, {right:?}] (weights={left_weight}:{right_weight}, reverse={reverse}, owned={owned})" @@ -1050,30 +932,18 @@ fn test_signed_zero_and_symmetric_centroids_survive_interpolation() { for magnitude in [0., f64::from_bits(1), 1., f64::MAX] { let heavy = 1_u64 << 54; for reverse in [false, true] { - let mut digest = deserialize_with_centroids( - 10, - -f64::MAX, - f64::MAX, - &[ + let mean = merged_centroid_mean( + [ (-f64::MAX, heavy), (-magnitude, 1), (magnitude, 1), (f64::MAX, heavy), ], + reverse, + false, ); - let mut bytes = digest.serialize(); - if reverse { - bytes[5] |= 1 << 2; - } - let mut digest = TDigestMut::deserialize(&bytes).unwrap(); - digest.update(-f64::MAX); - let bytes = digest.serialize(); - let merged = bytes[32..] - .chunks_exact(16) - .find(|c| u64::from_le_bytes(c[8..].try_into().unwrap()) == 2) - .expect("the equally weighted middle centroids should merge"); // Equal and opposite contributions cancel exactly, including subnormal values. - assert_eq!(f64::from_le_bytes(merged[..8].try_into().unwrap()), 0.); + assert_eq!(mean, 0.); } } @@ -1301,6 +1171,26 @@ fn test_quantile_uses_max_when_last_center_rounds_down() { ); } +#[test] +fn test_rank_max_uses_the_same_rounded_mass_as_last_centroid() { + // The last half-sample center can round up or down. At N = 2^53 - 2 it rounds to N, + // while 1 - 0.5/N rounds below 1 and would give max a lower rank than the last mean. + for (total_weight, expected) in [ + ((1_u64 << 53) - 2, 1.), + ((1_u64 << 53) - 1, 1_f64.next_down()), + (u64::MAX, 1.), + ] { + let mut digest = + deserialize_with_centroids(100, 0., 100., &[(0., total_weight - 1), (90., 1)]); + let frozen = digest.clone().freeze(); + let points = [90., 90_f64.next_up(), 95., 100_f64.next_down(), 100.]; + let cdf = digest.cdf(&points).unwrap(); + assert_eq!(&cdf[..points.len()], &[expected; 5]); + assert_eq!(cdf, frozen.cdf(&points).unwrap()); + assert!(frozen.pmf(&points).unwrap().iter().all(|&mass| mass >= 0.)); + } +} + #[test] fn test_rank_stays_monotonic_at_the_right_tail_boundary() { let digest = deserialize_with_centroids(100, 0., 300., &[(0., 1), (100., 29)]).freeze();