Line data Source code
1 : /*-------------------------------------------------------------------------
2 : *
3 : * hyperloglog.c
4 : * HyperLogLog cardinality estimator
5 : *
6 : * Portions Copyright (c) 2014-2017, 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 :
56 : static 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 : */
64 : void
65 138 : initHyperLogLog(hyperLogLogState *cState, uint8 bwidth)
66 : {
67 : double alpha;
68 :
69 138 : if (bwidth < 4 || bwidth > 16)
70 0 : elog(ERROR, "bit width must be between 4 and 16 inclusive");
71 :
72 138 : cState->registerWidth = bwidth;
73 138 : cState->nRegisters = (Size) 1 << bwidth;
74 138 : 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 138 : 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 138 : switch (cState->nRegisters)
89 : {
90 : case 16:
91 0 : alpha = 0.673;
92 0 : break;
93 : case 32:
94 0 : alpha = 0.697;
95 0 : break;
96 : case 64:
97 0 : alpha = 0.709;
98 0 : break;
99 : default:
100 138 : 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 138 : cState->alphaMM = alpha * cState->nRegisters * cState->nRegisters;
108 138 : }
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 : */
126 : void
127 0 : initHyperLogLogError(hyperLogLogState *cState, double error)
128 : {
129 0 : uint8 bwidth = 4;
130 :
131 0 : while (bwidth < 16)
132 : {
133 0 : double m = (Size) 1 << bwidth;
134 :
135 0 : if (1.04 / sqrt(m) < error)
136 0 : break;
137 0 : bwidth++;
138 : }
139 :
140 0 : initHyperLogLog(cState, bwidth);
141 0 : }
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 : */
149 : void
150 0 : freeHyperLogLog(hyperLogLogState *cState)
151 : {
152 0 : Assert(cState->hashesArr != NULL);
153 0 : pfree(cState->hashesArr);
154 0 : }
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 : */
165 : void
166 2219 : addHyperLogLog(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 2219 : index = hash >> (BITS_PER_BYTE * sizeof(uint32) - cState->registerWidth);
173 :
174 : /* Compute the rank of the remaining 32 - "k" (registerWidth) bits */
175 2219 : count = rho(hash << cState->registerWidth,
176 2219 : BITS_PER_BYTE * sizeof(uint32) - cState->registerWidth);
177 :
178 2219 : cState->hashesArr[index] = Max(count, cState->hashesArr[index]);
179 2219 : }
180 :
181 : /*
182 : * Estimates cardinality, based on elements added so far
183 : */
184 : double
185 6 : estimateHyperLogLog(hyperLogLogState *cState)
186 : {
187 : double result;
188 6 : double sum = 0.0;
189 : int i;
190 :
191 6150 : for (i = 0; i < cState->nRegisters; i++)
192 : {
193 6144 : sum += 1.0 / pow(2.0, cState->hashesArr[i]);
194 : }
195 :
196 : /* result set to "raw" HyperLogLog estimate (E in the HyperLogLog paper) */
197 6 : result = cState->alphaMM / sum;
198 :
199 6 : if (result <= (5.0 / 2.0) * cState->nRegisters)
200 : {
201 : /* Small range correction */
202 6 : int zero_count = 0;
203 :
204 6150 : for (i = 0; i < cState->nRegisters; i++)
205 : {
206 6144 : if (cState->hashesArr[i] == 0)
207 6127 : zero_count++;
208 : }
209 :
210 6 : if (zero_count != 0)
211 6 : result = cState->nRegisters * log((double) cState->nRegisters /
212 : zero_count);
213 : }
214 0 : else if (result > (1.0 / 30.0) * POW_2_32)
215 : {
216 : /* Large range correction */
217 0 : result = NEG_POW_2_32 * log(1.0 - (result / POW_2_32));
218 : }
219 :
220 6 : 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 : */
240 : static inline uint8
241 2219 : rho(uint32 x, uint8 b)
242 : {
243 2219 : uint8 j = 1;
244 :
245 5306 : while (j <= b && !(x & 0x80000000))
246 : {
247 868 : j++;
248 868 : x <<= 1;
249 : }
250 :
251 2219 : return j;
252 : }
|
