sketch_ams.nx source
↩ module page · 217 lines · 7386 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"
29import "nx_vecmath.nx"
30
31const NX_AMS_MIN_D: i64 = 3
32const NX_AMS_MAX_D: i64 = 32
33const NX_AMS_MIN_S: i64 = 4
34const NX_AMS_MAX_S: i64 = 1024
35
36struct AMS {
37 counters: *i64, // d * s
38 d: i64,
39 s: i64,
40 seed: i64,
41 total: i64,
42 ac: i64, // bits-up: precomputed A * C mod 2^64
43 bc: i64, // bits-up: precomputed B * C mod 2^64
44 // factors (key+j*A+k*B)*C as
45 // (key*C+seed) + j*ac + k*bc
46 // saves one mul per (j,k) cell in nx_ams_add.
47 scratch: *i64, // d-sized scratch hoisted from per-call sys_mmap.
48}
49
50// === construction =================================================
51
52func nx_ams_alloc(d: i64, s: i64, seed: i64) -> *AMS {
53 if d < NX_AMS_MIN_D { return 0 as *AMS }
54 if d > NX_AMS_MAX_D { return 0 as *AMS }
55 if s < NX_AMS_MIN_S { return 0 as *AMS }
56 if s > NX_AMS_MAX_S { return 0 as *AMS }
57 let raw: *u8 = sys_mmap(80)
58 let a: *AMS = raw as *AMS
59 let cells: i64 = d * s
60 a.counters = sys_mmap(cells * 8) as *i64
61 var i: i64 = 0
62 while i < cells {
63 a.counters[i] = 0
64 i = i + 1
65 }
66 a.d = d
67 a.s = s
68 a.seed = seed
69 a.total = 0
70 // Bits-up: precompute A*C and B*C (mod 2^64) so the per-cell sign
71 // computation drops from (key + j*A + k*B) * C + seed to
72 // (key*C + seed) + j*ac + k*bc -- one fewer mul per cell.
73 a.ac = (0x9E3779B9 * 0xC2B2AE3D27D4EB4F) & 0xFFFFFFFFFFFFFFFF
74 a.bc = (0xBF58476D1CE4E5B9 * 0xC2B2AE3D27D4EB4F) & 0xFFFFFFFFFFFFFFFF
75 a.scratch = sys_mmap(d * 8) as *i64
76 return a
77}
78
79// === sign hash function ===========================================
80//
81// For estimator (j, k) and item key, derive a deterministic ±1 sign.
82// Mix item key with row/col indices and seed; top bit -> sign.
83
84func nx_ams_sign(a: *AMS, key: i64, j: i64, k: i64) -> i64 {
85 // Bits-up: factored form using precomputed ac/bc. Mathematically
86 // identical to (key + j*A + k*B) * C + seed mod 2^64.
87 let mixed: i64 = ((key * 0xC2B2AE3D27D4EB4F + a.seed +
88 j * a.ac + k * a.bc)) & 0xFFFFFFFFFFFFFFFF
89 if (mixed & (1 << 63)) == 0 { return 1 }
90 return -1
91}
92
93func nx_ams_cell_idx(a: *AMS, j: i64, k: i64) -> i64 {
94 return j * a.s + k
95}
96
97// === add ==========================================================
98
99func nx_ams_add(a: *AMS, key: i64, count: i64) -> i64 {
100 if count == 0 { return 0 }
101 a.total = a.total + count
102 // Bits-up: precompute key*C + seed once. Per cell becomes
103 // key_mix + j*ac + k*bc (2 mul + 2 add, was 3 mul + 4 add).
104 let key_mix: i64 = (key * 0xC2B2AE3D27D4EB4F + a.seed) & 0xFFFFFFFFFFFFFFFF
105 var j: i64 = 0
106 while j < a.d {
107 let j_term: i64 = (j * a.ac) & 0xFFFFFFFFFFFFFFFF
108 let row_base: i64 = j * a.s
109 var k: i64 = 0
110 while k < a.s {
111 let mixed: i64 = (key_mix + j_term + k * a.bc) & 0xFFFFFFFFFFFFFFFF
112 var sign: i64 = 1
113 if (mixed & (1 << 63)) != 0 { sign = -1 }
114 let idx: i64 = row_base + k
115 a.counters[idx] = a.counters[idx] + sign * count
116 k = k + 1
117 }
118 j = j + 1
119 }
120 return 0
121}
122
123// === F_2 estimate =================================================
124//
125// Per estimator: F_2_jk = counter[j][k]²
126// Group average: F_2_j = (1/s) Σ_k counter[j][k]²
127// Median across j: F_2 = median(F_2_j)
128
129func nx_ams_isqrt(x: i64) -> i64 { return vm_isqrt(x) }
130
131func nx_ams_f2(a: *AMS) -> i64 {
132 // Bits-up: scratch hoisted to struct (was per-call sys_mmap).
133 let scratch: *i64 = a.scratch
134 var j: i64 = 0
135 while j < a.d {
136 var sum_sq: i64 = 0
137 var k: i64 = 0
138 while k < a.s {
139 let idx: i64 = nx_ams_cell_idx(a, j, k)
140 let c: i64 = a.counters[idx]
141 sum_sq = sum_sq + c * c
142 k = k + 1
143 }
144 scratch[j] = sum_sq / a.s
145 j = j + 1
146 }
147 // Insertion sort scratch[0..d).
148 var i: i64 = 1
149 while i < a.d {
150 let cur: i64 = scratch[i]
151 var p: i64 = i - 1
152 var done: i64 = 0
153 while done == 0 {
154 if p < 0 { done = 1 }
155 if done == 0 {
156 if scratch[p] <= cur { done = 1 }
157 if done == 0 {
158 scratch[p + 1] = scratch[p]
159 p = p - 1
160 }
161 }
162 }
163 scratch[p + 1] = cur
164 i = i + 1
165 }
166 // Median.
167 return scratch[a.d / 2]
168}
169
170// === typed envelope ===============================================
171
172func nx_ams_stderr_ppb(s: i64) -> i64 {
173 // Per-estimator relative stddev ~ 1/sqrt(s).
174 let isq: i64 = nx_ams_isqrt(s)
175 if isq == 0 { return 1000000000 }
176 return 1000000000 / isq
177}
178
179func nx_ams_conf_ppb(d: i64) -> i64 {
180 // (1 - 2^(-d/2)) confidence. d=5 -> 1 - 0.177 ~ 0.823.
181 if d <= 3 { return 500000000 }
182 if d <= 5 { return 823000000 }
183 if d <= 7 { return 875000000 }
184 if d <= 11 { return 968000000 }
185 return 992000000
186}
187
188func nx_ams_query_f2(a: *AMS) -> *ApproxI64 {
189 let f2: i64 = nx_ams_f2(a)
190 return nx_approx_new(f2, NX_ENV_REL_STDDEV, nx_ams_stderr_ppb(a.s),
191 nx_ams_conf_ppb(a.d),
192 NX_MATURITY_REFERENCE_IMPL,
193 NX_ADV_HONEST)
194}
195
196// === merge ========================================================
197//
198// Counters add element-wise. Matching d, s, seed required.
199
200func nx_ams_merge(a: *AMS, b: *AMS) -> *AMS {
201 if a.d != b.d { return 0 as *AMS }
202 if a.s != b.s { return 0 as *AMS }
203 if a.seed != b.seed { return 0 as *AMS }
204 let out: *AMS = nx_ams_alloc(a.d, a.s, a.seed)
205 let cells: i64 = a.d * a.s
206 var i: i64 = 0
207 while i < cells {
208 out.counters[i] = a.counters[i] + b.counters[i]
209 i = i + 1
210 }
211 out.total = a.total + b.total
212 return out
213}
214
215func nx_ams_memory_bytes(a: *AMS) -> i64 {
216 return 48 + a.d * a.s * 8
217}