nx_fd_f32.nx source
↩ module page · 425 lines · 17104 B
1// nx_fd_f32.nx -- FD v2: Frequent Directions in f32 (the cliff fix; the SIMD-substrate-serves-VRAM-math synthesis).
2//
3// module: nishi-core.ai.fd_f32 capability: VRAM_REDUCTION (low-rank KV-cache compression, precision-fixed)
4//
5// WHY (measured, 2026-06-16): the Q14 FD (sketch_freq_directions.nx) gives quality-preserving KV compression
6// only up to ~n=48 (12x), then COLLAPSES (n=64 -> 521% covariance error). Autonomous cliff-research confirmed
7// FD is numerically stable BY DESIGN in exact arithmetic (the shrink step bounds ||A^T A - B^T B||_2 <= eps||A||_F^2,
8// Liberty 1501.01711); the cliff is purely a FINITE-PRECISION artifact, and the precision-critical step is the
9// orthogonalization/deflation (Tropp 1902.08651: columns stay near-orthonormal only with enough precision).
10// FIX: do FD in f32 (24-bit mantissa vs Q14's ~14) -- arithmetic on the SSE hardware substrate (nx_f32_hw:
11// __f32_add/mul/div) + the one op hw lacks, nx_f32_sqrt. f32 compares done via plain bit-ops on the f32 bits.
12// This UNITES the SIMD/matmul perf arc with the VRAM-math arc -- the fast sovereign kernels SERVE the math.
13// Sovereign: nx_cc->nxasm SSE, no gcc. HONEST: still REFERENCE-SCALE (l<=4,d<=8 design); v3 = Jacobi/QR for d>8.
14import "syscalls.nx"
15import "nx_itoa_lib.nx" // shared MSB-first emitter (zero-alloc)
16import "nx_f32_hw.nx"
17import "nx_f32.nx"
18const FD32_MAGIC_1103515245: i64 = 1103515245
19const FD32_MAGIC_12345: i64 = 12345
20const FD32_MAGIC_7919: i64 = 7919
21const FD32_MAGIC_2048: i64 = 2048
22const FD32_MAGIC_4096: i64 = 4096
23
24const FD32_L_MAX: i64 = 4
25const FD32_D_MAX: i64 = 8
26const FD32_POWER_ITERS: i64 = 64
27const FD32_OK: i64 = 0
28const FD32_ERR_TOO_BIG: i64 = -1
29const FD32_ERR_DIM: i64 = -2
30
31struct FreqDir32 {
32 l: i64,
33 d: i64,
34 b: *i64, // l*d f32-bit entries, row-major
35 total_rows: i64,
36 scratch_v: *i64, // d
37 scratch_w: *i64, // d
38 scratch_bv: *i64, // l
39 scratch_basis: *i64, // l*d
40 scratch_s2: *i64, // l
41}
42
43// ---- f32 bit-level predicates (no extra deps; valid for IEEE binary32) ----
44func f32_is_zero_bits(x: i64) -> i64 { if (x & 0x7fffffff) == 0 { return 1 } return 0 }
45func f32_is_neg_bits(x: i64) -> i64 { if (x & 0x80000000) != 0 { if (x & 0x7fffffff) != 0 { return 1 } } return 0 }
46// non-negative f32 order == integer order of the bits.
47func f32_lt_nonneg(a: i64, b: i64) -> i64 { if a < b { return 1 } return 0 }
48
49// ---- construction ----
50func fd32_alloc(l: i64, d: i64) -> *FreqDir32 {
51 if l < 2 { return 0 as *FreqDir32 }
52 if l > FD32_L_MAX { return 0 as *FreqDir32 }
53 if d < 1 { return 0 as *FreqDir32 }
54 if d > FD32_D_MAX { return 0 as *FreqDir32 }
55 let fd: *FreqDir32 = sys_mmap(72) as *FreqDir32
56 fd.l = l
57 fd.d = d
58 fd.total_rows = 0
59 fd.b = sys_mmap(l * d * 8) as *i64
60 var i: i64 = 0
61 while i < l * d { fd.b[i] = 0; i = i + 1 }
62 fd.scratch_v = sys_mmap(d * 8) as *i64
63 fd.scratch_w = sys_mmap(d * 8) as *i64
64 fd.scratch_bv = sys_mmap(l * 8) as *i64
65 fd.scratch_basis = sys_mmap(l * d * 8) as *i64
66 fd.scratch_s2 = sys_mmap(l * 8) as *i64
67 return fd
68}
69
70func fd32_b_get(fd: *FreqDir32, r: i64, c: i64) -> i64 { return fd.b[r * fd.d + c] }
71func fd32_b_set(fd: *FreqDir32, r: i64, c: i64, v: i64) -> i64 { fd.b[r * fd.d + c] = v; return 0 }
72func fd32_basis_get(fd: *FreqDir32, idx: i64, c: i64) -> i64 { return fd.scratch_basis[idx * fd.d + c] }
73func fd32_basis_set(fd: *FreqDir32, idx: i64, c: i64, v: i64) -> i64 { fd.scratch_basis[idx * fd.d + c] = v; return 0 }
74
75// dot in f32
76func fd32_dot(a: *i64, b: *i64, n: i64) -> i64 {
77 var s: i64 = 0
78 var i: i64 = 0
79 while i < n { s = f32_add(s, f32_mul(a[i], b[i])); i = i + 1 }
80 return s
81}
82
83func fd32_row_norm_sq(fd: *FreqDir32, r: i64) -> i64 {
84 var s: i64 = 0
85 var c: i64 = 0
86 while c < fd.d { let v: i64 = fd32_b_get(fd, r, c); s = f32_add(s, f32_mul(v, v)); c = c + 1 }
87 return s
88}
89
90func fd32_empty_row_index(fd: *FreqDir32) -> i64 {
91 var r: i64 = 0
92 while r < fd.l { if f32_is_zero_bits(fd32_row_norm_sq(fd, r)) == 1 { return r } r = r + 1 }
93 return -1
94}
95
96func fd32_vscale(v: *i64, n: i64, factor: i64) -> i64 {
97 var i: i64 = 0
98 while i < n { v[i] = f32_mul(v[i], factor); i = i + 1 }
99 return 0
100}
101
102func fd32_vnormalize(v: *i64, n: i64) -> i64 {
103 let ns: i64 = fd32_dot(v, v, n)
104 if f32_is_zero_bits(ns) == 1 { return -1 }
105 let norm: i64 = nx_f32_sqrt(ns)
106 if f32_is_zero_bits(norm) == 1 { return -1 }
107 let inv: i64 = f32_div(f32_of(1), norm)
108 fd32_vscale(v, n, inv)
109 return 0
110}
111
112func fd32_matvec_Bv(fd: *FreqDir32, v: *i64, out: *i64) -> i64 {
113 var r: i64 = 0
114 while r < fd.l {
115 var s: i64 = 0
116 var c: i64 = 0
117 while c < fd.d { s = f32_add(s, f32_mul(fd32_b_get(fd, r, c), v[c])); c = c + 1 }
118 out[r] = s
119 r = r + 1
120 }
121 return 0
122}
123
124func fd32_matvec_Btw(fd: *FreqDir32, w: *i64, out: *i64) -> i64 {
125 var c: i64 = 0
126 while c < fd.d {
127 var s: i64 = 0
128 var r: i64 = 0
129 while r < fd.l { s = f32_add(s, f32_mul(fd32_b_get(fd, r, c), w[r])); r = r + 1 }
130 out[c] = s
131 c = c + 1
132 }
133 return 0
134}
135
136// orthogonalize v against the first `found` basis vectors (the precision-critical step).
137func fd32_deflate(fd: *FreqDir32, v: *i64, found: i64) -> i64 {
138 var j: i64 = 0
139 while j < found {
140 var s: i64 = 0
141 var c: i64 = 0
142 while c < fd.d { s = f32_add(s, f32_mul(v[c], fd32_basis_get(fd, j, c))); c = c + 1 }
143 c = 0
144 while c < fd.d { v[c] = f32_sub(v[c], f32_mul(s, fd32_basis_get(fd, j, c))); c = c + 1 }
145 j = j + 1
146 }
147 return 0
148}
149
150// REORTHOGONALIZE: classical Gram-Schmidt TWICE (DGKS) -- one pass leaves O(eps) leakage of the dominant
151// direction into later singular vectors under f32; a second pass drives it to O(eps^2). This is the fix for
152// the M1 cliff (naive single-deflate f32 was WORSE than Q14 because precision exposed the leak). [Tropp 1902.08651]
153func fd32_reorth(fd: *FreqDir32, v: *i64, found: i64) -> i64 {
154 fd32_deflate(fd, v, found)
155 fd32_deflate(fd, v, found)
156 return 0
157}
158
159// one singular vector via power iteration; sigma^2 -> scratch_s2[found_idx], vector -> basis[found_idx].
160func fd32_power_one(fd: *FreqDir32, found_idx: i64, init_idx: i64) -> i64 {
161 var c: i64 = 0
162 while c < fd.d { fd.scratch_v[c] = 0; c = c + 1 }
163 fd.scratch_v[init_idx] = f32_of(1)
164 fd32_reorth(fd, fd.scratch_v, found_idx)
165 if fd32_vnormalize(fd.scratch_v, fd.d) < 0 { fd.scratch_s2[found_idx] = 0; return 0 }
166 var iter: i64 = 0
167 while iter < FD32_POWER_ITERS {
168 fd32_matvec_Bv(fd, fd.scratch_v, fd.scratch_bv)
169 fd32_matvec_Btw(fd, fd.scratch_bv, fd.scratch_w)
170 fd32_reorth(fd, fd.scratch_w, found_idx)
171 if f32_is_zero_bits(fd32_dot(fd.scratch_w, fd.scratch_w, fd.d)) == 1 { fd.scratch_s2[found_idx] = 0; return 0 }
172 if fd32_vnormalize(fd.scratch_w, fd.d) < 0 { fd.scratch_s2[found_idx] = 0; return 0 }
173 c = 0
174 while c < fd.d { fd.scratch_v[c] = fd.scratch_w[c]; c = c + 1 }
175 iter = iter + 1
176 }
177 fd32_matvec_Bv(fd, fd.scratch_v, fd.scratch_bv)
178 fd.scratch_s2[found_idx] = fd32_dot(fd.scratch_bv, fd.scratch_bv, fd.l)
179 c = 0
180 while c < fd.d { fd32_basis_set(fd, found_idx, c, fd.scratch_v[c]); c = c + 1 }
181 return 0
182}
183
184// shrink: SVD via power-iter, subtract sigma^2_min, rebuild B = diag(sqrt(shrunk)) V^T.
185func fd32_shrink(fd: *FreqDir32) -> i64 {
186 var i: i64 = 0
187 while i < fd.l { fd32_power_one(fd, i, i % fd.d); i = i + 1 }
188 var s2_min: i64 = fd.scratch_s2[0]
189 i = 1
190 while i < fd.l { if f32_lt_nonneg(fd.scratch_s2[i], s2_min) == 1 { s2_min = fd.scratch_s2[i] } i = i + 1 }
191 i = 0
192 while i < fd.l {
193 var shrunk: i64 = f32_sub(fd.scratch_s2[i], s2_min)
194 if f32_is_neg_bits(shrunk) == 1 { shrunk = 0 }
195 let sigma: i64 = nx_f32_sqrt(shrunk)
196 var c: i64 = 0
197 while c < fd.d { fd32_b_set(fd, i, c, f32_mul(sigma, fd32_basis_get(fd, i, c))); c = c + 1 }
198 i = i + 1
199 }
200 return 0
201}
202
203func fd32_add_row(fd: *FreqDir32, row: *i64, n: i64) -> i64 {
204 if n != fd.d { return FD32_ERR_DIM }
205 var empty: i64 = fd32_empty_row_index(fd)
206 if empty < 0 {
207 fd32_shrink(fd)
208 empty = fd32_empty_row_index(fd)
209 if empty < 0 {
210 var min_idx: i64 = 0
211 var min_norm: i64 = fd32_row_norm_sq(fd, 0)
212 var i: i64 = 1
213 while i < fd.l {
214 let nrm: i64 = fd32_row_norm_sq(fd, i)
215 if f32_lt_nonneg(nrm, min_norm) == 1 { min_norm = nrm; min_idx = i }
216 i = i + 1
217 }
218 var c: i64 = 0
219 while c < fd.d { fd32_b_set(fd, min_idx, c, 0); c = c + 1 }
220 empty = min_idx
221 }
222 }
223 var c: i64 = 0
224 while c < fd.d { fd32_b_set(fd, empty, c, row[c]); c = c + 1 }
225 fd.total_rows = fd.total_rows + 1
226 return FD32_OK
227}
228
229// ================= gate / measurement =================
230
231func g_w(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 }
232// MIGRATED to the shared emitter (debt 1785563586). The old body mmapped a scratch buffer
233// per call and never freed it. At PAGE granularity that is 4096B leaked PER CALL -- the
234// defect that took 28.5GB of a 36GB host in nx_ts_lumadiff (2MB input, ~3.66M calls).
235// nxi_* is MSB-first, allocates NOTHING, and emits identical bytes including the sign.
236func g_wn(v: i64) -> i64 { nxi_out(v); return 0 }
237
238func f32frac(num: i64, den: i64) -> i64 { return f32_div(f32_of(num), f32_of(den)) }
239
240// rank-2 KV row builder (f32, values ~O(1))
241// pseudo-random stream: advancing LCG, HIGH bits only (low bits are non-random -- never use them).
242func g_lcg_next(st: *i64) -> i64 {
243 st[0] = st[0] * FD32_MAGIC_1103515245 + FD32_MAGIC_12345
244 return (st[0] >> 20) & 0xfff
245}
246
247// GENUINELY rank-r KV: K = C[n x r] @ B[r x d], C,B drawn from an INDEPENDENT pseudo-random stream
248// (non-separable -> true rank min(r,d)). The earlier basis*coeff form was separable (both factors were
249// products f(i)*g(j)) -> effective rank ~1; the rank-sweep correctly exposed that. This is the honest one.
250func g_build_rankr(K: *i64, n: i64, d: i64, r: i64) -> i64 {
251 let st: *i64 = sys_mmap(8) as *i64
252 st[0] = FD32_MAGIC_7919 + r * 131 + d * 17
253 let C: *i64 = sys_mmap(n * r * 8) as *i64
254 let B: *i64 = sys_mmap(r * d * 8) as *i64
255 var t: i64 = 0
256 while t < n * r { C[t] = f32frac(g_lcg_next(st) - FD32_MAGIC_2048, FD32_MAGIC_4096); t = t + 1 } // zero-mean -> flatter spectrum
257 t = 0
258 while t < r * d { B[t] = f32frac(g_lcg_next(st) - FD32_MAGIC_2048, FD32_MAGIC_4096); t = t + 1 }
259 var i: i64 = 0
260 while i < n {
261 var c: i64 = 0
262 while c < d {
263 var acc: i64 = 0
264 var j: i64 = 0
265 while j < r { acc = f32_add(acc, f32_mul(C[i * r + j], B[j * d + c])); j = j + 1 }
266 K[i * d + c] = acc
267 c = c + 1
268 }
269 i = i + 1
270 }
271 return 0
272}
273
274func g_build_lowrank(K: *i64, n: i64, d: i64) -> i64 { return g_build_rankr(K, n, d, 2) } // genuinely rank-2
275func g_build_fullrank(K: *i64, n: i64, d: i64) -> i64 { return g_build_rankr(K, n, d, d) } // genuinely full rank
276
277func g_gram_full(K: *i64, n: i64, d: i64, G: *i64) -> i64 {
278 var a: i64 = 0
279 while a < d {
280 var b: i64 = 0
281 while b < d {
282 var s: i64 = 0
283 var i: i64 = 0
284 while i < n { s = f32_add(s, f32_mul(K[i * d + a], K[i * d + b])); i = i + 1 }
285 G[a * d + b] = s
286 b = b + 1
287 }
288 a = a + 1
289 }
290 return 0
291}
292
293func g_gram_sketch(fd: *FreqDir32, d: i64, G: *i64) -> i64 {
294 var a: i64 = 0
295 while a < d {
296 var b: i64 = 0
297 while b < d {
298 var s: i64 = 0
299 var r: i64 = 0
300 while r < fd.l { s = f32_add(s, f32_mul(fd32_b_get(fd, r, a), fd32_b_get(fd, r, b))); r = r + 1 }
301 G[a * d + b] = s
302 b = b + 1
303 }
304 a = a + 1
305 }
306 return 0
307}
308
309// relative Frobenius covariance error in permille (f32 throughout)
310func g_rel_permille(Gf: *i64, Gs: *i64, d: i64) -> i64 {
311 var num: i64 = 0
312 var den: i64 = 0
313 var i: i64 = 0
314 while i < d * d {
315 let diff: i64 = f32_sub(Gf[i], Gs[i])
316 num = f32_add(num, f32_mul(diff, diff))
317 den = f32_add(den, f32_mul(Gf[i], Gf[i]))
318 i = i + 1
319 }
320 if f32_is_zero_bits(den) == 1 { return 0 }
321 let rel: i64 = nx_f32_sqrt(f32_div(num, den))
322 return f32_int(f32_mul(rel, f32_of(1000)))
323}
324
325func g_run(K: *i64, n: i64, d: i64, l: i64) -> i64 {
326 let fd: *FreqDir32 = fd32_alloc(l, d)
327 if fd == (0 as *FreqDir32) { return -1 }
328 let trow: *i64 = sys_mmap(d * 8) as *i64
329 var i: i64 = 0
330 while i < n {
331 var c: i64 = 0
332 while c < d { trow[c] = K[i * d + c]; c = c + 1 }
333 fd32_add_row(fd, trow, d)
334 i = i + 1
335 }
336 let Gf: *i64 = sys_mmap(d * d * 8) as *i64
337 let Gs: *i64 = sys_mmap(d * d * 8) as *i64
338 g_gram_full(K, n, d, Gf)
339 g_gram_sketch(fd, d, Gs)
340 return g_rel_permille(Gf, Gs, d)
341}
342
343func g_lowrank_err(n: i64, d: i64, l: i64) -> i64 {
344 let K: *i64 = sys_mmap(n * d * 8) as *i64
345 g_build_lowrank(K, n, d)
346 return g_run(K, n, d, l)
347}
348
349func g_rankr_err(n: i64, d: i64, l: i64, r: i64) -> i64 {
350 let K: *i64 = sys_mmap(n * d * 8) as *i64
351 g_build_rankr(K, n, d, r)
352 return g_run(K, n, d, l)
353}
354
355func main() -> i64 {
356 let d: i64 = 8
357 let l: i64 = 4
358 g_w("=== FD v2 (f32) low-rank KV compression: does precision kill the n=64 cliff? ===\n")
359 g_w(" d="); g_wn(d); g_w(" sketch_l="); g_wn(l); g_w(" arith=SSE f32 (nx_f32_hw)+nx_f32_sqrt (Q14 v1 cliffed at n=64=521%)\n\n")
360
361 let sketch_bytes: i64 = l * d * 8
362 let ns: *i64 = sys_mmap(8 * 8) as *i64
363 ns[0] = 8; ns[1] = 16; ns[2] = 32; ns[3] = 64; ns[4] = 128; ns[5] = 256
364 var ceiling: i64 = 0
365 var err64: i64 = 0
366 var idx: i64 = 0
367 while idx < 6 {
368 let n: i64 = ns[idx]
369 let e: i64 = g_lowrank_err(n, d, l)
370 let full: i64 = n * d * 8
371 g_w(" n="); g_wn(n); g_w(" full_KV="); g_wn(full); g_w("B sketch="); g_wn(sketch_bytes)
372 g_w("B x"); g_wn(full / sketch_bytes); g_w(" cov_err="); g_wn(e); g_w("permille")
373 if e < 150 { g_w(" [quality OK]"); ceiling = n } else { g_w(" [degraded]") }
374 g_w("\n")
375 if n == 64 { err64 = e }
376 idx = idx + 1
377 }
378
379 let Kf: *i64 = sys_mmap(64 * d * 8) as *i64
380 g_build_fullrank(Kf, 64, d)
381 let errFull: i64 = g_run(Kf, 64, d, l)
382 g_w(" full-rank rank8 n=64 (neg-control) cov_err="); g_wn(errFull); g_w("permille\n\n")
383
384 g_w(" => f32 quality-preserved ceiling: n<="); g_wn(ceiling)
385 g_w(" (Q14 v1 ceiling was n<=48; n=64 went 521%)\n\n")
386
387 // --- intrinsic-rank sweep @ n=64 (REALISM: 0permille above was the EASY exact-rank-2 case;
388 // real KV is approximately low-rank -> show graceful degradation as intrinsic rank exceeds sketch l) ---
389 g_w(" intrinsic-rank sweep @ n=64 (FD holds rank<=l-1="); g_wn(l - 1); g_w(" losslessly; shrink costs 1 dim):\n")
390 let rs: *i64 = sys_mmap(8 * 8) as *i64
391 rs[0] = 1; rs[1] = 2; rs[2] = 3; rs[3] = 4; rs[4] = 5; rs[5] = 6; rs[6] = 8
392 var err_r4: i64 = 0
393 var err_r8: i64 = 0
394 var ri: i64 = 0
395 while ri < 7 {
396 let r: i64 = rs[ri]
397 let er: i64 = g_rankr_err(64, d, l, r)
398 g_w(" rank="); g_wn(r); g_w(" cov_err="); g_wn(er); g_w("permille")
399 if r < l { g_w(" [<l: lossless]") } else { g_w(" [>=l: graceful loss]") }
400 g_w("\n")
401 if r == 4 { err_r4 = er }
402 if r == 8 { err_r8 = er }
403 ri = ri + 1
404 }
405 g_w("\n")
406
407 var pass: i64 = 0
408 var fail: i64 = 0
409 // T1: f32 HOLDS at n=64 where Q14 cliffed (the precision fix works)
410 if err64 < 150 { g_w(" T1 f32 holds at n=64 (cliff GONE, was 521%): PASS ("); g_wn(err64); g_w("permille)\n"); pass = pass + 1 }
411 else { g_w(" T1 f32 holds at n=64: FAIL ("); g_wn(err64); g_w("permille)\n"); fail = fail + 1 }
412 // T2: ceiling extended past the Q14 limit of 48
413 if ceiling > 48 { g_w(" T2 ceiling extended past Q14's n=48: PASS (n<="); g_wn(ceiling); g_w(")\n"); pass = pass + 1 }
414 else { g_w(" T2 ceiling extended: FAIL\n"); fail = fail + 1 }
415 // T3: neg-control still separates (metric real)
416 if err64 * 3 < errFull { g_w(" T3 low-rank err << full-rank err (real metric): PASS\n"); pass = pass + 1 }
417 else { g_w(" T3 neg-control separation: FAIL\n"); fail = fail + 1 }
418 // T4: graceful degradation -- rank<=l captured (low err), rank>l degrades (FD holds exactly l directions)
419 if err_r8 > err_r4 { g_w(" T4 graceful: rank8 err > rank4 err (holds l dirs, degrades beyond): PASS\n"); pass = pass + 1 }
420 else { g_w(" T4 graceful degradation: FAIL\n"); fail = fail + 1 }
421
422 g_w("\n PASS="); g_wn(pass); g_w("/4 ")
423 if fail == 0 { g_w("VERDICT=GREEN (f32 substrate fixes the cliff -- SIMD work SERVES the VRAM math)\n"); sys_exit(0); return 0 }
424 g_w("VERDICT=RED\n"); sys_exit(1); return 1
425}