nx_sketch_tdigest.nx source
↩ module page · 340 lines · 12025 B
1// sketch_tdigest.nx -- T-Digest (Dunning 2019) tail-tight quantile sketch.
2//
3// Maintain centroids (mean, weight) sorted by mean. Each centroid covers
4// a window of the cumulative-distribution domain [q_left, q_right]. The
5// MAX WEIGHT a centroid may absorb depends on its q-position via a scale
6// function that gives more room to centroids in the middle of the
7// distribution and TIGHTLY BOUNDS centroids near the tails -- the
8// headline T-Digest property.
9//
10// COMPLEMENTS KLL:
11// - KLL: provable uniform rank-error bound, lighter at middle
12// quantiles, looser at tails.
13// - T-Digest: empirically tail-tight (p99/p99.9 accuracy is the
14// headline), no global rank-error guarantee.
15// Together: both shipping means callers pick the right tool per query.
16//
17// DUNNING k1 vs SIMPLIFIED SCALE FUNCTION (v1 trade-off):
18// Dunning's k1(q) = (delta / 2π) · arcsin(2q - 1)
19// → w_max_norm(q) = k1_inv(k1(q) + 1) - q
20// exact: tail w_max ~ sqrt-like; middle ~ delta-bounded.
21// v1 here: w_max_norm(q) = 4·q·(1-q) · π / delta
22// The simplified scale UNDERESTIMATES tail capacity (more conservative
23// = MORE centroids near tails, better tail accuracy than spec). At
24// q=0.5 it gives 0.0314·N matching Dunning within 1%. Trade-off:
25// memory uses slightly more centroids than necessary; queried accuracy
26// is no worse than Dunning's bound, often better at tails.
27// v2 will swap in a tabulated arcsin-based scale.
28//
29// LOSSLESS-LANGUAGE DISCIPLINE (doc 20):
30// Query returns ApproxI64 with envelope_kind = NX_ENV_RANK_ERROR,
31// param_a = ~1% absolute rank error (conservative bound for delta=100),
32// maturity = ReferenceImpl, adv = Honest.
33
34// nx_safety_envelope:
35// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
36// sil_target: SIL1
37// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
38// verdict: NOT_YET_EVALUATED
39
40import "nx_syscalls.nx"
41import "nx_sketch_types.nx"
42
43const NX_TD_DELTA_DEFAULT: i64 = 100
44const NX_TD_BUFFER_CAP: i64 = 1000
45const NX_TD_MAX_CENTROIDS: i64 = 800 // delta * 8 headroom
46
47const NX_TD_PI_PPB: i64 = 3141592654 // π · 10^9
48
49// Centroid = (mean, weight) packed as 16 bytes.
50struct Centroid {
51 mean: i64,
52 weight: i64,
53}
54
55struct TDigest {
56 centroids: *Centroid,
57 n_centroids: i64,
58 buffer: *i64,
59 n_buf: i64,
60 total_weight: i64,
61 delta: i64,
62}
63
64// === construction =================================================
65
66func nx_tdigest_alloc(delta: i64) -> *TDigest {
67 if delta < 10 { return 0 as *TDigest }
68 if delta > 500 { return 0 as *TDigest }
69 let raw: *u8 = sys_mmap(48)
70 let td: *TDigest = raw as *TDigest
71 let cent_raw: *u8 = sys_mmap(NX_TD_MAX_CENTROIDS * 16)
72 td.centroids = cent_raw as *Centroid
73 td.n_centroids = 0
74 let buf_raw: *u8 = sys_mmap(NX_TD_BUFFER_CAP * 8)
75 td.buffer = buf_raw as *i64
76 td.n_buf = 0
77 td.total_weight = 0
78 td.delta = delta
79 return td
80}
81
82func nx_tdigest_centroid_at(td: *TDigest, i: i64) -> *Centroid {
83 return (td.centroids as i64 + i * 16) as *Centroid
84}
85
86// === scale function (simplified) ==================================
87//
88// w_max(q, N) returns the max weight a centroid covering q_left=q
89// may absorb in a sketch of total weight N.
90//
91// w_max_norm(q) = 4·q·(1-q) · π / delta (in ppb)
92// w_max(q, N) = w_max_norm(q) · N / 1e9
93//
94// Floor: 1 (else degenerate single-point centroids would multiply).
95
96func nx_tdigest_w_max(td: *TDigest, cum_w: i64, total: i64) -> i64 {
97 if total == 0 { return 1 }
98 // q in [0, 1] represented as (cum_w / total).
99 // 4q(1-q) = 4 · cum_w · (total - cum_w) / total^2
100 // Scaled to ppb: ((4 · cum_w · (total - cum_w)) · π_ppb) / (delta · total^2)
101 // For overflow safety with int64 arithmetic:
102 // max cum_w · (total-cum_w) ~ total^2/4
103 // then · π_ppb (3.14e9) could overflow at total^2/4 · 3.14e9 > 2^63 ≈ 9.2e18
104 // so total^2 < 1.17e10 → total < 108k. For larger streams we'd
105 // need a different fixed-point layout. At total > 100k we use
106 // a degraded but safe path: w_max_norm fixed at 0.01·N (1%).
107 if total > 100000 {
108 // Safe path: 1% of N per centroid.
109 let w_safe: i64 = total / 100
110 if w_safe < 1 { return 1 }
111 return w_safe
112 }
113 let qmul: i64 = cum_w * (total - cum_w)
114 if qmul == 0 { return 1 }
115 let num: i64 = 4 * qmul * NX_TD_PI_PPB
116 let denom: i64 = td.delta * total * total
117 if denom == 0 { return 1 }
118 let w_norm_ppb: i64 = num / denom
119 let w_max: i64 = (w_norm_ppb * total) / 1000000000
120 if w_max < 1 { return 1 }
121 return w_max
122}
123
124// === insertion sort over an i64 array =============================
125
126func nx_tdigest_sort_buffer(td: *TDigest) -> i64 {
127 var i: i64 = 1
128 while i < td.n_buf {
129 let cur: i64 = td.buffer[i]
130 var j: i64 = i - 1
131 var done: i64 = 0
132 while done == 0 {
133 if j < 0 { done = 1 }
134 if done == 0 {
135 if td.buffer[j] <= cur { done = 1 }
136 if done == 0 {
137 td.buffer[j + 1] = td.buffer[j]
138 j = j - 1
139 }
140 }
141 }
142 td.buffer[j + 1] = cur
143 i = i + 1
144 }
145 return 0
146}
147
148// === merge buffer + centroids =====================================
149//
150// Sort buffer; merge into the centroid stream. Walk both sorted
151// inputs in tandem, accumulating weight into a "current centroid"
152// until its weight would exceed w_max(q_left); then commit it and
153// start a new centroid. Output replaces the existing centroid list.
154
155func nx_tdigest_merge(td: *TDigest) -> i64 {
156 if td.n_buf == 0 { return 0 }
157 nx_tdigest_sort_buffer(td)
158 let new_total: i64 = td.total_weight + td.n_buf
159
160 // Scratch space for new centroids (output ordering changes as
161 // small centroids are merged into larger ones at the middle).
162 let scratch_raw: *u8 = sys_mmap(NX_TD_MAX_CENTROIDS * 16)
163 let scratch: *Centroid = scratch_raw as *Centroid
164
165 var n_scratch: i64 = 0
166 var bi: i64 = 0
167 var ci: i64 = 0
168 var cum_w: i64 = 0
169
170 // Current centroid being built. cur_started=0 means empty.
171 var cur_sum: i64 = 0
172 var cur_w: i64 = 0
173 var cur_q_left: i64 = 0
174 var cur_started: i64 = 0
175
176 var loop_done: i64 = 0
177 while loop_done == 0 {
178 let b_done: i64 = bi >= td.n_buf
179 let c_done: i64 = ci >= td.n_centroids
180 if b_done == 1 {
181 if c_done == 1 { loop_done = 1 }
182 }
183 if loop_done == 0 {
184 // Pick next (mean, weight) from the smaller of the two
185 // sorted heads.
186 var nm: i64 = 0
187 var nw: i64 = 0
188 if b_done == 1 {
189 let c: *Centroid = nx_tdigest_centroid_at(td, ci)
190 nm = c.mean
191 nw = c.weight
192 ci = ci + 1
193 }
194 if b_done == 0 {
195 if c_done == 1 {
196 nm = td.buffer[bi]
197 nw = 1
198 bi = bi + 1
199 }
200 if c_done == 0 {
201 let c: *Centroid = nx_tdigest_centroid_at(td, ci)
202 if td.buffer[bi] <= c.mean {
203 nm = td.buffer[bi]
204 nw = 1
205 bi = bi + 1
206 }
207 if td.buffer[bi] > c.mean {
208 nm = c.mean
209 nw = c.weight
210 ci = ci + 1
211 }
212 }
213 }
214
215 // Process the pulled item.
216 var just_init: i64 = 0
217 if cur_started == 0 {
218 cur_sum = nm * nw
219 cur_w = nw
220 cur_q_left = cum_w
221 cur_started = 1
222 just_init = 1
223 }
224 if just_init == 0 {
225 // Centroid was already running; decide absorb vs flush.
226 let w_after: i64 = cur_w + nw
227 let w_max: i64 = nx_tdigest_w_max(td, cur_q_left, new_total)
228 if w_after <= w_max {
229 cur_sum = cur_sum + nm * nw
230 cur_w = cur_w + nw
231 }
232 if w_after > w_max {
233 // Flush current centroid to scratch.
234 let cent: *Centroid = (scratch as i64 + n_scratch * 16) as *Centroid
235 cent.mean = cur_sum / cur_w
236 cent.weight = cur_w
237 n_scratch = n_scratch + 1
238 cum_w = cum_w + cur_w
239 // Start a new centroid with the pulled item.
240 cur_sum = nm * nw
241 cur_w = nw
242 cur_q_left = cum_w
243 }
244 }
245 }
246 }
247 // Flush trailing centroid.
248 if cur_started == 1 {
249 let cent: *Centroid = (scratch as i64 + n_scratch * 16) as *Centroid
250 cent.mean = cur_sum / cur_w
251 cent.weight = cur_w
252 n_scratch = n_scratch + 1
253 }
254 // Copy scratch back into td.centroids.
255 var k: i64 = 0
256 while k < n_scratch {
257 let src: *Centroid = (scratch as i64 + k * 16) as *Centroid
258 let dst: *Centroid = nx_tdigest_centroid_at(td, k)
259 dst.mean = src.mean
260 dst.weight = src.weight
261 k = k + 1
262 }
263 td.n_centroids = n_scratch
264 td.total_weight = new_total
265 td.n_buf = 0
266 return 0
267}
268
269// === add ==========================================================
270
271func nx_tdigest_add(td: *TDigest, value: i64) -> i64 {
272 td.buffer[td.n_buf] = value
273 td.n_buf = td.n_buf + 1
274 if td.n_buf >= NX_TD_BUFFER_CAP {
275 nx_tdigest_merge(td)
276 }
277 return 0
278}
279
280// === quantile =====================================================
281//
282// p_milli: 0..1000. Walk centroids accumulating weight; locate
283// the cumulative-weight target p · total / 1000; return the
284// centroid mean covering that point.
285
286func nx_tdigest_quantile(td: *TDigest, p_milli: i64) -> i64 {
287 // Always flush before querying so buffer is incorporated.
288 if td.n_buf > 0 { nx_tdigest_merge(td) }
289 if td.n_centroids == 0 { return 0 }
290 if td.total_weight == 0 { return 0 }
291 let target_w: i64 = (p_milli * td.total_weight) / 1000
292 var cum: i64 = 0
293 var i: i64 = 0
294 while i < td.n_centroids {
295 let c: *Centroid = nx_tdigest_centroid_at(td, i)
296 let next_cum: i64 = cum + c.weight
297 if next_cum >= target_w {
298 return c.mean
299 }
300 cum = next_cum
301 i = i + 1
302 }
303 // Past end: return last centroid mean.
304 let last: *Centroid = nx_tdigest_centroid_at(td, td.n_centroids - 1)
305 return last.mean
306}
307
308// === typed query ==================================================
309//
310// Rank-error envelope. For delta=100, expected absolute rank error
311// is ~1% at p=0.5, tightening to ~0.1% at p=0.99 -- the tail-tight
312// property. We declare the GLOBAL bound (~1%) as the envelope, but
313// callers querying tail quantiles get much better accuracy in
314// practice.
315
316func nx_tdigest_rank_error_ppb(delta: i64) -> i64 {
317 if delta >= 200 { return 5000000 } // 0.5%
318 if delta >= 100 { return 10000000 } // 1.0%
319 if delta >= 50 { return 20000000 } // 2.0%
320 return 30000000 // 3.0% for very-low delta
321}
322
323func nx_tdigest_query(td: *TDigest, p_milli: i64) -> *ApproxI64 {
324 let v: i64 = nx_tdigest_quantile(td, p_milli)
325 return nx_approx_new(v, NX_ENV_RANK_ERROR,
326 nx_tdigest_rank_error_ppb(td.delta),
327 950000000,
328 NX_MATURITY_REFERENCE_IMPL,
329 NX_ADV_HONEST)
330}
331
332// === introspection ================================================
333
334func nx_tdigest_memory_bytes(td: *TDigest) -> i64 {
335 return 48 + NX_TD_MAX_CENTROIDS * 16 + NX_TD_BUFFER_CAP * 8
336}
337
338func nx_tdigest_n_centroids(td: *TDigest) -> i64 {
339 return td.n_centroids
340}