nx_numeric_oracle.nx source
↩ module page · 429 lines · 17109 B
1// nx_numeric_oracle.nx -- the tool that verifies every future kernel.
2//
3// Per the "build the tools that build the systems" cardinal + user's
4// "use the most modern research so we can really do incredible things
5// with the math and performance" direction.
6//
7// Synthesises:
8//
9// * Csmith / EMI (Yang+ 2011, Le+ 2014) -- differential
10// testing: feed the same input to two
11// implementations of the same spec, refuse
12// to ship if they disagree.
13// * Kahan compensated (Kahan 1965) + Neumaier (1974) -- maintain
14// summation a separate running error term so naive
15// floating sum loss-of-precision is recovered.
16// * Metamorphic testing (Chen+ 1998) -- when no ground-truth oracle
17// exists, verify INVARIANTS the answer must
18// satisfy (matmul(A, I) == A, T(T(A)) == A,
19// sum-then-square vs square-then-sum, etc.)
20// * Property-based (Claessen+Hughes QuickCheck 2000) --
21// testing generate many random inputs from a typed
22// spec; failure produces a shrunken witness.
23// * ULP bounds (IEEE 754 standard, 2008/2019) -- when
24// comparing FP, count units-in-last-place
25// distance instead of absolute / relative
26// epsilon hacks.
27// * Reproducible BLAS (Demmel+ 2013) -- bit-identical reduction
28// (ReproBLAS) regardless of thread/order; achievable
29// via deterministic accumulators.
30// * Boldo+Melquiond formal FP verification framework that
31// (Flocq, 2011) informs which properties are decidable.
32// * Herbie (Panchekha+ 2015) -- given a FP expression,
33// find an equivalent expression with smaller
34// rounding error. Inverse direction useful
35// for kernel optimization audits.
36//
37// The oracle does NOT pick ONE comparison mode -- it offers the
38// sealed taxonomy and lets the kernel author pick the right one for
39// the op. Substrate refuses to call something WIN unless the oracle
40// has a structural verdict to back it.
41//
42// genealogy_id: csmith_2011 + kahan_1965 + metamorphic_testing_1998 +
43// ieee_754_ulp + reproblas_demmel_2013 + flocq_2011 +
44// herbie_panchekha_2015
45// lineage_id: substrate_numeric_oracle_v1
46
47// nx_safety_envelope:
48// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
49// sil_target: SIL1
50// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
51// verdict: NOT_YET_EVALUATED
52
53import "nx_syscalls.nx"
54import "nx_tier.nx"
55import "nx_tensor.nx"
56import "nx_prng.nx"
57const NX_MAGIC_1024: i64 = 1024
58
59// ===== Sealed-enum: ComparisonMode =================================
60//
61// What kind of equivalence are we asserting? Picked per-call by the
62// kernel author based on the op's semantics:
63//
64// BIT_EXACT -- two outputs must be byte-identical. Used for
65// deterministic integer ops (matmul on i64 with
66// fixed reduction order).
67// ULP_BOUNDED -- two FP outputs must agree to within N units-in-
68// last-place. Used for floating ops where the
69// reduction order may differ but error is bounded.
70// EPSILON_REL -- two outputs agree to relative epsilon (|a-b|
71// <= eps * max(|a|, |b|)). Used for nonlinear ops
72// like softmax or norm where exact accumulation
73// isn't realistic.
74// METAMORPHIC -- no ground truth; verify an invariant relation
75// the output must satisfy (matmul A*I == A, etc.).
76// STRUCTURAL -- compare types/shapes/dtypes without comparing
77// values; for catching shape-bug regressions cheap.
78
79const NX_NO_CMP_BIT_EXACT: nx_int = 0
80const NX_NO_CMP_ULP_BOUNDED: nx_int = 1
81const NX_NO_CMP_EPSILON_REL: nx_int = 2
82const NX_NO_CMP_METAMORPHIC: nx_int = 3
83const NX_NO_CMP_STRUCTURAL: nx_int = 4
84const NX_NO_CMP_N_KINDS: nx_int = 5
85
86func nx_no_cmp_kind_is_valid(k: nx_int) -> nx_int {
87 if k < 0 { return 0 }
88 if k >= NX_NO_CMP_N_KINDS { return 0 }
89 return 1
90}
91
92// ===== Sealed-enum: Verdict ========================================
93//
94// Closed taxonomy of oracle outcomes. Caller switches on this
95// without if/else chains that could grow over time.
96
97const NX_NO_VERDICT_EQUAL: nx_int = 0 // exact match in chosen mode
98const NX_NO_VERDICT_DIFFERS: nx_int = 1 // values differ beyond tolerance
99const NX_NO_VERDICT_SHAPE_MISMATCH: nx_int = 2 // shapes don't match
100const NX_NO_VERDICT_DTYPE_MISMATCH: nx_int = 3 // dtypes don't match
101const NX_NO_VERDICT_NOT_COMPARABLE: nx_int = 4 // dtype not yet supported in chosen mode
102const NX_NO_VERDICT_INVARIANT_BROKEN: nx_int = 5 // metamorphic check failed
103const NX_NO_VERDICT_N_KINDS: nx_int = 6
104
105func nx_no_verdict_is_valid(v: nx_int) -> nx_int {
106 if v < 0 { return 0 }
107 if v >= NX_NO_VERDICT_N_KINDS { return 0 }
108 return 1
109}
110
111// ===== Counterexample (witness) record =============================
112//
113// When DIFFERS or INVARIANT_BROKEN, the oracle captures the FIRST
114// failing position + the actual / expected values. Per Csmith's
115// "shrink to minimum" practice, the caller can use this witness as
116// a regression seed (re-run the failing input through the buggy
117// kernel after a fix).
118//
119// Flat-array layout (substrate convention):
120// 0 failing_index_0 first dim index of first differing element
121// 1 failing_index_1
122// 2 failing_index_2
123// 3 failing_index_3 (up to 4-D witness; sufficient for vision tensors)
124// 4 actual_i64 what the kernel produced
125// 5 expected_i64 what the reference says
126// 6 abs_diff |actual - expected|
127// 7 ulp_distance ULPs between them (0 for i64; meaningful when F32 lands)
128
129const NX_NO_WITNESS_F_IDX0: nx_int = 0
130const NX_NO_WITNESS_F_IDX1: nx_int = 1
131const NX_NO_WITNESS_F_IDX2: nx_int = 2
132const NX_NO_WITNESS_F_IDX3: nx_int = 3
133const NX_NO_WITNESS_F_ACTUAL: nx_int = 4
134const NX_NO_WITNESS_F_EXPECTED: nx_int = 5
135const NX_NO_WITNESS_F_ABS_DIFF: nx_int = 6
136const NX_NO_WITNESS_F_ULP_DISTANCE: nx_int = 7
137const NX_NO_WITNESS_FIELDS: nx_int = 8
138
139func nx_no_witness_clear(w: *i64) -> nx_int {
140 var i: nx_int = 0
141 while i < NX_NO_WITNESS_FIELDS {
142 w[i] = 0
143 i = i + 1
144 }
145 return 0
146}
147
148// ===== Structural pre-check (shape + dtype) =======================
149//
150// Always cheap; runs first. If shapes or dtypes don't match the
151// values comparison can't be meaningful.
152
153func _no_structural_check(a: *NxTensor, b: *NxTensor) -> nx_int {
154 if a.dtype != b.dtype { return NX_NO_VERDICT_DTYPE_MISMATCH }
155 if a.ndim != b.ndim { return NX_NO_VERDICT_SHAPE_MISMATCH }
156 var i: nx_int = 0
157 while i < a.ndim {
158 if a.shape[i] != b.shape[i] { return NX_NO_VERDICT_SHAPE_MISMATCH }
159 i = i + 1
160 }
161 return NX_NO_VERDICT_EQUAL
162}
163
164func nx_no_check_structural(a: *NxTensor, b: *NxTensor) -> nx_int {
165 return _no_structural_check(a, b)
166}
167
168// ===== BIT_EXACT comparator (i64 tensors) ==========================
169//
170// Walks the tensors element-by-element. First differing element
171// fills the witness and returns DIFFERS. All-equal returns EQUAL.
172
173func nx_no_check_bit_exact_i64(a: *NxTensor, b: *NxTensor,
174 witness_out: *i64) -> nx_int {
175 nx_no_witness_clear(witness_out)
176 let pre: nx_int = _no_structural_check(a, b)
177 if pre != NX_NO_VERDICT_EQUAL { return pre }
178 if a.dtype != NX_DT_I64 { return NX_NO_VERDICT_NOT_COMPARABLE }
179
180 // Iterate over numel; reconstruct multi-index by linear scan.
181 // For tensors up to NX_T_MAX_NDIM=8, this works on contiguous
182 // tensors via simple linear traversal. For non-contiguous,
183 // we'd need to walk by index vector; v1 routes through
184 // contiguity check.
185 if nx_t_is_contiguous(a) == 0 { return NX_NO_VERDICT_NOT_COMPARABLE }
186 if nx_t_is_contiguous(b) == 0 { return NX_NO_VERDICT_NOT_COMPARABLE }
187
188 let pa: *i64 = a.storage as *i64
189 let pb: *i64 = b.storage as *i64
190 var k: nx_int = 0
191 while k < a.numel {
192 if pa[k] != pb[k] {
193 // Decompose k into multi-index for the witness.
194 var lin: nx_int = k
195 var d: nx_int = a.ndim - 1
196 while d >= 0 {
197 let dim: nx_int = a.shape[d]
198 let idx_d: nx_int = lin - (lin / dim) * dim
199 lin = lin / dim
200 if d < 4 {
201 witness_out[d] = idx_d
202 }
203 d = d - 1
204 }
205 witness_out[NX_NO_WITNESS_F_ACTUAL] = pa[k]
206 witness_out[NX_NO_WITNESS_F_EXPECTED] = pb[k]
207 var diff: nx_int = pa[k] - pb[k]
208 if diff < 0 { diff = 0 - diff }
209 witness_out[NX_NO_WITNESS_F_ABS_DIFF] = diff
210 return NX_NO_VERDICT_DIFFERS
211 }
212 k = k + 1
213 }
214 return NX_NO_VERDICT_EQUAL
215}
216
217// ===== EPSILON_REL comparator (i64 with relative tolerance) ========
218//
219// |a - b| <= eps_q10 * max(|a|, |b|) / 1024
220// eps_q10 is Q10 fixed-point (1024 = 1.0). For Q10 substrate this
221// is the lingua franca tolerance.
222
223func nx_no_check_epsilon_rel_q10(a: *NxTensor, b: *NxTensor,
224 eps_q10: nx_int,
225 witness_out: *i64) -> nx_int {
226 nx_no_witness_clear(witness_out)
227 let pre: nx_int = _no_structural_check(a, b)
228 if pre != NX_NO_VERDICT_EQUAL { return pre }
229 if a.dtype != NX_DT_I64 { return NX_NO_VERDICT_NOT_COMPARABLE }
230 if nx_t_is_contiguous(a) == 0 { return NX_NO_VERDICT_NOT_COMPARABLE }
231 if nx_t_is_contiguous(b) == 0 { return NX_NO_VERDICT_NOT_COMPARABLE }
232
233 let pa: *i64 = a.storage as *i64
234 let pb: *i64 = b.storage as *i64
235 var k: nx_int = 0
236 while k < a.numel {
237 var va: nx_int = pa[k]
238 var vb: nx_int = pb[k]
239 var diff: nx_int = va - vb
240 if diff < 0 { diff = 0 - diff }
241 var abs_a: nx_int = va
242 if abs_a < 0 { abs_a = 0 - abs_a }
243 var abs_b: nx_int = vb
244 if abs_b < 0 { abs_b = 0 - abs_b }
245 var mx: nx_int = abs_a
246 if abs_b > mx { mx = abs_b }
247
248 let bound: nx_int = (eps_q10 * mx) / NX_MAGIC_1024
249 if diff > bound {
250 witness_out[NX_NO_WITNESS_F_ACTUAL] = va
251 witness_out[NX_NO_WITNESS_F_EXPECTED] = vb
252 witness_out[NX_NO_WITNESS_F_ABS_DIFF] = diff
253 // Decompose k
254 var lin: nx_int = k
255 var d: nx_int = a.ndim - 1
256 while d >= 0 {
257 let dim: nx_int = a.shape[d]
258 let idx_d: nx_int = lin - (lin / dim) * dim
259 lin = lin / dim
260 if d < 4 {
261 witness_out[d] = idx_d
262 }
263 d = d - 1
264 }
265 return NX_NO_VERDICT_DIFFERS
266 }
267 k = k + 1
268 }
269 return NX_NO_VERDICT_EQUAL
270}
271
272// ===== Metamorphic: transpose-of-transpose identity ================
273//
274// Verify that permute(permute(t, p1), p2) where p1=p2=reverse equals
275// the original. Returns EQUAL on success, INVARIANT_BROKEN on
276// failure. This is the simplest metamorphic check and exercises
277// the view sharing the permute primitive offers.
278
279func nx_no_check_double_transpose_id(t: *NxTensor, witness_out: *i64) -> nx_int {
280 nx_no_witness_clear(witness_out)
281 if t.dtype != NX_DT_I64 { return NX_NO_VERDICT_NOT_COMPARABLE }
282 if nx_t_is_contiguous(t) == 0 { return NX_NO_VERDICT_NOT_COMPARABLE }
283
284 // Reverse-permutation [ndim-1, ndim-2, ..., 0]
285 let perm: *i64 = (sys_mmap(NX_T_SHAPE_BYTES)) as *i64
286 var i: nx_int = 0
287 while i < t.ndim {
288 perm[i] = t.ndim - 1 - i
289 i = i + 1
290 }
291 let v1: *NxTensor = nx_t_permute(t, perm)
292 let v2: *NxTensor = nx_t_permute(v1, perm)
293 // v2 should index back into the same elements as t
294
295 let pa: *i64 = t.storage as *i64
296 var k: nx_int = 0
297 while k < t.numel {
298 // Build idx from k
299 let idx: *i64 = (sys_mmap(NX_T_SHAPE_BYTES)) as *i64
300 var lin: nx_int = k
301 var d: nx_int = t.ndim - 1
302 while d >= 0 {
303 let dim: nx_int = t.shape[d]
304 idx[d] = lin - (lin / dim) * dim
305 lin = lin / dim
306 d = d - 1
307 }
308 let orig: nx_int = nx_t_get_i64(t, idx)
309 let dt: nx_int = nx_t_get_i64(v2, idx)
310 if orig != dt {
311 witness_out[NX_NO_WITNESS_F_ACTUAL] = dt
312 witness_out[NX_NO_WITNESS_F_EXPECTED] = orig
313 return NX_NO_VERDICT_INVARIANT_BROKEN
314 }
315 k = k + 1
316 }
317 return NX_NO_VERDICT_EQUAL
318}
319
320// ===== Property-based generator: random i64 tensor =================
321//
322// Builds a tensor of given shape filled with prng-drawn i64s in
323// [lo, hi). Per QuickCheck convention the seed is caller-supplied so
324// failures are reproducible.
325
326func nx_no_gen_random_i64(shape: *i64, ndim: nx_int,
327 lo: nx_int, hi: nx_int,
328 prng_state: *i64,
329 err_out: *i64) -> *NxTensor {
330 let t: *NxTensor = nx_t_alloc(NX_DT_I64, shape, ndim, err_out)
331 if (t as nx_int) == 0 { return t }
332 let p: *i64 = t.storage as *i64
333 var k: nx_int = 0
334 while k < t.numel {
335 // nx_prng_range_lo_hi draws in [lo, hi); deterministic per seed.
336 p[k] = nx_prng_range_lo_hi(prng_state, lo, hi)
337 k = k + 1
338 }
339 return t
340}
341
342// ===== Differential test: run two ops, compare ====================
343//
344// Both implementations write to caller-supplied output tensors;
345// the oracle is invoked by name, not closure (NishiLang has no
346// function pointers v1; caller composes the diff harness). We
347// provide the COMPARISON layer; the caller wires the two calls.
348
349func nx_no_differential_verdict(out_a: *NxTensor, out_b: *NxTensor,
350 mode: nx_int, tolerance_q10: nx_int,
351 witness_out: *i64) -> nx_int {
352 if mode == NX_NO_CMP_STRUCTURAL {
353 return _no_structural_check(out_a, out_b)
354 }
355 if mode == NX_NO_CMP_BIT_EXACT {
356 return nx_no_check_bit_exact_i64(out_a, out_b, witness_out)
357 }
358 if mode == NX_NO_CMP_EPSILON_REL {
359 return nx_no_check_epsilon_rel_q10(out_a, out_b, tolerance_q10,
360 witness_out)
361 }
362 // ULP_BOUNDED + METAMORPHIC handled by dedicated helpers
363 return NX_NO_VERDICT_NOT_COMPARABLE
364}
365
366// ===== Reproducibility check: run same op twice ===================
367//
368// "Same input, same op, same output -- twice." Catches accidental
369// non-determinism (heap-allocator address-dependent ordering,
370// uninitialised buffer reuse, etc.) at the kernel layer.
371//
372// Caller computes a hash of the output via SHA-256 (or any
373// deterministic hash); oracle compares. Returns EQUAL if bytes
374// match across runs, DIFFERS otherwise.
375
376func nx_no_check_reproducible(hash_run_1: *u8, hash_run_2: *u8,
377 hash_bytes: nx_int) -> nx_int {
378 var i: nx_int = 0
379 while i < hash_bytes {
380 if hash_run_1[i] != hash_run_2[i] { return NX_NO_VERDICT_DIFFERS }
381 i = i + 1
382 }
383 return NX_NO_VERDICT_EQUAL
384}
385
386// ===== Kahan compensated summation (i64 placeholder for F32 path) =
387//
388// Documents the algorithm shape; for i64 there's no rounding error
389// to compensate. When F32 lands the same function signature applies
390// and the compensation kicks in.
391//
392// s = running sum
393// c = running compensation (accumulates the lost low-order bits)
394// y = next_value - c
395// t = s + y
396// c' = (t - s) - y
397// s' = t
398//
399// For i64 with no rounding loss this is just `s += x` with c always 0.
400// Keep the API in place so callers can audit "summation method:
401// KAHAN" up front and have the right behavior when F32 lands.
402
403func nx_no_kahan_sum_i64(values: *i64, n: nx_int) -> nx_int {
404 var s: nx_int = 0
405 var i: nx_int = 0
406 while i < n {
407 s = s + values[i]
408 i = i + 1
409 }
410 return s
411}
412
413// ===== Witness pretty-printer ====================================
414//
415// Caller passes a witness from a failed check + a fd to write to.
416// Useful for debugging output: substrate emits the witness as one
417// JSON line for log ingestion. We keep it deliberately minimal --
418// caller's nx_trace_emit will wrap this in a span.
419
420const NX_NO_WITNESS_JSON_CAP: nx_int = 256
421
422func nx_no_witness_to_buf(w: *i64, buf: *u8) -> nx_int {
423 // Layout: "{\"idx\":[I0,I1,I2,I3],\"actual\":A,\"expected\":E,
424 // \"abs_diff\":D,\"ulp\":U}"
425 // We don't emit the literal here -- caller writes via nx_trace_emit
426 // or nx_json_emit which already do this shape. Return the count
427 // of fields populated so the caller knows the witness shape.
428 return NX_NO_WITNESS_FIELDS
429}