nx_sketch_reservoir.nx source
↩ module page · 188 lines · 6450 B
1// sketch_reservoir.nx -- Vitter reservoir sampling + quantile.
2//
3// Vitter 1985 "Random Sampling with a Reservoir." Algorithm R:
4// maintain a cap-element reservoir. For the i-th input value
5// (1-indexed):
6// - if i <= cap: items[i-1] = value
7// - else: draw j uniformly from [1, i]; if j <= cap,
8// items[j-1] = value (else discard)
9//
10// At end of stream, items[] holds a UNIFORMLY RANDOM sample of
11// the input. Per-quantile estimation: sort items, look up the
12// floor(p * n_items)-th value. Hoeffding bound on quantile from
13// k samples: error ~ sqrt(ln(1/delta) / (2k)) ~ 1/sqrt(k) at fixed
14// delta. We declare 1/sqrt(cap) as the rank_error envelope.
15//
16// OPENS TWO ROADMAP AXES SIMULTANEOUSLY:
17// - Sampling family (Vitter / VarOpt -- DataSketches ships both;
18// this is the first axis primitive)
19// - Quantile family (rank_error envelope; T-Digest / KLL ship
20// this same family with tighter bounds; reservoir is the
21// SIMPLEST sample-based quantile)
22//
23// LCG-based PRNG for deterministic reproducibility. Per the
24// lossless-language discipline (doc 20): randomized sketches MUST
25// declare their RNG state in the typed envelope so two runs with
26// the same seed produce bit-identical output.
27
28// nx_safety_envelope:
29// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
30// sil_target: SIL1
31// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
32// verdict: NOT_YET_EVALUATED
33
34import "nx_syscalls.nx"
35import "nx_sketch_types.nx"
36
37const NX_RES_CAP_MIN: i64 = 4
38const NX_RES_CAP_MAX: i64 = 1000000
39
40// LCG constants -- glibc / Numerical Recipes lineage.
41const NX_RES_LCG_A: i64 = 1103515245
42const NX_RES_LCG_C: i64 = 12345
43const NX_RES_LCG_MOD: i64 = 0x7FFFFFFF
44
45struct Reservoir {
46 items: *i64,
47 cap: i64,
48 n_items: i64, // = min(total_seen, cap)
49 total_seen: i64,
50 rng_state: i64,
51}
52
53func nx_reservoir_alloc(cap: i64, seed: i64) -> *Reservoir {
54 if cap < NX_RES_CAP_MIN { return 0 as *Reservoir }
55 if cap > NX_RES_CAP_MAX { return 0 as *Reservoir }
56 let raw: *u8 = sys_mmap(40)
57 let r: *Reservoir = raw as *Reservoir
58 let items_raw: *u8 = sys_mmap(cap * 8)
59 r.items = items_raw as *i64
60 r.cap = cap
61 r.n_items = 0
62 r.total_seen = 0
63 r.rng_state = seed | 1 // never zero (LCG degenerate at 0)
64 return r
65}
66
67// LCG next: returns a 31-bit pseudorandom value. Mutates state.
68func nx_reservoir_rng_next(r: *Reservoir) -> i64 {
69 let next: i64 = ((r.rng_state * NX_RES_LCG_A) + NX_RES_LCG_C) & NX_RES_LCG_MOD
70 r.rng_state = next
71 return next
72}
73
74// Uniform integer in [0, n). n must be > 0.
75func nx_reservoir_rng_below(r: *Reservoir, n: i64) -> i64 {
76 let raw: i64 = nx_reservoir_rng_next(r)
77 return raw % n
78}
79
80func nx_reservoir_add(r: *Reservoir, value: i64) -> i64 {
81 r.total_seen = r.total_seen + 1
82 if r.n_items < r.cap {
83 r.items[r.n_items] = value
84 r.n_items = r.n_items + 1
85 return 0
86 }
87 // r.total_seen > cap: replace with probability cap / total_seen.
88 let j: i64 = nx_reservoir_rng_below(r, r.total_seen)
89 if j < r.cap {
90 r.items[j] = value
91 }
92 return 0
93}
94
95// === sort the sample (insertion sort; cap is small) ==============
96
97func nx_reservoir_sort(r: *Reservoir) -> i64 {
98 var i: i64 = 1
99 while i < r.n_items {
100 let cur: i64 = r.items[i]
101 var j: i64 = i - 1
102 var done: i64 = 0
103 while done == 0 {
104 if j < 0 { done = 1 }
105 if done == 0 {
106 let prev: i64 = r.items[j]
107 if prev <= cur { done = 1 }
108 if done == 0 {
109 r.items[j + 1] = prev
110 j = j - 1
111 }
112 }
113 }
114 r.items[j + 1] = cur
115 i = i + 1
116 }
117 return 0
118}
119
120// === quantile ====================================================
121//
122// p_milli: quantile in parts-per-thousand (500 = median, 990 = p99).
123// Returns the value at index floor(p_milli * n_items / 1000).
124// Caller responsible for calling nx_reservoir_sort first if the
125// reservoir hasn't been finalized -- we do it lazily here.
126
127func nx_reservoir_quantile(r: *Reservoir, p_milli: i64) -> i64 {
128 if r.n_items == 0 { return 0 }
129 nx_reservoir_sort(r)
130 var idx: i64 = (p_milli * r.n_items) / 1000
131 if idx < 0 { idx = 0 }
132 if idx >= r.n_items { idx = r.n_items - 1 }
133 return r.items[idx]
134}
135
136// === rank-of-value (inverse: how many items <= value, as fraction) ==
137
138func nx_reservoir_rank(r: *Reservoir, value: i64) -> i64 {
139 // Returns rank in parts-per-thousand. Sorts then binary-searches.
140 if r.n_items == 0 { return 0 }
141 nx_reservoir_sort(r)
142 // Count items <= value.
143 var lo: i64 = 0
144 var hi: i64 = r.n_items
145 while lo < hi {
146 let mid: i64 = (lo + hi) / 2
147 if r.items[mid] <= value {
148 lo = mid + 1
149 }
150 if r.items[mid] > value {
151 hi = mid
152 }
153 }
154 // lo now equals the number of items <= value.
155 return (lo * 1000) / r.n_items
156}
157
158// === error bound: rank_error ~ 1/sqrt(cap) ======================
159//
160// Hoeffding-style: P(|F_n - F| > eps) <= 2 * exp(-2 * n * eps^2)
161// For 95% confidence (delta=0.05): eps ~ sqrt(ln(40) / (2 * cap))
162// ~ 1.36 / sqrt(cap). We hardcode for typical cap values; for
163// other values, return a conservative bound based on cap class.
164
165func nx_reservoir_rank_error_ppb(cap: i64) -> i64 {
166 if cap <= 16 { return 340000000 } // 0.34
167 if cap <= 64 { return 170000000 } // 0.17
168 if cap <= 256 { return 85000000 } // 0.085
169 if cap <= 1024 { return 42500000 } // 0.0425
170 if cap <= 4096 { return 21250000 } // 0.02125
171 if cap <= 16384 { return 10625000 } // 0.0106
172 return 5000000 // 0.005 fallback
173}
174
175func nx_reservoir_query_quantile(r: *Reservoir, p_milli: i64) -> *ApproxI64 {
176 let v: i64 = nx_reservoir_quantile(r, p_milli)
177 return nx_approx_new(v, NX_ENV_RANK_ERROR,
178 nx_reservoir_rank_error_ppb(r.cap),
179 950000000, // 95% confidence
180 NX_MATURITY_REFERENCE_IMPL,
181 NX_ADV_HONEST)
182}
183
184// === memory introspection =======================================
185
186func nx_reservoir_memory_bytes(r: *Reservoir) -> i64 {
187 return 40 + r.cap * 8
188}