sketch_linear_counter.nx source
↩ module page · 223 lines · 7647 B
1// sketch_linear_counter.nx -- Linear Counting (Whang-Vander Zanden 1990).
2//
3// Pre-HLL cardinality primitive. Simpler than HLL; sometimes more
4// accurate at small cardinalities (n < m/10). Used by Google's
5// BigQuery as the small-cardinality fallback.
6//
7// ALGORITHM:
8// bitmap of m bits, all zero initially.
9// add(x): set bit at index = hash(x) % m
10// estimate(n) = -m * ln(zeros / m)
11//
12// where zeros = count of unset bits at query time. Derivation: if
13// items are hashed uniformly to m buckets and we observe `zeros`
14// empty buckets, then by occupancy theory the true count n satisfies
15// E[zeros] = m * (1 - 1/m)^n ≈ m * e^(-n/m)
16// inverting: n ≈ -m * ln(zeros / m).
17//
18// CAPABILITIES (vs HLL):
19// - LC: accurate at n < m (uses bitmap directly)
20// - HLL: accurate at n >> m (uses harmonic mean of register maxima)
21// Together: BigQuery + others run LC for n < threshold, HLL for n >=.
22//
23// INTEGER-FIXED-POINT IMPLEMENTATION:
24// ln(zeros/m) approximated via tabulated values keyed on (zeros*100)/m
25// (percentage of empty bits, in [0, 100]). Table values are
26// -ln(p/100) * 1_000_000 in PPM (negative of -log so positive).
27//
28// MEMORY: m bits = m/8 bytes. For m=8192: 1 KB.
29//
30// LOSSLESS-LANGUAGE DISCIPLINE:
31// nx_lc_query returns ApproxI64 with NX_ENV_REL_STDDEV. Standard
32// error per Whang 1990: sqrt(e^t - t - 1) / t where t = n/m (load).
33// Tabulated by load class.
34
35import "syscalls.nx"
36import "murmur3.nx"
37import "sketch_types.nx"
38
39const NX_LC_MIN_BITS: i64 = 64
40const NX_LC_MAX_BITS: i64 = 16777216 // 2 MB cap
41
42struct LinearCounter {
43 bits: *u8, // m / 8 bytes
44 m_bits: i64, // bitmap size (power of 2)
45 mask: i64, // m_bits - 1
46 seed: i64,
47}
48
49// === construction =================================================
50
51func nx_lc_is_pow2(n: i64) -> i64 {
52 if n < NX_LC_MIN_BITS { return 0 }
53 if (n & (n - 1)) != 0 { return 0 }
54 return 1
55}
56
57func nx_lc_alloc(m_bits: i64, seed: i64) -> *LinearCounter {
58 if nx_lc_is_pow2(m_bits) != 1 { return 0 as *LinearCounter }
59 if m_bits > NX_LC_MAX_BITS { return 0 as *LinearCounter }
60 let raw: *u8 = sys_mmap(40)
61 let lc: *LinearCounter = raw as *LinearCounter
62 let bytes: i64 = m_bits / 8
63 lc.bits = sys_mmap(bytes)
64 var i: i64 = 0
65 while i < bytes {
66 lc.bits[i] = 0
67 i = i + 1
68 }
69 lc.m_bits = m_bits
70 lc.mask = m_bits - 1
71 lc.seed = seed
72 return lc
73}
74
75// === bit helpers =================================================
76
77func nx_lc_set_bit(lc: *LinearCounter, bit_idx: i64) -> i64 {
78 let byte_idx: i64 = bit_idx >> 3
79 let bit_pos: i64 = bit_idx & 7
80 lc.bits[byte_idx] = lc.bits[byte_idx] | (1 << bit_pos)
81 return 0
82}
83
84func nx_lc_get_bit(lc: *LinearCounter, bit_idx: i64) -> i64 {
85 let byte_idx: i64 = bit_idx >> 3
86 let bit_pos: i64 = bit_idx & 7
87 return (lc.bits[byte_idx] >> bit_pos) & 1
88}
89
90// === add ==========================================================
91
92func nx_lc_add(lc: *LinearCounter, key: *u8, len: i64) -> i64 {
93 let h: i64 = murmur3_32(lc.seed, key, len) & 0xFFFFFFFF
94 let bit_idx: i64 = h & lc.mask
95 nx_lc_set_bit(lc, bit_idx)
96 return 0
97}
98
99// === zero count ==================================================
100//
101// Walks the bitmap counting unset bits. Worst-case linear in m/8.
102// Optimization: use popcount per byte then subtract from 8 * bytes.
103
104func nx_lc_popcount_byte(b: i64) -> i64 {
105 var n: i64 = 0
106 var t: i64 = b & 0xFF
107 while t > 0 {
108 n = n + (t & 1)
109 t = t >> 1
110 }
111 return n
112}
113
114func nx_lc_zeros(lc: *LinearCounter) -> i64 {
115 let bytes: i64 = lc.m_bits / 8
116 var set: i64 = 0
117 var i: i64 = 0
118 while i < bytes {
119 set = set + nx_lc_popcount_byte(lc.bits[i])
120 i = i + 1
121 }
122 return lc.m_bits - set
123}
124
125// === ln estimator ================================================
126//
127// Approximate -ln(zeros / m) via tabulation indexed by integer
128// percentage = (zeros * 100) / m. Table values are -ln(p/100) * 1e6.
129//
130// p=99: -ln(0.99) = 0.01005 -> 10050 PPM
131// p=90: -ln(0.90) = 0.10536 -> 105360 PPM
132// p=70: -ln(0.70) = 0.35667 -> 356670 PPM
133// p=50: -ln(0.50) = 0.69315 -> 693150 PPM
134// p=30: -ln(0.30) = 1.20397 -> 1203970 PPM
135// p=10: -ln(0.10) = 2.30259 -> 2302590 PPM
136// p=01: -ln(0.01) = 4.60517 -> 4605170 PPM
137// p=0.1: -ln(0.001)= 6.90776 -> 6907760 PPM
138//
139// For finer resolution we'd use a larger table. Coarse buckets are
140// fine: estimation error ~ 1% per bucket transition.
141
142func nx_lc_neg_ln_ppm(zeros_x_100: i64, m: i64) -> i64 {
143 let p: i64 = zeros_x_100 / m
144 if p >= 99 { return 10050 }
145 if p >= 95 { return 51290 }
146 if p >= 90 { return 105360 }
147 if p >= 80 { return 223140 }
148 if p >= 70 { return 356670 }
149 if p >= 60 { return 510830 }
150 if p >= 50 { return 693150 }
151 if p >= 40 { return 916290 }
152 if p >= 30 { return 1203970 }
153 if p >= 20 { return 1609440 }
154 if p >= 10 { return 2302590 }
155 if p >= 5 { return 2995730 }
156 if p >= 1 { return 4605170 }
157 return 6907760
158}
159
160// === estimate ====================================================
161//
162// estimate = -m * ln(zeros / m)
163//
164// If zeros == 0: bitmap saturated, estimate is much greater than m;
165// we return m as a conservative ceiling.
166
167func nx_lc_estimate(lc: *LinearCounter) -> i64 {
168 let zeros: i64 = nx_lc_zeros(lc)
169 if zeros == 0 { return lc.m_bits }
170 if zeros == lc.m_bits { return 0 }
171 let ln_ppm: i64 = nx_lc_neg_ln_ppm(zeros * 100, lc.m_bits)
172 // estimate = (m * ln_ppm) / 1e6
173 return (lc.m_bits * ln_ppm) / 1000000
174}
175
176// === typed envelope ==============================================
177//
178// Whang 1990 standard error: sqrt(e^t - t - 1) / t where t = n/m.
179// Tabulated by load class for the integer path.
180
181func nx_lc_stddev_rel_ppb(load_pct: i64) -> i64 {
182 if load_pct <= 1 { return 710000000 } // very small N, high relative variance
183 if load_pct <= 5 { return 226000000 }
184 if load_pct <= 10 { return 159000000 }
185 if load_pct <= 25 { return 95000000 }
186 if load_pct <= 50 { return 65000000 }
187 if load_pct <= 75 { return 53000000 }
188 if load_pct <= 90 { return 58000000 }
189 if load_pct <= 100 { return 80000000 }
190 return 300000000
191}
192
193func nx_lc_query(lc: *LinearCounter) -> *ApproxI64 {
194 let est: i64 = nx_lc_estimate(lc)
195 var load_pct: i64 = 0
196 if lc.m_bits > 0 { load_pct = (est * 100) / lc.m_bits }
197 return nx_approx_new(est, NX_ENV_REL_STDDEV,
198 nx_lc_stddev_rel_ppb(load_pct),
199 682700000,
200 NX_MATURITY_REFERENCE_IMPL,
201 NX_ADV_HONEST)
202}
203
204// === merge =======================================================
205//
206// Bit-wise OR. Both must have matching m and seed.
207
208func nx_lc_merge(a: *LinearCounter, b: *LinearCounter) -> *LinearCounter {
209 if a.m_bits != b.m_bits { return 0 as *LinearCounter }
210 if a.seed != b.seed { return 0 as *LinearCounter }
211 let out: *LinearCounter = nx_lc_alloc(a.m_bits, a.seed)
212 let bytes: i64 = a.m_bits / 8
213 var i: i64 = 0
214 while i < bytes {
215 out.bits[i] = a.bits[i] | b.bits[i]
216 i = i + 1
217 }
218 return out
219}
220
221func nx_lc_memory_bytes(lc: *LinearCounter) -> i64 {
222 return 40 + lc.m_bits / 8
223}