sketch_ddsketch.nx source
↩ module page · 264 lines · 8931 B
1// sketch_ddsketch.nx -- DDSketch (Masson et al, PVLDB 2019).
2//
3// Multiplicative-error quantile sketch. Different error model than
4// KLL (additive rank) and T-Digest (scale-function tail-tight):
5//
6// For query rank q, returned value v satisfies:
7// true_value(q) ∈ [v · (1-α), v · (1+α)]
8//
9// I.e., the returned VALUE is within α-relative-error of the true.
10// Perfect for log-distributed metrics (latency, sizes, durations).
11//
12// CAPABILITY: heavy-tail observability. Latency p99=100ms vs p99=1000ms
13// both reported with ±1% multiplicative error, regardless of magnitude.
14//
15// ALGORITHM:
16// bucket(x) = ⌈log_γ(x)⌉ where γ = (1+α)/(1-α)
17// For our integer implementation:
18// octave = bitlen(x) - 1 (floor log2)
19// mantissa = x - (1 << octave)
20// sub_bucket = mantissa * BPO / (1 << octave)
21// bucket_index = octave * BPO + sub_bucket
22//
23// BPO = buckets-per-octave; relative error ≈ 1/BPO.
24// BPO=64 -> ~1.5% multiplicative error
25// BPO=128 -> ~0.8%
26//
27// MEMORY: 64 * BPO * 8 bytes (~32KB for BPO=64, 64KB for BPO=128).
28// MERGE: count-wise addition; matching BPO required.
29//
30// LOSSLESS-LANGUAGE DISCIPLINE:
31// nx_dd_query_quantile returns NX_ENV_REL_STDDEV with param_a = 1/BPO
32// in PPB. Conf = 1e9 (deterministic). MaturityClass = ReferenceImpl.
33
34import "syscalls.nx"
35import "sketch_types.nx"
36
37const NX_DD_BPO_MIN: i64 = 16
38const NX_DD_BPO_MAX: i64 = 1024
39const NX_DD_MAX_OCTAVE: i64 = 62 // covers x up to 2^62
40
41struct DDSketch {
42 buckets: *i64, // size = (MAX_OCTAVE+1) * BPO
43 bpo: i64,
44 total: i64,
45 zero_count: i64, // x == 0 lands here
46 neg_count: i64, // x < 0 lands here (raw count; no neg buckets v1)
47 min_bucket: i64, // lowest populated bucket index, -1 = empty.
48 // bits-up path optimization 2026-05-20: lets
49 // nx_dd_quantile skip leading empty buckets.
50 max_bucket: i64, // highest populated bucket index, -1 = empty.
51 // lets quantile early-exit instead of walking
52 // all 63*BPO buckets.
53 lg_bpo: i64, // log2(bpo) if bpo is power of 2, else -1.
54 // bits-up math optimization: enables shift
55 // fast-path in nx_dd_bucket (saves a 64-bit
56 // divide per add on common bpo=64/128/256).
57}
58
59// === construction =================================================
60
61func nx_dd_alloc(bpo: i64) -> *DDSketch {
62 if bpo < NX_DD_BPO_MIN { return 0 as *DDSketch }
63 if bpo > NX_DD_BPO_MAX { return 0 as *DDSketch }
64 let raw: *u8 = sys_mmap(72)
65 let d: *DDSketch = raw as *DDSketch
66 let n_buckets: i64 = (NX_DD_MAX_OCTAVE + 1) * bpo
67 let buckets_raw: *u8 = sys_mmap(n_buckets * 8)
68 d.buckets = buckets_raw as *i64
69 var i: i64 = 0
70 while i < n_buckets {
71 d.buckets[i] = 0
72 i = i + 1
73 }
74 d.bpo = bpo
75 d.total = 0
76 d.zero_count = 0
77 d.neg_count = 0
78 d.min_bucket = -1
79 d.max_bucket = -1
80 // Bits-up: detect power-of-2 bpo for shift fast-path.
81 if (bpo & (bpo - 1)) == 0 {
82 var t: i64 = bpo
83 var lg: i64 = 0
84 while t > 1 {
85 t = t >> 1
86 lg = lg + 1
87 }
88 d.lg_bpo = lg
89 }
90 if (bpo & (bpo - 1)) != 0 {
91 d.lg_bpo = -1
92 }
93 return d
94}
95
96// === bit-length helper ============================================
97
98func nx_dd_bitlen(x: i64) -> i64 {
99 if x <= 0 { return 0 }
100 var n: i64 = 0
101 var t: i64 = x
102 while t > 0 {
103 t = t >> 1
104 n = n + 1
105 }
106 return n
107}
108
109// === bucket index =================================================
110//
111// For x > 0:
112// octave = bitlen(x) - 1
113// mantissa = x - (1 << octave)
114// sub_bucket = mantissa * bpo / (1 << octave)
115// bucket = octave * bpo + sub_bucket
116
117func nx_dd_bucket(d: *DDSketch, x: i64) -> i64 {
118 if x <= 0 { return 0 } // caller responsibility to not pass <= 0
119 let lg: i64 = nx_dd_bitlen(x)
120 let octave: i64 = lg - 1
121 if octave > NX_DD_MAX_OCTAVE { return (NX_DD_MAX_OCTAVE + 1) * d.bpo - 1 }
122 let base: i64 = 1 << octave
123 let mantissa: i64 = x - base
124 // Bits-up math: when bpo is power of 2, replace
125 // (mantissa * bpo) / base (one 64-bit divide)
126 // with
127 // (mantissa << lg_bpo) >> octave (two shifts)
128 var sub: i64 = 0
129 if d.lg_bpo >= 0 {
130 sub = (mantissa << d.lg_bpo) >> octave
131 }
132 if d.lg_bpo < 0 {
133 sub = (mantissa * d.bpo) / base
134 }
135 if sub < 0 { sub = 0 }
136 if sub >= d.bpo { sub = d.bpo - 1 }
137 return octave * d.bpo + sub
138}
139
140// === inverse: bucket -> representative value ======================
141//
142// For bucket index b in octave o, sub s:
143// o = b / bpo; s = b % bpo
144// value = (1 << o) + s * (1 << o) / bpo
145// = (1 << o) * (bpo + s) / bpo
146
147func nx_dd_value_at_bucket(d: *DDSketch, b: i64) -> i64 {
148 let octave: i64 = b / d.bpo
149 let sub: i64 = b - octave * d.bpo
150 let base: i64 = 1 << octave
151 return (base * (d.bpo + sub)) / d.bpo
152}
153
154// === add ==========================================================
155
156func nx_dd_add(d: *DDSketch, x: i64) -> i64 {
157 d.total = d.total + 1
158 if x == 0 {
159 d.zero_count = d.zero_count + 1
160 return 0
161 }
162 if x < 0 {
163 d.neg_count = d.neg_count + 1
164 return 0
165 }
166 let b: i64 = nx_dd_bucket(d, x)
167 d.buckets[b] = d.buckets[b] + 1
168 // Bits-up: track populated bucket range for fast quantile.
169 if d.min_bucket < 0 {
170 d.min_bucket = b
171 d.max_bucket = b
172 }
173 if d.min_bucket >= 0 {
174 if b < d.min_bucket { d.min_bucket = b }
175 if b > d.max_bucket { d.max_bucket = b }
176 }
177 return 0
178}
179
180// === quantile =====================================================
181//
182// Walk buckets accumulating count; return value-at-bucket for the
183// bucket containing the p·total rank.
184
185func nx_dd_quantile(d: *DDSketch, p_milli: i64) -> i64 {
186 if d.total == 0 { return 0 }
187 let target: i64 = (p_milli * d.total) / 1000
188 var cum: i64 = d.neg_count + d.zero_count
189 if cum >= target {
190 if target <= d.neg_count { return -1 } // negative quantile sentinel
191 return 0
192 }
193 // Bits-up: walk only [min_bucket, max_bucket] instead of full range.
194 // For typical heavy-tail streams populating <50% of buckets, this
195 // halves to tenths the query work. Sentinel -1 means no positive
196 // values inserted; the early returns above already handled zero/neg.
197 if d.min_bucket < 0 { return 0 }
198 var i: i64 = d.min_bucket
199 while i <= d.max_bucket {
200 let c: i64 = d.buckets[i]
201 if c > 0 {
202 cum = cum + c
203 if cum >= target {
204 return nx_dd_value_at_bucket(d, i)
205 }
206 }
207 i = i + 1
208 }
209 return 0
210}
211
212// === typed envelope ===============================================
213//
214// Relative error = 1/BPO. For BPO=64: 1.56% = 15_600_000 ppb.
215
216func nx_dd_query_quantile(d: *DDSketch, p_milli: i64) -> *ApproxI64 {
217 let v: i64 = nx_dd_quantile(d, p_milli)
218 let err_ppb: i64 = 1000000000 / d.bpo
219 return nx_approx_new(v, NX_ENV_REL_STDDEV, err_ppb,
220 1000000000,
221 NX_MATURITY_REFERENCE_IMPL,
222 NX_ADV_HONEST)
223}
224
225// === merge ========================================================
226
227func nx_dd_merge(a: *DDSketch, b: *DDSketch) -> *DDSketch {
228 if a.bpo != b.bpo { return 0 as *DDSketch }
229 let out: *DDSketch = nx_dd_alloc(a.bpo)
230 // Bits-up: walk only the union of populated bucket ranges of a, b
231 // instead of all 4032 buckets. alloc() already zeroed out.buckets,
232 // so unwritten cells stay zero.
233 var lo: i64 = a.min_bucket
234 var hi: i64 = a.max_bucket
235 if b.min_bucket >= 0 {
236 if lo < 0 { lo = b.min_bucket }
237 if b.min_bucket < lo { lo = b.min_bucket }
238 if b.max_bucket > hi { hi = b.max_bucket }
239 }
240 if lo >= 0 {
241 var i: i64 = lo
242 while i <= hi {
243 let sum: i64 = a.buckets[i] + b.buckets[i]
244 out.buckets[i] = sum
245 i = i + 1
246 }
247 out.min_bucket = lo
248 out.max_bucket = hi
249 }
250 out.total = a.total + b.total
251 out.zero_count = a.zero_count + b.zero_count
252 out.neg_count = a.neg_count + b.neg_count
253 return out
254}
255
256// === introspection ================================================
257
258func nx_dd_memory_bytes(d: *DDSketch) -> i64 {
259 return 48 + (NX_DD_MAX_OCTAVE + 1) * d.bpo * 8
260}
261
262func nx_dd_total(d: *DDSketch) -> i64 {
263 return d.total
264}