nx_matrix.nx source
↩ module page · 522 lines · 20257 B
1// nx_matrix.nx -- general matrix algebra (i64, row-major).
2//
3// Closes the gap that nx_kalman / nx_cam / nx_pc currently bridge with
4// inline matrix arithmetic. Provides:
5//
6// + alloc / set / get / copy
7// + identity / zeros / fill
8// + transpose
9// + multiply (C = A * B)
10// + 2x2 + 3x3 determinant
11// + 2x2 + 3x3 inverse (in Q14)
12//
13// Entries are pure i64. For fractional operations (inverse), gains
14// are emitted in Q14 fixed-point so subsequent multiplications can be
15// divided by NX_MATRIX_Q.
16//
17// genealogy_id: classical_linear_algebra + cramer_1750_determinant
18// lineage_id: row_major_matrix + cofactor_expansion
19
20// nx_safety_envelope:
21// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
22// sil_target: SIL1
23// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
24// verdict: NOT_YET_EVALUATED
25
26import "syscalls.nx"
27
28const NX_MATRIX_Q: i64 = 16384 // Q14
29
30struct Matrix {
31 rows: i64,
32 cols: i64,
33 data: *i64, // length rows*cols
34}
35
36// ===== Allocation + accessors ==========================================
37
38func nx_matrix_alloc(rows: i64, cols: i64) -> *Matrix {
39 let m: *Matrix = (sys_mmap(24)) as *Matrix
40 m.rows = rows
41 m.cols = cols
42 m.data = (sys_mmap(rows * cols * 8)) as *i64
43 var i: i64 = 0
44 while i < rows * cols {
45 m.data[i] = 0
46 i = i + 1
47 }
48 return m
49}
50
51func nx_matrix_set(m: *Matrix, r: i64, c: i64, v: i64) -> i64 {
52 m.data[r * m.cols + c] = v
53 return 0
54}
55
56func nx_matrix_get(m: *Matrix, r: i64, c: i64) -> i64 {
57 return m.data[r * m.cols + c]
58}
59
60func nx_matrix_fill(m: *Matrix, v: i64) -> i64 {
61 var i: i64 = 0
62 while i < m.rows * m.cols {
63 m.data[i] = v
64 i = i + 1
65 }
66 return 0
67}
68
69// Identity (square only). Diagonal entries get `unit_value` (use
70// NX_MATRIX_Q if you want a Q14 unit, else 1).
71
72func nx_matrix_identity(m: *Matrix, unit_value: i64) -> i64 {
73 if m.rows != m.cols { return -1 }
74 nx_matrix_fill(m, 0)
75 var i: i64 = 0
76 while i < m.rows {
77 nx_matrix_set(m, i, i, unit_value)
78 i = i + 1
79 }
80 return 0
81}
82
83// ===== Transpose ========================================================
84//
85// Returns a new Matrix with swapped dimensions.
86
87func nx_matrix_transpose(a: *Matrix) -> *Matrix {
88 let t: *Matrix = nx_matrix_alloc(a.cols, a.rows)
89 var r: i64 = 0
90 while r < a.rows {
91 var c: i64 = 0
92 while c < a.cols {
93 nx_matrix_set(t, c, r, nx_matrix_get(a, r, c))
94 c = c + 1
95 }
96 r = r + 1
97 }
98 return t
99}
100
101// ===== Multiply C = A * B ==============================================
102//
103// A: m x k, B: k x n, C: m x n. Returns 0 on success, -1 if shape
104// mismatch. Caller-provided output matrix; we use it directly.
105
106func nx_matrix_multiply(a: *Matrix, b: *Matrix, c: *Matrix) -> i64 {
107 if a.cols != b.rows { return -1 }
108 if c.rows != a.rows { return -1 }
109 if c.cols != b.cols { return -1 }
110 var i: i64 = 0
111 while i < a.rows {
112 var j: i64 = 0
113 while j < b.cols {
114 var sum: i64 = 0
115 var k: i64 = 0
116 while k < a.cols {
117 sum = sum + nx_matrix_get(a, i, k) * nx_matrix_get(b, k, j)
118 k = k + 1
119 }
120 nx_matrix_set(c, i, j, sum)
121 j = j + 1
122 }
123 i = i + 1
124 }
125 return 0
126}
127
128// ===== Determinant ======================================================
129
130// det of 2x2: ad - bc
131// NOTE: keep the expression on ONE line -- a continuation line starting with a binary op is parsed as a
132// SEPARATE statement (statement-boundary trap), which silently returned only the first product (ad) for
133// months until the run-and-compare bench executed this organ against the known answer (caught 2026-07-09).
134// ===== SYMMETRIC EIGENVALUES BY CYCLIC JACOBI ==========================
135//
136// WHY THIS EXISTS (measured 2026-09-03). nx_capsearch over 7,266 sources with corpus_complete=1 found NO
137// general eigensolver and NO SVD anywhere in this estate -- only ep_eig_dom/ep_eig_null, specialised to the
138// 3x3 essential matrix in the vision lane. /compare/appliedmath therefore graded this row ABSENT while all
139// four rivals grade Best, which was honest and is now one rung less true.
140//
141// WHY JACOBI RATHER THAN QR. Cyclic Jacobi is the algorithm whose every step is a plane rotation, so it
142// needs no shifts, no deflation heuristics and no pivoting -- and each rotation is exactly representable
143// in fixed point. It converges for ANY real symmetric matrix and its accuracy does not depend on the
144// eigenvalue spacing the way a shifted QR's does. For a fixed-point stack with no float anywhere, that
145// predictability is worth more than QR's asymptotic speed at the sizes this estate actually uses.
146//
147// THE ROTATION IS COMPUTED IN THE STABLE FORM. The textbook writes theta = (aqq-app)/(2*apq) and then
148// t = sign(theta)/(|theta|+sqrt(theta^2+1)) -- and theta OVERFLOWS in fixed point exactly when apq is
149// small, which is most of the time near convergence. The equivalent form below never forms theta at all:
150// r = sqrt(num^2 + den^2), t = den / (num +/- r)
151// so the quantity that gets small divides nothing.
152//
153// DECLARED IMPRECISION. This is Q14 fixed point, so the eigenvalues carry rounding from every rotation.
154// The gate does not assert exact equality; it asserts a RESIDUAL BOUND that the caller can read
155// (nx_matrix_residual), which is what the reference numerical libraries report and what a caller actually
156// needs in order to know what the answer is worth.
157import "nx_isqrt.nx"
158
159const NXM_EIG_OK: i64 = 0
160const NXM_EIG_ERR_NOT_SQUARE: i64 = 0 - 1
161const NXM_EIG_ERR_EMPTY: i64 = 0 - 2
162const NXM_EIG_ERR_ASYMMETRIC: i64 = 0 - 3 // Jacobi is only valid for symmetric input; refuse, never guess
163const NXM_EIG_SWEEPS_DEFAULT: i64 = 12
164// Off-diagonal magnitude at or below which a rotation is skipped. One Q14 unit is the smallest value the
165// representation can carry, so this skips only what is already at the noise floor.
166const NXM_EIG_TOL: i64 = 1
167
168func nxm_fmul(a: i64, b: i64) -> i64 { return (a * b) / NX_MATRIX_Q }
169
170func nxm_fdiv(a: i64, b: i64) -> i64 {
171 if b == 0 { return 0 }
172 return (a * NX_MATRIX_Q) / b
173}
174
175// sqrt of a Q14 value, in Q14. sqrt(x/Q)*Q == sqrt(x*Q), and the multiply happens BEFORE the root so the
176// low bits are not thrown away first.
177func nxm_fsqrt(x: i64) -> i64 {
178 if x <= 0 { return 0 }
179 return nx_isqrt(x * NX_MATRIX_Q)
180}
181
182func nxm_abs(v: i64) -> i64 { if v < 0 { return 0 - v } return v }
183
184// Is the matrix symmetric to the last representable unit? Asked BEFORE any work, because a Jacobi run on
185// an asymmetric matrix returns numbers that look like eigenvalues and are not.
186func nx_matrix_is_symmetric(m: *Matrix) -> i64 {
187 if m.rows != m.cols { return 0 }
188 var i: i64 = 0
189 while i < m.rows {
190 var j: i64 = i + 1
191 while j < m.cols {
192 if nxm_abs(nx_matrix_get(m, i, j) - nx_matrix_get(m, j, i)) > NXM_EIG_TOL { return 0 }
193 j = j + 1
194 }
195 i = i + 1
196 }
197 return 1
198}
199
200// Eigenvalues (into evals[0..n)) and eigenvectors (as COLUMNS of v) of a symmetric matrix.
201// The input matrix `a` is NOT modified: Jacobi destroys its working copy, so one is made here.
202func nx_matrix_eig_sym(a: *Matrix, evals: *i64, v: *Matrix, sweeps: i64) -> i64 {
203 if a.rows != a.cols { return NXM_EIG_ERR_NOT_SQUARE }
204 let n: i64 = a.rows
205 if n < 1 { return NXM_EIG_ERR_EMPTY }
206 if nx_matrix_is_symmetric(a) == 0 { return NXM_EIG_ERR_ASYMMETRIC }
207
208 let w: *Matrix = nx_matrix_alloc(n, n)
209 var i: i64 = 0
210 while i < n {
211 var j: i64 = 0
212 while j < n { nx_matrix_set(w, i, j, nx_matrix_get(a, i, j)); j = j + 1 }
213 i = i + 1
214 }
215 nx_matrix_identity(v, NX_MATRIX_Q)
216
217 var s_left: i64 = sweeps
218 if s_left <= 0 { s_left = NXM_EIG_SWEEPS_DEFAULT }
219 var sweep: i64 = 0
220 while sweep < s_left {
221 var p: i64 = 0
222 while p < n - 1 {
223 var q: i64 = p + 1
224 while q < n {
225 let apq: i64 = nx_matrix_get(w, p, q)
226 if nxm_abs(apq) > NXM_EIG_TOL {
227 let app: i64 = nx_matrix_get(w, p, p)
228 let aqq: i64 = nx_matrix_get(w, q, q)
229 let num: i64 = aqq - app
230 let den: i64 = 2 * apq
231 let r: i64 = nxm_fsqrt(nxm_fmul(num, num) + nxm_fmul(den, den))
232 var t: i64 = 0
233 if num >= 0 { t = nxm_fdiv(den, num + r) } else { t = nxm_fdiv(den, num - r) }
234 let c: i64 = nxm_fdiv(NX_MATRIX_Q, nxm_fsqrt(NX_MATRIX_Q + nxm_fmul(t, t)))
235 let s: i64 = nxm_fmul(t, c)
236
237 var k: i64 = 0
238 while k < n {
239 let akp: i64 = nx_matrix_get(w, k, p)
240 let akq: i64 = nx_matrix_get(w, k, q)
241 nx_matrix_set(w, k, p, nxm_fmul(c, akp) - nxm_fmul(s, akq))
242 nx_matrix_set(w, k, q, nxm_fmul(s, akp) + nxm_fmul(c, akq))
243 k = k + 1
244 }
245 k = 0
246 while k < n {
247 let apk: i64 = nx_matrix_get(w, p, k)
248 let aqk: i64 = nx_matrix_get(w, q, k)
249 nx_matrix_set(w, p, k, nxm_fmul(c, apk) - nxm_fmul(s, aqk))
250 nx_matrix_set(w, q, k, nxm_fmul(s, apk) + nxm_fmul(c, aqk))
251 k = k + 1
252 }
253 k = 0
254 while k < n {
255 let vkp: i64 = nx_matrix_get(v, k, p)
256 let vkq: i64 = nx_matrix_get(v, k, q)
257 nx_matrix_set(v, k, p, nxm_fmul(c, vkp) - nxm_fmul(s, vkq))
258 nx_matrix_set(v, k, q, nxm_fmul(s, vkp) + nxm_fmul(c, vkq))
259 k = k + 1
260 }
261 }
262 q = q + 1
263 }
264 p = p + 1
265 }
266 sweep = sweep + 1
267 }
268
269 var d: i64 = 0
270 while d < n { evals[d] = nx_matrix_get(w, d, d); d = d + 1 }
271 return NXM_EIG_OK
272}
273
274// THE RESIDUAL A CALLER CAN READ: max over all entries of |A*V - V*diag(evals)|, in Q14 units.
275// ★A NUMBER WITHOUT A STATED BOUND IS THE REAL HAZARD IN ANY FIXED-POINT STACK. The reference numerical
276// libraries return an error bound with every decomposition for exactly this reason; returning the factors
277// alone tells a caller what the answer is but never what it is worth.
278func nx_matrix_residual(a: *Matrix, v: *Matrix, evals: *i64) -> i64 {
279 let n: i64 = a.rows
280 var worst: i64 = 0
281 var i: i64 = 0
282 while i < n {
283 var j: i64 = 0
284 while j < n {
285 // (A*V)[i][j]
286 var av: i64 = 0
287 var k: i64 = 0
288 while k < n { av = av + nxm_fmul(nx_matrix_get(a, i, k), nx_matrix_get(v, k, j)); k = k + 1 }
289 // (V*Lambda)[i][j] = V[i][j] * evals[j]
290 let vl: i64 = nxm_fmul(nx_matrix_get(v, i, j), evals[j])
291 let e: i64 = nxm_abs(av - vl)
292 if e > worst { worst = e }
293 j = j + 1
294 }
295 i = i + 1
296 }
297 return worst
298}
299
300// ===== SINGULAR VALUE DECOMPOSITION ====================================
301//
302// A = U * diag(sv) * V^T for a general m x n matrix, built ON TOP of the symmetric eigensolver above
303// rather than as a second algorithm: V and sv come from the eigenproblem of A^T A, and U follows from
304// u_i = A v_i / sv_i. Composing the incumbent means there is exactly ONE rotation engine in this library
305// and a fix to it reaches both entry points.
306//
307// ⛔DECLARED IMPRECISION, BECAUSE IT IS THE KNOWN COST OF THIS METHOD AND HIDING IT WOULD BE DISHONEST:
308// forming A^T A SQUARES THE CONDITION NUMBER. For a well-conditioned matrix this is invisible; for a
309// nearly-rank-deficient one the small singular values lose roughly half their significant digits compared
310// with a Golub-Kahan bidiagonalisation, which never forms the product. That is the accepted trade for an
311// implementation that reuses one proven rotation engine and adds no new numerical machinery, and a caller
312// who needs the small singular values accurately should read the residual this returns rather than trust
313// the digits. nx_matrix_svd_residual exists precisely so that cost is MEASURABLE and not merely admitted.
314//
315// The singular values are returned in DESCENDING order with U and V permuted to match, because that is the
316// contract every caller of an SVD expects and an unsorted spectrum silently breaks rank-truncation.
317
318const NXM_SVD_OK: i64 = 0
319const NXM_SVD_ERR_EMPTY: i64 = 0 - 2
320const NXM_SVD_ERR_DIMS: i64 = 0 - 4 // u must be m x n and v must be n x n
321
322// out = A^T A, in Q14. Computed with the fixed-point multiply rather than nx_matrix_multiply because that
323// routine accumulates a PLAIN integer product and does not rescale -- the caller owns the fixed point, and
324// here the caller is this function.
325func nxm_ata(a: *Matrix, out: *Matrix) -> i64 {
326 let m: i64 = a.rows
327 let n: i64 = a.cols
328 var i: i64 = 0
329 while i < n {
330 var j: i64 = 0
331 while j < n {
332 var s: i64 = 0
333 var k: i64 = 0
334 while k < m { s = s + nxm_fmul(nx_matrix_get(a, k, i), nx_matrix_get(a, k, j)); k = k + 1 }
335 nx_matrix_set(out, i, j, s)
336 j = j + 1
337 }
338 i = i + 1
339 }
340 return 0
341}
342
343func nx_matrix_svd(a: *Matrix, u: *Matrix, sv: *i64, v: *Matrix, sweeps: i64) -> i64 {
344 let m: i64 = a.rows
345 let n: i64 = a.cols
346 if m < 1 { return NXM_SVD_ERR_EMPTY }
347 if n < 1 { return NXM_SVD_ERR_EMPTY }
348 if u.rows != m { return NXM_SVD_ERR_DIMS }
349 if u.cols != n { return NXM_SVD_ERR_DIMS }
350 if v.rows != n { return NXM_SVD_ERR_DIMS }
351 if v.cols != n { return NXM_SVD_ERR_DIMS }
352
353 let b: *Matrix = nx_matrix_alloc(n, n)
354 nxm_ata(a, b)
355 let lam: *i64 = (sys_mmap(8 * n)) as *i64
356 let rc: i64 = nx_matrix_eig_sym(b, lam, v, sweeps)
357 if rc != NXM_EIG_OK { return rc }
358
359 // Singular values are the square roots of the eigenvalues of A^T A. A tiny negative eigenvalue is
360 // rounding noise around zero, not a defect: clamp it to zero rather than taking a root of it.
361 var i: i64 = 0
362 while i < n {
363 var l: i64 = lam[i]
364 if l < 0 { l = 0 }
365 sv[i] = nxm_fsqrt(l)
366 i = i + 1
367 }
368
369 // DESCENDING ORDER, permuting V's columns with the values so the factorisation stays consistent.
370 var p: i64 = 0
371 while p < n - 1 {
372 var q: i64 = p + 1
373 while q < n {
374 if sv[q] > sv[p] {
375 let t: i64 = sv[p]; sv[p] = sv[q]; sv[q] = t
376 var r: i64 = 0
377 while r < n {
378 let vt: i64 = nx_matrix_get(v, r, p)
379 nx_matrix_set(v, r, p, nx_matrix_get(v, r, q))
380 nx_matrix_set(v, r, q, vt)
381 r = r + 1
382 }
383 }
384 q = q + 1
385 }
386 p = p + 1
387 }
388
389 // u_i = A v_i / sv_i. A zero singular value has no determined left vector, so that column is left
390 // ZERO rather than filled with a fabricated direction -- an invented column would reconstruct just as
391 // well and would be a lie about the null space.
392 var c: i64 = 0
393 while c < n {
394 var r: i64 = 0
395 while r < m {
396 var s: i64 = 0
397 var k: i64 = 0
398 while k < n { s = s + nxm_fmul(nx_matrix_get(a, r, k), nx_matrix_get(v, k, c)); k = k + 1 }
399 if sv[c] > 0 { nx_matrix_set(u, r, c, nxm_fdiv(s, sv[c])) } else { nx_matrix_set(u, r, c, 0) }
400 r = r + 1
401 }
402 c = c + 1
403 }
404 return NXM_SVD_OK
405}
406
407// max over all entries of |A - U*diag(sv)*V^T|, in Q14 units. THE NUMBER THAT SAYS WHAT THE FACTORISATION
408// IS WORTH -- and the one place the squared condition number of the A^T A method becomes visible instead
409// of being a footnote.
410func nx_matrix_svd_residual(a: *Matrix, u: *Matrix, sv: *i64, v: *Matrix) -> i64 {
411 let m: i64 = a.rows
412 let n: i64 = a.cols
413 var worst: i64 = 0
414 var r: i64 = 0
415 while r < m {
416 var c: i64 = 0
417 while c < n {
418 var acc: i64 = 0
419 var i: i64 = 0
420 while i < n {
421 acc = acc + nxm_fmul(nxm_fmul(nx_matrix_get(u, r, i), sv[i]), nx_matrix_get(v, c, i))
422 i = i + 1
423 }
424 let e: i64 = nxm_abs(nx_matrix_get(a, r, c) - acc)
425 if e > worst { worst = e }
426 c = c + 1
427 }
428 r = r + 1
429 }
430 return worst
431}
432
433func nx_matrix_det_2x2(m: *Matrix) -> i64 {
434 if m.rows != 2 { return 0 }
435 if m.cols != 2 { return 0 }
436 return nx_matrix_get(m, 0, 0) * nx_matrix_get(m, 1, 1) - nx_matrix_get(m, 0, 1) * nx_matrix_get(m, 1, 0)
437}
438
439// det of 3x3: cofactor expansion along first row.
440func nx_matrix_det_3x3(m: *Matrix) -> i64 {
441 if m.rows != 3 { return 0 }
442 if m.cols != 3 { return 0 }
443 let a: i64 = nx_matrix_get(m, 0, 0)
444 let b: i64 = nx_matrix_get(m, 0, 1)
445 let c: i64 = nx_matrix_get(m, 0, 2)
446 let d: i64 = nx_matrix_get(m, 1, 0)
447 let e: i64 = nx_matrix_get(m, 1, 1)
448 let f: i64 = nx_matrix_get(m, 1, 2)
449 let g: i64 = nx_matrix_get(m, 2, 0)
450 let h: i64 = nx_matrix_get(m, 2, 1)
451 let i: i64 = nx_matrix_get(m, 2, 2)
452 let cof_a: i64 = a * (e * i - f * h)
453 let cof_b: i64 = b * (d * i - f * g)
454 let cof_c: i64 = c * (d * h - e * g)
455 return cof_a - cof_b + cof_c
456}
457
458// ===== Inverse (output entries in Q14) ================================
459//
460// For exact integer matrices, the inverse generally has fractional
461// entries. We multiply by NX_MATRIX_Q before dividing by the
462// determinant so the result is a Q14 matrix.
463//
464// Caller is responsible for Q14 interpretation: when applying inv_M
465// to a vector, multiply then divide by NX_MATRIX_Q.
466
467func nx_matrix_inv_2x2(m: *Matrix, out: *Matrix) -> i64 {
468 if m.rows != 2 { return -1 }
469 if m.cols != 2 { return -1 }
470 if out.rows != 2 { return -1 }
471 if out.cols != 2 { return -1 }
472 let det: i64 = nx_matrix_det_2x2(m)
473 if det == 0 { return -1 }
474 let a: i64 = nx_matrix_get(m, 0, 0)
475 let b: i64 = nx_matrix_get(m, 0, 1)
476 let c: i64 = nx_matrix_get(m, 1, 0)
477 let d: i64 = nx_matrix_get(m, 1, 1)
478 nx_matrix_set(out, 0, 0, (d * NX_MATRIX_Q) / det)
479 nx_matrix_set(out, 0, 1, (-b * NX_MATRIX_Q) / det)
480 nx_matrix_set(out, 1, 0, (-c * NX_MATRIX_Q) / det)
481 nx_matrix_set(out, 1, 1, (a * NX_MATRIX_Q) / det)
482 return 0
483}
484
485// Inverse 3x3 via cofactor / adjugate / det. Output in Q14.
486func nx_matrix_inv_3x3(m: *Matrix, out: *Matrix) -> i64 {
487 if m.rows != 3 { return -1 }
488 if m.cols != 3 { return -1 }
489 if out.rows != 3 { return -1 }
490 if out.cols != 3 { return -1 }
491 let det: i64 = nx_matrix_det_3x3(m)
492 if det == 0 { return -1 }
493 let a: i64 = nx_matrix_get(m, 0, 0)
494 let b: i64 = nx_matrix_get(m, 0, 1)
495 let c: i64 = nx_matrix_get(m, 0, 2)
496 let d: i64 = nx_matrix_get(m, 1, 0)
497 let e: i64 = nx_matrix_get(m, 1, 1)
498 let f: i64 = nx_matrix_get(m, 1, 2)
499 let g: i64 = nx_matrix_get(m, 2, 0)
500 let h: i64 = nx_matrix_get(m, 2, 1)
501 let i: i64 = nx_matrix_get(m, 2, 2)
502 // Cofactor matrix transposed = adjugate / det.
503 let A: i64 = (e * i - f * h)
504 let B: i64 = -(d * i - f * g)
505 let C: i64 = (d * h - e * g)
506 let D: i64 = -(b * i - c * h)
507 let E: i64 = (a * i - c * g)
508 let F: i64 = -(a * h - b * g)
509 let G: i64 = (b * f - c * e)
510 let H: i64 = -(a * f - c * d)
511 let I: i64 = (a * e - b * d)
512 nx_matrix_set(out, 0, 0, (A * NX_MATRIX_Q) / det)
513 nx_matrix_set(out, 0, 1, (D * NX_MATRIX_Q) / det)
514 nx_matrix_set(out, 0, 2, (G * NX_MATRIX_Q) / det)
515 nx_matrix_set(out, 1, 0, (B * NX_MATRIX_Q) / det)
516 nx_matrix_set(out, 1, 1, (E * NX_MATRIX_Q) / det)
517 nx_matrix_set(out, 1, 2, (H * NX_MATRIX_Q) / det)
518 nx_matrix_set(out, 2, 0, (C * NX_MATRIX_Q) / det)
519 nx_matrix_set(out, 2, 1, (F * NX_MATRIX_Q) / det)
520 nx_matrix_set(out, 2, 2, (I * NX_MATRIX_Q) / det)
521 return 0
522}