quantile_sketch/

lib.rs

#![allow(rustdoc::bare_urls)]
#![doc = include_str!("../README.md")]
#![cfg_attr(not(any(feature = "std", test)), no_std)]

extern crate alloc;

use alloc::{boxed::Box, vec::Vec};
use core::{fmt, iter::Peekable, ptr};

#[cfg(feature = "loom")]
use loom::sync::atomic::{AtomicPtr, AtomicU64, Ordering};
#[cfg(not(feature = "loom"))]
use portable_atomic::{AtomicPtr, AtomicU64, Ordering};

#[cfg(all(test, not(feature = "loom")))]
mod concurrency_tests;
#[cfg(all(feature = "loom", test))]
mod loom_tests;
mod math;
#[cfg(all(test, not(feature = "loom")))]
mod property_tests;
#[cfg(feature = "serde")]
mod serde_impl;

use math::CubicMapping;

const BLOCK_SIZE: usize = 64;
const DEFAULT_ERROR: f64 = 0.01;

/// Error returned when sketches with different configurations are merged.
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct MergeError;

impl fmt::Display for MergeError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.write_str("cannot merge sketches with different configurations")
    }
}

#[cfg(feature = "std")]
impl std::error::Error for MergeError {}

/// Lazily allocated block
struct Block {
    counts: [AtomicU64; BLOCK_SIZE],
}

impl Block {
    #[inline]
    fn total(&self) -> u64 {
        self.counts.iter().map(|count| count.load(Ordering::Relaxed)).sum()
    }
}

impl Default for Block {
    fn default() -> Self {
        Self {
            counts: core::array::from_fn(|_| AtomicU64::new(0)),
        }
    }
}

#[inline]
fn bucket_at_rank(
    buckets: &mut Peekable<impl Iterator<Item = ((usize, usize), u64)>>,
    mut rank: u64,
) -> Option<(usize, usize)> {
    loop {
        let (coord, remaining) = buckets.peek_mut()?;
        if rank < *remaining {
            // Keep the selected sample as the origin for the next relative rank.
            *remaining -= rank;
            return Some(*coord);
        }
        rank -= *remaining;
        buckets.next();
    }
}

#[inline]
fn rank(q: f64, last_rank: u64) -> u64 {
    ((q * last_rank as f64) as u64).min(last_rank)
}

fn required_bucket_count(min_index: i64, max_index: i64) -> usize {
    max_index
        .checked_sub(min_index)
        .and_then(|count| count.checked_add(1))
        .and_then(|count| usize::try_from(count).ok())
        .expect("range requires too many buckets for this target")
}

/// DDSketch with atomic concurrent inserts.
pub struct ConcurrentDDSketch {
    blocks: Box<[AtomicPtr<Block>]>,
    mapping: CubicMapping,
    min: f64,
    max: f64,
    min_index: i64,
}

impl ConcurrentDDSketch {
    /// Creates a sketch with 1% relative error for values from zero through
    /// `f64::MAX`.
    pub fn new() -> Self {
        Self::with_err_and_range(DEFAULT_ERROR, 0.0, f64::MAX)
    }

    /// Creates a sketch with the given relative error for the full positive
    /// range.
    ///
    /// # Panics
    ///
    /// Panics if `err` is not strictly between zero and one or is too small to
    /// produce distinct buckets.
    pub fn with_err(err: f64) -> Self {
        Self::with_err_and_range(err, 0.0, f64::MAX)
    }

    /// Creates a sketch with 1% relative error for the inclusive range
    /// `min..=max`.
    ///
    /// # Panics
    ///
    /// Panics if the range is invalid or contains negative values.
    pub fn with_range(min: f64, max: f64) -> Self {
        Self::with_err_and_range(DEFAULT_ERROR, min, max)
    }

    /// Creates a sketch with the given relative error and inclusive range.
    ///
    /// # Panics
    ///
    /// Panics if `err` or the range is invalid.
    pub fn with_err_and_range(err: f64, min: f64, max: f64) -> Self {
        assert!(
            min == 0.0 || min >= f64::MIN_POSITIVE,
            "min must be zero or normal and positive"
        );
        assert!(min <= max, "min must not exceed max");
        assert!(max < f64::INFINITY, "max must be finite");

        let mapping = CubicMapping::new(err);
        let min_index = if min == 0.0 {
            mapping.index(f64::MIN_POSITIVE) - 1
        } else {
            mapping.index(min)
        };
        let max_index = if max < f64::MIN_POSITIVE {
            min_index
        } else {
            mapping.index(max)
        };
        let bucket_count = required_bucket_count(min_index, max_index);
        let block_count = bucket_count / BLOCK_SIZE + usize::from(bucket_count % BLOCK_SIZE != 0);
        let blocks = (0..block_count)
            .map(|_| AtomicPtr::new(ptr::null_mut()))
            .collect::<Vec<_>>();
        Self {
            blocks: blocks.into(),
            min_index,
            mapping,
            min,
            max,
        }
    }

