card_est_array/impls/
hyper_log_log.rs

1/*
2 * SPDX-FileCopyrightText: 2024 Matteo Dell'Acqua
3 * SPDX-FileCopyrightText: 2025 Sebastiano Vigna
4 *
5 * SPDX-License-Identifier: Apache-2.0 OR LGPL-2.1-or-later
6 */
7
8use anyhow::{Result, ensure};
9use common_traits::{CastableFrom, CastableInto, Number, UpcastableInto};
10use std::hash::*;
11use std::{borrow::Borrow, f64::consts::LN_2};
12use sux::{bits::BitFieldVec, traits::Word};
13use value_traits::slices::{SliceByValue, SliceByValueMut};
14
15use crate::traits::{EstimationLogic, MergeEstimationLogic, SliceEstimationLogic};
16
17use super::DefaultEstimator;
18
19/// The type returned by the hash function.
20type HashResult = u64;
21
22/// Estimator logic implementing the HyperLogLog algorithm.
23///
24/// Instances are built using [`HyperLogLogBuilder`], which provides convenient
25/// ways to set the internal parameters.
26///
27/// Note that `T` can be any type satisfying the [`Hash`] trait. The parameter
28/// `H` makes it possible to select a hashing algorithm, and `W` is the unsigned
29/// type used to store backends.
30///
31/// An important constraint is that `W` must be able to represent exactly the
32/// backend of an estimator. While usually `usize` will work (and it is the default
33/// type chosen by [`new`](HyperLogLogBuilder::new)), with odd register sizes
34/// and small number of registers it might be necessary to select a smaller
35/// type, resulting in slower merges. For example, using 16 5-bit registers one
36/// needs to use `u16`, whereas for 16 6-bit registers `u32` will be sufficient.
37#[derive(Debug, PartialEq)]
38pub struct HyperLogLog<T, H, W> {
39    build_hasher: H,
40    register_size: usize,
41    num_registers_minus_1: HashResult,
42    log_2_num_registers: usize,
43    sentinel_mask: HashResult,
44    num_registers: usize,
45    pub(super) words_per_estimator: usize,
46    alpha_m_m: f64,
47    msb_mask: Box<[W]>,
48    lsb_mask: Box<[W]>,
49    _marker: std::marker::PhantomData<T>,
50}
51
52// We implement Clone manually because we do not want to require that T is
53// Clone.
54impl<T, H: Clone, W: Clone> Clone for HyperLogLog<T, H, W> {
55    fn clone(&self) -> Self {
56        Self {
57            build_hasher: self.build_hasher.clone(),
58            register_size: self.register_size,
59            num_registers_minus_1: self.num_registers_minus_1,
60            log_2_num_registers: self.log_2_num_registers,
61            sentinel_mask: self.sentinel_mask,
62            num_registers: self.num_registers,
63            words_per_estimator: self.words_per_estimator,
64            alpha_m_m: self.alpha_m_m,
65            msb_mask: self.msb_mask.clone(),
66            lsb_mask: self.lsb_mask.clone(),
67            _marker: std::marker::PhantomData,
68        }
69    }
70}
71
72impl<T, H: Clone, W: Word> HyperLogLog<T, H, W> {
73    /// Returns the value contained in a register of a given backend.
74    #[inline(always)]
75    fn get_register_unchecked(&self, backend: impl AsRef<[W]>, index: usize) -> W {
76        let backend = backend.as_ref();
77        let bit_width = self.register_size;
78        let mask = W::MAX >> (W::BITS - bit_width);
79        let pos = index * bit_width;
80        let word_index = pos / W::BITS;
81        let bit_index = pos % W::BITS;
82
83        if bit_index + bit_width <= W::BITS {
84            (unsafe { *backend.get_unchecked(word_index) } >> bit_index) & mask
85        } else {
86            ((unsafe { *backend.get_unchecked(word_index) } >> bit_index)
87                | (unsafe { *backend.get_unchecked(word_index + 1) } << (W::BITS - bit_index)))
88                & mask
89        }
90    }
91
92    /// Sets the value contained in a register of a given backend.
93    #[inline(always)]
94    fn set_register_unchecked(&self, mut backend: impl AsMut<[W]>, index: usize, new_value: W) {
95        let backend = backend.as_mut();
96        let bit_width = self.register_size;
97        let mask = W::MAX >> (W::BITS - bit_width);
98        let pos = index * bit_width;
99        let word_index = pos / W::BITS;
100        let bit_index = pos % W::BITS;
101
102        if bit_index + bit_width <= W::BITS {
103            let mut word = unsafe { *backend.get_unchecked_mut(word_index) };
104            word &= !(mask << bit_index);
105            word |= new_value << bit_index;
106            unsafe { *backend.get_unchecked_mut(word_index) = word };
107        } else {
108            let mut word = unsafe { *backend.get_unchecked_mut(word_index) };
109            word &= (W::ONE << bit_index) - W::ONE;
110            word |= new_value << bit_index;
111            unsafe { *backend.get_unchecked_mut(word_index) = word };
112
113            let mut word = unsafe { *backend.get_unchecked_mut(word_index + 1) };
114            word &= !(mask >> (W::BITS - bit_index));
115            word |= new_value >> (W::BITS - bit_index);
116            unsafe { *backend.get_unchecked_mut(word_index + 1) = word };
117        }
118    }
119}
120
121impl<
122    T: Hash,
123    H: BuildHasher + Clone,
124    W: Word + UpcastableInto<HashResult> + CastableFrom<HashResult>,
125> SliceEstimationLogic<W> for HyperLogLog<T, H, W>
126{
127    fn backend_len(&self) -> usize {
128        self.words_per_estimator
129    }
130}
131
132impl<
133    T: Hash,
134    H: BuildHasher + Clone,
135    W: Word + UpcastableInto<HashResult> + CastableFrom<HashResult>,
136> EstimationLogic for HyperLogLog<T, H, W>
137{
138    type Item = T;
139    type Backend = [W];
140    type Estimator<'a>
141        = DefaultEstimator<Self, &'a Self, Box<[W]>>
142    where
143        T: 'a,
144        W: 'a,
145        H: 'a;
146
147    fn new_estimator(&self) -> Self::Estimator<'_> {
148        Self::Estimator::new(
149            self,
150            vec![W::ZERO; self.words_per_estimator].into_boxed_slice(),
151        )
152    }
153
154    fn add(&self, mut backend: &mut Self::Backend, element: impl Borrow<T>) {
155        let x = self.build_hasher.hash_one(element.borrow());
156        let j = x & self.num_registers_minus_1;
157        let r =
158            ((x >> self.log_2_num_registers) | self.sentinel_mask).trailing_zeros() as HashResult;
159        let register = j as usize;
160
161        debug_assert!(r < (1 << self.register_size) - 1);
162        debug_assert!(register < self.num_registers);
163
164        let current_value = self.get_register_unchecked(&mut backend, register);
165        let candidate_value = r + 1;
166        let new_value = std::cmp::max(current_value, candidate_value.cast());
167        if current_value != new_value {
168            self.set_register_unchecked(backend, register, new_value);
169        }
170    }
171
172    fn estimate(&self, backend: &[W]) -> f64 {
173        let mut harmonic_mean = 0.0;
174        let mut zeroes = 0;
175
176        for i in 0..self.num_registers {
177            let value: u64 = self.get_register_unchecked(backend, i).upcast();
178            if value == 0 {
179                zeroes += 1;
180            }
181            harmonic_mean += 1.0 / (1_u64 << value) as f64;
182        }
183
184        let mut estimate = self.alpha_m_m / harmonic_mean;
185        if zeroes != 0 && estimate < 2.5 * self.num_registers as f64 {
186            estimate = self.num_registers as f64 * (self.num_registers as f64 / zeroes as f64).ln();
187        }
188        estimate
189    }
190
191    fn clear(&self, backend: &mut [W]) {
192        backend.as_mut().fill(W::ZERO);
193    }
194
195    fn set(&self, dst: &mut [W], src: &[W]) {
196        debug_assert_eq!(dst.as_mut().len(), src.as_ref().len());
197        dst.as_mut().copy_from_slice(src.as_ref());
198    }
199}
200
201/// Helper for merge operations with [`HyperLogLog`] logic.
202pub struct HyperLogLogHelper<W> {
203    acc: Vec<W>,
204    mask: Vec<W>,
205}
206
207impl<
208    T: Hash,
209    H: BuildHasher + Clone,
210    W: Word + UpcastableInto<HashResult> + CastableFrom<HashResult>,
211> MergeEstimationLogic for HyperLogLog<T, H, W>
212{
213    type Helper = HyperLogLogHelper<W>;
214
215    fn new_helper(&self) -> Self::Helper {
216        HyperLogLogHelper {
217            acc: vec![W::ZERO; self.words_per_estimator],
218            mask: vec![W::ZERO; self.words_per_estimator],
219        }
220    }
221
222    fn merge_with_helper(&self, dst: &mut [W], src: &[W], helper: &mut Self::Helper) {
223        merge_hyperloglog_bitwise(
224            dst,
225            src,
226            self.msb_mask.as_ref(),
227            self.lsb_mask.as_ref(),
228            &mut helper.acc,
229            &mut helper.mask,
230            self.register_size,
231        );
232    }
233}
234
235/// Builds a [`HyperLogLog`] cardinality-estimator logic.
236#[derive(Debug, Clone)]
237pub struct HyperLogLogBuilder<H, W = usize> {
238    build_hasher: H,
239    log_2_num_registers: usize,
240    n: usize,
241    _marker: std::marker::PhantomData<(H, W)>,
242}
243
244impl HyperLogLogBuilder<BuildHasherDefault<DefaultHasher>, usize> {
245    /// Creates a new builder for a [`HyperLogLog`] logic with a default word
246    /// type of `usize`.
247    pub fn new(n: usize) -> Self {
248        Self {
249            build_hasher: BuildHasherDefault::default(),
250            log_2_num_registers: 4,
251            n,
252            _marker: std::marker::PhantomData,
253        }
254    }
255}
256
257fn min_alignment(bits: usize) -> String {
258    if bits % 128 == 0 {
259        "u128"
260    } else if bits % 64 == 0 {
261        "u64"
262    } else if bits % 32 == 0 {
263        "u32"
264    } else if bits % 16 == 0 {
265        "u16"
266    } else {
267        "u8"
268    }
269    .to_string()
270}
271
272impl HyperLogLog<(), (), ()> {
273    /// Returns the logarithm of the number of registers per estimator that are
274    /// necessary to attain a given relative standard deviation.
275    ///
276    /// # Arguments
277    /// * `rsd`: the relative standard deviation to be attained.
278    pub fn log_2_num_of_registers(rsd: f64) -> usize {
279        ((1.106 / rsd).pow(2.0)).log2().ceil() as usize
280    }
281
282    /// Returns the relative standard deviation corresponding to a given number
283    /// of registers per estimator.
284    ///
285    /// # Arguments
286    ///
287    /// * `log_2_num_registers`: the logarithm of the number of registers per
288    ///   estimator.
289    pub fn rel_std(log_2_num_registers: usize) -> f64 {
290        let tmp = match log_2_num_registers {
291            4 => 1.106,
292            5 => 1.070,
293            6 => 1.054,
294            7 => 1.046,
295            _ => 1.04,
296        };
297        tmp / ((1 << log_2_num_registers) as f64).sqrt()
298    }
299
300    /// Returns the register size in bits, given an upper bound on the number of
301    /// distinct elements.
302    ///
303    /// # Arguments
304    /// * `n`: an upper bound on the number of distinct elements.
305    pub fn register_size(n: usize) -> usize {
306        std::cmp::max(5, (((n as f64).ln() / LN_2) / LN_2).ln().ceil() as usize)
307    }
308}
309
310impl<H, W: Word> HyperLogLogBuilder<H, W> {
311    /// Sets the desired relative standard deviation.
312    ///
313    /// ## Note
314    ///
315    /// This is a high-level alternative to [`Self::log_2_num_reg`]. Calling one
316    /// after the other invalidates the work done by the first one.
317    ///
318    /// # Arguments
319    /// * `rsd`: the relative standard deviation to be attained.
320    pub fn rsd(self, rsd: f64) -> Self {
321        self.log_2_num_reg(HyperLogLog::log_2_num_of_registers(rsd))
322    }
323
324    /// Sets the base-2 logarithm of the number of register.
325    ///
326    /// ## Note
327    /// This is a low-level alternative to [`Self::rsd`]. Calling one after the
328    /// other invalidates the work done by the first one.
329    ///
330    /// # Arguments
331    /// * `log_2_num_registers`: the logarithm of the number of registers per
332    ///   estimator.
333    pub fn log_2_num_reg(mut self, log_2_num_registers: usize) -> Self {
334        self.log_2_num_registers = log_2_num_registers;
335        self
336    }
337
338    /// Sets the type `W` to use to represent backends.
339    ///
340    /// See the [`logic documentation`](HyperLogLog) for the limitations on the
341    /// choice of `W2`.
342    pub fn word_type<W2>(self) -> HyperLogLogBuilder<H, W2> {
343        HyperLogLogBuilder {
344            n: self.n,
345            build_hasher: self.build_hasher,
346            log_2_num_registers: self.log_2_num_registers,
347            _marker: std::marker::PhantomData,
348        }
349    }
350
351    /// Sets the upper bound on the number of elements.
352    pub fn num_elements(mut self, n: usize) -> Self {
353        self.n = n;
354        self
355    }
356
357    /// Sets the [`BuildHasher`] to use.
358    ///
359    /// Using this method you can select a specific hashed based on one or more
360    /// seeds.
361    pub fn build_hasher<H2>(self, build_hasher: H2) -> HyperLogLogBuilder<H2, W> {
362        HyperLogLogBuilder {
363            n: self.n,
364            log_2_num_registers: self.log_2_num_registers,
365            build_hasher,
366            _marker: std::marker::PhantomData,
367        }
368    }
369
370    /// Builds the logic.
371    ///
372    /// The type of objects the estimators keep track of is defined here by `T`,
373    /// but it is usually inferred by the compiler.
374    ///
375    /// # Errors
376    ///
377    /// Errors will be caused by consistency checks (at least 16 registers per
378    /// estimator, backend bits divisible exactly `W::BITS`)
379    pub fn build<T>(self) -> Result<HyperLogLog<T, H, W>> {
380        let log_2_num_registers = self.log_2_num_registers;
381        let num_elements = self.n;
382
383        // This ensures estimators are at least 16-bit-aligned.
384        ensure!(
385            log_2_num_registers >= 4,
386            "the logarithm of the number of registers per estimator should be at least 4; got {}",
387            log_2_num_registers
388        );
389
390        let number_of_registers = 1 << log_2_num_registers;
391        let register_size = HyperLogLog::register_size(num_elements);
392        let sentinel_mask = 1 << ((1 << register_size) - 2);
393        let alpha = match log_2_num_registers {
394            4 => 0.673,
395            5 => 0.697,
396            6 => 0.709,
397            _ => 0.7213 / (1.0 + 1.079 / number_of_registers as f64),
398        };
399        let num_registers_minus_1 = (number_of_registers - 1) as HashResult;
400
401        let est_size_in_bits = number_of_registers * register_size;
402
403        // This ensures estimators are always aligned to W
404        ensure!(
405            est_size_in_bits % W::BITS == 0,
406            "W should allow estimator backends to be aligned. Use {} or smaller unsigned integer types",
407            min_alignment(est_size_in_bits)
408        );
409        let est_size_in_words = est_size_in_bits / W::BITS;
410
411        let mut msb = BitFieldVec::new(register_size, number_of_registers);
412        let mut lsb = BitFieldVec::new(register_size, number_of_registers);
413        let msb_w = W::ONE << (register_size - 1);
414        let lsb_w = W::ONE;
415        for i in 0..number_of_registers {
416            msb.set_value(i, msb_w);
417            lsb.set_value(i, lsb_w);
418        }
419
420        Ok(HyperLogLog {
421            num_registers: number_of_registers,
422            num_registers_minus_1,
423            log_2_num_registers,
424            register_size,
425            alpha_m_m: alpha * (number_of_registers as f64).powi(2),
426            sentinel_mask,
427            build_hasher: self.build_hasher,
428            msb_mask: msb.as_slice().into(),
429            lsb_mask: lsb.as_slice().into(),
430            words_per_estimator: est_size_in_words,
431            _marker: std::marker::PhantomData,
432        })
433    }
434}
435
436impl<T, H, W> std::fmt::Display for HyperLogLog<T, H, W> {
437    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
438        write!(
439            f,
440            "HyperLogLog with relative standard deviation: {}% ({} registers/estimator, {} bits/register, {} bytes/estimator)",
441            100.0 * HyperLogLog::rel_std(self.log_2_num_registers),
442            self.num_registers,
443            self.register_size,
444            (self.num_registers * self.register_size) / 8
445        )
446    }
447}
448
449/// Performs a multiple precision subtraction, leaving the result in the first operand.
450/// The operands MUST have the same length.
451///
452/// # Arguments
453/// * `x`: the first operand. This will contain the final result.
454/// * `y`: the second operand that will be subtracted from `x`.
455#[inline(always)]
456pub(super) fn subtract<W: Word>(x: &mut [W], y: &[W]) {
457    debug_assert_eq!(x.len(), y.len());
458    let mut borrow = false;
459
460    for (x_word, &y) in x.iter_mut().zip(y.iter()) {
461        let mut x = *x_word;
462        if !borrow {
463            borrow = x < y;
464        } else if x != W::ZERO {
465            x = x.wrapping_sub(W::ONE);
466            borrow = x < y;
467        } else {
468            x = x.wrapping_sub(W::ONE);
469        }
470        *x_word = x.wrapping_sub(y);
471    }
472}
473
474fn merge_hyperloglog_bitwise<W: Word>(
475    mut x: impl AsMut<[W]>,
476    y: impl AsRef<[W]>,
477    msb_mask: impl AsRef<[W]>,
478    lsb_mask: impl AsRef<[W]>,
479    acc: &mut Vec<W>,
480    mask: &mut Vec<W>,
481    register_size: usize,
482) {
483    let x = x.as_mut();
484    let y = y.as_ref();
485    let msb_mask = msb_mask.as_ref();
486    let lsb_mask = lsb_mask.as_ref();
487
488    debug_assert_eq!(x.len(), y.len());
489    debug_assert_eq!(x.len(), msb_mask.len());
490    debug_assert_eq!(x.len(), lsb_mask.len());
491
492    let register_size_minus_1 = register_size - 1;
493    let num_words_minus_1 = x.len() - 1;
494    let shift_register_size_minus_1 = W::BITS - register_size_minus_1;
495
496    acc.clear();
497    mask.clear();
498
499    /* We work in two phases. Let H_r (msb_mask) be the mask with the
500     * highest bit of each register (of size r) set, and L_r (lsb_mask)
501     * be the mask with the lowest bit of each register set.
502     * We describe the algorithm on a single word.
503     *
504     * In the first phase we perform an unsigned strict register-by-register
505     * comparison of x and y, using the formula
506     *
507     * z = ((((y | H_r) - (x & !H_r)) | (y ^ x)) ^ (y | !x)) & H_r
508     *
509     * Then, we generate a register-by-register mask of all ones or
510     * all zeroes, depending on the result of the comparison, using the
511     * formula
512     *
513     * (((z >> r-1 | H_r) - L_r) | H_r) ^ z
514     *
515     * At that point, it is trivial to select from x and y the right values.
516     */
517
518    // We load y | H_r into the accumulator.
519    acc.extend(
520        y.iter()
521            .zip(msb_mask)
522            .map(|(&y_word, &msb_word)| y_word | msb_word),
523    );
524
525    // We load x & !H_r into mask as temporary storage.
526    mask.extend(
527        x.iter()
528            .zip(msb_mask)
529            .map(|(&x_word, &msb_word)| x_word & !msb_word),
530    );
531
532    // We subtract x & !H_r, using mask as temporary storage
533    subtract(acc, mask);
534
535    // We OR with y ^ x, XOR with (y | !x), and finally AND with H_r.
536    acc.iter_mut()
537        .zip(x.iter())
538        .zip(y.iter())
539        .zip(msb_mask.iter())
540        .for_each(|(((acc_word, &x_word), &y_word), &msb_word)| {
541            *acc_word = ((*acc_word | (y_word ^ x_word)) ^ (y_word | !x_word)) & msb_word
542        });
543
544    // We shift by register_size - 1 places and put the result into mask.
545    {
546        let (mask_last, mask_slice) = mask.split_last_mut().unwrap();
547        let (&msb_last, msb_slice) = msb_mask.split_last().unwrap();
548        mask_slice
549            .iter_mut()
550            .zip(acc[0..num_words_minus_1].iter())
551            .zip(acc[1..].iter())
552            .zip(msb_slice.iter())
553            .rev()
554            .for_each(|(((mask_word, &acc_word), &next_acc_word), &msb_word)| {
555                // W is always unsigned so the shift is always with a 0
556                *mask_word = (acc_word >> register_size_minus_1)
557                    | (next_acc_word << shift_register_size_minus_1)
558                    | msb_word
559            });
560        *mask_last = (acc[num_words_minus_1] >> register_size_minus_1) | msb_last;
561    }
562
563    // We subtract L_r from mask.
564    subtract(mask, lsb_mask);
565
566    // We OR with H_r and XOR with the accumulator.
567    mask.iter_mut()
568        .zip(msb_mask.iter())
569        .zip(acc.iter())
570        .for_each(|((mask_word, &msb_word), &acc_word)| {
571            *mask_word = (*mask_word | msb_word) ^ acc_word
572        });
573
574    // Finally, we use mask to select the right bits from x and y and store the result.
575    x.iter_mut()
576        .zip(y.iter())
577        .zip(mask.iter())
578        .for_each(|((x_word, &y_word), &mask_word)| {
579            *x_word = *x_word ^ ((*x_word ^ y_word) & mask_word);
580        });
581}
card_est_array/impls/hyper_log_log.rs

card_est_array/impls/
hyper_log_log.rs
