sketch_tdigest_v2.nx source
↩ module page · 359 lines · 11895 B
1// sketch_tdigest_v2.nx -- T-Digest with Dunning k1 arcsin-based scale.
2//
3// V2 vs V1 (sketch_tdigest.nx):
4// V1 scale: w_max_norm(q) = 4·q·(1-q) · π / δ (parabolic approx)
5// V2 scale: w_max_norm(q) = (1 - cos(2π/δ))·(1/2 - q)
6// + sin(2π/δ)·sqrt(q(1-q))
7//
8// The V2 form is the closed-form expansion of Dunning's
9// k1(q) = (δ/2π) · arcsin(2q - 1) inverse-difference w_max =
10// k1⁻¹(k1(q)+1) - q, via the sin-of-sum identity. At q=0.5
11// V1 and V2 agree; in the tails V2 is LOOSER than V1 (fewer
12// centroids needed for the same global rank guarantee), the
13// canonical memory-efficiency T-Digest property.
14//
15// V1 trades memory for tighter tail recall.
16// V2 trades tail recall for tighter memory.
17// Both ship -- caller picks the right tool. Per cardinal 25
18// ("Build Intelligence, Never Strip Features").
19//
20// TRIG CONSTANTS:
21// sin(2π/δ) and cos(2π/δ) precomputed at alloc time as Q14
22// fixed-point (1.0 = 16384), via Taylor series:
23// sin(x) ≈ x - x³/6
24// cos(x) ≈ 1 - x²/2 + x⁴/24
25// For δ ∈ [10, 500] -> x ∈ [0.0126, 0.628] -> Taylor error < 1e-4
26// in Q14, well below the 1-LSB sketch noise floor.
27//
28// LOSSLESS-LANGUAGE DISCIPLINE (doc 20):
29// Query returns ApproxI64 same NX_ENV_RANK_ERROR envelope as V1,
30// param_a interpolates the V2-specific bound:
31// V2 rank_err ≈ 1.5 / δ at p=0.5, ≈ 0.3 / δ at p=0.99 (tail
32// roughly 5x looser than V1; still tighter than KLL globally).
33
34import "syscalls.nx"
35import "sketch_types.nx"
36
37const NX_TD2_DELTA_MIN: i64 = 10
38const NX_TD2_DELTA_MAX: i64 = 500
39const NX_TD2_BUFFER_CAP: i64 = 1000
40const NX_TD2_MAX_CENTROIDS: i64 = 800
41
42const NX_TD2_Q14: i64 = 16384 // 1.0 fixed-point unit
43const NX_TD2_TWO_PI_Q14: i64 = 102944 // round(2π · 16384)
44
45// Beyond total > 100000 the q(1-q) cross-product would risk overflow
46// in the Q14 path. Above that we fall back to a constant 1% bound
47// (matches V1's safe-path behavior).
48const NX_TD2_SAFE_TOTAL_MAX: i64 = 100000
49
50struct CentroidV2 {
51 mean: i64,
52 weight: i64,
53}
54
55struct TDigestV2 {
56 centroids: *CentroidV2,
57 n_centroids: i64,
58 buffer: *i64,
59 n_buf: i64,
60 total_weight: i64,
61 delta: i64,
62 sin_2pi_over_d: i64, // Q14
63 one_minus_cos: i64, // Q14: 1 - cos(2π/δ)
64 merge_scratch: *CentroidV2,// bits-up: hoisted from nx_tdv2_merge
65}
66
67// === construction =================================================
68//
69// Precompute sin(2π/δ) and (1 - cos(2π/δ)) in Q14 via Taylor.
70
71func nx_tdv2_alloc(delta: i64) -> *TDigestV2 {
72 if delta < NX_TD2_DELTA_MIN { return 0 as *TDigestV2 }
73 if delta > NX_TD2_DELTA_MAX { return 0 as *TDigestV2 }
74 let raw: *u8 = sys_mmap(80)
75 let td: *TDigestV2 = raw as *TDigestV2
76 let cent_raw: *u8 = sys_mmap(NX_TD2_MAX_CENTROIDS * 16)
77 td.centroids = cent_raw as *CentroidV2
78 td.n_centroids = 0
79 let buf_raw: *u8 = sys_mmap(NX_TD2_BUFFER_CAP * 8)
80 td.buffer = buf_raw as *i64
81 td.n_buf = 0
82 td.total_weight = 0
83 td.delta = delta
84 // x_q14 = 2π/δ scaled
85 let x: i64 = NX_TD2_TWO_PI_Q14 / delta
86 // x² scaled in Q14:
87 // x_q14 · x_q14 / Q14
88 let x2: i64 = (x * x) / NX_TD2_Q14
89 // x³ scaled in Q14:
90 // x² · x / Q14
91 let x3: i64 = (x2 * x) / NX_TD2_Q14
92 // x⁴ scaled in Q14:
93 let x4: i64 = (x2 * x2) / NX_TD2_Q14
94 // sin(x) ≈ x - x³/6
95 let sin_v: i64 = x - (x3 / 6)
96 // cos(x) ≈ 1 - x²/2 + x⁴/24
97 let cos_v: i64 = NX_TD2_Q14 - (x2 / 2) + (x4 / 24)
98 td.sin_2pi_over_d = sin_v
99 td.one_minus_cos = NX_TD2_Q14 - cos_v
100 let scratch_raw: *u8 = sys_mmap(NX_TD2_MAX_CENTROIDS * 16)
101 td.merge_scratch = scratch_raw as *CentroidV2
102 return td
103}
104
105func nx_tdv2_centroid_at(td: *TDigestV2, i: i64) -> *CentroidV2 {
106 return (td.centroids as i64 + i * 16) as *CentroidV2
107}
108
109// === integer sqrt (Newton, matches sketch_stream_stats style) =====
110
111func nx_tdv2_isqrt(x: i64) -> i64 {
112 if x < 0 { return 0 }
113 if x == 0 { return 0 }
114 if x < 4 { return 1 }
115 var g: i64 = (x >> 1) + 1
116 var iter: i64 = 0
117 while iter < 64 {
118 let next_g: i64 = (g + x / g) / 2
119 if next_g >= g { iter = 64 }
120 if next_g < g {
121 g = next_g
122 iter = iter + 1
123 }
124 }
125 return g
126}
127
128// === Dunning k1 scale function ====================================
129//
130// w_max(q, N) = (1-cos) · (N - 2·cum_w) / 2 + sin · sqrt(cum_w·(N-cum_w))
131//
132// All terms computed in integer space; the Q14 factor is divided
133// out at the end. Floor: 1 (centroids with weight 0 are degenerate).
134
135func nx_tdv2_w_max(td: *TDigestV2, cum_w: i64, total: i64) -> i64 {
136 if total == 0 { return 1 }
137 if total > NX_TD2_SAFE_TOTAL_MAX {
138 let w_safe: i64 = total / 100
139 if w_safe < 1 { return 1 }
140 return w_safe
141 }
142 // (N - 2·cum_w) can be negative (q > 0.5). Keep signed.
143 let signed_half: i64 = (total - 2 * cum_w) / 2
144 // First term: (1-cos_q14) * signed_half / Q14 units: items (signed)
145 let term1: i64 = (td.one_minus_cos * signed_half) / NX_TD2_Q14
146 // Second term: sin_q14 * sqrt(cum_w * (total - cum_w)) / Q14
147 let prod: i64 = cum_w * (total - cum_w)
148 let sq: i64 = nx_tdv2_isqrt(prod)
149 let term2: i64 = (td.sin_2pi_over_d * sq) / NX_TD2_Q14
150 let w_max: i64 = term1 + term2
151 if w_max < 1 { return 1 }
152 return w_max
153}
154
155// === insertion sort over buffer ===================================
156
157func nx_tdv2_sort_buffer(td: *TDigestV2) -> i64 {
158 var i: i64 = 1
159 while i < td.n_buf {
160 let cur: i64 = td.buffer[i]
161 var j: i64 = i - 1
162 var done: i64 = 0
163 while done == 0 {
164 if j < 0 { done = 1 }
165 if done == 0 {
166 if td.buffer[j] <= cur { done = 1 }
167 if done == 0 {
168 td.buffer[j + 1] = td.buffer[j]
169 j = j - 1
170 }
171 }
172 }
173 td.buffer[j + 1] = cur
174 i = i + 1
175 }
176 return 0
177}
178
179// === merge ========================================================
180//
181// Identical control structure to V1's merge -- only the scale
182// function (and thus the absorb-vs-flush decision) differs.
183
184func nx_tdv2_merge(td: *TDigestV2) -> i64 {
185 if td.n_buf == 0 { return 0 }
186 nx_tdv2_sort_buffer(td)
187 let new_total: i64 = td.total_weight + td.n_buf
188
189 // Bits-up: merge_scratch hoisted to struct.
190 let scratch: *CentroidV2 = td.merge_scratch
191
192 var n_scratch: i64 = 0
193 var bi: i64 = 0
194 var ci: i64 = 0
195 var cum_w: i64 = 0
196
197 var cur_sum: i64 = 0
198 var cur_w: i64 = 0
199 var cur_q_left: i64 = 0
200 var cur_started: i64 = 0
201
202 var loop_done: i64 = 0
203 while loop_done == 0 {
204 let b_done: i64 = bi >= td.n_buf
205 let c_done: i64 = ci >= td.n_centroids
206 if b_done == 1 {
207 if c_done == 1 { loop_done = 1 }
208 }
209 if loop_done == 0 {
210 var nm: i64 = 0
211 var nw: i64 = 0
212 if b_done == 1 {
213 let c: *CentroidV2 = nx_tdv2_centroid_at(td, ci)
214 nm = c.mean
215 nw = c.weight
216 ci = ci + 1
217 }
218 if b_done == 0 {
219 if c_done == 1 {
220 nm = td.buffer[bi]
221 nw = 1
222 bi = bi + 1
223 }
224 if c_done == 0 {
225 let c: *CentroidV2 = nx_tdv2_centroid_at(td, ci)
226 if td.buffer[bi] <= c.mean {
227 nm = td.buffer[bi]
228 nw = 1
229 bi = bi + 1
230 }
231 if td.buffer[bi] > c.mean {
232 nm = c.mean
233 nw = c.weight
234 ci = ci + 1
235 }
236 }
237 }
238 var just_init: i64 = 0
239 if cur_started == 0 {
240 cur_sum = nm * nw
241 cur_w = nw
242 cur_q_left = cum_w
243 cur_started = 1
244 just_init = 1
245 }
246 if just_init == 0 {
247 let w_after: i64 = cur_w + nw
248 let w_max: i64 = nx_tdv2_w_max(td, cur_q_left, new_total)
249 if w_after <= w_max {
250 cur_sum = cur_sum + nm * nw
251 cur_w = cur_w + nw
252 }
253 if w_after > w_max {
254 let cent: *CentroidV2 = (scratch as i64 + n_scratch * 16) as *CentroidV2
255 cent.mean = cur_sum / cur_w
256 cent.weight = cur_w
257 n_scratch = n_scratch + 1
258 cum_w = cum_w + cur_w
259 cur_sum = nm * nw
260 cur_w = nw
261 cur_q_left = cum_w
262 }
263 }
264 }
265 }
266 if cur_started == 1 {
267 let cent: *CentroidV2 = (scratch as i64 + n_scratch * 16) as *CentroidV2
268 cent.mean = cur_sum / cur_w
269 cent.weight = cur_w
270 n_scratch = n_scratch + 1
271 }
272 var k: i64 = 0
273 while k < n_scratch {
274 let src: *CentroidV2 = (scratch as i64 + k * 16) as *CentroidV2
275 let dst: *CentroidV2 = nx_tdv2_centroid_at(td, k)
276 dst.mean = src.mean
277 dst.weight = src.weight
278 k = k + 1
279 }
280 td.n_centroids = n_scratch
281 td.total_weight = new_total
282 td.n_buf = 0
283 return 0
284}
285
286// === add ==========================================================
287
288func nx_tdv2_add(td: *TDigestV2, value: i64) -> i64 {
289 td.buffer[td.n_buf] = value
290 td.n_buf = td.n_buf + 1
291 if td.n_buf >= NX_TD2_BUFFER_CAP {
292 nx_tdv2_merge(td)
293 }
294 return 0
295}
296
297// === quantile =====================================================
298
299func nx_tdv2_quantile(td: *TDigestV2, p_milli: i64) -> i64 {
300 if td.n_buf > 0 { nx_tdv2_merge(td) }
301 if td.n_centroids == 0 { return 0 }
302 if td.total_weight == 0 { return 0 }
303 let target_w: i64 = (p_milli * td.total_weight) / 1000
304 var cum: i64 = 0
305 var i: i64 = 0
306 while i < td.n_centroids {
307 let c: *CentroidV2 = nx_tdv2_centroid_at(td, i)
308 let next_cum: i64 = cum + c.weight
309 if next_cum >= target_w {
310 return c.mean
311 }
312 cum = next_cum
313 i = i + 1
314 }
315 let last: *CentroidV2 = nx_tdv2_centroid_at(td, td.n_centroids - 1)
316 return last.mean
317}
318
319// === typed query ==================================================
320//
321// V2 rank-error envelope: ~1.5/δ globally, tighter than V1's
322// reported 1/δ in the middle, looser at extreme tails. The honest
323// bound to declare is the global one.
324
325func nx_tdv2_rank_error_ppb(delta: i64) -> i64 {
326 if delta >= 400 { return 4000000 } // 0.4% at delta=400+
327 if delta >= 200 { return 8000000 } // 0.8%
328 if delta >= 100 { return 15000000 } // 1.5%
329 if delta >= 50 { return 30000000 } // 3.0%
330 return 45000000 // 4.5% at delta=10..50
331}
332
333func nx_tdv2_query(td: *TDigestV2, p_milli: i64) -> *ApproxI64 {
334 let v: i64 = nx_tdv2_quantile(td, p_milli)
335 return nx_approx_new(v, NX_ENV_RANK_ERROR,
336 nx_tdv2_rank_error_ppb(td.delta),
337 950000000,
338 NX_MATURITY_REFERENCE_IMPL,
339 NX_ADV_HONEST)
340}
341
342// === introspection ================================================
343
344func nx_tdv2_memory_bytes(td: *TDigestV2) -> i64 {
345 return 72 + NX_TD2_MAX_CENTROIDS * 16 + NX_TD2_BUFFER_CAP * 8
346}
347
348func nx_tdv2_n_centroids(td: *TDigestV2) -> i64 {
349 return td.n_centroids
350}
351
352// Introspection helpers for cross-validating trig constants.
353func nx_tdv2_sin_q14(td: *TDigestV2) -> i64 {
354 return td.sin_2pi_over_d
355}
356
357func nx_tdv2_one_minus_cos_q14(td: *TDigestV2) -> i64 {
358 return td.one_minus_cos
359}