HyperLogLog: Harmonic Mean & Stochastic Averaging in 1.5 KB

From Flajolet-Martin to HyperLogLog Harmonic Mean

Given a 64-bit uniform hash $h(x)$ and precision parameter $p \in [4, 18]$ (so $m = 2^p$ registers $M[0..m-1]$ initialized to $0$):

1. Substream Bucketing (Stochastic Averaging): Extract the top $p$ bits of $h(x)$ as the register index $j = h(x) \gg (64 - p) \in [0, m-1]$.
2. Leading-Zero Rank $\rho(w)$: Let $w$ be the remaining $64 - p$ bits of $h(x)$. Define $\rho(w) = \text{clz}(w) + 1 \in [1, 64 - p + 1]$ as the 1-based position of the most significant 1-bit in $w$. Update register $M[j] \gets \max(M[j], \rho(w))$.
3. Harmonic Mean Estimator: Why did earlier LogLog (2003) use the geometric mean ($1.30 / \sqrt{m}$ error) while HyperLogLog (2007) uses the harmonic mean ($1.04 / \sqrt{m}$ error)? Because $2^{M[j]}$ has a heavy right-skewed tail ($\mathbb{E}[2^{M[j]}] = \infty$), an arithmetic mean explodes if a single register hits $\rho = 40$. The harmonic mean averages $2^{-M[j]}$, completely neutralizing large outlier registers:

   $$E_{\text{raw}} = \alpha_m \cdot m^2 \cdot \left( \sum_{j=0}^{m-1} 2^{-M[j]} \right)^{-1}$$

   where $\alpha_m = \left(m \int_0^\infty \left(\log_2\left(\frac{2+u}{1+u}\right)\right)^m du\right)^{-1} \approx \frac{0.7213}{1 + 1.079/m}$ corrects asymptotic multiplicative bias.
4. Small-Range LinearCounting & Sparse Representations: When $n$ is small and $V > 0$ registers are still $0$, and $E_{\text{raw}} \le \frac{5}{2} m$, HLL switches to LinearCounting (Whang et al., 1990) based on the coupon-collector occupancy rate: $E = m \ln(m / V)$.

Visualizing 64-Bit Hash Splitting ($p = 14$) & Harmonic Mean Aggregation

Production C++20 Implementation: 64-Bit HyperLogLog with LinearCounting & Merge

// C++20: HyperLogLog (64-bit hash, p in [4..18]) with Sentinel CLZ & LinearCounting
#include <vector>
#include <cstdint>
#include <bit>
#include <cmath>
#include <algorithm>

class HyperLogLog {
    uint32_t p_;
    uint32_t m_;
    double alphaMM_;
    std::vector<uint8_t> regs_;

    static uint64_t hash64(uint64_t x) noexcept {
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL;
        x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL;
        return x ^ (x >> 31);
    }

public:
    explicit HyperLogLog(uint32_t p = 14)
        : p_(std::clamp(p, 4U, 18U)), m_(1U << p_), regs_(m_, 0) {
        double alpha = 0.7213 / (1.0 + 1.079 / static_cast<double>(m_));
        if (m_ == 16) alpha = 0.673;
        else if (m_ == 32) alpha = 0.697;
        else if (m_ == 64) alpha = 0.709;
        alphaMM_ = alpha * static_cast<double>(m_) * static_cast<double>(m_);
    }

    void add(uint64_t item) noexcept {
        uint64_t h = hash64(item);
        uint32_t j = static_cast<uint32_t>(h >> (64 - p_));
        // Shift remaining (64 - p) bits to MSB and set sentinel bit at (p - 1)
        uint64_t w = (h << p_) | (1ULL << (p_ - 1));
        uint8_t rho = static_cast<uint8_t>(std::countl_zero(w) + 1);
        if (rho > regs_[j]) regs_[j] = rho;
    }

    double count() const noexcept {
        double sumInv = 0.0;
        uint32_t zeroRegs = 0;
        for (uint8_t r : regs_) {
            sumInv += std::ldexp(1.0, -static_cast<int>(r)); // 2^{-r}
            if (r == 0) ++zeroRegs;
        }
        double est = alphaMM_ / sumInv;
        if (est <= 2.5 * static_cast<double>(m_) && zeroRegs > 0) {
            return static_cast<double>(m_) * std::log(static_cast<double>(m_) / zeroRegs);
        }
        return est;
    }

    // Lossless Commutative Semilattice Merge (Redis PFMERGE)
    void merge(const HyperLogLog& other) noexcept {
        for (uint32_t j = 0; j < m_; ++j) {
            regs_[j] = std::max(regs_[j], other.regs_[j]);
        }
    }
};