sketch_ams.nx source
↩ module page · 231 lines · 7695 B
1// sketch_ams.nx -- AMS sketch (Alon-Matias-Szegedy 1996).
2//
3// Estimates F_2 = Σ f_i² (second frequency moment) of a stream over
4// implicit-keyed items. Used for:
5// - self-join size estimation in databases
6// - query-plan cardinality estimation
7// - skew detection (high F_2 = skewed distribution)
8// - L2 norm of frequency vector
9//
10// ALGORITHM:
11// For each of d * s estimators, a random ±1 sign function ξ_jk.
12// On (item x, count c): counter[j][k] += ξ_jk(x) * c for all (j,k).
13// F_2 estimate per estimator: counter[j][k]²
14// Within-group AVERAGE: F_2_j = mean of counter[j][k]² over k.
15// Across-group MEDIAN: F_2 ≈ median(F_2_j) across j.
16//
17// Variance: average reduces variance by 1/s. Median over d
18// independent estimates boosts confidence to 1 - 2^(-d/2).
19//
20// MEMORY: d * s * 8 bytes. d=5, s=64 -> 2560 bytes for F_2 with
21// 12.5% relative error at 87.5% confidence.
22//
23// LOSSLESS-LANGUAGE DISCIPLINE: nx_ams_query returns ApproxI64 with
24// NX_ENV_REL_STDDEV = 1/sqrt(s) per estimator (further tightened by
25// median-of-d boost in practice). conf_ppb tracks (1 - 2^(-d/2)).
26
27import "syscalls.nx"
28import "sketch_types.nx"
29
30const NX_AMS_MIN_D: i64 = 3
31const NX_AMS_MAX_D: i64 = 32
32const NX_AMS_MIN_S: i64 = 4
33const NX_AMS_MAX_S: i64 = 1024
34
35struct AMS {
36 counters: *i64, // d * s
37 d: i64,
38 s: i64,
39 seed: i64,
40 total: i64,
41 ac: i64, // bits-up: precomputed A * C mod 2^64
42 bc: i64, // bits-up: precomputed B * C mod 2^64
43 // factors (key+j*A+k*B)*C as
44 // (key*C+seed) + j*ac + k*bc
45 // saves one mul per (j,k) cell in nx_ams_add.
46 scratch: *i64, // d-sized scratch hoisted from per-call sys_mmap.
47}
48
49// === construction =================================================
50
51func nx_ams_alloc(d: i64, s: i64, seed: i64) -> *AMS {
52 if d < NX_AMS_MIN_D { return 0 as *AMS }
53 if d > NX_AMS_MAX_D { return 0 as *AMS }
54 if s < NX_AMS_MIN_S { return 0 as *AMS }
55 if s > NX_AMS_MAX_S { return 0 as *AMS }
56 let raw: *u8 = sys_mmap(80)
57 let a: *AMS = raw as *AMS
58 let cells: i64 = d * s
59 a.counters = sys_mmap(cells * 8) as *i64
60 var i: i64 = 0
61 while i < cells {
62 a.counters[i] = 0
63 i = i + 1
64 }
65 a.d = d
66 a.s = s
67 a.seed = seed
68 a.total = 0
69 // Bits-up: precompute A*C and B*C (mod 2^64) so the per-cell sign
70 // computation drops from (key + j*A + k*B) * C + seed to
71 // (key*C + seed) + j*ac + k*bc -- one fewer mul per cell.
72 a.ac = (0x9E3779B9 * 0xC2B2AE3D27D4EB4F) & 0xFFFFFFFFFFFFFFFF
73 a.bc = (0xBF58476D1CE4E5B9 * 0xC2B2AE3D27D4EB4F) & 0xFFFFFFFFFFFFFFFF
74 a.scratch = sys_mmap(d * 8) as *i64
75 return a
76}
77
78// === sign hash function ===========================================
79//
80// For estimator (j, k) and item key, derive a deterministic ±1 sign.
81// Mix item key with row/col indices and seed; top bit -> sign.
82
83func nx_ams_sign(a: *AMS, key: i64, j: i64, k: i64) -> i64 {
84 // Bits-up: factored form using precomputed ac/bc. Mathematically
85 // identical to (key + j*A + k*B) * C + seed mod 2^64.
86 let mixed: i64 = ((key * 0xC2B2AE3D27D4EB4F + a.seed +
87 j * a.ac + k * a.bc)) & 0xFFFFFFFFFFFFFFFF
88 if (mixed & (1 << 63)) == 0 { return 1 }
89 return -1
90}
91
92func nx_ams_cell_idx(a: *AMS, j: i64, k: i64) -> i64 {
93 return j * a.s + k
94}
95
96// === add ==========================================================
97
98func nx_ams_add(a: *AMS, key: i64, count: i64) -> i64 {
99 if count == 0 { return 0 }
100 a.total = a.total + count
101 // Bits-up: precompute key*C + seed once. Per cell becomes
102 // key_mix + j*ac + k*bc (2 mul + 2 add, was 3 mul + 4 add).
103 let key_mix: i64 = (key * 0xC2B2AE3D27D4EB4F + a.seed) & 0xFFFFFFFFFFFFFFFF
104 var j: i64 = 0
105 while j < a.d {
106 let j_term: i64 = (j * a.ac) & 0xFFFFFFFFFFFFFFFF
107 let row_base: i64 = j * a.s
108 var k: i64 = 0
109 while k < a.s {
110 let mixed: i64 = (key_mix + j_term + k * a.bc) & 0xFFFFFFFFFFFFFFFF
111 var sign: i64 = 1
112 if (mixed & (1 << 63)) != 0 { sign = -1 }
113 let idx: i64 = row_base + k
114 a.counters[idx] = a.counters[idx] + sign * count
115 k = k + 1
116 }
117 j = j + 1
118 }
119 return 0
120}
121
122// === F_2 estimate =================================================
123//
124// Per estimator: F_2_jk = counter[j][k]²
125// Group average: F_2_j = (1/s) Σ_k counter[j][k]²
126// Median across j: F_2 = median(F_2_j)
127
128func nx_ams_isqrt(x: i64) -> i64 {
129 if x < 0 { return 0 }
130 if x == 0 { return 0 }
131 if x < 4 { return 1 }
132 var g: i64 = (x >> 1) + 1
133 var iter: i64 = 0
134 while iter < 64 {
135 let next_g: i64 = (g + x / g) / 2
136 if next_g >= g { iter = 64 }
137 if next_g < g {
138 g = next_g
139 iter = iter + 1
140 }
141 }
142 return g
143}
144
145func nx_ams_f2(a: *AMS) -> i64 {
146 // Bits-up: scratch hoisted to struct (was per-call sys_mmap).
147 let scratch: *i64 = a.scratch
148 var j: i64 = 0
149 while j < a.d {
150 var sum_sq: i64 = 0
151 var k: i64 = 0
152 while k < a.s {
153 let idx: i64 = nx_ams_cell_idx(a, j, k)
154 let c: i64 = a.counters[idx]
155 sum_sq = sum_sq + c * c
156 k = k + 1
157 }
158 scratch[j] = sum_sq / a.s
159 j = j + 1
160 }
161 // Insertion sort scratch[0..d).
162 var i: i64 = 1
163 while i < a.d {
164 let cur: i64 = scratch[i]
165 var p: i64 = i - 1
166 var done: i64 = 0
167 while done == 0 {
168 if p < 0 { done = 1 }
169 if done == 0 {
170 if scratch[p] <= cur { done = 1 }
171 if done == 0 {
172 scratch[p + 1] = scratch[p]
173 p = p - 1
174 }
175 }
176 }
177 scratch[p + 1] = cur
178 i = i + 1
179 }
180 // Median.
181 return scratch[a.d / 2]
182}
183
184// === typed envelope ===============================================
185
186func nx_ams_stderr_ppb(s: i64) -> i64 {
187 // Per-estimator relative stddev ~ 1/sqrt(s).
188 let isq: i64 = nx_ams_isqrt(s)
189 if isq == 0 { return 1000000000 }
190 return 1000000000 / isq
191}
192
193func nx_ams_conf_ppb(d: i64) -> i64 {
194 // (1 - 2^(-d/2)) confidence. d=5 -> 1 - 0.177 ~ 0.823.
195 if d <= 3 { return 500000000 }
196 if d <= 5 { return 823000000 }
197 if d <= 7 { return 875000000 }
198 if d <= 11 { return 968000000 }
199 return 992000000
200}
201
202func nx_ams_query_f2(a: *AMS) -> *ApproxI64 {
203 let f2: i64 = nx_ams_f2(a)
204 return nx_approx_new(f2, NX_ENV_REL_STDDEV, nx_ams_stderr_ppb(a.s),
205 nx_ams_conf_ppb(a.d),
206 NX_MATURITY_REFERENCE_IMPL,
207 NX_ADV_HONEST)
208}
209
210// === merge ========================================================
211//
212// Counters add element-wise. Matching d, s, seed required.
213
214func nx_ams_merge(a: *AMS, b: *AMS) -> *AMS {
215 if a.d != b.d { return 0 as *AMS }
216 if a.s != b.s { return 0 as *AMS }
217 if a.seed != b.seed { return 0 as *AMS }
218 let out: *AMS = nx_ams_alloc(a.d, a.s, a.seed)
219 let cells: i64 = a.d * a.s
220 var i: i64 = 0
221 while i < cells {
222 out.counters[i] = a.counters[i] + b.counters[i]
223 i = i + 1
224 }
225 out.total = a.total + b.total
226 return out
227}
228
229func nx_ams_memory_bytes(a: *AMS) -> i64 {
230 return 48 + a.d * a.s * 8
231}