    #[inline]
    fn coord(&self, val: f64) -> (usize, usize) {
        assert!(val >= self.min && val <= self.max, "value out of range");
        let index = if val < f64::MIN_POSITIVE {
            self.min_index
        } else {
            self.mapping.index(val)
        };
        let offset = (index - self.min_index) as usize;
        (offset / BLOCK_SIZE, offset % BLOCK_SIZE)
    }

    #[inline]
    fn allocated_blocks(&self) -> impl DoubleEndedIterator<Item = (usize, &Block)> {
        self.blocks.iter().enumerate().filter_map(|(index, slot)| {
            let block = slot.load(Ordering::Acquire);
            if block.is_null() {
                None
            } else {
                // SAFETY: Every non-null pointer was created by
                // `Box::into_raw` in `init_block`. Published pointers are never
                // replaced and remain allocated until `self` is dropped.
                Some((index, unsafe { &*block }))
            }
        })
    }

    #[inline]
    fn buckets(&self) -> impl DoubleEndedIterator<Item = ((usize, usize), u64)> + '_ {
        self.allocated_blocks().flat_map(|(block_index, block)| {
            block
                .counts
                .iter()
                .enumerate()
                .map(move |(count_index, count)| ((block_index, count_index), count.load(Ordering::Relaxed)))
        })
    }

    #[inline]
    fn last_rank(&self) -> Option<u64> {
        self.allocated_blocks()
            .map(|(_, block)| block.total())
            .sum::<u64>()
            .checked_sub(1)
    }

    #[inline]
    fn value_at_coord(&self, block_index: usize, count_index: usize) -> f64 {
        let index = self.min_index + (block_index * BLOCK_SIZE + count_index) as i64;
        if self.min == 0.0 && index == self.min_index {
            0.0
        } else {
            self.mapping.value(index).clamp(self.min, self.max)
        }
    }

    fn write_ranked_quantiles(
        &self,
        ranks: impl Iterator<Item = (usize, u64)>,
        buckets: impl Iterator<Item = ((usize, usize), u64)>,
        output: &mut [Option<f64>],
    ) {
        let mut buckets = buckets.peekable();
        let mut prev_rank = 0;
        for (output_index, rank) in ranks {
            let (block_index, count_index) =
                bucket_at_rank(&mut buckets, rank - prev_rank).expect("bucket totals must match bucket counts");
            output[output_index] = Some(self.value_at_coord(block_index, count_index));
            prev_rank = rank;
        }
    }

    fn write_quantiles(&self, quantiles: &[f64], last_rank: u64, reverse: bool, output: &mut [Option<f64>]) {
        if reverse {
            let r = quantiles.iter().map(|&q| last_rank - rank(q, last_rank)).enumerate();
            self.write_ranked_quantiles(r.rev(), self.buckets().rev(), output);
        } else {
            let r = quantiles.iter().map(|&q| rank(q, last_rank)).enumerate();
            self.write_ranked_quantiles(r, self.buckets(), output);
        }
    }

    #[cold]
    #[inline(never)]
    fn init_block(slot: &AtomicPtr<Block>) -> *mut Block {
        let new = Box::into_raw(Box::<Block>::default());
        match slot.compare_exchange(ptr::null_mut(), new, Ordering::Release, Ordering::Acquire) {
            Ok(_) => new,
            Err(existing) => {
                // Another thread initialized this block first.
                // SAFETY: `new` came from `Box::into_raw` above and was not
                // published because the compare-exchange failed.
                unsafe { drop(Box::from_raw(new)) };
                existing
            }
        }
    }

    #[inline]
    fn get_or_init_block(slot: &AtomicPtr<Block>) -> &Block {
        let block = slot.load(Ordering::Acquire);
        let block = if block.is_null() { Self::init_block(slot) } else { block };
        // SAFETY: `block` is non-null and points to a published `Block` that is
        // not reclaimed until the sketch is dropped. The returned reference is
        // tied to the borrow of the slot within that sketch.
        unsafe { &*block }
    }

    /// Adds a value to the sketch.
    ///
    /// # Panics
    ///
    /// Panics if `val` is NaN or outside the configured range.
    #[inline]
    pub fn insert(&self, val: f64) {
        let (block_index, count_index) = self.coord(val);
        Self::get_or_init_block(&self.blocks[block_index]).counts[count_index].fetch_add(1, Ordering::Relaxed);
    }

    /// Adds the counts from a sketch with the same error and supported range.
    ///
    /// This is not an atomic snapshot: each source counter is read and added
    /// once, so concurrent inserts may or may not be included and concurrent
    /// queries may observe a partial merge.
    pub fn merge(&self, other: &Self) -> Result<(), MergeError> {
        if self.mapping.error() != other.mapping.error() || self.min != other.min || self.max != other.max {
            return Err(MergeError);
        }

        for (block_index, source) in other.allocated_blocks() {
            let target = Self::get_or_init_block(&self.blocks[block_index]);
            for (target, source) in target.counts.iter().zip(&source.counts) {
                let count = source.load(Ordering::Relaxed);
                if count != 0 {
                    target.fetch_add(count, Ordering::Relaxed);
                }
            }
        }
        Ok(())
    }

    /// Returns the configured inclusive minimum supported value.
    pub fn min(&self) -> f64 {
        self.min
    }

    /// Returns the configured inclusive maximum supported value.
    pub fn max(&self) -> f64 {
        self.max
    }

    /// Returns the lower `q`-quantile, or `None` when the sketch is empty.
    /// Concurrent inserts may or may not be reflected in the result.
    ///
    /// # Panics
    ///
    /// Panics if `q` is NaN or outside `0.0..=1.0`.
    pub fn quantile(&self, q: f64) -> Option<f64> {
        assert!((0.0..=1.0).contains(&q), "quantile must be between 0 and 1");

        let last_rank = self.last_rank()?;
        let rank = rank(q, last_rank);

        let found = if rank <= last_rank - rank {
            bucket_at_rank(&mut self.buckets().peekable(), rank)
        } else {
            bucket_at_rank(&mut self.buckets().rev().peekable(), last_rank - rank)
        };

        let (block_index, count_index) = found.expect("bucket totals must match bucket counts");
        Some(self.value_at_coord(block_index, count_index))
    }

    /// Writes lower estimates for sorted `quantiles` to `output` without
    /// allocating.
    ///
    /// Results are nondecreasing but are not an atomic snapshot of concurrent
    /// inserts. Separate [`Self::quantile`] calls can return p50 > p99 if
    /// inserts occur between them. An empty sketch fills `output` with `None`.
    ///
    /// # Panics
    ///
    /// Panics if the lengths differ or a quantile is unsorted, NaN, or outside
    /// `0.0..=1.0`. Validation failures leave `output` unchanged.
    pub fn quantiles(&self, quantiles: &[f64], output: &mut [Option<f64>]) {
        assert_eq!(output.len(), quantiles.len(), "output length must match quantiles");
        for &q in quantiles {
            assert!((0.0..=1.0).contains(&q), "quantile must be between 0 and 1");
        }
        assert!(
            quantiles.windows(2).all(|pair| pair[0] <= pair[1]),
            "quantiles must be sorted"
        );

        output.fill(None);
        let Some((&first, &last)) = quantiles.first().zip(quantiles.last()) else {
            return;
        };
        let Some(last_rank) = self.last_rank() else {
            return;
        };

        let reverse = rank(last, last_rank) > last_rank - rank(first, last_rank);
        self.write_quantiles(quantiles, last_rank, reverse, output);
    }

    /// Returns the estimated fraction of samples at or below `val`.
    ///
    /// The result is in `0.0..=1.0`, or `None` when the sketch is empty. Since
    /// values in the same bucket are indistinguishable, this includes the
    /// entire bucket containing `val`. Concurrent inserts may or may not be
    /// reflected in the result.
    ///
    /// # Panics
    ///
    /// Panics if `val` is NaN or outside the configured range.
    pub fn percentile_rank(&self, val: f64) -> Option<f64> {
        let (target_block, target_count) = self.coord(val);
        let mut cumulative = 0_u64;
        let mut total = 0_u64;

        for (block_index, block) in self.allocated_blocks() {
            if block_index == target_block {
                for (count_index, count) in block.counts.iter().enumerate() {
                    let count = count.load(Ordering::Relaxed);
                    total += count;
                    if count_index <= target_count {
                        cumulative += count;
                    }
                }
                continue;
            }

            let count = block.total();
            total += count;
            if block_index < target_block {
                cumulative += count;
            }
        }

        (total != 0).then_some(cumulative as f64 / total as f64)
    }
}

