quantile-sketch 0.1.0

A fast, concurrent DDSketch for relative-error quantiles.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130

//! Datadog's cubically interpolated logarithmic mapping.
//!
//! See the [DDSketch paper](https://arxiv.org/abs/1908.10693) and Datadog's
//! [reference implementation](https://github.com/DataDog/sketches-go/blob/master/ddsketch/mapping/cubically_interpolated_mapping.go).

const MANTISSA_MASK: u64 = (1_u64 << 52) - 1;
const ONE_BITS: u64 = 1023_u64 << 52;
const A: f64 = 6.0 / 35.0;
const B: f64 = -3.0 / 5.0;
const C: f64 = 10.0 / 7.0;

pub(crate) struct CubicMapping {
    multiplier: f64,
    error: f64,
}

impl CubicMapping {
    pub(crate) fn is_valid_error(error: f64) -> bool {
        if !(error > 0.0 && error < 1.0) {
            return false;
        }
        let gamma = (1.0 + error) / (1.0 - error);
        gamma > 1.0 && gamma.is_finite()
    }

    pub(crate) fn new(error: f64) -> Self {
        assert!(
            Self::is_valid_error(error),
            "relative error must be representable and between 0 and 1"
        );
        let gamma = (1.0 + error) / (1.0 - error);
        let gamma = pow(gamma, 10.0 * core::f64::consts::LN_2 / 7.0);
        Self {
            multiplier: 1.0 / log2(gamma),
            error,
        }
    }

    pub(crate) fn error(&self) -> f64 {
        self.error
    }

    #[inline]
    pub(crate) fn index(&self, value: f64) -> i64 {
        let bits = value.to_bits();
        let exponent = ((bits >> 52) & 0x7ff) as i64 - 1023;
        let s = f64::from_bits((bits & MANTISSA_MASK) | ONE_BITS) - 1.0;
        let log2 = ((A * s + B) * s + C) * s + exponent as f64;
        floor(log2 * self.multiplier) as i64
    }

    #[inline]
    pub(crate) fn value(&self, index: i64) -> f64 {
        let x = index as f64 / self.multiplier;
        let exponent = floor(x);
        let d0 = B * B - 3.0 * A * C;
        let d1 = 2.0 * B * B * B - 9.0 * A * B * C - 27.0 * A * A * (x - exponent);
        let p = cbrt((d1 - sqrt(d1 * d1 - 4.0 * d0 * d0 * d0)) / 2.0);
        let significand = -(B + p + d0 / p) / (3.0 * A) + 1.0;
        build_f64(exponent as i64, significand) * (1.0 + self.error)
    }
}

#[inline]
fn build_f64(exponent: i64, significand: f64) -> f64 {
    let fraction = ((significand - 1.0) * (1_u64 << 52) as f64) as u64;
    f64::from_bits((((exponent + 1023) as u64) << 52) | (fraction & MANTISSA_MASK))
}

#[inline]
fn floor(value: f64) -> f64 {
    #[cfg(feature = "std")]
    return value.floor();
    #[cfg(not(feature = "std"))]
    return libm::floor(value);
}

#[inline]
fn pow(value: f64, exponent: f64) -> f64 {
    #[cfg(feature = "std")]
    return value.powf(exponent);
    #[cfg(not(feature = "std"))]
    return libm::pow(value, exponent);
}

#[inline]
fn log2(value: f64) -> f64 {
    #[cfg(feature = "std")]
    return value.log2();
    #[cfg(not(feature = "std"))]
    return libm::log2(value);
}

#[inline]
fn sqrt(value: f64) -> f64 {
    #[cfg(feature = "std")]
    return value.sqrt();
    #[cfg(not(feature = "std"))]
    return libm::sqrt(value);
}

#[inline]
fn cbrt(value: f64) -> f64 {
    #[cfg(feature = "std")]
    return value.cbrt();
    #[cfg(not(feature = "std"))]
    return libm::cbrt(value);
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn mapping_has_expected_relative_error() {
        const ERROR: f64 = 0.01;
        #[cfg(miri)]
        const SAMPLES: usize = 8_000;
        #[cfg(not(miri))]
        const SAMPLES: usize = 100_000;

        let mapping = CubicMapping::new(ERROR);
        for sample in 0..=SAMPLES {
            let value = 10.0_f64.powf(-9.0 + 18.0 * sample as f64 / SAMPLES as f64);
            let estimate = mapping.value(mapping.index(value));
            let error = (estimate - value).abs() / value;
            assert!(error <= ERROR * 1.000_001, "{value}: {error}");
        }
    }
}