1/*-------------------------------------------------------------------------
2 *
3 * hyperloglog.c
4 * HyperLogLog cardinality estimator
5 *
6 * Portions Copyright (c) 2014-2019, PostgreSQL Global Development Group
7 *
8 * Based on Hideaki Ohno's C++ implementation. This is probably not ideally
9 * suited to estimating the cardinality of very large sets; in particular, we
10 * have not attempted to further optimize the implementation as described in
11 * the Heule, Nunkesser and Hall paper "HyperLogLog in Practice: Algorithmic
12 * Engineering of a State of The Art Cardinality Estimation Algorithm".
13 *
14 * A sparse representation of HyperLogLog state is used, with fixed space
15 * overhead.
16 *
17 * The copyright terms of Ohno's original version (the MIT license) follow.
18 *
19 * IDENTIFICATION
20 * src/backend/lib/hyperloglog.c
21 *
22 *-------------------------------------------------------------------------
23 */
24
25/*
26 * Copyright (c) 2013 Hideaki Ohno <hide.o.j55{at}gmail.com>
27 *
28 * Permission is hereby granted, free of charge, to any person obtaining a copy
29 * of this software and associated documentation files (the 'Software'), to
30 * deal in the Software without restriction, including without limitation the
31 * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
32 * sell copies of the Software, and to permit persons to whom the Software is
33 * furnished to do so, subject to the following conditions:
34 *
35 * The above copyright notice and this permission notice shall be included in
36 * all copies or substantial portions of the Software.
37 *
38 * THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
39 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
40 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
41 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
42 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
43 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
44 * IN THE SOFTWARE.
45 */
46
47#include "postgres.h"
48
49#include <math.h>
50
51#include "lib/hyperloglog.h"
52
53#define POW_2_32 (4294967296.0)
54#define NEG_POW_2_32 (-4294967296.0)
55
56static inline uint8 rho(uint32 x, uint8 b);
57
58/*
59 * Initialize HyperLogLog track state, by bit width
60 *
61 * bwidth is bit width (so register size will be 2 to the power of bwidth).
62 * Must be between 4 and 16 inclusive.
63 */
64void
65initHyperLogLog(hyperLogLogState *cState, uint8 bwidth)
66{
67 double alpha;
68
69 if (bwidth < 4 || bwidth > 16)
70 elog(ERROR, "bit width must be between 4 and 16 inclusive");
71
72 cState->registerWidth = bwidth;
73 cState->nRegisters = (Size) 1 << bwidth;
74 cState->arrSize = sizeof(uint8) * cState->nRegisters + 1;
75
76 /*
77 * Initialize hashes array to zero, not negative infinity, per discussion
78 * of the coupon collector problem in the HyperLogLog paper
79 */
80 cState->hashesArr = palloc0(cState->arrSize);
81
82 /*
83 * "alpha" is a value that for each possible number of registers (m) is
84 * used to correct a systematic multiplicative bias present in m ^ 2 Z (Z
85 * is "the indicator function" through which we finally compute E,
86 * estimated cardinality).
87 */
88 switch (cState->nRegisters)
89 {
90 case 16:
91 alpha = 0.673;
92 break;
93 case 32:
94 alpha = 0.697;
95 break;
96 case 64:
97 alpha = 0.709;
98 break;
99 default:
100 alpha = 0.7213 / (1.0 + 1.079 / cState->nRegisters);
101 }
102
103 /*
104 * Precalculate alpha m ^ 2, later used to generate "raw" HyperLogLog
105 * estimate E
106 */
107 cState->alphaMM = alpha * cState->nRegisters * cState->nRegisters;
108}
109
110/*
111 * Initialize HyperLogLog track state, by error rate
112 *
113 * Instead of specifying bwidth (number of bits used for addressing the
114 * register), this method allows sizing the counter for particular error
115 * rate using a simple formula from the paper:
116 *
117 * e = 1.04 / sqrt(m)
118 *
119 * where 'm' is the number of registers, i.e. (2^bwidth). The method
120 * finds the lowest bwidth with 'e' below the requested error rate, and
121 * then uses it to initialize the counter.
122 *
123 * As bwidth has to be between 4 and 16, the worst possible error rate
124 * is between ~25% (bwidth=4) and 0.4% (bwidth=16).
125 */
126void
127initHyperLogLogError(hyperLogLogState *cState, double error)
128{
129 uint8 bwidth = 4;
130
131 while (bwidth < 16)
132 {
133 double m = (Size) 1 << bwidth;
134
135 if (1.04 / sqrt(m) < error)
136 break;
137 bwidth++;
138 }
139
140 initHyperLogLog(cState, bwidth);
141}
142
143/*
144 * Free HyperLogLog track state
145 *
146 * Releases allocated resources, but not the state itself (in case it's not
147 * allocated by palloc).
148 */
149void
150freeHyperLogLog(hyperLogLogState *cState)
151{
152 Assert(cState->hashesArr != NULL);
153 pfree(cState->hashesArr);
154}
155
156/*
157 * Adds element to the estimator, from caller-supplied hash.
158 *
159 * It is critical that the hash value passed be an actual hash value, typically
160 * generated using hash_any(). The algorithm relies on a specific bit-pattern
161 * observable in conjunction with stochastic averaging. There must be a
162 * uniform distribution of bits in hash values for each distinct original value
163 * observed.
164 */
165void
166addHyperLogLog(hyperLogLogState *cState, uint32 hash)
167{
168 uint8 count;
169 uint32 index;
170
171 /* Use the first "k" (registerWidth) bits as a zero based index */
172 index = hash >> (BITS_PER_BYTE * sizeof(uint32) - cState->registerWidth);
173
174 /* Compute the rank of the remaining 32 - "k" (registerWidth) bits */
175 count = rho(hash << cState->registerWidth,
176 BITS_PER_BYTE * sizeof(uint32) - cState->registerWidth);
177
178 cState->hashesArr[index] = Max(count, cState->hashesArr[index]);
179}
180
181/*
182 * Estimates cardinality, based on elements added so far
183 */
184double
185estimateHyperLogLog(hyperLogLogState *cState)
186{
187 double result;
188 double sum = 0.0;
189 int i;
190
191 for (i = 0; i < cState->nRegisters; i++)
192 {
193 sum += 1.0 / pow(2.0, cState->hashesArr[i]);
194 }
195
196 /* result set to "raw" HyperLogLog estimate (E in the HyperLogLog paper) */
197 result = cState->alphaMM / sum;
198
199 if (result <= (5.0 / 2.0) * cState->nRegisters)
200 {
201 /* Small range correction */
202 int zero_count = 0;
203
204 for (i = 0; i < cState->nRegisters; i++)
205 {
206 if (cState->hashesArr[i] == 0)
207 zero_count++;
208 }
209
210 if (zero_count != 0)
211 result = cState->nRegisters * log((double) cState->nRegisters /
212 zero_count);
213 }
214 else if (result > (1.0 / 30.0) * POW_2_32)
215 {
216 /* Large range correction */
217 result = NEG_POW_2_32 * log(1.0 - (result / POW_2_32));
218 }
219
220 return result;
221}
222
223/*
224 * Worker for addHyperLogLog().
225 *
226 * Calculates the position of the first set bit in first b bits of x argument
227 * starting from the first, reading from most significant to least significant
228 * bits.
229 *
230 * Example (when considering fist 10 bits of x):
231 *
232 * rho(x = 0b1000000000) returns 1
233 * rho(x = 0b0010000000) returns 3
234 * rho(x = 0b0000000000) returns b + 1
235 *
236 * "The binary address determined by the first b bits of x"
237 *
238 * Return value "j" used to index bit pattern to watch.
239 */
240static inline uint8
241rho(uint32 x, uint8 b)
242{
243 uint8 j = 1;
244
245 while (j <= b && !(x & 0x80000000))
246 {
247 j++;
248 x <<= 1;
249 }
250
251 return j;
252}
253