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