nx_nofloat_autograd.nx source
↩ module page · 775 lines · 38884 B
1// nx_nofloat_autograd.nx -- NO-FLOAT (integer Q16 fixed-point) reverse-mode TENSOR autograd: the genuine
2// MISSING GENERATION in the no-float lineage (genealogy DeepMind ladder, operator first law = don't lie).
3//
4// HONEST landscape this fills (verified by reading the organs, 2026-06-22):
5// - GENERAL reverse-mode autograd ALREADY EXISTS, but on the SOFTWARE-FLOAT tower: nx_autograd /
6// nx_autograd_tensor / nx_tgrad_core all build on nx_f32 (IEEE-754 emulated). Float, not no-float.
7// - NO-FLOAT (integer Q16) training ALREADY EXISTS, but only nx_nn_train = a SINGLE linear layer with a
8// HAND-CODED analytic gradient. No general graph, no composition.
9// - What was genuinely missing = a GENERAL reverse-mode autograd in the SAME integer Q16 tower the
10// no-float Qwen inference pipeline actually uses (nx_nofloat_llm: Q16=65536, qmul=(a*b)>>16,
11// accumulate-then-shift matmul). This file is that bridge: a tape of arbitrary ops + ONE reverse
12// sweep, so an MLP (matvec->vadd->relu->matvec->vadd->mse) trains by the SAME loop that will train the
13// transformer -- all deterministic integer arithmetic (no nx_f32, no FPU, no GPU, no third-party AD).
14//
15// TAPE: stride-7 nodes {op, ai, bi, rows, cols, valp, gradp}; valp/gradp index a per-build BUMP ARENA of
16// Q16 cells (st[1] counter; no allocator -> deterministic). Forward EAGER (Q16), backward ONE reverse sweep
17// (nodes built in topological order, so reverse index order visits parents before children).
18// Backward identities (Q16) = the proven float identities with nx_f32_mul->qmul, nx_f32_add->integer +:
19// matvec y=W*x: dW[i][j] += qmul(gy[i], x[j]); dx[j] += qmul(W[i][j], gy[i])
20// vadd y=a+b: ga += gy ; gb += gy
21// relu: gx += gy iff the INPUT cell was > 0
22// mse L=(1/n)sum(p-t)^2: dp_i += qmul(gL, (2*(p_i-t_i))/n)
23// Determinism EXCEED-AXIS: integer add is EXACTLY associative -> bit-identical training independent of
24// order (the f32 tower is only bit-identical because it pins ONE order; here it is structural).
25// genealogy_id: linnainmaa_1970_reverse_mode_ad + rumelhart_1986_backprop, ported to fixed-point Q16
26// lineage_id: sovereign_nofloat_tensor_tape_autograd_v1
27// license_tier: ORIGINAL No `main` (pure library; proof lives in nx_nofloat_autograd_gate.nx).
28import "nx_syscalls.nx"
29import "nx_vecmath.nx"
30const NFA_MAGIC_5040: i64 = 5040
31const NFA_MAGIC_362880: i64 = 362880
32const NFA_MAGIC_40320: i64 = 40320
33const NFA_MAGIC_45426: i64 = 45426
34
35const NFA_LEAF: i64 = 0
36const NFA_MATVEC: i64 = 1
37const NFA_VADD: i64 = 2
38const NFA_RELU: i64 = 3
39const NFA_MSE: i64 = 4
40const NFA_SOFTMAX: i64 = 5 // attention nonlinearity
41const NFA_SILU: i64 = 6 // FFN activation (x*sigmoid(x))
42const NFA_RMSNORM: i64 = 7 // transformer normalization
43const NFA_MATMUL: i64 = 8 // C[m,p] = A[m,k] . B[k,p]
44const NFA_MATMUL_NT: i64 = 9 // S[m,p] = A[m,k] . B[p,k]^T (the Q.K^T attention contraction)
45const NFA_CMUL: i64 = 10 // elementwise scale by a Q16 constant (bi holds the constant, not a node)
46const NFA_SOFTMAX_ROWS: i64 = 11 // per-row softmax (bi=1 -> CAUSAL: row i normalizes over j<=i only)
47const NFA_ROPE: i64 = 12 // rotary position embedding on [T,hd] (parameter-free; pos = row index)
48const NFA_HADAMARD: i64 = 13 // elementwise vector multiply a (*) b (the SwiGLU gate)
49const NFA_RMSNORM_ROWS: i64 = 14 // PER-ROW (per-token) RMSNorm over the cols of an [T,d] node
50const NFA_EMBED: i64 = 15 // gather rows of E[V,dm] by integer token ids -> X[T,dm] (bi = ids ptr)
51const NFA_SOFTCE_ROWS: i64 = 16 // fused per-row softmax cross-entropy over targets (bi = target ids ptr)
52const NFA_SLICE_COLS: i64 = 17 // extract cols [c0,c0+w) of an [T,C] node -> [T,w] (bi = c0); for multi-head
53const NFA_CONCAT_COLS: i64 = 18 // concat two [T,*] nodes column-wise -> [T, wa+wb]; for multi-head
54const NFA_Q16: i64 = 65536 // 1.0 in Q16 fixed-point (matches nx_nofloat_llm)
55// Q16 transcendental constants -- copied VERBATIM from the canonical nx_nofloat_llm so the TRAINING forward
56// uses the SAME fixed-point math as inference (DRY tension noted: a future shared nx_nofloat_math imported by
57// both is the clean resolution; kept inline so this lib still imports only nx_syscalls = no-float purity is
58// trivially auditable).
59const NFA_LOG2E: i64 = 94548
60const NFA_PC0: i64 = 65536
61const NFA_PC1: i64 = 45426
62const NFA_PC2: i64 = 15743
63const NFA_PC3: i64 = 4367
64// RoPE: Q16 angle constants + Qwen2.5 rope base ln(1e6) (verbatim from nx_nofloat_llm)
65const NFA_HALF_PI: i64 = 102944
66const NFA_PI: i64 = 205887
67const NFA_3HALF_PI: i64 = 308831
68const NFA_TWO_PI: i64 = 411775
69const NFA_LN_BASE: i64 = 905421
70
71// Q16 multiply: (a*b)>>16. Accumulate in full i64 precision, ONE shift -- the no-float matmul convention.
72func nfa_qmul(a: i64, b: i64) -> i64 { return (a * b) >> 16 }
73
74// ---- Q16 fixed-point math primitives (mirror nx_nofloat_llm; pure integer, no float) ----
75func nfa_isqrt(v: i64) -> i64 { return vm_isqrt(v) }
76// exp for arg <= 0 (softmax/sigmoid call it only after a max-subtract -> non-positive arg); returns Q16.
77func nfa_fxexp(x: i64) -> i64 { var xm: i64=0-x; if x>0 { xm=0 } let ym: i64=(xm*NFA_LOG2E)>>16; let yi: i64=ym>>16; let yf: i64=ym-(yi<<16); let g: i64=NFA_Q16-yf; var t: i64=NFA_PC3; t=NFA_PC2+((g*t)>>16); t=NFA_PC1+((g*t)>>16); t=NFA_PC0+((g*t)>>16); t=t>>1; if yi>=31 { return 0 } return t>>yi }
78func nfa_sigmoid(x: i64) -> i64 { if x>=0 { let ex: i64=nfa_fxexp(0-x); return (NFA_Q16*NFA_Q16)/(NFA_Q16+ex) } let ex: i64=nfa_fxexp(x); let sp: i64=(NFA_Q16*NFA_Q16)/(NFA_Q16+ex); return NFA_Q16-sp }
79func nfa_siluf(x: i64) -> i64 { return nfa_qmul(x, nfa_sigmoid(x)) } // silu(x) = x*sigmoid(x)
80func nfa_silud(x: i64) -> i64 { let s: i64=nfa_sigmoid(x); return s + nfa_qmul(x, nfa_qmul(s, NFA_Q16 - s)) } // silu'(x) = s + x*s*(1-s)
81// Q16 sin/cos (Taylor + quadrant reduction), verbatim from nx_nofloat_llm -- for RoPE rotation angles.
82func nfa_sinq(x: i64) -> i64 { let x2: i64=nfa_qmul(x,x); let x3: i64=nfa_qmul(x2,x); let x5: i64=nfa_qmul(x3,x2); let x7: i64=nfa_qmul(x5,x2); let x9: i64=nfa_qmul(x7,x2); return x - x3/6 + x5/120 - x7/NFA_MAGIC_5040 + x9/NFA_MAGIC_362880 }
83func nfa_cosq(x: i64) -> i64 { let x2: i64=nfa_qmul(x,x); let x4: i64=nfa_qmul(x2,x2); let x6: i64=nfa_qmul(x4,x2); let x8: i64=nfa_qmul(x6,x2); return NFA_Q16 - x2/2 + x4/24 - x6/720 + x8/NFA_MAGIC_40320 }
84func nfa_reduce2pi(a: i64) -> i64 { var t: i64=a; while t<0 { t=t+NFA_TWO_PI } while t>=NFA_TWO_PI { t=t-NFA_TWO_PI } return t }
85func nfa_sinf(a: i64) -> i64 { let t: i64=nfa_reduce2pi(a); if t<NFA_HALF_PI { return nfa_sinq(t) } if t<NFA_PI { return nfa_sinq(NFA_PI-t) } if t<NFA_3HALF_PI { return 0-nfa_sinq(t-NFA_PI) } return 0-nfa_sinq(NFA_TWO_PI-t) }
86func nfa_cosf(a: i64) -> i64 { let t: i64=nfa_reduce2pi(a); if t<NFA_HALF_PI { return nfa_cosq(t) } if t<NFA_PI { return 0-nfa_cosq(NFA_PI-t) } if t<NFA_3HALF_PI { return 0-nfa_cosq(t-NFA_PI) } return nfa_cosq(NFA_TWO_PI-t) }
87// Q16 natural log for x>0: x = m*2^e with m in [1,2); ln(x)=e*ln2 + 2*atanh((m-1)/(m+1)). The atanh series
88// converges fast because (m-1)/(m+1) in [0,1/3]. For the cross-entropy VALUE only (backward uses the exact
89// softmax identity, no log). 45426 = ln(2) in Q16.
90func nfa_ln(x: i64) -> i64 {
91 if x <= 0 { return 0 }
92 var e: i64 = 0; var m: i64 = x
93 while m >= 2*NFA_Q16 { m = m >> 1; e = e + 1 }
94 while m < NFA_Q16 { m = m << 1; e = e - 1 }
95 let y: i64 = ((m - NFA_Q16) << 16) / (m + NFA_Q16)
96 let y2: i64 = nfa_qmul(y, y)
97 var term: i64 = y; var sum: i64 = y
98 term = nfa_qmul(term, y2); sum = sum + term/3
99 term = nfa_qmul(term, y2); sum = sum + term/5
100 term = nfa_qmul(term, y2); sum = sum + term/7
101 term = nfa_qmul(term, y2); sum = sum + term/9
102 return e * NFA_MAGIC_45426 + 2*sum
103}
104
105// st[0] = next node index, st[1] = next free arena cell. Allocate a node, return its index.
106func nfa_new(tape: *i64, st: *i64, op: i64, ai: i64, bi: i64, rows: i64, cols: i64) -> i64 {
107 let k: i64 = st[0]
108 let off: i64 = st[1]
109 tape[7*k+0] = op; tape[7*k+1] = ai; tape[7*k+2] = bi
110 tape[7*k+3] = rows; tape[7*k+4] = cols
111 tape[7*k+5] = off; tape[7*k+6] = off
112 st[0] = k + 1
113 st[1] = off + rows * cols
114 return k
115}
116
117// leaf holding rows*cols Q16 cells copied from src[soff..].
118func nfa_leaf(tape: *i64, vals: *i64, st: *i64, rows: i64, cols: i64, src: *i64, soff: i64) -> i64 {
119 let k: i64 = nfa_new(tape, st, NFA_LEAF, 0 - 1, 0 - 1, rows, cols)
120 let off: i64 = tape[7*k+5]
121 let sz: i64 = rows * cols
122 var i: i64 = 0
123 while i < sz { vals[off + i] = src[soff + i]; i = i + 1 }
124 return k
125}
126
127// y = W * x (W is r*c, x is c*1, y is r*1): full-precision accumulate, ONE shift (no per-term underflow).
128func nfa_matvec(tape: *i64, vals: *i64, st: *i64, aW: i64, bx: i64) -> i64 {
129 let r: i64 = tape[7*aW+3]
130 let c: i64 = tape[7*aW+4]
131 let k: i64 = nfa_new(tape, st, NFA_MATVEC, aW, bx, r, 1)
132 let offy: i64 = tape[7*k+5]
133 let offW: i64 = tape[7*aW+5]
134 let offx: i64 = tape[7*bx+5]
135 var i: i64 = 0
136 while i < r {
137 var acc: i64 = 0
138 var j: i64 = 0
139 while j < c { acc = acc + vals[offW + i*c + j] * vals[offx + j]; j = j + 1 }
140 vals[offy + i] = acc >> 16
141 i = i + 1
142 }
143 return k
144}
145
146// y = a + b (same shape)
147func nfa_vadd(tape: *i64, vals: *i64, st: *i64, a: i64, b: i64) -> i64 {
148 let r: i64 = tape[7*a+3]
149 let c: i64 = tape[7*a+4]
150 let k: i64 = nfa_new(tape, st, NFA_VADD, a, b, r, c)
151 let offy: i64 = tape[7*k+5]
152 let offa: i64 = tape[7*a+5]
153 let offb: i64 = tape[7*b+5]
154 let sz: i64 = r * c
155 var i: i64 = 0
156 while i < sz { vals[offy + i] = vals[offa + i] + vals[offb + i]; i = i + 1 }
157 return k
158}
159
160// y = relu(a) elementwise
161func nfa_relu(tape: *i64, vals: *i64, st: *i64, a: i64) -> i64 {
162 let r: i64 = tape[7*a+3]
163 let c: i64 = tape[7*a+4]
164 let k: i64 = nfa_new(tape, st, NFA_RELU, a, 0 - 1, r, c)
165 let offy: i64 = tape[7*k+5]
166 let offa: i64 = tape[7*a+5]
167 let sz: i64 = r * c
168 var i: i64 = 0
169 while i < sz {
170 var v: i64 = vals[offa + i]
171 if v < 0 { v = 0 }
172 vals[offy + i] = v
173 i = i + 1
174 }
175 return k
176}
177
178// L = (1/n) sum_i (pred_i - target_i)^2 (scalar, 1 cell). Q16: di=p_q-t_q is Q16, di*di is Q32,
179// L_q = (sum di*di) / (Q16 * n) -> back to Q16 of the mean squared error.
180func nfa_mse(tape: *i64, vals: *i64, st: *i64, pred: i64, target: i64) -> i64 {
181 let np: i64 = tape[7*pred+3] * tape[7*pred+4]
182 let k: i64 = nfa_new(tape, st, NFA_MSE, pred, target, 1, 1)
183 let offL: i64 = tape[7*k+5]
184 let offp: i64 = tape[7*pred+5]
185 let offt: i64 = tape[7*target+5]
186 var acc: i64 = 0
187 var i: i64 = 0
188 while i < np {
189 let di: i64 = vals[offp + i] - vals[offt + i]
190 acc = acc + di * di
191 i = i + 1
192 }
193 vals[offL] = acc / (NFA_Q16 * np)
194 return k
195}
196
197// ---- transformer-defining ops (the Qwen sublayer nonlinearities), each with a Q16 reverse-mode backward ----
198// y = softmax(a) over the node's n=rows*cols cells (one vector). Stores y for the backward Jacobian-vector product.
199func nfa_softmax(tape: *i64, vals: *i64, st: *i64, a: i64) -> i64 {
200 let r: i64 = tape[7*a+3]; let c: i64 = tape[7*a+4]; let n: i64 = r*c
201 let k: i64 = nfa_new(tape, st, NFA_SOFTMAX, a, 0 - 1, r, c)
202 let offy: i64 = tape[7*k+5]; let offa: i64 = tape[7*a+5]
203 var mx: i64 = vals[offa]; var i: i64 = 1
204 while i < n { if vals[offa+i] > mx { mx = vals[offa+i] } i = i + 1 }
205 var sum: i64 = 0; i = 0
206 while i < n { let e: i64 = nfa_fxexp(vals[offa+i] - mx); vals[offy+i] = e; sum = sum + e; i = i + 1 }
207 if sum <= 0 { sum = 1 }
208 i = 0
209 while i < n { vals[offy+i] = (vals[offy+i] << 16) / sum; i = i + 1 }
210 return k
211}
212// y = silu(a) elementwise = a*sigmoid(a) (the FFN activation)
213func nfa_silu(tape: *i64, vals: *i64, st: *i64, a: i64) -> i64 {
214 let r: i64 = tape[7*a+3]; let c: i64 = tape[7*a+4]; let n: i64 = r*c
215 let k: i64 = nfa_new(tape, st, NFA_SILU, a, 0 - 1, r, c)
216 let offy: i64 = tape[7*k+5]; let offa: i64 = tape[7*a+5]
217 var i: i64 = 0
218 while i < n { vals[offy+i] = nfa_siluf(vals[offa+i]); i = i + 1 }
219 return k
220}
221// y = rmsnorm(a) (no gamma): y_i = a_i / sqrt(mean(a^2)). The transformer normalization.
222func nfa_rmsnorm(tape: *i64, vals: *i64, st: *i64, a: i64) -> i64 {
223 let r: i64 = tape[7*a+3]; let c: i64 = tape[7*a+4]; let n: i64 = r*c
224 let k: i64 = nfa_new(tape, st, NFA_RMSNORM, a, 0 - 1, r, c)
225 let offy: i64 = tape[7*k+5]; let offa: i64 = tape[7*a+5]
226 var ss: i64 = 0; var i: i64 = 0
227 while i < n { ss = ss + nfa_qmul(vals[offa+i], vals[offa+i]); i = i + 1 }
228 let ms: i64 = ss / n
229 let sd: i64 = nfa_isqrt((ms + 1) << 16)
230 i = 0
231 while i < n { if sd > 0 { vals[offy+i] = (vals[offa+i] << 16) / sd } else { vals[offy+i] = 0 } i = i + 1 }
232 return k
233}
234
235// ---- attention-contraction ops (the QK^T / AV matmuls + causal row-softmax), each with Q16 backward ----
236// C = A . B (A is [m,k], B is [k,p] -> C is [m,p]): full-precision accumulate, ONE shift.
237func nfa_matmul(tape: *i64, vals: *i64, st: *i64, a: i64, b: i64) -> i64 {
238 let m: i64 = tape[7*a+3]; let kk: i64 = tape[7*a+4]; let p: i64 = tape[7*b+4]
239 let nd: i64 = nfa_new(tape, st, NFA_MATMUL, a, b, m, p)
240 let offC: i64 = tape[7*nd+5]; let offA: i64 = tape[7*a+5]; let offB: i64 = tape[7*b+5]
241 var i: i64 = 0
242 while i < m {
243 var j: i64 = 0
244 while j < p {
245 var acc: i64 = 0; var l: i64 = 0
246 while l < kk { acc = acc + vals[offA + i*kk + l] * vals[offB + l*p + j]; l = l + 1 }
247 vals[offC + i*p + j] = acc >> 16
248 j = j + 1
249 }
250 i = i + 1
251 }
252 return nd
253}
254// S = A . B^T (A is [m,k], B is [p,k] -> S is [m,p]): S[i][j] = sum_l A[i][l] B[j][l] (the Q.K^T form).
255func nfa_matmul_nt(tape: *i64, vals: *i64, st: *i64, a: i64, b: i64) -> i64 {
256 let m: i64 = tape[7*a+3]; let kk: i64 = tape[7*a+4]; let p: i64 = tape[7*b+3]
257 let nd: i64 = nfa_new(tape, st, NFA_MATMUL_NT, a, b, m, p)
258 let offS: i64 = tape[7*nd+5]; let offA: i64 = tape[7*a+5]; let offB: i64 = tape[7*b+5]
259 var i: i64 = 0
260 while i < m {
261 var j: i64 = 0
262 while j < p {
263 var acc: i64 = 0; var l: i64 = 0
264 while l < kk { acc = acc + vals[offA + i*kk + l] * vals[offB + j*kk + l]; l = l + 1 }
265 vals[offS + i*p + j] = acc >> 16
266 j = j + 1
267 }
268 i = i + 1
269 }
270 return nd
271}
272// y = a * c elementwise, c a Q16 constant (stored in bi). The 1/sqrt(d) attention scale.
273func nfa_cmul(tape: *i64, vals: *i64, st: *i64, a: i64, c_q: i64) -> i64 {
274 let r: i64 = tape[7*a+3]; let c: i64 = tape[7*a+4]; let n: i64 = r*c
275 let nd: i64 = nfa_new(tape, st, NFA_CMUL, a, c_q, r, c)
276 let offy: i64 = tape[7*nd+5]; let offa: i64 = tape[7*a+5]
277 var i: i64 = 0
278 while i < n { vals[offy+i] = nfa_qmul(vals[offa+i], c_q); i = i + 1 }
279 return nd
280}
281// y = per-row softmax(a); causal=1 -> row i normalizes over columns j<=i only (the rest are 0).
282func nfa_softmax_rows(tape: *i64, vals: *i64, st: *i64, a: i64, causal: i64) -> i64 {
283 let r: i64 = tape[7*a+3]; let c: i64 = tape[7*a+4]
284 let nd: i64 = nfa_new(tape, st, NFA_SOFTMAX_ROWS, a, causal, r, c)
285 let offy: i64 = tape[7*nd+5]; let offa: i64 = tape[7*a+5]
286 var i: i64 = 0
287 while i < r {
288 var lim: i64 = c
289 if causal == 1 { lim = i + 1 }
290 let base: i64 = i * c
291 var mx: i64 = vals[offa+base]; var j: i64 = 1
292 while j < lim { if vals[offa+base+j] > mx { mx = vals[offa+base+j] } j = j + 1 }
293 var sum: i64 = 0; j = 0
294 while j < lim { let e: i64 = nfa_fxexp(vals[offa+base+j] - mx); vals[offy+base+j] = e; sum = sum + e; j = j + 1 }
295 j = lim
296 while j < c { vals[offy+base+j] = 0; j = j + 1 }
297 if sum <= 0 { sum = 1 }
298 j = 0
299 while j < lim { vals[offy+base+j] = (vals[offy+base+j] << 16) / sum; j = j + 1 }
300 i = i + 1
301 }
302 return nd
303}
304// y = RoPE(a): rotary position embedding on [T, hd] (hd even). Row t (= position t) rotates each (2i,2i+1)
305// pair by angle t*theta_i, theta_i = base^(-2i/hd). Parameter-free; backward = rotate by -angle (orthogonal).
306func nfa_rope(tape: *i64, vals: *i64, st: *i64, a: i64) -> i64 {
307 let T: i64 = tape[7*a+3]; let hd: i64 = tape[7*a+4]; let np: i64 = hd/2
308 let nd: i64 = nfa_new(tape, st, NFA_ROPE, a, 0 - 1, T, hd)
309 let offy: i64 = tape[7*nd+5]; let offa: i64 = tape[7*a+5]
310 var t: i64 = 0
311 while t < T {
312 var i: i64 = 0
313 while i < np {
314 let theta: i64 = nfa_fxexp(0 - (i*NFA_LN_BASE)/np)
315 let ang: i64 = t * theta
316 let c: i64 = nfa_cosf(ang); let s: i64 = nfa_sinf(ang)
317 let av: i64 = vals[offa + t*hd + 2*i]; let bv: i64 = vals[offa + t*hd + 2*i + 1]
318 vals[offy + t*hd + 2*i] = nfa_qmul(av,c) - nfa_qmul(bv,s)
319 vals[offy + t*hd + 2*i + 1] = nfa_qmul(av,s) + nfa_qmul(bv,c)
320 i = i + 1
321 }
322 t = t + 1
323 }
324 return nd
325}
326// y = a (*) b elementwise (Hadamard); the SwiGLU gate: silu(gate) (*) up.
327func nfa_hadamard(tape: *i64, vals: *i64, st: *i64, a: i64, b: i64) -> i64 {
328 let r: i64 = tape[7*a+3]; let c: i64 = tape[7*a+4]; let n: i64 = r*c
329 let nd: i64 = nfa_new(tape, st, NFA_HADAMARD, a, b, r, c)
330 let offy: i64 = tape[7*nd+5]; let offa: i64 = tape[7*a+5]; let offb: i64 = tape[7*b+5]
331 var i: i64 = 0
332 while i < n { vals[offy+i] = nfa_qmul(vals[offa+i], vals[offb+i]); i = i + 1 }
333 return nd
334}
335// y = rmsnorm of EACH ROW (per-token) of an [r,c] node: y[i][j] = x[i][j]/sqrt(mean_j(x[i]^2)). Transformer pre-norm.
336func nfa_rmsnorm_rows(tape: *i64, vals: *i64, st: *i64, a: i64) -> i64 {
337 let r: i64 = tape[7*a+3]; let c: i64 = tape[7*a+4]
338 let nd: i64 = nfa_new(tape, st, NFA_RMSNORM_ROWS, a, 0 - 1, r, c)
339 let offy: i64 = tape[7*nd+5]; let offa: i64 = tape[7*a+5]
340 var i: i64 = 0
341 while i < r {
342 let base: i64 = i * c
343 var ss: i64 = 0; var j: i64 = 0
344 while j < c { ss = ss + nfa_qmul(vals[offa+base+j], vals[offa+base+j]); j = j + 1 }
345 let ms: i64 = ss / c
346 let sd: i64 = nfa_isqrt((ms + 1) << 16)
347 j = 0
348 while j < c { if sd > 0 { vals[offy+base+j] = (vals[offa+base+j] << 16) / sd } else { vals[offy+base+j] = 0 } j = j + 1 }
349 i = i + 1
350 }
351 return nd
352}
353// ---- LM head ops: embedding gather + fused softmax cross-entropy (ids are RAW ints, not Q16; no grad to ids) ----
354// X = embed(E, ids): E is [V, dm], ids is *i64 of T token ids -> X[T, dm], X[t] = E[ids[t]]. ids ptr stashed in bi.
355func nfa_embed(tape: *i64, vals: *i64, st: *i64, e_node: i64, ids: *i64, T: i64) -> i64 {
356 let dm: i64 = tape[7*e_node+4]
357 let nd: i64 = nfa_new(tape, st, NFA_EMBED, e_node, ids as i64, T, dm)
358 let offx: i64 = tape[7*nd+5]; let offE: i64 = tape[7*e_node+5]
359 var t: i64 = 0
360 while t < T {
361 let id: i64 = ids[t]
362 var j: i64 = 0
363 while j < dm { vals[offx + t*dm + j] = vals[offE + id*dm + j]; j = j + 1 }
364 t = t + 1
365 }
366 return nd
367}
368// L = (1/T) sum_t [ logsumexp(logits[t]) - logits[t][id_t] ]; logits is [T,V], tgt is *i64 of T target ids.
369// Stored value = the scalar Q16 loss. Backward (in nfa_backward) uses the exact identity dlogit=(softmax-onehot)/T.
370func nfa_softce_rows(tape: *i64, vals: *i64, st: *i64, logits: i64, tgt: *i64) -> i64 {
371 let T: i64 = tape[7*logits+3]; let V: i64 = tape[7*logits+4]
372 let nd: i64 = nfa_new(tape, st, NFA_SOFTCE_ROWS, logits, tgt as i64, 1, 1)
373 let offL: i64 = tape[7*nd+5]; let offp: i64 = tape[7*logits+5]
374 var total: i64 = 0; var t: i64 = 0
375 while t < T {
376 let base: i64 = t*V
377 var mx: i64 = vals[offp+base]; var j: i64 = 1
378 while j < V { if vals[offp+base+j] > mx { mx = vals[offp+base+j] } j = j + 1 }
379 var sum: i64 = 0; j = 0
380 while j < V { sum = sum + nfa_fxexp(vals[offp+base+j] - mx); j = j + 1 }
381 let lse: i64 = mx + nfa_ln(sum)
382 total = total + (lse - vals[offp+base+tgt[t]])
383 t = t + 1
384 }
385 vals[offL] = total / T
386 return nd
387}
388// ---- column slice / concat (the head-split bookkeeping for MULTI-HEAD attention) ----
389// y = a[:, c0:c0+w] (a is [T,C] -> y is [T,w]); c0 stored in bi. backward scatters gy into a's columns.
390func nfa_slice_cols(tape: *i64, vals: *i64, st: *i64, a: i64, c0: i64, w: i64) -> i64 {
391 let T: i64 = tape[7*a+3]; let C: i64 = tape[7*a+4]
392 let nd: i64 = nfa_new(tape, st, NFA_SLICE_COLS, a, c0, T, w)
393 let offy: i64 = tape[7*nd+5]; let offa: i64 = tape[7*a+5]
394 var t: i64 = 0
395 while t < T { var j: i64 = 0; while j < w { vals[offy + t*w + j] = vals[offa + t*C + c0 + j]; j = j + 1 } t = t + 1 }
396 return nd
397}
398// y = [a | b] (a is [T,wa], b is [T,wb] -> y is [T,wa+wb]). backward splits gy back to a and b.
399func nfa_concat_cols(tape: *i64, vals: *i64, st: *i64, a: i64, b: i64) -> i64 {
400 let T: i64 = tape[7*a+3]; let wa: i64 = tape[7*a+4]; let wb: i64 = tape[7*b+4]; let wc: i64 = wa + wb
401 let nd: i64 = nfa_new(tape, st, NFA_CONCAT_COLS, a, b, T, wc)
402 let offy: i64 = tape[7*nd+5]; let offa: i64 = tape[7*a+5]; let offb: i64 = tape[7*b+5]
403 var t: i64 = 0
404 while t < T {
405 var j: i64 = 0
406 while j < wa { vals[offy + t*wc + j] = vals[offa + t*wa + j]; j = j + 1 }
407 j = 0
408 while j < wb { vals[offy + t*wc + wa + j] = vals[offb + t*wb + j]; j = j + 1 }
409 t = t + 1
410 }
411 return nd
412}
413
414func nfa_val(tape: *i64, vals: *i64, k: i64, c: i64) -> i64 { return vals[tape[7*k+5] + c] }
415func nfa_grad(tape: *i64, grads: *i64, k: i64, c: i64) -> i64 { return grads[tape[7*k+6] + c] }
416
417// reverse sweep: zero all grads, seed grad[root]=Q16(1.0), accumulate the identities backward.
418func nfa_backward(tape: *i64, vals: *i64, grads: *i64, n: i64, root: i64) -> i64 {
419 var k: i64 = 0
420 while k < n {
421 let off: i64 = tape[7*k+5]
422 let sz: i64 = tape[7*k+3] * tape[7*k+4]
423 var i: i64 = 0
424 while i < sz { grads[off + i] = 0; i = i + 1 }
425 k = k + 1
426 }
427 grads[tape[7*root+6] + 0] = NFA_Q16
428 k = n - 1
429 while k >= 0 {
430 let op: i64 = tape[7*k+0]
431 if op == NFA_MATVEC {
432 let aW: i64 = tape[7*k+1]
433 let bx: i64 = tape[7*k+2]
434 let r: i64 = tape[7*k+3]
435 let c: i64 = tape[7*aW+4]
436 let offgy: i64 = tape[7*k+6]
437 let offW: i64 = tape[7*aW+5]
438 let gW: i64 = tape[7*aW+6]
439 let offx: i64 = tape[7*bx+5]
440 let gx: i64 = tape[7*bx+6]
441 var i: i64 = 0
442 while i < r {
443 let gyi: i64 = grads[offgy + i]
444 var j: i64 = 0
445 while j < c {
446 grads[gW + i*c + j] = grads[gW + i*c + j] + nfa_qmul(gyi, vals[offx + j])
447 grads[gx + j] = grads[gx + j] + nfa_qmul(vals[offW + i*c + j], gyi)
448 j = j + 1
449 }
450 i = i + 1
451 }
452 }
453 if op == NFA_VADD {
454 let a: i64 = tape[7*k+1]
455 let b: i64 = tape[7*k+2]
456 let sz: i64 = tape[7*k+3] * tape[7*k+4]
457 let offgy: i64 = tape[7*k+6]
458 let ga: i64 = tape[7*a+6]
459 let gb: i64 = tape[7*b+6]
460 var i: i64 = 0
461 while i < sz {
462 let g: i64 = grads[offgy + i]
463 grads[ga + i] = grads[ga + i] + g
464 grads[gb + i] = grads[gb + i] + g
465 i = i + 1
466 }
467 }
468 if op == NFA_RELU {
469 let a: i64 = tape[7*k+1]
470 let sz: i64 = tape[7*k+3] * tape[7*k+4]
471 let offgy: i64 = tape[7*k+6]
472 let ga: i64 = tape[7*a+6]
473 let offa: i64 = tape[7*a+5]
474 var i: i64 = 0
475 while i < sz {
476 if vals[offa + i] > 0 { grads[ga + i] = grads[ga + i] + grads[offgy + i] }
477 i = i + 1
478 }
479 }
480 if op == NFA_MSE {
481 let a: i64 = tape[7*k+1]
482 let b: i64 = tape[7*k+2]
483 let np: i64 = tape[7*a+3] * tape[7*a+4]
484 let offp: i64 = tape[7*a+5]
485 let gp: i64 = tape[7*a+6]
486 let offt: i64 = tape[7*b+5]
487 let gL: i64 = grads[tape[7*k+6] + 0]
488 var i: i64 = 0
489 while i < np {
490 let di: i64 = vals[offp + i] - vals[offt + i]
491 grads[gp + i] = grads[gp + i] + nfa_qmul(gL, (2 * di) / np)
492 i = i + 1
493 }
494 }
495 if op == NFA_SOFTMAX {
496 // dx_i = y_i * (gy_i - sum_j y_j gy_j) -- the softmax Jacobian-vector product
497 let a: i64 = tape[7*k+1]
498 let n: i64 = tape[7*k+3] * tape[7*k+4]
499 let offy: i64 = tape[7*k+5]
500 let offgy: i64 = tape[7*k+6]
501 let ga: i64 = tape[7*a+6]
502 var dot: i64 = 0; var j: i64 = 0
503 while j < n { dot = dot + nfa_qmul(vals[offy+j], grads[offgy+j]); j = j + 1 }
504 var i: i64 = 0
505 while i < n { grads[ga+i] = grads[ga+i] + nfa_qmul(vals[offy+i], grads[offgy+i] - dot); i = i + 1 }
506 }
507 if op == NFA_SILU {
508 // dx_i = gy_i * silu'(x_i)
509 let a: i64 = tape[7*k+1]
510 let n: i64 = tape[7*k+3] * tape[7*k+4]
511 let offgy: i64 = tape[7*k+6]
512 let offa: i64 = tape[7*a+5]
513 let ga: i64 = tape[7*a+6]
514 var i: i64 = 0
515 while i < n { grads[ga+i] = grads[ga+i] + nfa_qmul(grads[offgy+i], nfa_silud(vals[offa+i])); i = i + 1 }
516 }
517 if op == NFA_RMSNORM {
518 // dx_j = gy_j/r - y_j*dot/(n*ms), r=rms, ms=mean(x^2), dot=sum_i gy_i x_i (Q16; recompute fwd stats)
519 let a: i64 = tape[7*k+1]
520 let n: i64 = tape[7*k+3] * tape[7*k+4]
521 let offy: i64 = tape[7*k+5]
522 let offgy: i64 = tape[7*k+6]
523 let offa: i64 = tape[7*a+5]
524 let ga: i64 = tape[7*a+6]
525 var ss: i64 = 0; var i: i64 = 0
526 while i < n { ss = ss + nfa_qmul(vals[offa+i], vals[offa+i]); i = i + 1 }
527 let ms: i64 = ss / n
528 let sd: i64 = nfa_isqrt((ms + 1) << 16)
529 var dot: i64 = 0; i = 0
530 while i < n { dot = dot + nfa_qmul(grads[offgy+i], vals[offa+i]); i = i + 1 }
531 if sd > 0 { if ms > 0 {
532 i = 0
533 while i < n {
534 let t1: i64 = (grads[offgy+i] << 16) / sd
535 let t2: i64 = (vals[offy+i] * dot) / (n * ms)
536 grads[ga+i] = grads[ga+i] + t1 - t2
537 i = i + 1
538 }
539 } }
540 }
541 if op == NFA_MATMUL {
542 // C=A.B: dA[i][l] += sum_j gy[i][j] B[l][j]; dB[l][j] += sum_i A[i][l] gy[i][j]
543 let a: i64 = tape[7*k+1]; let b: i64 = tape[7*k+2]
544 let m: i64 = tape[7*k+3]; let p: i64 = tape[7*k+4]; let kk: i64 = tape[7*a+4]
545 let offgy: i64 = tape[7*k+6]; let offA: i64 = tape[7*a+5]; let gA: i64 = tape[7*a+6]; let offB: i64 = tape[7*b+5]; let gB: i64 = tape[7*b+6]
546 var i: i64 = 0
547 while i < m {
548 var j: i64 = 0
549 while j < p {
550 let g: i64 = grads[offgy + i*p + j]; var l: i64 = 0
551 while l < kk {
552 grads[gA + i*kk + l] = grads[gA + i*kk + l] + nfa_qmul(g, vals[offB + l*p + j])
553 grads[gB + l*p + j] = grads[gB + l*p + j] + nfa_qmul(vals[offA + i*kk + l], g)
554 l = l + 1
555 }
556 j = j + 1
557 }
558 i = i + 1
559 }
560 }
561 if op == NFA_MATMUL_NT {
562 // S=A.B^T: dA[i][l] += sum_j gy[i][j] B[j][l]; dB[j][l] += sum_i gy[i][j] A[i][l]
563 let a: i64 = tape[7*k+1]; let b: i64 = tape[7*k+2]
564 let m: i64 = tape[7*k+3]; let p: i64 = tape[7*k+4]; let kk: i64 = tape[7*a+4]
565 let offgy: i64 = tape[7*k+6]; let offA: i64 = tape[7*a+5]; let gA: i64 = tape[7*a+6]; let offB: i64 = tape[7*b+5]; let gB: i64 = tape[7*b+6]
566 var i: i64 = 0
567 while i < m {
568 var j: i64 = 0
569 while j < p {
570 let g: i64 = grads[offgy + i*p + j]; var l: i64 = 0
571 while l < kk {
572 grads[gA + i*kk + l] = grads[gA + i*kk + l] + nfa_qmul(g, vals[offB + j*kk + l])
573 grads[gB + j*kk + l] = grads[gB + j*kk + l] + nfa_qmul(g, vals[offA + i*kk + l])
574 l = l + 1
575 }
576 j = j + 1
577 }
578 i = i + 1
579 }
580 }
581 if op == NFA_CMUL {
582 let a: i64 = tape[7*k+1]; let c_q: i64 = tape[7*k+2]
583 let n: i64 = tape[7*k+3] * tape[7*k+4]
584 let offgy: i64 = tape[7*k+6]; let ga: i64 = tape[7*a+6]
585 var i: i64 = 0
586 while i < n { grads[ga+i] = grads[ga+i] + nfa_qmul(grads[offgy+i], c_q); i = i + 1 }
587 }
588 if op == NFA_SOFTMAX_ROWS {
589 // per-row Jacobian-vector product over j<lim (causal: lim=i+1); masked entries had 0 effect -> 0 grad
590 let a: i64 = tape[7*k+1]; let causal: i64 = tape[7*k+2]
591 let r: i64 = tape[7*k+3]; let c: i64 = tape[7*k+4]
592 let offy: i64 = tape[7*k+5]; let offgy: i64 = tape[7*k+6]; let ga: i64 = tape[7*a+6]
593 var i: i64 = 0
594 while i < r {
595 var lim: i64 = c
596 if causal == 1 { lim = i + 1 }
597 let base: i64 = i * c
598 var dot: i64 = 0; var j: i64 = 0
599 while j < lim { dot = dot + nfa_qmul(vals[offy+base+j], grads[offgy+base+j]); j = j + 1 }
600 j = 0
601 while j < lim { grads[ga+base+j] = grads[ga+base+j] + nfa_qmul(vals[offy+base+j], grads[offgy+base+j] - dot); j = j + 1 }
602 i = i + 1
603 }
604 }
605 if op == NFA_ROPE {
606 // backward = rotate gy by -ang (orthogonal Jacobian): gv[2i]=g0*c+g1*s; gv[2i+1]=g1*c-g0*s
607 let a: i64 = tape[7*k+1]; let T: i64 = tape[7*k+3]; let hd: i64 = tape[7*k+4]; let np: i64 = hd/2
608 let offgy: i64 = tape[7*k+6]; let ga: i64 = tape[7*a+6]
609 var t: i64 = 0
610 while t < T {
611 var i: i64 = 0
612 while i < np {
613 let theta: i64 = nfa_fxexp(0 - (i*NFA_LN_BASE)/np)
614 let ang: i64 = t * theta
615 let c: i64 = nfa_cosf(ang); let s: i64 = nfa_sinf(ang)
616 let g0: i64 = grads[offgy + t*hd + 2*i]; let g1: i64 = grads[offgy + t*hd + 2*i + 1]
617 grads[ga + t*hd + 2*i] = grads[ga + t*hd + 2*i] + nfa_qmul(g0,c) + nfa_qmul(g1,s)
618 grads[ga + t*hd + 2*i + 1] = grads[ga + t*hd + 2*i + 1] + nfa_qmul(g1,c) - nfa_qmul(g0,s)
619 i = i + 1
620 }
621 t = t + 1
622 }
623 }
624 if op == NFA_HADAMARD {
625 let a: i64 = tape[7*k+1]; let b: i64 = tape[7*k+2]
626 let n: i64 = tape[7*k+3] * tape[7*k+4]
627 let offgy: i64 = tape[7*k+6]; let offa: i64 = tape[7*a+5]; let offb: i64 = tape[7*b+5]
628 let ga: i64 = tape[7*a+6]; let gb: i64 = tape[7*b+6]
629 var i: i64 = 0
630 while i < n {
631 let g: i64 = grads[offgy+i]
632 grads[ga+i] = grads[ga+i] + nfa_qmul(g, vals[offb+i])
633 grads[gb+i] = grads[gb+i] + nfa_qmul(g, vals[offa+i])
634 i = i + 1
635 }
636 }
637 if op == NFA_RMSNORM_ROWS {
638 // per-row RMSNorm backward: dx_j = gy_j/r - y_j*dot/(c*ms) per row (dot = sum_j gy_j x_j), Q16
639 let a: i64 = tape[7*k+1]; let r: i64 = tape[7*k+3]; let c: i64 = tape[7*k+4]
640 let offy: i64 = tape[7*k+5]; let offgy: i64 = tape[7*k+6]; let offa: i64 = tape[7*a+5]; let ga: i64 = tape[7*a+6]
641 var i: i64 = 0
642 while i < r {
643 let base: i64 = i * c
644 var ss: i64 = 0; var j: i64 = 0
645 while j < c { ss = ss + nfa_qmul(vals[offa+base+j], vals[offa+base+j]); j = j + 1 }
646 let ms: i64 = ss / c
647 let sd: i64 = nfa_isqrt((ms + 1) << 16)
648 var dot: i64 = 0; j = 0
649 while j < c { dot = dot + nfa_qmul(grads[offgy+base+j], vals[offa+base+j]); j = j + 1 }
650 if sd > 0 { if ms > 0 {
651 j = 0
652 while j < c {
653 let t1: i64 = (grads[offgy+base+j] << 16) / sd
654 let t2: i64 = (vals[offy+base+j] * dot) / (c * ms)
655 grads[ga+base+j] = grads[ga+base+j] + t1 - t2
656 j = j + 1
657 }
658 } }
659 i = i + 1
660 }
661 }
662 if op == NFA_EMBED {
663 // scatter-ADD gy rows into dE[ids[t]] (a repeated token id accumulates grads); no grad to ids
664 let e_node: i64 = tape[7*k+1]; let ids: *i64 = tape[7*k+2] as *i64
665 let T: i64 = tape[7*k+3]; let dm: i64 = tape[7*k+4]
666 let offgy: i64 = tape[7*k+6]; let gE: i64 = tape[7*e_node+6]
667 var t: i64 = 0
668 while t < T {
669 let id: i64 = ids[t]
670 var j: i64 = 0
671 while j < dm { grads[gE + id*dm + j] = grads[gE + id*dm + j] + grads[offgy + t*dm + j]; j = j + 1 }
672 t = t + 1
673 }
674 }
675 if op == NFA_SOFTCE_ROWS {
676 // dlogit[t][j] = qmul(gL, (softmax[t][j] - onehot[t][j]) / T) -- the exact fused CE identity (no log)
677 let logits: i64 = tape[7*k+1]; let tgt: *i64 = tape[7*k+2] as *i64
678 let T: i64 = tape[7*logits+3]; let V: i64 = tape[7*logits+4]
679 let offp: i64 = tape[7*logits+5]; let gp: i64 = tape[7*logits+6]
680 let gL: i64 = grads[tape[7*k+6] + 0]
681 var t: i64 = 0
682 while t < T {
683 let base: i64 = t*V
684 var mx: i64 = vals[offp+base]; var j: i64 = 1
685 while j < V { if vals[offp+base+j] > mx { mx = vals[offp+base+j] } j = j + 1 }
686 var sum: i64 = 0; j = 0
687 while j < V { sum = sum + nfa_fxexp(vals[offp+base+j] - mx); j = j + 1 }
688 if sum <= 0 { sum = 1 }
689 let id: i64 = tgt[t]
690 j = 0
691 while j < V {
692 let sm: i64 = (nfa_fxexp(vals[offp+base+j] - mx) << 16) / sum
693 var oh: i64 = 0
694 if j == id { oh = NFA_Q16 }
695 grads[gp+base+j] = grads[gp+base+j] + nfa_qmul(gL, (sm - oh) / T)
696 j = j + 1
697 }
698 t = t + 1
699 }
700 }
701 if op == NFA_SLICE_COLS {
702 let a: i64 = tape[7*k+1]; let c0: i64 = tape[7*k+2]; let T: i64 = tape[7*k+3]; let w: i64 = tape[7*k+4]; let C: i64 = tape[7*a+4]
703 let offgy: i64 = tape[7*k+6]; let ga: i64 = tape[7*a+6]
704 var t: i64 = 0
705 while t < T { var j: i64 = 0; while j < w { grads[ga + t*C + c0 + j] = grads[ga + t*C + c0 + j] + grads[offgy + t*w + j]; j = j + 1 } t = t + 1 }
706 }
707 if op == NFA_CONCAT_COLS {
708 let a: i64 = tape[7*k+1]; let b: i64 = tape[7*k+2]; let T: i64 = tape[7*k+3]; let wc: i64 = tape[7*k+4]
709 let wa: i64 = tape[7*a+4]; let wb: i64 = tape[7*b+4]
710 let offgy: i64 = tape[7*k+6]; let ga: i64 = tape[7*a+6]; let gb: i64 = tape[7*b+6]
711 var t: i64 = 0
712 while t < T {
713 var j: i64 = 0
714 while j < wa { grads[ga + t*wa + j] = grads[ga + t*wa + j] + grads[offgy + t*wc + j]; j = j + 1 }
715 j = 0
716 while j < wb { grads[gb + t*wb + j] = grads[gb + t*wb + j] + grads[offgy + t*wc + wa + j]; j = j + 1 }
717 t = t + 1
718 }
719 }
720 k = k - 1
721 }
722 return 0
723}
724
725// SGD optimizer step on a parameter vector: w[i] -= lr * g[i], lr in Q16. Step = (lr_q*g_q)/Q16 with
726// integer division truncating TOWARD ZERO -> sign-symmetric (no arithmetic-shift floor bias). General
727// (any param vector, any Q16 rate); the same step that scales from this 2-weight demo to the full model.
728func nfa_sgd(w: *i64, g: *i64, n: i64, lr_q: i64) -> i64 {
729 var i: i64 = 0
730 while i < n { w[i] = w[i] - (lr_q * g[i]) / NFA_Q16; i = i + 1 }
731 return 0
732}
733
734// ERROR-FEEDBACK SGD: like nfa_sgd but sub-quantum step remainders ACCUMULATE in resid[n] (persistent,
735// zero-init) instead of being truncated away -- no gradient is ever lost to the /Q16 quantization. This is
736// the integer-training lesson the neural-image arc proved (small systematic gradients otherwise quantize to
737// ZERO updates, and models freeze in degenerate states, e.g. the reader's uniform-collapse absorbing state).
738// acc = resid + lr*g (raw Q32-ish); apply the whole-quantum part; keep the remainder (sign-symmetric).
739func nfa_sgd_ef(w: *i64, g: *i64, resid: *i64, n: i64, lr_q: i64) -> i64 {
740 var i: i64 = 0
741 while i < n {
742 let acc: i64 = resid[i] + lr_q * g[i]
743 let step: i64 = acc / NFA_Q16
744 w[i] = w[i] - step
745 resid[i] = acc - step * NFA_Q16
746 i = i + 1
747 }
748 return 0
749}
750
751// AdamW optimizer step (Q16 integer; port of nx_tgrad_core's ad_step to fixed-point). Per-param adaptive
752// moments + bias-correction -> normalizes each parameter's step by its own gradient scale, so it converges
753// where plain SGD (one global rate) struggles on ill-conditioned losses. m/v are persistent Q16 state (zero-init);
754// t = 1-based step. All integer -> bit-exact/deterministic. sqrt via nfa_isqrt. b1/b2/eps/wd/lr all Q16.
755func nfa_adamw(w: *i64, g: *i64, m: *i64, v: *i64, n: i64, lr: i64, b1: i64, b2: i64, eps: i64, wd: i64, t: i64) -> i64 {
756 var c1: i64 = NFA_Q16; var c2: i64 = NFA_Q16
757 var k: i64 = 0
758 while k < t { c1 = nfa_qmul(c1, b1); c2 = nfa_qmul(c2, b2); k = k + 1 } // b1^t, b2^t
759 let bc1: i64 = NFA_Q16 - c1
760 let bc2: i64 = NFA_Q16 - c2
761 var i: i64 = 0
762 while i < n {
763 m[i] = nfa_qmul(b1, m[i]) + nfa_qmul(NFA_Q16 - b1, g[i])
764 v[i] = nfa_qmul(b2, v[i]) + nfa_qmul(NFA_Q16 - b2, nfa_qmul(g[i], g[i]))
765 var mh: i64 = m[i]; if bc1 > 0 { mh = (m[i] << 16) / bc1 }
766 var vh: i64 = v[i]; if bc2 > 0 { vh = (v[i] << 16) / bc2 }
767 if vh < 0 { vh = 0 }
768 let sq: i64 = nfa_isqrt(vh << 16) // sqrt(vh) in Q16
769 let denom: i64 = sq + eps
770 var upd: i64 = 0; if denom > 0 { upd = (mh << 16) / denom } // mhat/(sqrt(vhat)+eps) in Q16
771 w[i] = w[i] - nfa_qmul(lr, upd) - nfa_qmul(nfa_qmul(lr, wd), w[i])
772 i = i + 1
773 }
774 return 0
775}