impl Default for ConcurrentDDSketch {
    fn default() -> Self {
        Self::new()
    }
}

impl fmt::Debug for ConcurrentDDSketch {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.debug_struct("ConcurrentDDSketch")
            .field("error", &self.mapping.error())
            .field("min", &self.min)
            .field("max", &self.max)
            .finish_non_exhaustive()
    }
}

impl Drop for ConcurrentDDSketch {
    fn drop(&mut self) {
        for slot in self.blocks.iter_mut() {
            let block = slot.load(Ordering::Relaxed);
            if !block.is_null() {
                // SAFETY: Drop has exclusive access to the sketch. Each
                // published pointer came from `Box::into_raw`, is unique to its
                // slot, and has not previously been reclaimed.
                unsafe { drop(Box::from_raw(block)) };
            }
        }
    }
}

#[cfg(test)]
fn assert_estimate_eq(actual: Option<f64>, expected: Option<f64>) {
    #[cfg(not(miri))]
    assert_eq!(actual, expected);

    #[cfg(miri)]
    match (actual, expected) {
        (Some(actual), Some(expected)) => {
            // Miri injects permitted error into imprecise float operations, so
            // equivalent estimates are not always bit-for-bit equal.
            let tolerance = 8_192.0 * f64::EPSILON * expected.abs().max(1.0);
            assert!((actual - expected).abs() <= tolerance, "{actual} != {expected}");
        }
        _ => assert_eq!(actual, expected),
    }
}

