nx_linalg.nx source
↩ module page · 315 lines · 12378 B
1// nx_linalg.nx -- vector + matrix engine.
2//
3// Native NishiLang linear algebra over integer (Q-format) and real
4// elements. Vectors and matrices stored as flat arrays for tight
5// memory + cache behaviour. Engine-style: dimensions are runtime
6// parameters, not types -- one set of functions handles any size.
7
8// nx_safety_envelope:
9// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
10// sil_target: SIL1
11// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
12// verdict: NOT_YET_EVALUATED
13
14import "nx_kernel_v2.nx"
15import "nx_fixq30_lib.nx"
16
17// ===== Vector =======================================================
18// flat array of nx_int (Q10 fixed point or raw int -- caller decides
19// the meaning; the engine just does arithmetic).
20struct Vec {
21 data: *nx_int,
22 n: nx_int,
23}
24const NX_VEC_BYTES: nx_int = 16
25
26func nx_vec_new(n: nx_int) -> *Vec {
27 let v: *Vec = (sys_mmap(NX_VEC_BYTES as i64)) as *Vec
28 v.data = (sys_mmap((n * 8) as i64)) as *nx_int
29 v.n = n
30 return v
31}
32
33func nx_vec_set(v: *Vec, i: nx_int, x: nx_int) -> nx_int {
34 let p: *nx_int = ((v.data as nx_int) + (i * 8)) as *nx_int
35 p[0] = x
36 return x
37}
38
39func nx_vec_get(v: *Vec, i: nx_int) -> nx_int {
40 let p: *nx_int = ((v.data as nx_int) + (i * 8)) as *nx_int
41 return p[0]
42}
43
44func nx_vec_add(a: *Vec, b: *Vec) -> *Vec {
45 let n: nx_int = a.n
46 let r: *Vec = nx_vec_new(n)
47 var i: nx_int = 0
48 while i < n {
49 let _s: nx_int = nx_vec_set(r, i, nx_vec_get(a, i) + nx_vec_get(b, i))
50 i = i + 1
51 }
52 return r
53}
54
55func nx_vec_scale(a: *Vec, k: nx_int) -> *Vec {
56 let r: *Vec = nx_vec_new(a.n)
57 var i: nx_int = 0
58 while i < a.n {
59 let _s: nx_int = nx_vec_set(r, i, nx_vec_get(a, i) * k)
60 i = i + 1
61 }
62 return r
63}
64
65func nx_vec_dot(a: *Vec, b: *Vec) -> nx_int {
66 var sum: nx_int = 0
67 var i: nx_int = 0
68 while i < a.n {
69 sum = sum + (nx_vec_get(a, i) * nx_vec_get(b, i))
70 i = i + 1
71 }
72 return sum
73}
74
75func nx_vec_norm_sq(a: *Vec) -> nx_int {
76 return nx_vec_dot(a, a)
77}
78
79// ===== Matrix (row-major) ==========================================
80struct Mat {
81 data: *nx_int,
82 rows: nx_int,
83 cols: nx_int,
84}
85const NX_MAT_BYTES: nx_int = 24
86
87func nx_mat_new(rows: nx_int, cols: nx_int) -> *Mat {
88 let m: *Mat = (sys_mmap(NX_MAT_BYTES as i64)) as *Mat
89 m.data = (sys_mmap((rows * cols * 8) as i64)) as *nx_int
90 m.rows = rows; m.cols = cols
91 return m
92}
93
94func nx_mat_set(m: *Mat, r: nx_int, c: nx_int, x: nx_int) -> nx_int {
95 let p: *nx_int = ((m.data as nx_int) + ((r * m.cols + c) * 8)) as *nx_int
96 p[0] = x
97 return x
98}
99
100func nx_mat_get(m: *Mat, r: nx_int, c: nx_int) -> nx_int {
101 let p: *nx_int = ((m.data as nx_int) + ((r * m.cols + c) * 8)) as *nx_int
102 return p[0]
103}
104
105// matmul: A is (n x m), B is (m x p) -> result is (n x p).
106func nx_mat_mul(a: *Mat, b: *Mat) -> *Mat {
107 let r: *Mat = nx_mat_new(a.rows, b.cols)
108 var i: nx_int = 0
109 while i < a.rows {
110 var j: nx_int = 0
111 while j < b.cols {
112 var sum: nx_int = 0
113 var k: nx_int = 0
114 while k < a.cols {
115 sum = sum + (nx_mat_get(a, i, k) * nx_mat_get(b, k, j))
116 k = k + 1
117 }
118 let _s: nx_int = nx_mat_set(r, i, j, sum)
119 j = j + 1
120 }
121 i = i + 1
122 }
123 return r
124}
125
126// 2x2 determinant -- exact integer arithmetic.
127func nx_mat_det2(m: *Mat) -> nx_int {
128 return nx_mat_get(m, 0, 0) * nx_mat_get(m, 1, 1) - nx_mat_get(m, 0, 1) * nx_mat_get(m, 1, 0)
129}
130
131// 3x3 determinant via cofactor expansion along row 0.
132func nx_mat_det3(m: *Mat) -> nx_int {
133 let a: nx_int = nx_mat_get(m, 0, 0)
134 let b: nx_int = nx_mat_get(m, 0, 1)
135 let c: nx_int = nx_mat_get(m, 0, 2)
136 let d: nx_int = nx_mat_get(m, 1, 0) * nx_mat_get(m, 2, 1) - nx_mat_get(m, 1, 1) * nx_mat_get(m, 2, 0)
137 let e: nx_int = nx_mat_get(m, 1, 0) * nx_mat_get(m, 2, 2) - nx_mat_get(m, 1, 2) * nx_mat_get(m, 2, 0)
138 let f: nx_int = nx_mat_get(m, 1, 1) * nx_mat_get(m, 2, 2) - nx_mat_get(m, 1, 2) * nx_mat_get(m, 2, 1)
139 return a * f - b * e + c * d
140}
141
142// Transpose of an (r x c) matrix gives (c x r).
143func nx_mat_transpose(m: *Mat) -> *Mat {
144 let t: *Mat = nx_mat_new(m.cols, m.rows)
145 var i: nx_int = 0
146 while i < m.rows {
147 var j: nx_int = 0
148 while j < m.cols {
149 let _s: nx_int = nx_mat_set(t, j, i, nx_mat_get(m, i, j))
150 j = j + 1
151 }
152 i = i + 1
153 }
154 return t
155}
156
157// ===== Q30 CHOLESKY (graphics GR47) ================================================================
158//
159// WHY: the correlated identity draw (GR24 bg_identity_space) needs a full-covariance sample rather than
160// independent per-channel sampling, and a covariance is sampled by factoring it. Measured 2026-09-03, this
161// estate had determinant and inverse to 3x3 and NO decomposition anywhere by name, so that draw had no
162// primitive. Independent per-channel sampling is precisely what makes a generated face read as an
163// implausible combination of plausible parts -- the axis a finite asset library cannot enter.
164//
165// A = L * L^T with L lower-triangular, in Q30 over nx_fixq30_lib arithmetic (fq_mul / fq_div / fq_sqrt).
166// The recurrence is the standard Cholesky-Banachiewicz one, taken column by column:
167// L[j][j] = sqrt( A[j][j] - sum_{k<j} L[j][k]^2 )
168// L[i][j] = ( A[i][j] - sum_{k<j} L[i][k] L[j][k] ) / L[j][j] for i > j
169//
170// NON-SPD IS REFUSED BY NAME, NEVER APPROXIMATED. The factorisation itself is the test: if any diagonal
171// radicand is <= 0 the matrix is not positive definite and there is no real L, so this returns
172// LA_CHOL_NOT_PD naming the failure rather than clamping to zero and handing back a matrix that looks
173// like a factor and is not one. Non-square and non-symmetric are separate named refusals, because a
174// caller that passed a transposed or ragged buffer needs to know WHICH mistake it made.
175//
176// THE BAND IS MEASURED, NOT ASSERTED. Every fq_mul, fq_div and fq_sqrt truncates at Q30, and those errors
177// compound through the division by L[j][j] and again through the n-term reconstruction sum, so the residual
178// could grow with dimension and with conditioning. la_chol_max_resid_q30 IS the ruler: it recomputes
179// L * L^T and returns the worst absolute entry error in Q30 units, so a caller (and the gate) can read the
180// real number instead of trusting a claim in this comment. LA_CHOL_BAND_Q30 is the ratchet the gate holds
181// it under, and the gate measures at the dimension cap rather than only at convenient sizes.
182//
183// DIMENSION IS CAPPED AND THE REASON IS THE BAND, NOT TASTE: past LA_CHOL_MAX_N the residual growth has not
184// been measured here, so the organ REFUSES rather than returning an unvalidated factorisation.
185
186const LA_CHOL_OK: nx_int = 0
187const LA_CHOL_NOT_SQUARE: nx_int = 1
188const LA_CHOL_NOT_SYMMETRIC: nx_int = 2
189const LA_CHOL_NOT_PD: nx_int = 3
190const LA_CHOL_BAD_DIM: nx_int = 4
191const LA_CHOL_N: nx_int = 5
192
193const LA_CHOL_MAX_N: nx_int = 32 // the largest dimension whose residual band has been measured here
194// RATCHET, SET FROM MEASUREMENT AND NOT FROM COMFORT. Worst |A - L L^T| over the gate's fixtures is 4 Q30
195// units (identity 0, textbook 3x3 0, irrational-root 3x3 3, dim-6 4, dim-32 4) -- notably FLAT in dimension,
196// so the error is dominated by per-operation truncation rather than by accumulation. 64 is 16x that worst
197// case: enough headroom for worse-conditioned SPD inputs than these fixtures, and still ~6e-8 relative to
198// FQ_ONE. It was first written at 4096 and the gate's own ratchet tooth REFUSED it as decorative -- a bar
199// three orders above the worst observation would never have caught a regression.
200const LA_CHOL_BAND_Q30: nx_int = 64 // ratchet: worst |A - L L^T| entry the gate permits, in Q30 units
201
202func la_chol_status_name(s: nx_int) -> *u8 {
203 if s == LA_CHOL_OK { return "OK" as *u8 }
204 if s == LA_CHOL_NOT_SQUARE { return "NOT-SQUARE" as *u8 }
205 if s == LA_CHOL_NOT_SYMMETRIC { return "NOT-SYMMETRIC" as *u8 }
206 if s == LA_CHOL_NOT_PD { return "NOT-POSITIVE-DEFINITE" as *u8 }
207 if s == LA_CHOL_BAD_DIM { return "DIMENSION-OUT-OF-MEASURED-RANGE" as *u8 }
208 return "UNKNOWN" as *u8
209}
210
211// exact symmetry: entries are Q30 integers, so a symmetric input matches bit for bit and anything else is
212// a different matrix. Approximate symmetry is deliberately NOT accepted -- it would silently factor
213// whatever the caller happened to put in the upper triangle.
214func la_is_symmetric(a: *Mat) -> nx_int {
215 if a.rows != a.cols { return 0 }
216 var i: nx_int = 0
217 while i < a.rows {
218 var j: nx_int = 0
219 while j < i {
220 if nx_mat_get(a, i, j) != nx_mat_get(a, j, i) { return 0 }
221 j = j + 1
222 }
223 i = i + 1
224 }
225 return 1
226}
227
228// Factor A into lower-triangular L, both n x n Q30. L is fully written (upper triangle zeroed) so a caller
229// can hand it straight to nx_mat_mul with its transpose. Returns LA_CHOL_*; on any refusal L is left zeroed
230// rather than half-written, because a partially-filled factor is the shape most likely to be used by mistake.
231func la_cholesky_q30(a: *Mat, l: *Mat) -> nx_int {
232 if a.rows != a.cols { return LA_CHOL_NOT_SQUARE }
233 if l.rows != a.rows { return LA_CHOL_NOT_SQUARE }
234 if l.cols != a.cols { return LA_CHOL_NOT_SQUARE }
235 let n: nx_int = a.rows
236 if n < 1 { return LA_CHOL_BAD_DIM }
237 if n > LA_CHOL_MAX_N { return LA_CHOL_BAD_DIM }
238 var zi: nx_int = 0
239 while zi < n {
240 var zj: nx_int = 0
241 while zj < n { let _z: nx_int = nx_mat_set(l, zi, zj, 0); zj = zj + 1 }
242 zi = zi + 1
243 }
244 if la_is_symmetric(a) == 0 { return LA_CHOL_NOT_SYMMETRIC }
245 var j: nx_int = 0
246 while j < n {
247 var d: nx_int = nx_mat_get(a, j, j)
248 var k: nx_int = 0
249 while k < j {
250 let ljk: nx_int = nx_mat_get(l, j, k)
251 d = d - fq_mul(ljk, ljk)
252 k = k + 1
253 }
254 // THE FACTORISATION IS ITS OWN POSITIVE-DEFINITENESS TEST: a non-positive radicand means no real
255 // factor exists. Zeroing L first means the caller cannot mistake a partial result for a factor.
256 if d <= 0 {
257 var ci: nx_int = 0
258 while ci < n {
259 var cj: nx_int = 0
260 while cj < n { let _c: nx_int = nx_mat_set(l, ci, cj, 0); cj = cj + 1 }
261 ci = ci + 1
262 }
263 return LA_CHOL_NOT_PD
264 }
265 let djj: nx_int = fq_sqrt(d)
266 if djj <= 0 { return LA_CHOL_NOT_PD }
267 let _sd: nx_int = nx_mat_set(l, j, j, djj)
268 var i: nx_int = j + 1
269 while i < n {
270 var s: nx_int = nx_mat_get(a, i, j)
271 var kk: nx_int = 0
272 while kk < j {
273 s = s - fq_mul(nx_mat_get(l, i, kk), nx_mat_get(l, j, kk))
274 kk = kk + 1
275 }
276 let _si: nx_int = nx_mat_set(l, i, j, fq_div(s, djj))
277 i = i + 1
278 }
279 j = j + 1
280 }
281 return LA_CHOL_OK
282}
283
284// THE RULER FOR THE DONE-RULE: worst absolute entry of (A - L L^T), in Q30 units. Recomputed from L rather
285// than tracked during the factorisation, so it is an independent check of the result and not a running
286// total the same arithmetic produced. -1 when the shapes do not match.
287func la_chol_max_resid_q30(a: *Mat, l: *Mat) -> nx_int {
288 if a.rows != a.cols { return 0 - 1 }
289 if l.rows != a.rows { return 0 - 1 }
290 if l.cols != a.cols { return 0 - 1 }
291 let n: nx_int = a.rows
292 var worst: nx_int = 0
293 var i: nx_int = 0
294 while i < n {
295 var j: nx_int = 0
296 while j < n {
297 var s: nx_int = 0
298 var k: nx_int = 0
299 // L is lower-triangular, so the sum runs only to min(i,j); the zeroed upper triangle makes the
300 // full loop correct too, and this bound simply avoids the wasted terms.
301 var lim: nx_int = i
302 if j < i { lim = j }
303 while k <= lim {
304 s = s + fq_mul(nx_mat_get(l, i, k), nx_mat_get(l, j, k))
305 k = k + 1
306 }
307 var e: nx_int = nx_mat_get(a, i, j) - s
308 if e < 0 { e = 0 - e }
309 if e > worst { worst = e }
310 j = j + 1
311 }
312 i = i + 1
313 }
314 return worst
315}