1#include <math/fixed.h>
2#include <math/hash.h>
3#include <math/hyperloglog.h>
4#include <mem/alloc.h>
5
6/* Notes on naming: all HLL_ and hll_ prefixed functions are INTERNAL,
7 * all hyperloglog_ and HYPERLOGLOG_ prefixed functions are EXTERNAL */
8
9#define HLL_ALPHA_16 FX(0.673)
10#define HLL_ALPHA_32 FX(0.697)
11#define HLL_ALPHA_64 FX(0.709)
12
13struct hyperloglog {
14 uint8_t register_width;
15 fx32_32_t alpha_mm;
16 size_t nr_registers;
17 uint8_t *registers cc_counted_by(nr_registers);
18};
19
20static inline fx32_32_t hll_alpha_128_plus(size_t nr_registers) {
21 fx32_32_t numerator = FX(0.7213);
22 fx32_32_t denominator = fx_div(FX(2.079), b: fx_from_int(x: nr_registers));
23 return fx_div(a: numerator, b: denominator);
24}
25
26struct hyperloglog *hyperloglog_create(uint8_t register_width) {
27 size_t nr_registers = 1 << register_width;
28 uint8_t *reg_storage = vnews(reg_storage, nr_registers);
29 struct hyperloglog *hll_ret = vnew(hll_ret);
30
31 if (!reg_storage || !hll_ret) {
32 vfree(reg_storage);
33 vfree(hll_ret);
34 return NULL;
35 }
36
37 hll_ret->registers = reg_storage;
38 hll_ret->nr_registers = nr_registers;
39 hll_ret->register_width = register_width;
40
41 fx32_32_t alpha;
42 switch (nr_registers) {
43 case 16: alpha = HLL_ALPHA_16; break;
44 case 32: alpha = HLL_ALPHA_32; break;
45 case 64: alpha = HLL_ALPHA_64; break;
46 default: alpha = hll_alpha_128_plus(nr_registers); break;
47 }
48
49 hll_ret->alpha_mm = fx_mul(a: fx_mul(a: alpha, b: fx_from_int(x: nr_registers)),
50 b: fx_from_int(x: nr_registers));
51 return hll_ret;
52}
53
54void hyperloglog_destroy(struct hyperloglog *hll) {
55 vfree(hll->registers);
56 vfree(hll);
57}
58
59void hyperloglog_add(struct hyperloglog *hll, uint32_t hash) {
60 uint8_t rw = hll->register_width;
61
62 uint32_t idx = hash >> (32U - rw);
63 uint8_t rank = MIN(32U - hll->register_width, cw_clz(hash << rw)) + 1;
64
65 hll->registers[idx] = MAX(rank, hll->registers[idx]);
66}
67
68fx32_32_t hyperloglog_estimate(struct hyperloglog *hll) {
69 fx32_32_t sum = FX(0.0);
70 size_t zero_registers = 0;
71
72 for (size_t i = 0; i < hll->nr_registers; i++) {
73 uint8_t val = hll->registers[i];
74 if (val == 0)
75 zero_registers++;
76
77 /* 1.0 / 2^val */
78 sum += fx_div(FX(1.0), b: fx_from_int(x: 1ULL << val));
79 }
80
81 /* estimate E = alpha_m * m^2 / sum */
82 fx32_32_t estimate = fx_div(a: hll->alpha_mm, b: sum);
83
84 /* range correction: E <= (5/2) * m */
85 fx32_32_t threshold_small = fx_mul(FX(2.5), b: fx_from_int(x: hll->nr_registers));
86 if (estimate <= threshold_small) {
87 if (zero_registers > 0) {
88 /* linear counting: m * ln(m / V) */
89 fx32_32_t m = fx_from_int(x: hll->nr_registers);
90 fx32_32_t v = fx_from_int(x: zero_registers);
91 estimate = fx_mul(a: m, b: fx_ln(x: fx_div(a: m, b: v)));
92 }
93 return estimate;
94 }
95
96 /* Large range correction: E > (1/30) * 2^32 */
97 fx32_32_t threshold_large =
98 FX(143165576.533); /* (1.0 / 30.0) * (1ULL << 32) */
99
100 if (estimate > threshold_large) {
101 /* -2^32 * ln(1 - E / 2^32) */
102 fx32_32_t two_32 = fx_from_int(x: 1ULL << 32);
103 fx32_32_t ratio = fx_div(a: estimate, b: two_32);
104
105 if (ratio < FX(1.0))
106 estimate = fx_mul(a: -two_32, b: fx_ln(FX(1.0) - ratio));
107 }
108
109 return estimate;
110}
111