#[cfg(all(test, not(feature = "loom")))]
mod tests {
    use super::*;

    fn assert_panics(f: impl FnOnce()) {
        assert!(std::panic::catch_unwind(std::panic::AssertUnwindSafe(f)).is_err());
    }

    fn nonzero_bins(sketch: &ConcurrentDDSketch) -> Vec<(usize, u64)> {
        sketch
            .allocated_blocks()
            .flat_map(|(block_index, block)| {
                block.counts.iter().enumerate().filter_map(move |(count_index, count)| {
                    let count = count.load(Ordering::Relaxed);
                    (count != 0).then_some((block_index * BLOCK_SIZE + count_index, count))
                })
            })
            .collect()
    }

    #[cfg(feature = "serde")]
    #[derive(Debug, PartialEq)]
    struct SketchState {
        min: f64,
        max: f64,
        quantiles: [Option<f64>; 3],
    }

    #[cfg(feature = "serde")]
    impl<'de> serde::Deserialize<'de> for SketchState {
        fn deserialize<D: serde::Deserializer<'de>>(deserializer: D) -> Result<Self, D::Error> {
            let sketch = ConcurrentDDSketch::deserialize(deserializer)?;
            Ok(Self {
                min: sketch.min(),
                max: sketch.max(),
                quantiles: [sketch.quantile(0.0), sketch.quantile(0.5), sketch.quantile(1.0)],
            })
        }
    }

    #[cfg(feature = "serde")]
    fn serde_tokens(max_offset: u64) -> Vec<serde_test::Token> {
        use serde_test::Token;

        vec![
            Token::Struct { name: "Repr", len: 5 },
            Token::Str("version"),
            Token::U8(1),
            Token::Str("error"),
            Token::F64(0.02),
            Token::Str("min"),
            Token::F64(0.0),
            Token::Str("max"),
            Token::F64(10.0),
            Token::Str("bins"),
            Token::Seq { len: Some(2) },
            Token::Tuple { len: 2 },
            Token::U64(0),
            Token::U64(2),
            Token::TupleEnd,
            Token::Tuple { len: 2 },
            Token::U64(max_offset),
            Token::U64(1),
            Token::TupleEnd,
            Token::SeqEnd,
            Token::StructEnd,
        ]
    }

    #[test]
    fn computes_number_of_blocks() {
        for (min, max, expected) in [(1.0, 1.0, 1), (0.001, 1_000.0, 11), (1.0e-9, 1.0e9, 33)] {
            assert_eq!(ConcurrentDDSketch::with_range(min, max).blocks.len(), expected);
        }
    }

    #[test]
    fn constructors_set_expected_ranges() {
        let default = ConcurrentDDSketch::new();
        assert_eq!((default.min, default.max), (0.0, f64::MAX));
        assert_eq!(default.mapping.index(default.max) - default.min_index + 1, 71_609);
        assert_eq!(default.blocks.len(), 1_119);
        assert!(default.blocks.iter().all(|slot| slot.load(Ordering::Relaxed).is_null()));

        for (sketch, error, range) in [
            (ConcurrentDDSketch::with_err(0.02), 0.02, (0.0, f64::MAX)),
            (ConcurrentDDSketch::with_range(1.0, 10.0), DEFAULT_ERROR, (1.0, 10.0)),
            (
                ConcurrentDDSketch::with_err_and_range(0.02, 1.0, 10.0),
                0.02,
                (1.0, 10.0),
            ),
        ] {
            assert_eq!(sketch.mapping.index(10.0), CubicMapping::new(error).index(10.0));
            assert_eq!((sketch.min, sketch.max), range);
        }
    }

    #[test]
    fn rejects_invalid_configurations() {
        for (error, min, max) in [
            (0.0, 0.0, 1.0),
            (1.0, 0.0, 1.0),
            (f64::NAN, 0.0, 1.0),
            (f64::EPSILON / 4.0, 1.0, 1.0),
            (DEFAULT_ERROR, -1.0, 1.0),
            (DEFAULT_ERROR, f64::from_bits(1), 1.0),
            (DEFAULT_ERROR, 2.0, 1.0),
            (DEFAULT_ERROR, 0.0, f64::INFINITY),
            (DEFAULT_ERROR, 0.0, f64::NAN),
        ] {
            assert_panics(|| drop(ConcurrentDDSketch::with_err_and_range(error, min, max)));
        }
    }

    #[test]
    #[should_panic(expected = "range requires too many buckets for this target")]
    fn rejects_bucket_counts_that_cannot_be_indexed() {
        required_bucket_count(i64::MIN, i64::MAX);
    }

    #[test]
    fn handles_single_value_and_extreme_ranges() {
        let zero = ConcurrentDDSketch::with_range(0.0, 0.0);
        zero.insert(0.0);
        assert_eq!(
            [zero.quantile(0.0), zero.quantile(0.5), zero.quantile(1.0)],
            [Some(0.0); 3]
        );

        let single = ConcurrentDDSketch::with_range(1.0, 1.0);
        single.insert(1.0);
        assert_eq!(
            [single.quantile(0.0), single.quantile(0.5), single.quantile(1.0)],
            [Some(1.0); 3]
        );

        let full = ConcurrentDDSketch::new();
        for value in [0.0, f64::MIN_POSITIVE, f64::MAX] {
            full.insert(value);
        }
        assert_eq!(full.quantile(0.0), Some(0.0));
        let estimate = full.quantile(1.0).unwrap();
        assert!((estimate - f64::MAX).abs() / f64::MAX <= DEFAULT_ERROR);
        assert_eq!(full.percentile_rank(f64::MAX), Some(1.0));
    }

    #[test]
    fn zero_and_underflow_have_their_own_bucket() {
        let sketch = ConcurrentDDSketch::new();
        sketch.insert(0.0);
        sketch.insert(-0.0);
        sketch.insert(f64::from_bits(1));
        sketch.insert(f64::MIN_POSITIVE);
        sketch.insert(1.0);
        assert_eq!(sketch.quantile(0.0), Some(0.0));
        assert_eq!(sketch.quantile(0.5), Some(0.0));
        assert!(sketch.quantile(0.75).unwrap() > 0.0);
    }

    #[test]
    fn inserts_into_the_expected_bucket() {
        let sketch = ConcurrentDDSketch::with_err_and_range(0.01, 0.001, 1_000.0);
        let (block_index, count_index) = sketch.coord(1.0);

        sketch.insert(1.0);
        sketch.insert(1.0);

        let block = sketch.blocks[block_index].load(Ordering::Acquire);
        let block = unsafe { &*block };
        assert_eq!(block.counts[count_index].load(Ordering::Relaxed), 2);
    }

    #[test]
    fn computes_percentile_ranks() {
        let sketch = ConcurrentDDSketch::with_range(1.0, 4.0);
        assert_eq!(sketch.percentile_rank(2.0), None);

        for value in [1.0, 2.0, 2.0, 4.0] {
            sketch.insert(value);
        }

        assert_eq!(sketch.percentile_rank(1.0), Some(0.25));
        assert_eq!(sketch.percentile_rank(2.0), Some(0.75));
        assert_eq!(sketch.percentile_rank(3.0), Some(0.75));
        assert_eq!(sketch.percentile_rank(4.0), Some(1.0));
    }

    #[test]
    fn computes_multiple_sorted_quantiles() {
        let sketch = ConcurrentDDSketch::with_range(0.0, 4.0);
        sketch.quantiles(&[], &mut []);
        let mut estimates = [Some(1.0); 3];
        sketch.quantiles(&[0.0, 0.5, 1.0], &mut estimates);
        assert_eq!(estimates, [None; 3]);

        for value in [0.0, 1.0, 2.0, 3.0, 4.0] {
            sketch.insert(value);
        }

        let quantiles = [0.0, 0.5, 0.5, 0.75, 1.0];
        let mut estimates = [None; 5];
        sketch.quantiles(&quantiles, &mut estimates);
        for (&q, estimate) in quantiles.iter().zip(estimates) {
            assert_estimate_eq(estimate, sketch.quantile(q));
        }

        let high_quantiles = [0.75, 0.99, 0.99, 1.0];
        let mut high_estimates = [None; 4];
        sketch.quantiles(&high_quantiles, &mut high_estimates);
        for (&q, estimate) in high_quantiles.iter().zip(high_estimates) {
            assert_estimate_eq(estimate, sketch.quantile(q));
        }

        let mut ordered = [None; 4];
        sketch.quantiles(&[0.0, 0.5, 0.99, 1.0], &mut ordered);
        assert!(ordered.windows(2).all(|pair| pair[0] <= pair[1]));
    }

    #[test]
    fn rejects_out_of_range_inputs() {
        let sketch = ConcurrentDDSketch::with_range(1.0, 10.0);
        for value in [0.0, 11.0, f64::NAN] {
            assert_panics(|| sketch.insert(value));
            assert_panics(|| {
                sketch.percentile_rank(value);
            });
        }
        for quantile in [-f64::EPSILON, 1.0 + f64::EPSILON, f64::NAN] {
            assert_panics(|| {
                sketch.quantile(quantile);
            });
            assert_panics(|| sketch.quantiles(&[0.5, quantile], &mut [None; 2]));
        }
        assert_panics(|| sketch.quantiles(&[0.9, 0.5], &mut [None; 2]));
        assert_panics(|| sketch.quantiles(&[0.5], &mut [None; 2]));
        assert_eq!(sketch.quantile(0.5), None);
    }

    #[test]
    fn reports_configured_min_and_max() {
        let default = ConcurrentDDSketch::new();
        assert_eq!(default.min(), 0.0);
        assert_eq!(default.max(), f64::MAX);

        let bounded = ConcurrentDDSketch::with_range(1.0, 10.0);
        assert_eq!(bounded.min(), 1.0);
        assert_eq!(bounded.max(), 10.0);
    }

    #[test]
    fn merges_compatible_sketches() {
        let merged = ConcurrentDDSketch::with_range(0.0, 1_000.0);
        let empty = ConcurrentDDSketch::with_range(0.0, 1_000.0);
        merged.merge(&empty).unwrap();
        assert!(merged.blocks.iter().all(|slot| slot.load(Ordering::Relaxed).is_null()));

        merged.insert(1.0);
        let other = ConcurrentDDSketch::with_range(0.0, 1_000.0);
        for value in [0.0, f64::from_bits(1), 10.0, 1_000.0] {
            other.insert(value);
        }

        let expected = ConcurrentDDSketch::with_range(0.0, 1_000.0);
        for value in [0.0, f64::from_bits(1), 1.0, 10.0, 1_000.0] {
            expected.insert(value);
        }

        merged.merge(&other).unwrap();
        assert_eq!(nonzero_bins(&merged), nonzero_bins(&expected));
        for q in [0.0, 0.25, 0.5, 0.75, 1.0] {
            assert_estimate_eq(merged.quantile(q), expected.quantile(q));
        }
        assert_eq!(other.percentile_rank(0.0), Some(0.5));

        let self_merge = ConcurrentDDSketch::with_range(1.0, 2.0);
        self_merge.insert(1.0);
        self_merge.insert(2.0);
        self_merge.merge(&self_merge).unwrap();
        assert_eq!(self_merge.percentile_rank(1.0), Some(0.5));
    }

    #[test]
    fn rejects_incompatible_merges() {
        let sketch = ConcurrentDDSketch::with_range(1.0, 10.0);
        for other in [
            ConcurrentDDSketch::with_err_and_range(0.02, 1.0, 10.0),
            ConcurrentDDSketch::with_range(0.0, 10.0),
            ConcurrentDDSketch::with_range(1.0, 20.0),
        ] {
            assert_eq!(sketch.merge(&other), Err(MergeError));
        }
        assert_eq!(sketch.quantile(0.5), None);
    }

    #[test]
    fn debug_reports_configuration() {
        let sketch = ConcurrentDDSketch::with_err_and_range(0.02, 1.0, 10.0);
        assert_eq!(
            format!("{sketch:?}"),
            "ConcurrentDDSketch { error: 0.02, min: 1.0, max: 10.0, .. }"
        );
    }

    #[cfg(feature = "serde")]
    #[test]
    fn serde_preserves_state() {
        use serde_test::{assert_de_tokens, assert_ser_tokens};

        let sketch = ConcurrentDDSketch::with_err_and_range(0.02, 0.0, 10.0);
        for value in [0.0, 0.0, 10.0] {
            sketch.insert(value);
        }
        let max_offset = (sketch.mapping.index(10.0) - sketch.min_index) as u64;
        let tokens = serde_tokens(max_offset);
        let state = SketchState {
            min: sketch.min(),
            max: sketch.max(),
            quantiles: [sketch.quantile(0.0), sketch.quantile(0.5), sketch.quantile(1.0)],
        };

        assert_ser_tokens(&sketch, &tokens);
        assert_de_tokens(&state, &tokens);
    }

    #[cfg(feature = "serde")]
    #[test]
    fn serde_rejects_invalid_state() {
        use serde_test::{assert_de_tokens_error, Token};

        let sketch = ConcurrentDDSketch::with_err_and_range(0.02, 0.0, 10.0);
        let max_offset = (sketch.mapping.index(10.0) - sketch.min_index) as u64;
        let tokens = serde_tokens(max_offset);
        for (index, token, error) in [
            (2, Token::U8(2), "unsupported quantile-sketch version"),
            (
                4,
                Token::F64(0.0),
                "relative error must be representable and between 0 and 1",
            ),
            (
                4,
                Token::F64(f64::EPSILON / 4.0),
                "relative error must be representable and between 0 and 1",
            ),
            (6, Token::F64(f64::from_bits(1)), "invalid sketch range"),
            (8, Token::F64(f64::INFINITY), "invalid sketch range"),
            (13, Token::U64(0), "invalid sketch bin"),
            (16, Token::U64(max_offset + 1), "invalid sketch bin"),
            (16, Token::U64(0), "duplicate sketch bin"),
        ] {
            let mut invalid = tokens.clone();
            invalid[index] = token;
            assert_de_tokens_error::<ConcurrentDDSketch>(&invalid, error);
        }
    }

    #[test]
    fn quantiles_have_expected_relative_error() {
        const ERROR: f64 = 0.01;
        #[cfg(miri)]
        const SAMPLES: usize = 8_000;
        #[cfg(not(miri))]
        const SAMPLES: usize = 100_000;
        #[cfg(miri)]
        const CHECKS: usize = 801;
        #[cfg(not(miri))]
        const CHECKS: usize = 1_001;

        let sketch = ConcurrentDDSketch::with_err_and_range(ERROR, 1.0, SAMPLES as f64);
        assert_eq!(sketch.quantile(0.5), None);
        for value in 1..=SAMPLES {
            sketch.insert(value as f64);
        }

        let mut error_sum = 0.0;
        let mut max_error = 0.0_f64;
        for sample in 0..CHECKS {
            let q = sample as f64 / (CHECKS - 1) as f64;
            let exact = 1.0 + (q * (SAMPLES - 1) as f64).floor();
            let estimate = sketch.quantile(q).unwrap();
            let error = (estimate - exact).abs() / exact;
            error_sum += error;
            max_error = max_error.max(error);
        }

        let average_error = error_sum / CHECKS as f64;
        assert!(average_error <= ERROR * 0.55, "average error: {average_error}");
        assert!(max_error <= ERROR * 1.001, "maximum error: {max_error}");
    }
}