nx_matrix.nx source
↩ module page · 223 lines · 7388 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).
134func nx_matrix_det_2x2(m: *Matrix) -> i64 {
135 if m.rows != 2 { return 0 }
136 if m.cols != 2 { return 0 }
137 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)
138}
139
140// det of 3x3: cofactor expansion along first row.
141func nx_matrix_det_3x3(m: *Matrix) -> i64 {
142 if m.rows != 3 { return 0 }
143 if m.cols != 3 { return 0 }
144 let a: i64 = nx_matrix_get(m, 0, 0)
145 let b: i64 = nx_matrix_get(m, 0, 1)
146 let c: i64 = nx_matrix_get(m, 0, 2)
147 let d: i64 = nx_matrix_get(m, 1, 0)
148 let e: i64 = nx_matrix_get(m, 1, 1)
149 let f: i64 = nx_matrix_get(m, 1, 2)
150 let g: i64 = nx_matrix_get(m, 2, 0)
151 let h: i64 = nx_matrix_get(m, 2, 1)
152 let i: i64 = nx_matrix_get(m, 2, 2)
153 let cof_a: i64 = a * (e * i - f * h)
154 let cof_b: i64 = b * (d * i - f * g)
155 let cof_c: i64 = c * (d * h - e * g)
156 return cof_a - cof_b + cof_c
157}
158
159// ===== Inverse (output entries in Q14) ================================
160//
161// For exact integer matrices, the inverse generally has fractional
162// entries. We multiply by NX_MATRIX_Q before dividing by the
163// determinant so the result is a Q14 matrix.
164//
165// Caller is responsible for Q14 interpretation: when applying inv_M
166// to a vector, multiply then divide by NX_MATRIX_Q.
167
168func nx_matrix_inv_2x2(m: *Matrix, out: *Matrix) -> i64 {
169 if m.rows != 2 { return -1 }
170 if m.cols != 2 { return -1 }
171 if out.rows != 2 { return -1 }
172 if out.cols != 2 { return -1 }
173 let det: i64 = nx_matrix_det_2x2(m)
174 if det == 0 { return -1 }
175 let a: i64 = nx_matrix_get(m, 0, 0)
176 let b: i64 = nx_matrix_get(m, 0, 1)
177 let c: i64 = nx_matrix_get(m, 1, 0)
178 let d: i64 = nx_matrix_get(m, 1, 1)
179 nx_matrix_set(out, 0, 0, (d * NX_MATRIX_Q) / det)
180 nx_matrix_set(out, 0, 1, (-b * NX_MATRIX_Q) / det)
181 nx_matrix_set(out, 1, 0, (-c * NX_MATRIX_Q) / det)
182 nx_matrix_set(out, 1, 1, (a * NX_MATRIX_Q) / det)
183 return 0
184}
185
186// Inverse 3x3 via cofactor / adjugate / det. Output in Q14.
187func nx_matrix_inv_3x3(m: *Matrix, out: *Matrix) -> i64 {
188 if m.rows != 3 { return -1 }
189 if m.cols != 3 { return -1 }
190 if out.rows != 3 { return -1 }
191 if out.cols != 3 { return -1 }
192 let det: i64 = nx_matrix_det_3x3(m)
193 if det == 0 { return -1 }
194 let a: i64 = nx_matrix_get(m, 0, 0)
195 let b: i64 = nx_matrix_get(m, 0, 1)
196 let c: i64 = nx_matrix_get(m, 0, 2)
197 let d: i64 = nx_matrix_get(m, 1, 0)
198 let e: i64 = nx_matrix_get(m, 1, 1)
199 let f: i64 = nx_matrix_get(m, 1, 2)
200 let g: i64 = nx_matrix_get(m, 2, 0)
201 let h: i64 = nx_matrix_get(m, 2, 1)
202 let i: i64 = nx_matrix_get(m, 2, 2)
203 // Cofactor matrix transposed = adjugate / det.
204 let A: i64 = (e * i - f * h)
205 let B: i64 = -(d * i - f * g)
206 let C: i64 = (d * h - e * g)
207 let D: i64 = -(b * i - c * h)
208 let E: i64 = (a * i - c * g)
209 let F: i64 = -(a * h - b * g)
210 let G: i64 = (b * f - c * e)
211 let H: i64 = -(a * f - c * d)
212 let I: i64 = (a * e - b * d)
213 nx_matrix_set(out, 0, 0, (A * NX_MATRIX_Q) / det)
214 nx_matrix_set(out, 0, 1, (D * NX_MATRIX_Q) / det)
215 nx_matrix_set(out, 0, 2, (G * NX_MATRIX_Q) / det)
216 nx_matrix_set(out, 1, 0, (B * NX_MATRIX_Q) / det)
217 nx_matrix_set(out, 1, 1, (E * NX_MATRIX_Q) / det)
218 nx_matrix_set(out, 1, 2, (H * NX_MATRIX_Q) / det)
219 nx_matrix_set(out, 2, 0, (C * NX_MATRIX_Q) / det)
220 nx_matrix_set(out, 2, 1, (F * NX_MATRIX_Q) / det)
221 nx_matrix_set(out, 2, 2, (I * NX_MATRIX_Q) / det)
222 return 0
223}