1	# bend-ml-tensor-array: tensors over a flat Array<F32>, with the shape in the type.
2	#
3	#   import bend-ml-tensor-array@0.1.11.0/main.bend as TA
4	#
5	# Mat<r, c> stores r*c numbers in an Array<F32>, row by row (index i*c + j). The
6	# dimensions are erased type parameters: the same guarantee as bend-ml-tensor
7	# (a product with wrong dimensions does not compile), but ~50x faster than lists:
8	# Array.get/set by index, and products in parallel blocks of rows.
9	#
10	# An Array is affine (a single owner): every operation that reads an operand returns it
11	# together with the result (as Array.get does), in records MM, AR, ... Use Mat.clone
12	# when you need two copies.
13	#
14	# What is NOT proved: the invariant "the Array has capacity >= r*c" holds because
15	# the constructors (Mat.zeros, Mat.of_list) allocate the right size; and F32 is
16	# validated by tests against PyTorch (reference/test_tensor_array.py), not by proof.
17	
18	import Base
19	
20	# ---------------------------------------------------------------------
21	# helper types
22	# ---------------------------------------------------------------------
23	
24	type R is Type:
25	  R{a: Array<F32>, b: Array<F32>, x: F32}
26	
27	type G is Type:
28	  G{a: Array<F32>, b: Array<F32>, c: Array<F32>}
29	
30	# A dot product with four partial sums. Each step adds its product to a0 and rotates the sums
31	# (a1, a2, a3, a0), so each sum gets every fourth product and the processor can work on the next
32	# addition before the last one has finished; at the end the sums are added in pairs. The order of
33	# the additions is fixed, so the result is the same on every run and every thread count, but it is
34	# not the order of a single running sum (NOTES.md, exp. 17).
35	def dotr(k: Nat, +ia: U32, +ib: U32, +sa: U32, +sb: U32, a0: F32, a1: F32, a2: F32, a3: F32, ra: Array<F32> & F32, rb: Array<F32> & F32) -> R:
36	  match k ra rb:
37	    case 0n Tuple{a, x} Tuple{b, y}:
38	      R{a, b, ((a0 + a1 : F32) + (a2 + a3 : F32) : F32)}
39	    case 1n+p Tuple{a, +x} Tuple{b, +y}:
40	      dotr(p, (ia + sa : U32), (ib + sb : U32), sa, sb, a1, a2, a3, (a0 + (x * y : F32) : F32), Array.get(F32, a, (ia + sa : U32)), Array.get(F32, b, (ib + sb : U32)))
41	
42	def dotg(k: Nat, +ia: U32, +ib: U32, +sa: U32, +sb: U32, acc: F32, ra: Array<F32> & F32, rb: Array<F32> & F32) -> R:
43	  dotr(k, ia, ib, sa, sb, acc, 0.0, 0.0, 0.0, ra, rb)
44	
45	# writes one cell of C: c[idx] = x, or c[idx] = c[idx] + x when ac (an accumulating product, 0.1.10.0). x is computed
46	# exactly as without ac and added once, so C + A·B has the bits of adding a separately computed A·B to C.
47	
48	def cwr_add(+idx: U32, +x: F32, r: Array<F32> & F32) -> Array<F32>:
49	  match r:
50	    case Tuple{c, +v}:
51	      Array.set(F32, c, idx, (v + x : F32))
52	
53	def cwr(ac: Bool, +idx: U32, +x: F32, c: Array<F32>) -> Array<F32>:
54	  match ac:
55	    case True{}:
56	      cwr_add(idx, x, Array.get(F32, c, idx))
57	    case False{}:
58	      Array.set(F32, c, idx, x)
59	
60	# row i of C, columns j.. ; r is the product of the current column (already computed)
61	
62	def gcols(+ac: Bool, mleft: Nat, +kn: Nat, +ia0: U32, +ib0: U32, +sa: U32, +sb: U32, +bcol: U32, +cidx: U32, c: Array<F32>, r: R) -> G:
63	  match mleft r:
64	    case 0n R{a, b, +x}:
65	      G{a, b, cwr(ac, cidx, x, c)}
66	    case 1n+p R{a, b, +x}:
67	      gcols(ac, p, kn, ia0, (ib0 + bcol : U32), sa, sb, bcol, (cidx + 1 : U32), cwr(ac, cidx, x, c), dotg(kn, ia0, (ib0 + bcol : U32), sa, sb, 0.0, Array.get(F32, a, ia0), Array.get(F32, b, (ib0 + bcol : U32))))
68	
69	def grows(+ac: Bool, nleft: Nat, +mm1: Nat, +kn: Nat, +i: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +mu: U32, g: G) -> G:
70	  match nleft g:
71	    case 0n G{a, b, c}:
72	      G{a, b, c}
73	    case 1n+p G{a, b, c}:
74	      grows(ac, p, mm1, kn, (i + 1 : U32), arow, sa, sb, bcol, mu, gcols(ac, mm1, kn, (i * arow : U32), 0, sa, sb, bcol, (i * mu : U32), c, dotg(kn, (i * arow : U32), 0, sa, sb, 0.0, Array.get(F32, a, (i * arow : U32)), Array.get(F32, b, 0))))
75	
76	# C (n x m) = A · B with k steps, or C + A · B when ac; c is the output array (with capacity n*m)
77	
78	def gemm_ac(+ac: Bool, +n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
79	  grows(ac, n, Nat.sub(m, 1n), k, 0, arow, sa, sb, bcol, U32.from_nat(m), G{a, b, c})
80	
81	def gemm(+n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
82	  gemm_ac(False{}, n, m, k, arow, sa, sb, bcol, a, b, c)
83	
84	
85	# ---------------------------------------------------------------------
86	# parallel gemm: 2^d blocks of rows of C, each with its own copy of A and B,
87	# returning lists that are written into C at the end
88	# ---------------------------------------------------------------------
89	
90	type S3 is Type:
91	  S3{a: Array<F32>, b: Array<F32>, out: List<&2, F32>}
92	
93	def lcols(mleft: Nat, +kn: Nat, +ia0: U32, +ib0: U32, +sa: U32, +sb: U32, +bcol: U32, out: List<&2, F32>, r: R) -> S3:
94	  match mleft r:
95	    case 0n R{a, b, +x}:
96	      S3{a, b, x <> out}
97	    case 1n+p R{a, b, +x}:
98	      lcols(p, kn, ia0, (ib0 + bcol : U32), sa, sb, bcol, x <> out, dotg(kn, ia0, (ib0 + bcol : U32), sa, sb, 0.0, Array.get(F32, a, ia0), Array.get(F32, b, (ib0 + bcol : U32))))
99	
100	def lrows(nleft: Nat, +mm1: Nat, +kn: Nat, +i: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, st: S3) -> S3:
101	  match nleft st:
102	    case 0n S3{a, b, out}:
103	      S3{a, b, out}
104	    case 1n+p S3{a, b, out}:
105	      lrows(p, mm1, kn, (i + 1 : U32), arow, sa, sb, bcol, lcols(mm1, kn, (i * arow : U32), 0, sa, sb, bcol, out, dotg(kn, (i * arow : U32), 0, sa, sb, 0.0, Array.get(F32, a, (i * arow : U32)), Array.get(F32, b, 0))))
106	
107	def lblock_fin(st: S3) -> List<&2, F32>:
108	  match st:
109	    case S3{a, b, out}:
110	      List.reverse(&2, F32, out)
111	
112	def lblock(cnt: Nat, +mm1: Nat, +kn: Nat, +i: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>) -> List<&2, F32>:
113	  lblock_fin(lrows(cnt, mm1, kn, i, arow, sa, sb, bcol, S3{a, b, Nil{}}))
114	
115	def lpar(d: Nat, +cnt: Nat, +mm1: Nat, +kn: Nat, +i: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, ra: Array<F32> & Array<F32>, rb: Array<F32> & Array<F32>) -> List<&2, F32>:
116	  match d ra rb:
117	    case 0n Tuple{a, a2} Tuple{b, b2}:
118	      lblock(cnt, mm1, kn, i, arow, sa, sb, bcol, a, b)
119	    case 1n++q Tuple{a, a2} Tuple{b, b2}:
120	      x y = lpar(q, Nat.div(cnt, 2n), mm1, kn, i, arow, sa, sb, bcol, Array.clone(F32, a), Array.clone(F32, b)) lpar(q, Nat.sub(cnt, Nat.div(cnt, 2n)), mm1, kn, (i + U32.from_nat(Nat.div(cnt, 2n)) : U32), arow, sa, sb, bcol, Array.clone(F32, a2), Array.clone(F32, b2))
121	      List.append(&2, F32, x, y)
122	
123	def write_l(xs: List<&2, F32>, +idx: U32, c: Array<F32>) -> Array<F32>:
124	  match xs:
125	    case Nil{}:
126	      c
127	    case Con{h, t}:
128	      write_l(t, (idx + 1 : U32), Array.set(F32, c, idx, h))
129	
130	# c[idx..] = list, or c[idx..] += list when ac
131	
132	def write_la(+ac: Bool, xs: List<&2, F32>, +idx: U32, c: Array<F32>) -> Array<F32>:
133	  match xs:
134	    case Nil{}:
135	      c
136	    case Con{+h, t}:
137	      write_la(ac, t, (idx + 1 : U32), cwr(ac, idx, h, c))
138	
139	def gp3(+ac: Bool, a: Array<F32>, b: Array<F32>, c: Array<F32>, xs: List<&2, F32>) -> G:
140	  G{a, b, write_la(ac, xs, 0, c)}
141	
142	def gp2(d: Nat, +n: Nat, +mm1: Nat, +kn: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, c: Array<F32>, pa: Array<F32> & Array<F32>, pb: Array<F32> & Array<F32>) -> G:
143	  match pa pb:
144	    case Tuple{a1, a2} Tuple{b1, b2}:
145	      gp3(False{}, a1, b1, c, lpar(d, n, mm1, kn, 0, arow, sa, sb, bcol, Array.clone(F32, a2), Array.clone(F32, b2)))
146	
147	# C (n x m) = A · B, with 2^d blocks of rows in parallel; returns A, B (intact) and C
148	
149	def gemm_par(d: Nat, +n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
150	  gp2(d, n, Nat.sub(m, 1n), k, arow, sa, sb, bcol, c, Array.clone(F32, a), Array.clone(F32, b))
151	
152	
153	# ---------------------------------------------------------------------
154	# parallel gemm by COLUMNS, for n = 1 (matrix · vector): each block computes a
155	# range of columns of C, with its own copy of the operands
156	# ---------------------------------------------------------------------
157	
158	def cblock_fin(st: S3) -> List<&2, F32>:
159	  match st:
160	    case S3{a, b, out}:
161	      List.reverse(&2, F32, out)
162	
163	# cnt columns starting at j0 (cnt >= 1); row 0 of A is at ia0
164	
165	def cblock_go(cnt1: Nat, +kn: Nat, +ia0: U32, +j0: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>) -> List<&2, F32>:
166	  cblock_fin(lcols(cnt1, kn, ia0, (j0 * bcol : U32), sa, sb, bcol, Nil{}, dotg(kn, ia0, (j0 * bcol : U32), sa, sb, 0.0, Array.get(F32, a, ia0), Array.get(F32, b, (j0 * bcol : U32)))))
167	
168	def cblock(cnt: Nat, +kn: Nat, +ia0: U32, +j0: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>) -> List<&2, F32>:
169	  match cnt:
170	    case 0n:
171	      Nil{}
172	    case 1n+p:
173	      cblock_go(p, kn, ia0, j0, sa, sb, bcol, a, b)
174	
175	def cpar(d: Nat, +cnt: Nat, +kn: Nat, +ia0: U32, +j0: U32, +sa: U32, +sb: U32, +bcol: U32, ra: Array<F32> & Array<F32>, rb: Array<F32> & Array<F32>) -> List<&2, F32>:
176	  match d ra rb:
177	    case 0n Tuple{a, a2} Tuple{b, b2}:
178	      cblock(cnt, kn, ia0, j0, sa, sb, bcol, a, b)
179	    case 1n++q Tuple{a, a2} Tuple{b, b2}:
180	      x y = cpar(q, Nat.div(cnt, 2n), kn, ia0, j0, sa, sb, bcol, Array.clone(F32, a), Array.clone(F32, b)) cpar(q, Nat.sub(cnt, Nat.div(cnt, 2n)), kn, ia0, (j0 + U32.from_nat(Nat.div(cnt, 2n)) : U32), sa, sb, bcol, Array.clone(F32, a2), Array.clone(F32, b2))
181	      List.append(&2, F32, x, y)
182	
183	def cp2(+ac: Bool, d: Nat, +m: Nat, +kn: Nat, +sa: U32, +sb: U32, +bcol: U32, c: Array<F32>, pa: Array<F32> & Array<F32>, pb: Array<F32> & Array<F32>) -> G:
184	  match pa pb:
185	    case Tuple{a1, a2} Tuple{b1, b2}:
186	      gp3(ac, a1, b1, c, cpar(d, m, kn, 0, 0, sa, sb, bcol, Array.clone(F32, a2), Array.clone(F32, b2)))
187	
188	# C (1 x m) = A (1 x k) · B, with 2^d blocks of columns in parallel
189	
190	def gemm_par_cols(+ac: Bool, d: Nat, +m: Nat, +k: Nat, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
191	  cp2(ac, d, m, k, sa, sb, bcol, c, Array.clone(F32, a), Array.clone(F32, b))
192	
193	# ---------------------------------------------------------------------
194	# element-wise operations (in place, by index)
195	# ---------------------------------------------------------------------
196	
197	# ---------------------------------------------------------------------
198	# gemm: C[i][j] = sum_p A[i*arow + p*sa] * B[p*sb + j*bcol]
199	# (nn: arow=K, sa=1, sb=M, bcol=1 | nt: B stored m x k: sb=1, bcol=K | tn: A stored k x n: arow=1, sa=N)
200	# ---------------------------------------------------------------------
201	
202	# acc += A[ia + t*sa] * B[ib + t*sb] for k steps; ra, rb already carry the first pair
203	
204	type P2 is Type:
205	  P2{a: Array<F32>, b: Array<F32>}
206	
207	type AL is Type:
208	  AL{c: Array<F32>, l: List<&2, F32>}
209	
210	# list of m elements starting at idx (reads with Array.get, returning the array)
211	
212	def read_l(mleft: Nat, +idx: U32, acc: List<&2, F32>, rc: Array<F32> & F32) -> AL:
213	  match mleft rc:
214	    case 0n Tuple{c, +v}:
215	      AL{c, List.reverse(&2, F32, v <> acc)}
216	    case 1n+p Tuple{c, +v}:
217	      read_l(p, (idx + 1 : U32), v <> acc, Array.get(F32, c, (idx + 1 : U32)))
218	
219	# writes the list starting at idx
220	
221	def addl(xs: List<&2, F32>, +idx: U32, rc: Array<F32> & F32) -> Array<F32>:
222	  match xs rc:
223	    case Nil{} Tuple{c, v}:
224	      c
225	    case Con{+h, t} Tuple{c, +v}:
226	      addl(t, (idx + 1 : U32), Array.get(F32, Array.set(F32, c, idx, (v + h : F32)), (idx + 1 : U32)))
227	
228	# adds a bias (a list of m) to each of the n rows of c (n x m)
229	
230	def bias_rows(nleft: Nat, +i: U32, +mu: U32, +bl: List<&2, F32>, c: Array<F32>) -> Array<F32>:
231	  match nleft:
232	    case 0n:
233	      c
234	    case 1n+p:
235	      bias_rows(p, (i + 1 : U32), mu, bl, addl(bl, (i * mu : U32), Array.get(F32, c, (i * mu : U32))))
236	
237	# in-place relu on the first n elements
238	
239	def relu_l(nleft: Nat, +t: U32, rc: Array<F32> & F32) -> Array<F32>:
240	  match nleft rc:
241	    case 0n Tuple{c, v}:
242	      c
243	    case 1n+p Tuple{c, +v}:
244	      relu_l(p, (t + 1 : U32), Array.get(F32, Array.set(F32, c, t, F32.max(v, 0.0)), (t + 1 : U32)))
245	
246	def pos2(b: Bool) -> F32:
247	  match b:
248	    case True{}:
249	      1.0
250	    case False{}:
251	      0.0
252	
253	def pos(+y: F32) -> F32:
254	  pos2(F32.is_gt(y, 0.0))
255	
256	# d[t] *= (h[t] > 0): the relu gradient; returns both arrays
257	
258	def mask_l(nleft: Nat, +t: U32, rd: Array<F32> & F32, rh: Array<F32> & F32) -> P2:
259	  match nleft rd rh:
260	    case 0n Tuple{d, x} Tuple{h, y}:
261	      P2{d, h}
262	    case 1n+p Tuple{d, +x} Tuple{h, +y}:
263	      mask_l(p, (t + 1 : U32), Array.get(F32, Array.set(F32, d, t, (x * pos(y) : F32)), (t + 1 : U32)), Array.get(F32, h, (t + 1 : U32)))
264	
265	# w[t] -= lr * dw[t]
266	
267	def maxl2(xs: List<&2, F32>, cur: F32) -> F32:
268	  match xs:
269	    case Nil{}:
270	      cur
271	    case Con{+h, t}:
272	      maxl2(t, F32.max(cur, h))
273	
274	def maxl(xs: List<&2, F32>) -> F32:
275	  match xs:
276	    case Nil{}:
277	      0.0
278	    case Con{+h, t}:
279	      maxl2(t, h)
280	
281	def expl(xs: List<&2, F32>, +m: F32) -> List<&2, F32>:
282	  match xs:
283	    case Nil{}:
284	      Nil{}
285	    case Con{+h, t}:
286	      F32.exp((h - m : F32)) <> expl(t, m)
287	
288	def suml(xs: List<&2, F32>, acc: F32) -> F32:
289	  match xs:
290	    case Nil{}:
291	      acc
292	    case Con{+h, t}:
293	      suml(t, (acc + h : F32))
294	
295	def divl(xs: List<&2, F32>, +s: F32) -> List<&2, F32>:
296	  match xs:
297	    case Nil{}:
298	      Nil{}
299	    case Con{+h, t}:
300	      (h / s : F32) <> divl(t, s)
301	
302	def softmax2(+es: List<&2, F32>) -> List<&2, F32>:
303	  divl(es, suml(es, 0.0))
304	
305	def softmax(+xs: List<&2, F32>) -> List<&2, F32>:
306	  softmax2(expl(xs, maxl(xs)))
307	
308	# p[label]
309	
310	def sgd_l(nleft: Nat, +t: U32, +lr: F32, rw: Array<F32> & F32, rd: Array<F32> & F32) -> P2:
311	  match nleft rw rd:
312	    case 0n Tuple{w, x} Tuple{d, y}:
313	      P2{w, d}
314	    case 1n+p Tuple{w, +x} Tuple{d, +y}:
315	      sgd_l(p, (t + 1 : U32), lr, Array.get(F32, Array.set(F32, w, t, (x - (lr * y : F32) : F32)), (t + 1 : U32)), Array.get(F32, d, (t + 1 : U32)))
316	
317	# ---------------------------------------------------------------------
318	# softmax + cross-entropy per row (lists of 10 elements: negligible cost)
319	# ---------------------------------------------------------------------
320	
321	def pick(ps: List<&2, F32>, n: Nat) -> F32:
322	  match ps n:
323	    case Con{+h, t} 0n:
324	      h
325	    case Con{h, t} 1n+q:
326	      pick(t, q)
327	    case Nil{} _:
328	      0.0
329	
330	# (p - onehot(label)) / n
331	
332	def hit(b: Bool) -> F32:
333	  match b:
334	    case True{}:
335	      1.0
336	    case False{}:
337	      0.0
338	
339	def grad_l(ps: List<&2, F32>, +label: Nat, +j: Nat, +nf: F32) -> List<&2, F32>:
340	  match ps:
341	    case Nil{}:
342	      Nil{}
343	    case Con{+h, t}:
344	      (((h - hit(Nat.is_eq(label, j)) : F32) / nf) : F32) <> grad_l(t, label, 1n+j, nf)
345	
346	type L2 is Type:
347	  L2{loss: F32, c: Array<F32>}
348	
349	def ce_row2(+ps: List<&2, F32>, +label: Nat, +nf: F32, +idx: U32, loss: F32, c: Array<F32>) -> L2:
350	  L2{(loss - F32.log(F32.max(pick(ps, label), 0.000000000001)) : F32), write_l(grad_l(ps, label, 0n, nf), idx, c)}
351	
352	def ce_row(al: AL, +label: Nat, +nf: F32, +idx: U32, loss: F32) -> L2:
353	  match al:
354	    case AL{c, l}:
355	      ce_row2(softmax(l), label, nf, idx, loss, c)
356	
357	# walks the rows of z2 (n x m): r carries the accumulated loss and the array; at the end the array has become the gradient
358	
359	type Mat<-r: Nat, -c: Nat> is Type:
360	  Mat{d: Array<F32>}
361	
362	def ce_rows(labels: List<&2, Nat>, +mm1: Nat, +mu: U32, +nf: F32, +i: U32, r: L2) -> L2:
363	  match labels r:
364	    case Nil{} L2{loss, c}:
365	      L2{loss, c}
366	    case Con{+y, t} L2{loss, c}:
367	      ce_rows(t, mm1, mu, nf, (i + 1 : U32), ce_row(read_l(mm1, (i * mu : U32), Nil{}, Array.get(F32, c, (i * mu : U32))), y, nf, (i * mu : U32), loss))
368	
369	# ---- helper: hits per row ----
370	
371	type CR is Type:
372	  CR{c: Array<F32>, n: Nat}
373	
374	def gt_head(xs: List<&2, F32>, +v: F32) -> Bool:
375	  match xs:
376	    case Nil{}:
377	      False{}
378	    case Con{+h, t}:
379	      F32.is_gt(h, v)
380	
381	def argmax2(xs: List<&2, F32>, g: Bool, +i: Nat, +bi: Nat, +bv: F32) -> Nat:
382	  match xs g:
383	    case Nil{} _:
384	      bi
385	    case Con{+h, +t} True{}:
386	      argmax2(t, gt_head(t, h), 1n+i, i, h)
387	    case Con{h, +t} False{}:
388	      argmax2(t, gt_head(t, bv), 1n+i, bi, bv)
389	
390	def argmax(xs: List<&2, F32>) -> Nat:
391	  match xs:
392	    case Nil{}:
393	      0n
394	    case Con{+h, +t}:
395	      argmax2(t, gt_head(t, h), 1n, 0n, h)
396	
397	def one(b: Bool) -> Nat:
398	  match b:
399	    case True{}:
400	      1n
401	    case False{}:
402	      0n
403	
404	def count_rows3(+y: Nat, n: Nat, al: AL) -> CR:
405	  match al:
406	    case AL{c, l}:
407	      CR{c, Nat.add(n, one(Nat.is_eq(argmax(l), y)))}
408	
409	def count_rows(labels: List<&2, Nat>, +mm1: Nat, +mu: U32, +i: U32, r: CR) -> CR:
410	  match labels r:
411	    case Nil{} CR{c, n}:
412	      CR{c, n}
413	    case Con{+y, t} CR{c, n}:
414	      count_rows(t, mm1, mu, (i + 1 : U32), count_rows3(y, n, read_l(mm1, (i * mu : U32), Nil{}, Array.get(F32, c, (i * mu : U32)))))
415	
416	
417	# =====================================================================
418	# Typed API: Mat<r, c> (the dimensions enter the type)
419	# =====================================================================
420	
421	# =====================================================================
422	# Array capacity: a PROVED depth for the number of slots
423	# =====================================================================
424	#
425	# cap_depth(n) is the d such that an Array with 2^d slots holds n elements. It is computed by
426	# structural recursion on n (not with U32.log2) so that the law cap_ok below can prove
427	# `n <= 2^cap_depth(n)`. What stays trusted is only that Array.new(T, d, v) allocates 2^d slots.
428	
429	# a type chosen by a Bool: rewriting through it refutes True == False
430	def BD(b: Bool, t: Type, f: Type) -> Type:
431	  match b:
432	    case True{}:
433	      t
434	    case False{}:
435	      f
436	
437	# a <= b, by direct recursion
438	def leb(a: Nat, b: Nat) -> Bool:
439	  match a b:
440	    case 0n _:
441	      True{}
442	    case 1n+p 0n:
443	      False{}
444	    case 1n+p 1n+q:
445	      leb(p, q)
446	
447	def pow2(d: Nat) -> Nat:
448	  match d:
449	    case 0n:
450	      1n
451	    case 1n++q:
452	      Nat.add(pow2(q), pow2(q))
453	
454	# smallest power of two >= n, by structural recursion on n: returns (depth, 2^depth)
455	def cap_pick(ok: Bool, d: Nat, +p: Nat) -> Nat & Nat:
456	  match ok:
457	    case True{}:
458	      (d, p)
459	    case False{}:
460	      (1n+d, Nat.add(p, p))
461	
462	def cap_step(r: Nat & Nat, +n: Nat) -> Nat & Nat:
463	  match r:
464	    case Tuple{+d, +p}:
465	      cap_pick(leb(n, p), d, p)
466	
467	def cap_pair(n: Nat) -> Nat & Nat:
468	  match n:
469	    case 0n:
470	      (0n, 1n)
471	    case 1n++m:
472	      cap_step(cap_pair(m), 1n+m)
473	
474	def fst_nat(r: Nat & Nat) -> Nat:
475	  match r:
476	    case Tuple{d, p}:
477	      d
478	
479	def snd_nat(r: Nat & Nat) -> Nat:
480	  match r:
481	    case Tuple{d, p}:
482	      p
483	
484	# the same pair by doubling (0.1.10.0): p doubles until n <= p, log2(n) steps in a tail loop, where cap_pair takes n
485	# steps of a recursion that is not a tail call, each with a comparison of up to n steps. Its fuel is n, which always
486	# suffices (2^n >= n); if it ran out, the answer would be cap_pair's, so the law below needs no arithmetic about how
487	# many doublings there are. ok is leb(n, p), computed by the caller (match only inspects a parameter).
488	def cap_dbl(f: Nat, ok: Bool, +n: Nat, +d: Nat, +p: Nat) -> Nat & Nat:
489	  match f ok:
490	    case _ True{}:
491	      (d, p)
492	    case 0n False{}:
493	      cap_pair(n)
494	    case 1n+g False{}:
495	      cap_dbl(g, leb(n, Nat.add(p, p)), n, 1n+d, Nat.add(p, p))
496	
497	# the depth d such that an Array with 2^d slots holds n elements (the smallest such d)
498	def cap_depth(+n: Nat) -> Nat:
499	  fst_nat(cap_dbl(n, leb(n, 1n), n, 0n, 1n))
500	
501	# ---- lemmas ----
502	
503	# a <= b  ->  a <= 1 + b
504	def leb_succ(a: Nat, b: Nat, e: {True{} == leb(a, b) : Bool}) -> {True{} == leb(a, 1n+b) : Bool}:
505	  match a b:
506	    case 0n _:
507	      {==}
508	    case 1n+p 0n:
509	      %e : BD(_, Unit, {True{} == leb(1n+p, 1n) : Bool})
510	      Unit{}
511	    case 1n+p 1n+q:
512	      leb_succ(p, q, e)
513	
514	# a <= b  ->  a <= c + b
515	def leb_up_l(+a: Nat, +b: Nat, c: Nat, +e: {True{} == leb(a, b) : Bool}) -> {True{} == leb(a, Nat.add(c, b)) : Bool}:
516	  match c:
517	    case 0n:
518	      e
519	    case 1n++c2:
520	      leb_succ(a, Nat.add(c2, b), leb_up_l(a, b, c2, e))
521	
522	# a <= b  ->  a <= b + c
523	def leb_up(a: Nat, b: Nat, c: Nat, e: {True{} == leb(a, b) : Bool}) -> {True{} == leb(a, Nat.add(b, c)) : Bool}:
524	  match a b:
525	    case 0n _:
526	      {==}
527	    case 1n+p 0n:
528	      %e : BD(_, Unit, {True{} == leb(1n+p, Nat.add(0n, c)) : Bool})
529	      Unit{}
530	    case 1n+p 1n+q:
531	      leb_up(p, q, c, e)
532	
533	# a <= b  and  c <= d  ->  a + c <= b + d
534	def leb_add(a: Nat, b: Nat, c: Nat, d: Nat, e1: {True{} == leb(a, b) : Bool}, e2: {True{} == leb(c, d) : Bool}) -> {True{} == leb(Nat.add(a, c), Nat.add(b, d)) : Bool}:
535	  match a b:
536	    case 0n _:
537	      leb_up_l(c, d, b, e2)
538	    case 1n+p 0n:
539	      %e1 : BD(_, Unit, {True{} == leb(Nat.add(1n+p, c), Nat.add(0n, d)) : Bool})
540	      Unit{}
541	    case 1n+p 1n+q:
542	      leb_add(p, q, c, d, e1, e2)
543	
544	# 1 <= 2^d
545	def pow2_pos(d: Nat) -> {True{} == leb(1n, pow2(d)) : Bool}:
546	  match d:
547	    case 0n:
548	      {==}
549	    case 1n++q:
550	      leb_up(1n, pow2(q), pow2(q), pow2_pos(q))
551	
552	# 1 <= p when p = 2^d
553	def pos_p(d: Nat, p: Nat, i1: {p == pow2(d) : Nat}) -> {True{} == leb(1n, p) : Bool}:
554	  %Equal.sym(Nat, p, pow2(d), i1) : {True{} == leb(1n, _) : Bool}
555	  pow2_pos(d)
556	
557	# p + p == 2^(1+d) when p == 2^d
558	def step_eq(d: Nat, p: Nat, i1: {p == pow2(d) : Nat}) -> {Nat.add(p, p) == pow2(1n+d) : Nat}:
559	  %i1 : {Nat.add(p, p) == Nat.add(_, _) : Nat}
560	  {==}
561	
562	# one step: given the invariant for m, the invariant for 1+m holds whichever way the comparison goes
563	def step_ok(ok: Bool, +d: Nat, +p: Nat, +m: Nat, +e: {ok == leb(1n+m, p) : Bool}, +i1: {p == pow2(d) : Nat}, +i2: {True{} == leb(m, p) : Bool}) -> {snd_nat(cap_pick(ok, d, p)) == pow2(fst_nat(cap_pick(ok, d, p))) : Nat} & {True{} == leb(1n+m, snd_nat(cap_pick(ok, d, p))) : Bool}:
564	  match ok:
565	    case True{}:
566	      (i1, e)
567	    case False{}:
568	      (step_eq(d, p, i1), leb_add(1n, p, m, p, pos_p(d, p, i1), i2))
569	
570	
571	# the invariant for cap_step on a pair, opened
572	def step_pair_ok(+m: Nat, r: Nat & Nat, i1: {snd_nat(r) == pow2(fst_nat(r)) : Nat}, i2: {True{} == leb(m, snd_nat(r)) : Bool}) -> {snd_nat(cap_step(r, 1n+m)) == pow2(fst_nat(cap_step(r, 1n+m))) : Nat} & {True{} == leb(1n+m, snd_nat(cap_step(r, 1n+m))) : Bool}:
573	  match r:
574	    case Tuple{+d, +p}:
575	      step_ok(leb(1n+m, p), d, p, m, {==}, i1, i2)
576	
577	def cap_succ(+m: Nat, ih: {snd_nat(cap_pair(m)) == pow2(fst_nat(cap_pair(m))) : Nat} & {True{} == leb(m, snd_nat(cap_pair(m))) : Bool}) -> {snd_nat(cap_pair(1n+m)) == pow2(fst_nat(cap_pair(1n+m))) : Nat} & {True{} == leb(1n+m, snd_nat(cap_pair(1n+m))) : Bool}:
578	  match ih:
579	    case Tuple{i1, i2}:
580	      step_pair_ok(m, cap_pair(m), i1, i2)
581	
582	# the invariant of cap_pair: the second component is 2^(first) and is >= n
583	def cap_inv(n: Nat) -> {snd_nat(cap_pair(n)) == pow2(fst_nat(cap_pair(n))) : Nat} & {True{} == leb(n, snd_nat(cap_pair(n))) : Bool}:
584	  match n:
585	    case 0n:
586	      ({==}, {==})
587	    case 1n++m:
588	      cap_succ(m, cap_inv(m))
589	
590	# the invariant of cap_dbl: whichever way it ends (n <= p found, or the fuel out and cap_pair's answer), the second
591	# component is 2^(first) and is >= n
592	def dbl_inv(f: Nat, ok: Bool, +n: Nat, +d: Nat, +p: Nat, e: {ok == leb(n, p) : Bool}, i1: {p == pow2(d) : Nat}) -> {snd_nat(cap_dbl(f, ok, n, d, p)) == pow2(fst_nat(cap_dbl(f, ok, n, d, p))) : Nat} & {True{} == leb(n, snd_nat(cap_dbl(f, ok, n, d, p))) : Bool}:
593	  match f ok:
594	    case _ True{}:
595	      (i1, e)
596	    case 0n False{}:
597	      cap_inv(n)
598	    case 1n+g False{}:
599	      dbl_inv(g, leb(n, Nat.add(p, p)), n, 1n+d, Nat.add(p, p), {==}, step_eq(d, p, i1))
600	
601	def dbl_final(n: Nat, pr: {snd_nat(cap_dbl(n, leb(n, 1n), n, 0n, 1n)) == pow2(fst_nat(cap_dbl(n, leb(n, 1n), n, 0n, 1n))) : Nat} & {True{} == leb(n, snd_nat(cap_dbl(n, leb(n, 1n), n, 0n, 1n))) : Bool}) -> {True{} == leb(n, pow2(cap_depth(n))) : Bool}:
602	  match pr:
603	    case Tuple{i1, i2}:
604	      %i1 : {True{} == leb(n, _) : Bool}
605	      i2
606	
607	# LAW (cap_ok): an Array with 2^cap_depth(n) slots always has room for n elements.
608	law cap_ok:
609	  for +n: Nat
610	  {True{} == leb(n, pow2(cap_depth(n))) : Bool}
611	
612	def cap_ok(n):
613	  dbl_final(n, dbl_inv(n, leb(n, 1n), n, 0n, 1n, {==}, {==}))
614	
615	# r x c of zeros, with the exact capacity (smallest power of 2)
616	
617	def Mat.zeros(+r: Nat, +c: Nat) -> Mat<r, c>:
618	  Mat{Array.new(F32, cap_depth(Nat.mul(r, c)), 0.0)}
619	
620	def Mat.fill(+r: Nat, +c: Nat, +v: F32) -> Mat<r, c>:
621	  Mat{Array.new(F32, cap_depth(Nat.mul(r, c)), v)}
622	
623	def fill_list(xs: List<&2, F32>, +i: U32, a: Array<F32>) -> Array<F32>:
624	  match xs:
625	    case Nil{}:
626	      a
627	    case Con{h, t}:
628	      fill_list(t, (i + 1 : U32), Array.set(F32, a, i, h))
629	
630	def len_is(xs: List<&2, F32>, n: Nat) -> Bool:
631	  match xs n:
632	    case Nil{} 0n:
633	      True{}
634	    case Con{h, t} 1n+p:
635	      len_is(t, p)
636	    case _ _:
637	      False{}
638	
639	def Mat.of_list2(-r: Nat, -c: Nat, +rr: Nat, +cc: Nat, ok: Bool, xs: List<&2, F32>) -> Maybe<&1, Mat<r, c>>:
640	  match ok:
641	    case True{}:
642	      Some{Mat{fill_list(xs, 0, Array.new(F32, cap_depth(Nat.mul(rr, cc)), 0.0))}}
643	    case False{}:
644	      None{}
645	
646	# Mat<r, c> from r*c numbers in row order; None if the size does not match
647	
648	def Mat.of_list(+r: Nat, +c: Nat, +xs: List<&2, F32>) -> Maybe<&1, Mat<r, c>>:
649	  Mat.of_list2(r, c, r, c, len_is(xs, Nat.mul(r, c)), xs)
650	
651	# reads the first n numbers, returning the matrix
652	
653	# Mat<r, c> without checking the list size (a short list leaves the rest at zero, a long one is cut
654	# by the capacity): only for when the size is known by construction. Prefer Mat.of_list.
655	
656	def Mat.from_list(+r: Nat, +c: Nat, xs: List<&2, F32>) -> Mat<r, c>:
657	  Mat{fill_list(xs, 0, Array.new(F32, cap_depth(Nat.mul(r, c)), 0.0))}
658	
659	# writes the numbers of xs into m starting at flat index i (row-major); numbers past the
660	# capacity are dropped. For loading a large matrix block by block without building one big list.
661	def Mat.fill_at(-r: Nat, -c: Nat, +i: U32, xs: List<&2, F32>, m: Mat<r, c>) -> Mat<r, c>:
662	  match m:
663	    case Mat{d}:
664	      Mat{fill_list(xs, i, d)}
665	
666	type ML<-r: Nat, -c: Nat> is Type:
667	  ML{m: Mat<r, c>, l: List<&2, F32>}
668	
669	def Mat.to_list2(-r: Nat, -c: Nat, al: AL) -> ML<r, c>:
670	  match al:
671	    case AL{d, l}:
672	      ML{Mat{d}, l}
673	
674	# two independent copies
675	
676	def Mat.to_list(+r: Nat, +c: Nat, m: Mat<r, c>) -> ML<r, c>:
677	  match m:
678	    case Mat{d}:
679	      Mat.to_list2(r, c, read_l(Nat.sub(Nat.mul(r, c), 1n), 0, Nil{}, Array.get(F32, d, 0)))
680	
681	type MM2<-r: Nat, -c: Nat> is Type:
682	  MM2{a: Mat<r, c>, b: Mat<r, c>}
683	
684	def Mat.clone2(-r: Nat, -c: Nat, p: Array<F32> & Array<F32>) -> MM2<r, c>:
685	  match p:
686	    case Tuple{x, y}:
687	      MM2{Mat{x}, Mat{y}}
688	
689	# ---- products: C (n x m) = A · B, with 2^par blocks of rows in parallel ----
690	
691	# C = A(n x k) · B(k x m)
692	
693	def Mat.clone(-r: Nat, -c: Nat, m: Mat<r, c>) -> MM2<r, c>:
694	  match m:
695	    case Mat{d}:
696	      Mat.clone2(r, c, Array.clone(F32, d))
697	
698	type MMul<-n: Nat, -k: Nat, -m: Nat> is Type:
699	  MMul{a: Mat<n, k>, b: Mat<k, m>, c: Mat<n, m>}
700	
701	# C = A(n x k) · Bᵀ, with B stored m x k
702	
703	type MMulNT<-n: Nat, -k: Nat, -m: Nat> is Type:
704	  MMulNT{a: Mat<n, k>, b: Mat<m, k>, c: Mat<n, m>}
705	
706	# C = Aᵀ · B, with A stored k x n
707	
708	type MMulTN<-n: Nat, -k: Nat, -m: Nat> is Type:
709	  MMulTN{a: Mat<k, n>, b: Mat<k, m>, c: Mat<n, m>}
710	
711	# ---------------------------------------------------------------------
712	# gemm in tiles of 8 (0.1.8.0): for products with 2 or more rows, A's rows and B's columns are packed once into cells of
713	# 8 F32 (a row of k numbers is ceil(k/8) cells, zeros after the last), so the inner loop reads 8 numbers per array
714	# access instead of one. The order of the additions is that of `dotr`: product p of a dot goes to partial sum p mod 4,
715	# in order, and the four sums are combined as dotr combines them after k steps (its rotation depends on k mod 4). The
716	# padding adds +0 to a sum that starts at +0, which changes nothing, so the results are the same bits as the scalar
717	# loop's. Products with one row (matrix · vector) keep the scalar loop: packing B would cost as much as the product.
718	# ---------------------------------------------------------------------
719	
720	type T8 is Data:
721	  T8{v0: F32, v1: F32, v2: F32, v3: F32, v4: F32, v5: F32, v6: F32, v7: F32}
722	
723	def t8_zero() -> T8:
724	  T8{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}
725	
726	def pick0(b: Bool, +v: F32) -> F32:
727	  match b:
728	    case True{}:
729	      v
730	    case False{}:
731	      0.0
732	
733	type PK is Type:
734	  PK{a: Array<F32>, t: Array<T8>}
735	
736	# one cell: the 8 numbers of a at idx, idx + st, ..., those at positions q0 + q >= lim (past the end of the line) as 0
737	def pk_g7(+v0: F32, +v1: F32, +v2: F32, +v3: F32, +v4: F32, +v5: F32, +v6: F32, +q0: U32, +lim: U32, +idx: U32, +st: U32, +o: U32, t: Array<T8>, r: Array<F32> & F32) -> PK:
738	  match r:
739	    case Tuple{a, +v}:
740	      PK{a, Array.set(T8, t, o, T8{v0, v1, v2, v3, v4, v5, v6, pick0(U32.is_lt((q0 + 7 : U32), lim), v)})}
741	
742	def pk_g6(+v0: F32, +v1: F32, +v2: F32, +v3: F32, +v4: F32, +v5: F32, +q0: U32, +lim: U32, +idx: U32, +st: U32, +o: U32, t: Array<T8>, r: Array<F32> & F32) -> PK:
743	  match r:
744	    case Tuple{a, +v}:
745	      pk_g7(v0, v1, v2, v3, v4, v5, pick0(U32.is_lt((q0 + 6 : U32), lim), v), q0, lim, idx, st, o, t, Array.get(F32, a, (idx + (st * 7 : U32) : U32)))
746	
747	def pk_g5(+v0: F32, +v1: F32, +v2: F32, +v3: F32, +v4: F32, +q0: U32, +lim: U32, +idx: U32, +st: U32, +o: U32, t: Array<T8>, r: Array<F32> & F32) -> PK:
748	  match r:
749	    case Tuple{a, +v}:
750	      pk_g6(v0, v1, v2, v3, v4, pick0(U32.is_lt((q0 + 5 : U32), lim), v), q0, lim, idx, st, o, t, Array.get(F32, a, (idx + (st * 6 : U32) : U32)))
751	
752	def pk_g4(+v0: F32, +v1: F32, +v2: F32, +v3: F32, +q0: U32, +lim: U32, +idx: U32, +st: U32, +o: U32, t: Array<T8>, r: Array<F32> & F32) -> PK:
753	  match r:
754	    case Tuple{a, +v}:
755	      pk_g5(v0, v1, v2, v3, pick0(U32.is_lt((q0 + 4 : U32), lim), v), q0, lim, idx, st, o, t, Array.get(F32, a, (idx + (st * 5 : U32) : U32)))
756	
757	def pk_g3(+v0: F32, +v1: F32, +v2: F32, +q0: U32, +lim: U32, +idx: U32, +st: U32, +o: U32, t: Array<T8>, r: Array<F32> & F32) -> PK:
758	  match r:
759	    case Tuple{a, +v}:
760	      pk_g4(v0, v1, v2, pick0(U32.is_lt((q0 + 3 : U32), lim), v), q0, lim, idx, st, o, t, Array.get(F32, a, (idx + (st * 4 : U32) : U32)))
761	
762	def pk_g2(+v0: F32, +v1: F32, +q0: U32, +lim: U32, +idx: U32, +st: U32, +o: U32, t: Array<T8>, r: Array<F32> & F32) -> PK:
763	  match r:
764	    case Tuple{a, +v}:
765	      pk_g3(v0, v1, pick0(U32.is_lt((q0 + 2 : U32), lim), v), q0, lim, idx, st, o, t, Array.get(F32, a, (idx + (st * 3 : U32) : U32)))
766	
767	def pk_g1(+v0: F32, +q0: U32, +lim: U32, +idx: U32, +st: U32, +o: U32, t: Array<T8>, r: Array<F32> & F32) -> PK:
768	  match r:
769	    case Tuple{a, +v}:
770	      pk_g2(v0, pick0(U32.is_lt((q0 + 1 : U32), lim), v), q0, lim, idx, st, o, t, Array.get(F32, a, (idx + (st * 2 : U32) : U32)))
771	
772	def pk_g0(+q0: U32, +lim: U32, +idx: U32, +st: U32, +o: U32, t: Array<T8>, r: Array<F32> & F32) -> PK:
773	  match r:
774	    case Tuple{a, +v}:
775	      pk_g1(pick0(U32.is_lt((q0 + 0 : U32), lim), v), q0, lim, idx, st, o, t, Array.get(F32, a, (idx + (st * 1 : U32) : U32)))
776	
777	# the cells of one line, cell o.. of t; q0 counts the line's numbers before the cell
778	def pk_line(tl: Nat, +q0: U32, +lim: U32, +base: U32, +st: U32, +o: U32, r: PK) -> PK:
779	  match tl r:
780	    case 0n r:
781	      r
782	    case 1n+q PK{a, t}:
783	      pk_line(q, (q0 + 8 : U32), lim, (base + (st * 8 : U32) : U32), st, (o + 1 : U32), pk_g0(q0, lim, base, st, o, t, Array.get(F32, a, base)))
784	
785	# lines count lines; source line i starts at i * lstep (in a), and goes to cell row o (cells o * tku .. in t)
786	def pk_lines(cnt: Nat, +tk: Nat, +tku: U32, +ku: U32, +i: U32, +o: U32, +lstep: U32, +st: U32, r: PK) -> PK:
787	  match cnt:
788	    case 0n:
789	      r
790	    case 1n+p:
791	      pk_lines(p, tk, tku, ku, (i + 1 : U32), (o + 1 : U32), lstep, st, pk_line(tk, 0, ku, (i * lstep : U32), st, (o * tku : U32), r))
792	
793	def pk_pack(+cnt: Nat, +tk: Nat, +k: Nat, +lstep: U32, +st: U32, a: Array<F32>) -> PK:
794	  pk_lines(cnt, tk, U32.from_nat(tk), U32.from_nat(k), 0, 0, lstep, st, PK{a, Array.new(T8, cap_depth(Nat.mul(cnt, tk)), t8_zero())})
795	
796	# the dot of two lines of cells: tk cells from ia (in a) and ib (in b); ra, rb hold the arrays and the cells at ia, ib,
797	# already read. Product q of a cell goes to sum q mod 4, the first four before the last four, as dotr adds them.
798	def t8_dot(tk: Nat, +ia: U32, +ib: U32, s0: F32, s1: F32, s2: F32, s3: F32, ra: Array<T8> & T8, rb: Array<T8> & T8) -> TD:
799	  match tk ra rb:
800	    case 0n Tuple{a, x} Tuple{b, y}:
801	      TD{a, b, s0, s1, s2, s3}
802	    case 1n+p Tuple{a, T8{x0, x1, x2, x3, x4, x5, x6, x7}} Tuple{b, T8{y0, y1, y2, y3, y4, y5, y6, y7}}:
803	      t8_dot(p, (ia + 1 : U32), (ib + 1 : U32), ((s0 + (x0 * y0 : F32) : F32) + (x4 * y4 : F32) : F32), ((s1 + (x1 * y1 : F32) : F32) + (x5 * y5 : F32) : F32), ((s2 + (x2 * y2 : F32) : F32) + (x6 * y6 : F32) : F32), ((s3 + (x3 * y3 : F32) : F32) + (x7 * y7 : F32) : F32), Array.get(T8, a, (ia + 1 : U32)), Array.get(T8, b, (ib + 1 : U32)))
804	
805	type TD is Type:
806	  TD{a: Array<T8>, b: Array<T8>, s0: F32, s1: F32, s2: F32, s3: F32}
807	
808	# dotr's final sum after k products: the sums sit rotated by k mod 4 when its loop ends
809	def t8_comb(+rot: Nat, +s0: F32, +s1: F32, +s2: F32, +s3: F32) -> F32:
810	  match rot:
811	    case 0n:
812	      ((s0 + s1 : F32) + (s2 + s3 : F32) : F32)
813	    case 1n:
814	      ((s1 + s2 : F32) + (s3 + s0 : F32) : F32)
815	    case 2n:
816	      ((s2 + s3 : F32) + (s0 + s1 : F32) : F32)
817	    case _:
818	      ((s3 + s0 : F32) + (s1 + s2 : F32) : F32)
819	
820	type TR is Type:
821	  TR{a: Array<T8>, b: Array<T8>, x: F32}
822	
823	def t8_fin(+rot: Nat, d: TD) -> TR:
824	  match d:
825	    case TD{a, b, +s0, +s1, +s2, +s3}:
826	      TR{a, b, t8_comb(rot, s0, s1, s2, s3)}
827	
828	def t8_cell(+tk: Nat, +tku: U32, +rot: Nat, +i: U32, +j: U32, a: Array<T8>, b: Array<T8>) -> TR:
829	  t8_fin(rot, t8_dot(tk, (i * tku : U32), (j * tku : U32), 0.0, 0.0, 0.0, 0.0, Array.get(T8, a, (i * tku : U32)), Array.get(T8, b, (j * tku : U32))))
830	
831	type TG is Type:
832	  TG{a: Array<T8>, b: Array<T8>, c: Array<F32>}
833	
834	def t8_set(+ac: Bool, +idx: U32, c: Array<F32>, r: TR) -> TG:
835	  match r:
836	    case TR{a, b, +x}:
837	      TG{a, b, cwr(ac, idx, x, c)}
838	
839	def t8_cols(+ac: Bool, mleft: Nat, +tk: Nat, +tku: U32, +rot: Nat, +i: U32, +j: U32, +cidx: U32, g: TG) -> TG:
840	  match mleft g:
841	    case 0n g:
842	      g
843	    case 1n+p TG{a, b, c}:
844	      t8_cols(ac, p, tk, tku, rot, i, (j + 1 : U32), (cidx + 1 : U32), t8_set(ac, cidx, c, t8_cell(tk, tku, rot, i, j, a, b)))
845	
846	def t8_rows(+ac: Bool, nleft: Nat, +m: Nat, +tk: Nat, +tku: U32, +rot: Nat, +i: U32, +mu: U32, g: TG) -> TG:
847	  match nleft:
848	    case 0n:
849	      g
850	    case 1n+p:
851	      t8_rows(ac, p, m, tk, tku, rot, (i + 1 : U32), mu, t8_cols(ac, m, tk, tku, rot, i, 0, (i * mu : U32), g))
852	
853	# rows i0 .. i0+cnt-1 of C as a list, in order (the parallel version's leaves)
854	type TL is Type:
855	  TL{a: Array<T8>, b: Array<T8>, out: List<&2, F32>}
856	
857	# cells j.. of row i, m - j of them; r is cell j, already computed
858	def t8_lcols(mleft: Nat, +tk: Nat, +tku: U32, +rot: Nat, +i: U32, +j: U32, out: List<&2, F32>, r: TR) -> TL:
859	  match mleft r:
860	    case 0n TR{a, b, +x}:
861	      TL{a, b, x <> out}
862	    case 1n+p TR{a, b, +x}:
863	      t8_lcols(p, tk, tku, rot, i, (j + 1 : U32), x <> out, t8_cell(tk, tku, rot, i, (j + 1 : U32), a, b))
864	
865	def t8_lrows(nleft: Nat, +m: Nat, +tk: Nat, +tku: U32, +rot: Nat, +i: U32, l: TL) -> TL:
866	  match nleft l:
867	    case 0n l:
868	      l
869	    case 1n+p TL{a, b, out}:
870	      t8_lrows(p, m, tk, tku, rot, (i + 1 : U32), t8_lcols(Nat.sub(m, 1n), tk, tku, rot, i, 0, out, t8_cell(tk, tku, rot, i, 0, a, b)))
871	
872	def t8_lfin(l: TL) -> List<&2, F32>:
873	  match l:
874	    case TL{a, b, out}:
875	      List.reverse(&2, F32, out)
876	
877	type PKK is Type:
878	  PKK{a: Array<F32>, b: Array<F32>, ta: Array<T8>, tb: Array<T8>}
879	
880	def t8_pk2(pa: PK, pb: PK) -> PKK:
881	  match pa pb:
882	    case PK{a, ta} PK{b, tb}:
883	      PKK{a, b, ta, tb}
884	
885	def t8_g(a: Array<F32>, b: Array<F32>, g: TG) -> G:
886	  match g:
887	    case TG{ta, tb, c}:
888	      G{a, b, c}
889	
890	def t8_out_seq(+ac: Bool, +n: Nat, +m: Nat, +tk: Nat, +tku: U32, +rot: Nat, c: Array<F32>, p: PKK) -> G:
891	  match p:
892	    case PKK{a, b, ta, tb}:
893	      t8_g(a, b, t8_rows(ac, n, m, tk, tku, rot, 0, U32.from_nat(m), TG{ta, tb, c}))
894	
895	# ---- the parallel path (0.1.9.0): the larger operand is packed into one piece per block, so a block owns its piece;
896	# the other operand is copied once per node of the fork tree (2^d - 1 copies), instead of both at every node ----
897	
898	type PT is Type:
899	  PLeaf{t: Array<T8>}
900	  PNode{l: PT, r: PT}
901	
902	type PKT is Type:
903	  PKT{a: Array<F32>, t: PT}
904	
905	def pkt_leaf(r: PK) -> PKT:
906	  match r:
907	    case PK{a, t}:
908	      PKT{a, PLeaf{t}}
909	
910	def pkt_then(r: PKT, k: Array<F32> -> PT -> PKT) -> PKT:
911	  match r:
912	    case PKT{a, t}:
913	      k(a, t)
914	
915	# lines i .. i+cnt-1 of a, packed into 2^d pieces (halves of the line count, as the fork tree splits them)
916	def pk_tree(+d: Nat, +cnt: Nat, +i: U32, +tk: Nat, +tku: U32, +ku: U32, +lstep: U32, +st: U32, a: Array<F32>) -> PKT:
917	  match d:
918	    case 0n:
919	      pkt_leaf(pk_lines(cnt, tk, tku, ku, i, 0, lstep, st, PK{a, Array.new(T8, cap_depth(Nat.mul(cnt, tk)), t8_zero())}))
920	    case 1n+q:
921	      pkt_then(pk_tree(q, Nat.div(cnt, 2n), i, tk, tku, ku, lstep, st, a), a2 => l => pkt_then(pk_tree(q, Nat.sub(cnt, Nat.div(cnt, 2n)), (i + U32.from_nat(Nat.div(cnt, 2n)) : U32), tk, tku, ku, lstep, st, a2), a3 => r => PKT{a3, PNode{l, r}}))
922	
923	def t8_with2(-R: Type, p: Array<T8> & Array<T8>, k: Array<T8> -> Array<T8> -> R) -> R:
924	  match p:
925	    case Tuple{x, y}:
926	      k(x, y)
927	
928	# split by rows of C: the pieces are A's rows, B is shared; the result is C's rows in order
929	def t8_rpar(t: PT, +cnt: Nat, +m: Nat, +tk: Nat, +tku: U32, +rot: Nat, b: Array<T8>) -> List<&2, F32>:
930	  match t:
931	    case PLeaf{a}:
932	      t8_lfin(t8_lrows(cnt, m, tk, tku, rot, 0, TL{a, b, Nil{}}))
933	    case PNode{l, r}:
934	      t8_with2(List<&2, F32>, Array.clone(T8, b), b1 => b2 => (
935	        x y = t8_rpar(l, Nat.div(cnt, 2n), m, tk, tku, rot, b1) t8_rpar(r, Nat.sub(cnt, Nat.div(cnt, 2n)), m, tk, tku, rot, b2)
936	        List.append(&2, F32, x, y)))
937	
938	# a block of C's columns j0 .. j0+w-1 for every row, row by row
939	type CB is Data:
940	  CB{j0: U32, w: Nat, xs: List<&2, F32>}
941	
942	# split by columns of C: the pieces are B's columns, A is shared
943	def t8_cpar(t: PT, +cnt: Nat, +j0: U32, +n: Nat, +tk: Nat, +tku: U32, +rot: Nat, a: Array<T8>) -> List<&2, CB>:
944	  match t:
945	    case PLeaf{b}:
946	      [CB{j0, cnt, t8_lfin(t8_lrows(n, cnt, tk, tku, rot, 0, TL{a, b, Nil{}}))}]
947	    case PNode{l, r}:
948	      t8_with2(List<&2, CB>, Array.clone(T8, a), a1 => a2 => (
949	        x y = t8_cpar(l, Nat.div(cnt, 2n), j0, n, tk, tku, rot, a1) t8_cpar(r, Nat.sub(cnt, Nat.div(cnt, 2n)), (j0 + U32.from_nat(Nat.div(cnt, 2n)) : U32), n, tk, tku, rot, a2)
950	        List.append(&2, CB, x, y)))
951	
952	type RW is Type:
953	  RW{c: Array<F32>, xs: List<&2, F32>}
954	
955	# w numbers of xs into c from idx; what is left of xs comes back
956	def cb_row(+ac: Bool, w: Nat, +idx: U32, c: Array<F32>, xs: List<&2, F32>) -> RW:
957	  match w xs:
958	    case 1n+p Con{+h, t}:
959	      cb_row(ac, p, (idx + 1 : U32), cwr(ac, idx, h, c), t)
960	    case _ _:
961	      RW{c, xs}
962	
963	def cb_rows(+ac: Bool, n: Nat, +i: U32, +j0: U32, +w: Nat, +mu: U32, r: RW) -> Array<F32>:
964	  match n r:
965	    case 0n RW{c, xs}:
966	      c
967	    case 1n+p RW{c, xs}:
968	      cb_rows(ac, p, (i + 1 : U32), j0, w, mu, cb_row(ac, w, ((i * mu : U32) + j0 : U32), c, xs))
969	
970	def cb_write(+ac: Bool, bs: List<&2, CB>, +n: Nat, +mu: U32, c: Array<F32>) -> Array<F32>:
971	  match bs:
972	    case Nil{}:
973	      c
974	    case Con{CB{+j0, +w, xs}, t}:
975	      cb_write(ac, t, n, mu, cb_rows(ac, n, 0, j0, w, mu, RW{c, xs}))
976	
977	def t8_pr2(+ac: Bool, +n: Nat, +m: Nat, +tk: Nat, +tku: U32, +rot: Nat, c: Array<F32>, ra: PKT, rb: PK) -> G:
978	  match ra rb:
979	    case PKT{a, t} PK{b, tb}:
980	      G{a, b, write_la(ac, t8_rpar(t, n, m, tk, tku, rot, tb), 0, c)}
981	
982	def t8_pc2(+ac: Bool, +n: Nat, +m: Nat, +tk: Nat, +tku: U32, +rot: Nat, c: Array<F32>, ra: PK, rb: PKT) -> G:
983	  match ra rb:
984	    case PK{a, ta} PKT{b, t}:
985	      G{a, b, cb_write(ac, t8_cpar(t, m, 0, n, tk, tku, rot, ta), n, U32.from_nat(m), c)}
986	
987	# 2^d blocks along the larger of n (rows of C, pieces of A) and m (columns of C, pieces of B)
988	def t8_par2(+ac: Bool, byrows: Bool, d: Nat, +n: Nat, +m: Nat, +k: Nat, +tk: Nat, +tku: U32, +rot: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
989	  match byrows:
990	    case True{}:
991	      t8_pr2(ac, n, m, tk, tku, rot, c, pk_tree(d, n, 0, tk, tku, U32.from_nat(k), arow, sa, a), pk_pack(m, tk, k, bcol, sb, b))
992	    case False{}:
993	      t8_pc2(ac, n, m, tk, tku, rot, c, pk_pack(n, tk, k, arow, sa, a), pk_tree(d, m, 0, tk, tku, U32.from_nat(k), bcol, sb, b))
994	
995	def gemm_t8b(+ac: Bool, d: Nat, +n: Nat, +m: Nat, +k: Nat, +tk: Nat, +tku: U32, +rot: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
996	  match d:
997	    case 0n:
998	      t8_out_seq(ac, n, m, tk, tku, rot, c, t8_pk2(pk_pack(n, tk, k, arow, sa, a), pk_pack(m, tk, k, bcol, sb, b)))
999	    case 1n+q:
1000	      t8_par2(ac, Nat.is_le(m, n), 1n+q, n, m, k, tk, tku, rot, arow, sa, sb, bcol, a, b, c)
1001	
1002	# C (n x m) = A · B (or C + A · B when ac) in cells of 8, sequential (d = 0) or with 2^d blocks of rows in parallel; A
1003	# and B come back intact
1004	def gemm_t8_ac(+ac: Bool, d: Nat, +n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
1005	  gemm_t8b(ac, d, n, m, k, Nat.div(Nat.add(k, 7n), 8n), U32.from_nat(Nat.div(Nat.add(k, 7n), 8n)), Nat.mod(k, 4n), arow, sa, sb, bcol, a, b, c)
1006	
1007	def gemm_t8(d: Nat, +n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
1008	  gemm_t8_ac(False{}, d, n, m, k, arow, sa, sb, bcol, a, b, c)
1009	
1010	# one row (matrix · vector): the scalar loop, by columns in parallel when d > 0; two or more rows: cells of 8
1011	def gemm_any3(+ac: Bool, one: Bool, d: Nat, +n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
1012	  match one d:
1013	    case True{} 0n:
1014	      gemm_ac(ac, n, m, k, arow, sa, sb, bcol, a, b, c)
1015	    case True{} 1n+q:
1016	      gemm_par_cols(ac, 1n+q, m, k, sa, sb, bcol, a, b, c)
1017	    case False{} d:
1018	      gemm_t8_ac(ac, d, n, m, k, arow, sa, sb, bcol, a, b, c)
1019	
1020	def gemm_any2(+ac: Bool, d: Nat, +n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
1021	  gemm_any3(ac, Nat.is_le(n, 1n), d, n, m, k, arow, sa, sb, bcol, a, b, c)
1022	
1023	def gemm_any(+ac: Bool, d: Nat, +n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array<F32>, b: Array<F32>, c: Array<F32>) -> G:
1024	  gemm_any2(ac, d, n, m, k, arow, sa, sb, bcol, a, b, c)
1025	
1026	def mmul_fin(-n: Nat, -k: Nat, -m: Nat, g: G) -> MMul<n, k, m>:
1027	  match g:
1028	    case G{a, b, c}:
1029	      MMul{Mat{a}, Mat{b}, Mat{c}}
1030	
1031	def mmul_nt_fin(-n: Nat, -k: Nat, -m: Nat, g: G) -> MMulNT<n, k, m>:
1032	  match g:
1033	    case G{a, b, c}:
1034	      MMulNT{Mat{a}, Mat{b}, Mat{c}}
1035	
1036	def mmul_tn_fin(-n: Nat, -k: Nat, -m: Nat, g: G) -> MMulTN<n, k, m>:
1037	  match g:
1038	    case G{a, b, c}:
1039	      MMulTN{Mat{a}, Mat{b}, Mat{c}}
1040	
1041	# ---- the split depth of a product (0.1.11.0): par is an upper bound, and a product too small to gain splits 0 times ----
1042	#
1043	# A product of n x k by k x m splits into 2^par blocks only when it does at least 2^23 multiply-adds and each cell of C
1044	# sums at least 64 of them; otherwise it runs sequentially. Measured on 16 threads (bend2-lab, B8 sweep): below 2^23 the
1045	# split costs more than it saves, and with a short inner dimension (a weight gradient over a few tokens) every split
1046	# was 2-3.6x slower, because each block hands its rows of C back as a list that is then written into C, a cost per cell
1047	# of C that k multiply-adds per cell do not pay for. The result does not depend on the depth: the same bits at every
1048	# depth (tests/t8.bend, tests/acc.bend, tests/split.bend).
1049	def split_pick(big: Bool, d: Nat) -> Nat:
1050	  match big:
1051	    case True{}:
1052	      d
1053	    case False{}:
1054	      0n
1055	
1056	def split_depth(d: Nat, +n: Nat, +m: Nat, +k: Nat) -> Nat:
1057	  split_pick(Bool.and(Nat.is_le(64n, k), Nat.is_le(8388608n, Nat.mul(Nat.mul(n, m), k))), d)
1058	
1059	def mm_a(+n: Nat, +k: Nat, +m: Nat, d: Nat, a: Array<F32>, b: Array<F32>) -> MMul<n, k, m>:
1060	  mmul_fin(n, k, m, gemm_any(False{}, d, n, m, k, U32.from_nat(k), 1, U32.from_nat(m), 1, a, b, Array.new(F32, cap_depth(Nat.mul(n, m)), 0.0)))
1061	
1062	# C = A · B   (n x k) · (k x m); par = log2 of the most parallel blocks (0 = sequential); split_depth picks 0 or par
1063	
1064	def Mat.matmul(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: Mat<n, k>, b: Mat<k, m>) -> MMul<n, k, m>:
1065	  match a b:
1066	    case Mat{x} Mat{y}:
1067	      mm_a(n, k, m, split_depth(par, n, m, k), x, y)
1068	
1069	# C + A · B, into C (n x m); A and B come back with the result, as Mat.matmul returns them. Each cell of A · B is
1070	# computed as Mat.matmul computes it and added once to C's cell: the same bits as Mat.add(c, Mat.matmul(a, b)'s c),
1071	# without allocating the product or walking it a second time.
1072	
1073	def Mat.matmul_acc(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: Mat<n, k>, b: Mat<k, m>, c: Mat<n, m>) -> MMul<n, k, m>:
1074	  match a b c:
1075	    case Mat{x} Mat{y} Mat{z}:
1076	      mmul_fin(n, k, m, gemm_any(True{}, split_depth(par, n, m, k), n, m, k, U32.from_nat(k), 1, U32.from_nat(m), 1, x, y, z))
1077	
1078	def mm_nt(+n: Nat, +k: Nat, +m: Nat, d: Nat, a: Array<F32>, b: Array<F32>) -> MMulNT<n, k, m>:
1079	  mmul_nt_fin(n, k, m, gemm_any(False{}, d, n, m, k, U32.from_nat(k), 1, 1, U32.from_nat(k), a, b, Array.new(F32, cap_depth(Nat.mul(n, m)), 0.0)))
1080	
1081	# C = A · Bᵀ   (n x k) · (m x k)ᵀ: for weights stored with one row per output
1082	
1083	def Mat.matmul_nt(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: Mat<n, k>, b: Mat<m, k>) -> MMulNT<n, k, m>:
1084	  match a b:
1085	    case Mat{x} Mat{y}:
1086	      mm_nt(n, k, m, split_depth(par, n, m, k), x, y)
1087	
1088	# C + A · Bᵀ, into C (n x m)
1089	
1090	def Mat.matmul_nt_acc(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: Mat<n, k>, b: Mat<m, k>, c: Mat<n, m>) -> MMulNT<n, k, m>:
1091	  match a b c:
1092	    case Mat{x} Mat{y} Mat{z}:
1093	      mmul_nt_fin(n, k, m, gemm_any(True{}, split_depth(par, n, m, k), n, m, k, U32.from_nat(k), 1, 1, U32.from_nat(k), x, y, z))
1094	
1095	def mm_tn(+n: Nat, +k: Nat, +m: Nat, d: Nat, a: Array<F32>, b: Array<F32>) -> MMulTN<n, k, m>:
1096	  mmul_tn_fin(n, k, m, gemm_any(False{}, d, n, m, k, 1, U32.from_nat(n), U32.from_nat(m), 1, a, b, Array.new(F32, cap_depth(Nat.mul(n, m)), 0.0)))
1097	
1098	# C = Aᵀ · B   (k x n)ᵀ · (k x m): weight gradient, X^T · dY
1099	
1100	def Mat.matmul_tn(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: Mat<k, n>, b: Mat<k, m>) -> MMulTN<n, k, m>:
1101	  match a b:
1102	    case Mat{x} Mat{y}:
1103	      mm_tn(n, k, m, split_depth(par, n, m, k), x, y)
1104	
1105	# C + Aᵀ · B, into C (n x m): a weight gradient added into the gradient it accumulates, dW += Xᵀ · dY
1106	
1107	def Mat.matmul_tn_acc(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: Mat<k, n>, b: Mat<k, m>, c: Mat<n, m>) -> MMulTN<n, k, m>:
1108	  match a b c:
1109	    case Mat{x} Mat{y} Mat{z}:
1110	      mmul_tn_fin(n, k, m, gemm_any(True{}, split_depth(par, n, m, k), n, m, k, 1, U32.from_nat(n), U32.from_nat(m), 1, x, y, z))
1111	
1112	# ---- element-wise and training operations ----
1113	
1114	# adds the bias (one row, m numbers) to each row of C (n x m)
1115	
1116	type MBias<-n: Nat, -m: Nat> is Type:
1117	  MBias{c: Mat<n, m>, b: Mat<1n, m>}
1118	
1119	def add_row2(-n: Nat, -m: Nat, +nn: Nat, +mm: Nat, cd: Array<F32>, al: AL) -> MBias<n, m>:
1120	  match al:
1121	    case AL{bd, bl}:
1122	      MBias{Mat{bias_rows(nn, 0, U32.from_nat(mm), bl, cd)}, Mat{bd}}
1123	
1124	def Mat.add_row(+n: Nat, +m: Nat, c: Mat<n, m>, b: Mat<1n, m>) -> MBias<n, m>:
1125	  match c b:
1126	    case Mat{cd} Mat{bd}:
1127	      add_row2(n, m, n, m, cd, read_l(Nat.sub(m, 1n), 0, Nil{}, Array.get(F32, bd, 0)))
1128	
1129	# pair of matrices of the same shape (the result of an operation that reads two operands)
1130	
1131	type MPair<-r: Nat, -c: Nat> is Type:
1132	  MPair{a: Mat<r, c>, b: Mat<r, c>}
1133	
1134	def mpair_fin(-r: Nat, -c: Nat, p: P2) -> MPair<r, c>:
1135	  match p:
1136	    case P2{a, b}:
1137	      MPair{Mat{a}, Mat{b}}
1138	
1139	def relu_fin(-r: Nat, -c: Nat, d: Array<F32>) -> Mat<r, c>:
1140	  Mat{d}
1141	
1142	def Mat.relu(+r: Nat, +c: Nat, a: Mat<r, c>) -> Mat<r, c>:
1143	  match a:
1144	    case Mat{d}:
1145	      relu_fin(r, c, relu_l(Nat.mul(r, c), 0, Array.get(F32, d, 0)))
1146	
1147	# (dy * (h > 0), h): the gradient that passes through the relu
1148	
1149	def Mat.relu_bwd(+r: Nat, +c: Nat, dy: Mat<r, c>, h: Mat<r, c>) -> MPair<r, c>:
1150	  match dy h:
1151	    case Mat{d} Mat{x}:
1152	      mpair_fin(r, c, mask_l(Nat.mul(r, c), 0, Array.get(F32, d, 0), Array.get(F32, x, 0)))
1153	
1154	# (w - lr * dw, dw): one SGD step
1155	
1156	def Mat.sgd(+r: Nat, +c: Nat, +lr: F32, w: Mat<r, c>, dw: Mat<r, c>) -> MPair<r, c>:
1157	  match w dw:
1158	    case Mat{x} Mat{y}:
1159	      mpair_fin(r, c, sgd_l(Nat.mul(r, c), 0, lr, Array.get(F32, x, 0), Array.get(F32, y, 0)))
1160	
1161	# (a + b, b): element-wise sum (a - (-1)*b is exactly a + b)
1162	
1163	def Mat.add(+r: Nat, +c: Nat, a: Mat<r, c>, b: Mat<r, c>) -> MPair<r, c>:
1164	  match a b:
1165	    case Mat{x} Mat{y}:
1166	      mpair_fin(r, c, sgd_l(Nat.mul(r, c), 0, (0.0 - 1.0 : F32), Array.get(F32, x, 0), Array.get(F32, y, 0)))
1167	
1168	# sum of each column (n x m -> 1 x m), returning the matrix
1169	
1170	type MColSum<-n: Nat, -m: Nat> is Type:
1171	  MColSum{a: Mat<n, m>, s: Mat<1n, m>}
1172	
1173	def colsum_fin(-n: Nat, -m: Nat, g: G) -> MColSum<n, m>:
1174	  match g:
1175	    case G{ones, a, s}:
1176	      MColSum{Mat{a}, Mat{s}}
1177	
1178	def Mat.col_sums(+n: Nat, +m: Nat, a: Mat<n, m>) -> MColSum<n, m>:
1179	  match a:
1180	    case Mat{d}:
1181	      colsum_fin(n, m, gemm(1n, m, n, 0, 1, U32.from_nat(m), 1, Array.new(F32, cap_depth(n), 1.0), d, Array.new(F32, cap_depth(m), 0.0)))
1182	
1183	# s + the sum of each column, into s (1 x m): a bias gradient added into its accumulator, with the bits of
1184	# Mat.add(s, Mat.col_sums(a)'s s)
1185	
1186	def Mat.col_sums_acc(+n: Nat, +m: Nat, a: Mat<n, m>, s: Mat<1n, m>) -> MColSum<n, m>:
1187	  match a s:
1188	    case Mat{d} Mat{z}:
1189	      colsum_fin(n, m, gemm_ac(True{}, 1n, m, n, 0, 1, U32.from_nat(m), 1, Array.new(F32, cap_depth(n), 1.0), d, z))
1190	
1191	# mean cross-entropy and its gradient on the logits: (softmax - one-hot) / n.
1192	# labels has one entry per row (n); the gradient replaces the logits.
1193	
1194	type MCE<-n: Nat, -c: Nat> is Type:
1195	  MCE{loss: F32, grad: Mat<n, c>}
1196	
1197	def ce_fin(-n: Nat, -c: Nat, +nf: F32, r: L2) -> MCE<n, c>:
1198	  match r:
1199	    case L2{loss, d}:
1200	      MCE{(loss / nf : F32), Mat{d}}
1201	
1202	def Mat.softmax_ce(+n: Nat, +c: Nat, z: Mat<n, c>, labels: List<&2, Nat>) -> MCE<n, c>:
1203	  match z:
1204	    case Mat{d}:
1205	      ce_fin(n, c, F32.from_nat(n), ce_rows(labels, Nat.sub(c, 1n), U32.from_nat(c), F32.from_nat(n), 0, L2{0.0, d}))
1206	
1207	# argmax hits per row against labels
1208	
1209	type MHits<-n: Nat, -c: Nat> is Type:
1210	  MHits{z: Mat<n, c>, hits: Nat}
1211	
1212	def hits_fin(-n: Nat, -c: Nat, r: CR) -> MHits<n, c>:
1213	  match r:
1214	    case CR{d, k}:
1215	      MHits{Mat{d}, k}
1216	
1217	def Mat.count_correct(+n: Nat, +c: Nat, z: Mat<n, c>, labels: List<&2, Nat>) -> MHits<n, c>:
1218	  match z:
1219	    case Mat{d}:
1220	      hits_fin(n, c, count_rows(labels, Nat.sub(c, 1n), U32.from_nat(c), 0, CR{d, 0n}))
1221	
1222	# a row as a list (returning the matrix) and the other way around
1223	
1224	type MRow<-r: Nat, -c: Nat> is Type:
1225	  MRow{m: Mat<r, c>, row: List<&2, F32>}
1226	
1227	def row_fin(-r: Nat, -c: Nat, al: AL) -> MRow<r, c>:
1228	  match al:
1229	    case AL{d, l}:
1230	      MRow{Mat{d}, l}
1231	
1232	def Mat.read_row(+r: Nat, +c: Nat, m: Mat<r, c>, +i: U32) -> MRow<r, c>:
1233	  match m:
1234	    case Mat{d}:
1235	      row_fin(r, c, read_l(Nat.sub(c, 1n), (i * U32.from_nat(c) : U32), Nil{}, Array.get(F32, d, (i * U32.from_nat(c) : U32))))
1236	
1237	def Mat.write_row(-r: Nat, +c: Nat, m: Mat<r, c>, +i: U32, xs: List<&2, F32>) -> Mat<r, c>:
1238	  match m:
1239	    case Mat{d}:
1240	      Mat{write_l(xs, (i * U32.from_nat(c) : U32), d)}
1241	
1242	# ---- variants that check the number of labels ----
1243	
1244	# does the list of labels have exactly n entries (one per row)?
1245	def labels_len_is(xs: List<&2, Nat>, n: Nat) -> Bool:
1246	  match xs n:
1247	    case Nil{} 0n:
1248	      True{}
1249	    case Con{h, t} 1n+p:
1250	      labels_len_is(t, p)
1251	    case _ _:
1252	      False{}
1253	
1254	def ce_checked(ok: Bool, +n: Nat, +c: Nat, z: Mat<n, c>, labels: List<&2, Nat>) -> Maybe<&1, MCE<n, c>>:
1255	  match ok:
1256	    case True{}:
1257	      Some{Mat.softmax_ce(n, c, z, labels)}
1258	    case False{}:
1259	      None{}
1260	
1261	# like Mat.softmax_ce, but returns None unless labels has exactly n entries
1262	def Mat.softmax_ce_checked(+n: Nat, +c: Nat, z: Mat<n, c>, +labels: List<&2, Nat>) -> Maybe<&1, MCE<n, c>>:
1263	  ce_checked(labels_len_is(labels, n), n, c, z, labels)
1264	
1265	def hits_checked(ok: Bool, +n: Nat, +c: Nat, z: Mat<n, c>, labels: List<&2, Nat>) -> Maybe<&1, MHits<n, c>>:
1266	  match ok:
1267	    case True{}:
1268	      Some{Mat.count_correct(n, c, z, labels)}
1269	    case False{}:
1270	      None{}
1271	
1272	# like Mat.count_correct, but returns None unless labels has exactly n entries
1273	def Mat.count_correct_checked(+n: Nat, +c: Nat, z: Mat<n, c>, +labels: List<&2, Nat>) -> Maybe<&1, MHits<n, c>>:
1274	  hits_checked(labels_len_is(labels, n), n, c, z, labels)
1275	
1276	# ---------------------------------------------------------------------
1277	# Bands<r, c>: a matrix cut into bands of rows, for matrix · vector in parallel
1278	# without copying the weights.
1279	#
1280	# Mat.matmul_nt with n = 1 can only run in parallel by cloning the whole weight
1281	# matrix for every task, and in matrix · vector each weight is read once, so the
1282	# copies cost more than the products (NOTES.md, exp. 9). Bands keeps each band in
1283	# its own Array: a task takes its band and nothing is copied except the input vector.
1284	#
1285	# A node with r rows has children with half(r) and r - half(r) rows. The split is
1286	# part of the type, so the row count of every band follows from r, and the law
1287	# half_cover says that the two halves always cover the r rows exactly.
1288	# ---------------------------------------------------------------------
1289	
1290	# n / 2 rounded down. A tail loop with an accumulator: a non-tail recursion here would put a
1291	# continuation on every node of the band tree and keep the bands from running in parallel.
1292	def half.go(n: Nat, acc: Nat) -> Nat:
1293	  match n:
1294	    case 0n:
1295	      acc
1296	    case 1n+p:
1297	      match p:
1298	        case 0n:
1299	          acc
1300	        case 1n+q:
1301	          half.go(q, 1n+acc)
1302	
1303	def half(n: Nat) -> Nat:
1304	  half.go(n, 0n)
1305	
1306	# a <= a
1307	def leb_refl(+a: Nat) -> {True{} == leb(a, a) : Bool}:
1308	  match a:
1309	    case 0n:
1310	      {==}
1311	    case 1n+p:
1312	      leb_refl(p)
1313	
1314	# 1 + (a + b) == a + (1 + b)
1315	def add_succ_r(+a: Nat, -b: Nat) -> {1n+Nat.add(a, b) == Nat.add(a, 1n+b) : Nat}:
1316	  match a:
1317	    case 0n:
1318	      {==}
1319	    case 1n+p:
1320	      %add_succ_r(p, b) : {2n+Nat.add(p, b) == 1n+_ : Nat}
1321	      {==}
1322	
1323	# a + 0 == a
1324	def add_zero_r(+a: Nat) -> {Nat.add(a, 0n) == a : Nat}:
1325	  match a:
1326	    case 0n:
1327	      {==}
1328	    case 1n+p:
1329	      %add_zero_r(p) : {1n+Nat.add(p, 0n) == 1n+_ : Nat}
1330	      {==}
1331	
1332	# half.go(n, acc) <= n + acc
1333	def half_go_le(+n: Nat, +acc: Nat) -> {True{} == leb(half.go(n, acc), Nat.add(n, acc)) : Bool}:
1334	  match n:
1335	    case 0n:
1336	      leb_refl(acc)
1337	    case 1n+p:
1338	      match p:
1339	        case 0n:
1340	          leb_succ(acc, acc, leb_refl(acc))
1341	        case 1n+q:
1342	          %Equal.sym(Nat, 1n+Nat.add(q, acc), Nat.add(q, 1n+acc), add_succ_r(q, acc)) : {True{} == leb(half.go(q, 1n+acc), 1n+_) : Bool}
1343	          leb_succ(half.go(q, 1n+acc), Nat.add(q, 1n+acc), half_go_le(q, 1n+acc))
1344	
1345	# half(n) <= n
1346	def half_le(+n: Nat) -> {True{} == leb(half(n), n) : Bool}:
1347	  %add_zero_r(n) : {True{} == leb(half(n), _) : Bool}
1348	  half_go_le(n, 0n)
1349	
1350	# a <= r  ->  r == a + (r - a)
1351	def add_sub(+a: Nat, +r: Nat, e: {True{} == leb(a, r) : Bool}) -> {r == Nat.add(a, Nat.sub(r, a)) : Nat}:
1352	  match a r:
1353	    case 0n 0n:
1354	      {==}
1355	    case 0n 1n+q:
1356	      {==}
1357	    case 1n+p 0n:
1358	      %e : BD(_, Unit, {0n == Nat.add(1n+p, Nat.sub(0n, 1n+p)) : Nat})
1359	      Unit{}
1360	    case 1n+p 1n+q:
1361	      %add_sub(p, q, e) : {1n+q == 1n+_ : Nat}
1362	      {==}
1363	
1364	# LAW (half_cover): the two bands of a node, half(r) rows and r - half(r) rows,
1365	# add up to exactly r rows: no row is lost and none is counted twice.
1366	law half_cover:
1367	  for +n: Nat
1368	  {n == Nat.add(half(n), Nat.sub(n, half(n))) : Nat}
1369	
1370	def half_cover(n):
1371	  add_sub(half(n), n, half_le(n))
1372	
1373	type Bands<-r: Nat, -c: Nat> is Type:
1374	  BLeaf{m: Mat<r, c>}
1375	  BNode{x: Bands<half(r), c>, y: Bands<Nat.sub(r, half(r)), c>}
1376	
1377	# r x c of zeros in 2^d bands (d = 0: a single band, the same as a Mat)
1378	def Bands.zeros(+d: Nat, +r: Nat, +c: Nat) -> Bands<r, c>:
1379	  match d:
1380	    case 0n:
1381	      BLeaf{Mat.zeros(r, c)}
1382	    case 1n+q:
1383	      BNode{Bands.zeros(q, half(r), c), Bands.zeros(q, Nat.sub(r, half(r)), c)}
1384	
1385	# the row counts of the two bands follow the row count of the node: if rr == r at run
1386	# time, the same holds for both halves (the bands carry r only in their type)
1387	def eq_half(-a: Nat, -b: Nat, e: {a == b : Nat}) -> {half(a) == half(b) : Nat}:
1388	  %e : {half(a) == half(_) : Nat}
1389	  {==}
1390	
1391	def eq_rest(-a: Nat, -b: Nat, e: {a == b : Nat}) -> {Nat.sub(a, half(a)) == Nat.sub(b, half(b)) : Nat}:
1392	  %e : {Nat.sub(a, half(a)) == Nat.sub(_, half(_)) : Nat}
1393	  {==}
1394	
1395	# does a block of n numbers starting at flat index i end before the second band of a
1396	# node with r rows? does it start in the second band?
1397	def bfill_l(+r: Nat, +c: Nat, +i: Nat, +n: Nat) -> Bool:
1398	  Nat.is_le(Nat.add(i, n), Nat.mul(half(r), c))
1399	
1400	def bfill_r(+r: Nat, +c: Nat, +i: Nat) -> Bool:
1401	  Nat.is_le(Nat.mul(half(r), c), i)
1402	
1403	# writes the first n numbers of xs at index i (and nothing past them)
1404	def fill_n(n: Nat, xs: List<&2, F32>, +i: U32, a: Array<F32>) -> Array<F32>:
1405	  match n xs:
1406	    case 1n+p Con{h, t}:
1407	      fill_n(p, t, (i + 1 : U32), Array.set(F32, a, i, h))
1408	    case _ _:
1409	      a
1410	
1411	def fill_band(-r: Nat, -c: Nat, n: Nat, xs: List<&2, F32>, +i: U32, m: Mat<r, c>) -> Mat<r, c>:
1412	  match m:
1413	    case Mat{d}:
1414	      Mat{fill_n(n, xs, i, d)}
1415	
1416	# writes the first n numbers of xs at flat index i; l and rt say which bands they touch
1417	# (bfill_l, bfill_r); rr is the row count r at run time. A block that crosses into the second
1418	# band is not copied: the first band writes only its part, the second gets the rest by List.drop.
1419	def bfill(+c: Nat, -r: Nat, b: Bands<r, c>, l: Bool, rt: Bool, +rr: Nat, -e: {rr == r : Nat}, +i: Nat, +n: Nat, +xs: List<&2, F32>) -> Bands<r, c>:
1420	  match b l rt:
1421	    case BLeaf{m} _ _:
1422	      BLeaf{fill_band(r, c, n, xs, U32.from_nat(i), m)}
1423	    case BNode{x, y} True{} _:
1424	      BNode{bfill(c, half(r), x, bfill_l(half(rr), c, i, n), bfill_r(half(rr), c, i), half(rr), eq_half(rr, r, e), i, n, xs), y}
1425	    case BNode{x, y} False{} True{}:
1426	      BNode{x, bfill(c, Nat.sub(r, half(r)), y, bfill_l(Nat.sub(rr, half(rr)), c, Nat.sub(i, Nat.mul(half(rr), c)), n), bfill_r(Nat.sub(rr, half(rr)), c, Nat.sub(i, Nat.mul(half(rr), c))), Nat.sub(rr, half(rr)), eq_rest(rr, r, e), Nat.sub(i, Nat.mul(half(rr), c)), n, xs)}
1427	    case BNode{x, y} False{} False{}:
1428	      BNode{bfill(c, half(r), x, bfill_l(half(rr), c, i, Nat.sub(Nat.mul(half(rr), c), i)), bfill_r(half(rr), c, i), half(rr), eq_half(rr, r, e), i, Nat.sub(Nat.mul(half(rr), c), i), xs), bfill(c, Nat.sub(r, half(r)), y, bfill_l(Nat.sub(rr, half(rr)), c, 0n, Nat.sub(n, Nat.sub(Nat.mul(half(rr), c), i))), bfill_r(Nat.sub(rr, half(rr)), c, 0n), Nat.sub(rr, half(rr)), eq_rest(rr, r, e), 0n, Nat.sub(n, Nat.sub(Nat.mul(half(rr), c), i)), List.drop(&2, F32, xs, Nat.sub(Nat.mul(half(rr), c), i)))}
1429	
1430	def bfill_top(ok: Bool, +r: Nat, +c: Nat, b: Bands<r, c>, +i: Nat, +n: Nat, +xs: List<&2, F32>) -> Bands<r, c>:
1431	  match ok:
1432	    case True{}:
1433	      bfill(c, r, b, bfill_l(r, c, i, n), bfill_r(r, c, i), r, {==}, i, n, xs)
1434	    case False{}:
1435	      bfill(c, r, b, bfill_l(r, c, i, Nat.sub(Nat.mul(r, c), i)), bfill_r(r, c, i), r, {==}, i, Nat.sub(Nat.mul(r, c), i), xs)
1436	
1437	# writes the n numbers of xs into b starting at flat index i (row-major, as Mat.fill_at);
1438	# numbers past r*c are dropped. For loading a large matrix block by block.
1439	def Bands.fill_at(+r: Nat, +c: Nat, +i: Nat, +n: Nat, +xs: List<&2, F32>, b: Bands<r, c>) -> Bands<r, c>:
1440	  bfill_top(Nat.is_le(Nat.add(i, n), Nat.mul(r, c)), r, c, b, i, n, xs)
1441	
1442	# ---- y = B · x for B with r rows of c numbers (a weight matrix stored out x in) ----
1443	
1444	# B (intact) and y = B · x, typed like Mat.matmul_nt: x is 1 x c, y is 1 x r
1445	type BV<-r: Nat, -c: Nat> is Type:
1446	  BV{b: Bands<r, c>, y: Mat<1n, r>}
1447	
1448	# the bands and y as a list in row order (the result of Bands.matvec_l)
1449	type BL<-r: Nat, -c: Nat> is Type:
1450	  BL{b: Bands<r, c>, y: List<&2, F32>}
1451	
1452	def bv_g2(-r: Nat, -c: Nat, w: Array<F32>, al: AL) -> BL<r, c>:
1453	  match al:
1454	    case AL{d, l}:
1455	      BL{BLeaf{Mat{w}}, l}
1456	
1457	# the band's 1 + m rows came out in an Array (one number per row); read them back as a list
1458	def bv_g(-r: Nat, -c: Nat, +m: Nat, g: G) -> BL<r, c>:
1459	  match g:
1460	    case G{a, b, y}:
1461	      bv_g2(r, c, b, read_l(m, 0, Nil{}, Array.get(F32, y, 0)))
1462	
1463	# the rows of one band, one dot product each, with the gemm of Mat.matmul_nt (n = 1), so the
1464	# results are identical. Writing the results into an Array and reading them back is faster than
1465	# building the list row by row (NOTES.md, exp. 14).
1466	def bv_gemm(-r: Nat, +c: Nat, rows: Nat, +rr: Nat, x: Array<F32>, w: Array<F32>) -> BL<r, c>:
1467	  match rows:
1468	    case 0n:
1469	      BL{BLeaf{Mat{w}}, Nil{}}
1470	    case 1n+p:
1471	      bv_g(r, c, p, gemm(1n, rr, c, U32.from_nat(c), 1, 1, U32.from_nat(c), x, w, Array.new(F32, cap_depth(rr), 0.0)))
1472	
1473	def bv_mat(-r: Nat, +c: Nat, +rr: Nat, x: Array<F32>, m: Mat<r, c>) -> BL<r, c>:
1474	  match m:
1475	    case Mat{w}:
1476	      bv_gemm(r, c, rr, rr, x, w)
1477	
1478	def bv_join(-r: Nat, -c: Nat, p: BL<half(r), c>, q: BL<Nat.sub(r, half(r)), c>) -> BL<r, c>:
1479	  match p q:
1480	    case BL{bx, yx} BL{by, yy}:
1481	      BL{BNode{bx, by}, List.append(&2, F32, yx, yy)}
1482	
1483	# the two bands of a node run in parallel, each with its own copy of x (c numbers);
1484	# xp is a pair of copies of x (a leaf uses one). xp comes after the erased e: up to Bend 2.0.36 an erased
1485	# parameter at the end of the list turned the parallel let into two sequential calls (NOTES.md, exp. 11; fixed in
1486	# 2.0.37 by bendlang/bend#1377). The order is kept, and harmless.
1487	def bmv(+c: Nat, -r: Nat, b: Bands<r, c>, +rr: Nat, -e: {rr == r : Nat}, xp: Array<F32> & Array<F32>) -> BL<r, c>:
1488	  match b xp:
1489	    case BLeaf{m} Tuple{x, x2}:
1490	      bv_mat(r, c, rr, x, m)
1491	    case BNode{bx, by} Tuple{x1, x2}:
1492	      p q = bmv(c, half(r), bx, half(rr), eq_half(rr, r, e), Array.clone(F32, x1)) bmv(c, Nat.sub(r, half(r)), by, Nat.sub(rr, half(rr)), eq_rest(rr, r, e), Array.clone(F32, x2))
1493	      bv_join(r, c, p, q)
1494	
1495	def bv_fin(+r: Nat, -c: Nat, bl: BL<r, c>) -> BV<r, c>:
1496	  match bl:
1497	    case BL{b, y}:
1498	      BV{b, Mat.from_list(1n, r, y)}
1499	
1500	# y = B · x for B (r x c, one row per output) and x (1 x c); returns B intact and y (1 x r).
1501	# The same product as Mat.matmul_nt(1n, c, r, ..), and the same numbers, bit for bit.
1502	def Bands.matvec(+r: Nat, +c: Nat, b: Bands<r, c>, x: Mat<1n, c>) -> BV<r, c>:
1503	  match x:
1504	    case Mat{xa}:
1505	      bv_fin(r, c, bmv(c, r, b, r, {==}, Array.clone(F32, xa)))
1506	
1507	# the same product with x and y as lists (x: c numbers, y: r numbers in row order), without the
1508	# Mat conversions of Bands.matvec: for callers that already hold the vector as a list
1509	def Bands.matvec_l(+r: Nat, +c: Nat, b: Bands<r, c>, +x: List<&2, F32>) -> BL<r, c>:
1510	  bmv(c, r, b, r, {==}, Array.clone(F32, fill_list(x, 0, Array.new(F32, cap_depth(c), 0.0))))
1511	
1512	# ---- reading rows and the whole matrix ----
1513	
1514	type BRow<-r: Nat, -c: Nat> is Type:
1515	  BRow{b: Bands<r, c>, row: List<&2, F32>}
1516	
1517	def brow_fin(-r: Nat, -c: Nat, al: AL) -> BRow<r, c>:
1518	  match al:
1519	    case AL{d, l}:
1520	      BRow{BLeaf{Mat{d}}, l}
1521	
1522	# n numbers starting at flat index i of one band (a band may have no rows)
1523	def brow_leaf2(-r: Nat, -c: Nat, n: Nat, +i: U32, d: Array<F32>) -> BRow<r, c>:
1524	  match n:
1525	    case 0n:
1526	      BRow{BLeaf{Mat{d}}, Nil{}}
1527	    case 1n+p:
1528	      brow_fin(r, c, read_l(p, i, Nil{}, Array.get(F32, d, i)))
1529	
1530	def brow_leaf(-r: Nat, -c: Nat, m: Mat<r, c>, +i: U32, n: Nat) -> BRow<r, c>:
1531	  match m:
1532	    case Mat{d}:
1533	      brow_leaf2(r, c, n, i, d)
1534	
1535	def brow_l(-r: Nat, -c: Nat, p: BRow<half(r), c>, y: Bands<Nat.sub(r, half(r)), c>) -> BRow<r, c>:
1536	  match p:
1537	    case BRow{bx, row}:
1538	      BRow{BNode{bx, y}, row}
1539	
1540	def brow_r(-r: Nat, -c: Nat, x: Bands<half(r), c>, q: BRow<Nat.sub(r, half(r)), c>) -> BRow<r, c>:
1541	  match q:
1542	    case BRow{by, row}:
1543	      BRow{BNode{x, by}, row}
1544	
1545	# row i; left says whether i is in the first band of the node (i < half(rr))
1546	def brow(+c: Nat, -r: Nat, b: Bands<r, c>, left: Bool, +rr: Nat, -e: {rr == r : Nat}, +i: Nat) -> BRow<r, c>:
1547	  match b left:
1548	    case BLeaf{m} _:
1549	      brow_leaf(r, c, m, U32.from_nat(Nat.mul(i, c)), c)
1550	    case BNode{x, y} True{}:
1551	      brow_l(r, c, brow(c, half(r), x, Nat.is_lt(i, half(half(rr))), half(rr), eq_half(rr, r, e), i), y)
1552	    case BNode{x, y} False{}:
1553	      brow_r(r, c, x, brow(c, Nat.sub(r, half(r)), y, Nat.is_lt(Nat.sub(i, half(rr)), half(Nat.sub(rr, half(rr)))), Nat.sub(rr, half(rr)), eq_rest(rr, r, e), Nat.sub(i, half(rr))))
1554	
1555	# row i (c numbers), for example an embedding lookup
1556	def Bands.read_row(+r: Nat, +c: Nat, b: Bands<r, c>, +i: Nat) -> BRow<r, c>:
1557	  brow(c, r, b, Nat.is_lt(i, half(r)), r, {==}, i)
1558	
1559	def bl_join(-r: Nat, -c: Nat, p: BRow<half(r), c>, q: BRow<Nat.sub(r, half(r)), c>) -> BRow<r, c>:
1560	  match p q:
1561	    case BRow{bx, xs} BRow{by, ys}:
1562	      BRow{BNode{bx, by}, List.append(&2, F32, xs, ys)}
1563	
1564	def blist(+c: Nat, -r: Nat, b: Bands<r, c>, +rr: Nat, -e: {rr == r : Nat}) -> BRow<r, c>:
1565	  match b:
1566	    case BLeaf{m}:
1567	      brow_leaf(r, c, m, 0, Nat.mul(rr, c))
1568	    case BNode{x, y}:
1569	      bl_join(r, c, blist(c, half(r), x, half(rr), eq_half(rr, r, e)), blist(c, Nat.sub(r, half(r)), y, Nat.sub(rr, half(rr)), eq_rest(rr, r, e)))
1570	
1571	# all r*c numbers in row order (in BRow.row)
1572	def Bands.to_list(+r: Nat, +c: Nat, b: Bands<r, c>) -> BRow<r, c>:
1573	  blist(c, r, b, r, {==})
1574	
1575	# ---- the list path returns exactly r numbers ----
1576	#
1577	# The numbers do not matter here, only how many there are: each lemma below counts the list a
1578	# function builds, following the function's own recursion.
1579	
1580	def bl_len(-r: Nat, -c: Nat, bl: BL<r, c>) -> Nat:
1581	  match bl:
1582	    case BL{b, y}:
1583	      List.length(&2, F32, y)
1584	
1585	# s == a + (1 + b)  ->  s == 1 + (a + b)
1586	def succ_out(-s: Nat, +a: Nat, +b: Nat, e: {s == Nat.add(a, 1n+b) : Nat}) -> {s == 1n+Nat.add(a, b) : Nat}:
1587	  %Equal.sym(Nat, 1n+Nat.add(a, b), Nat.add(a, 1n+b), add_succ_r(a, b)) : {s == _ : Nat}
1588	  e
1589	
1590	# s == 1 + (p + (1 + l))  ->  s == 2 + (p + l)
1591	def len_step(-s: Nat, +p: Nat, +l: Nat, e: {s == 1n+Nat.add(p, 1n+l) : Nat}) -> {s == 2n+Nat.add(p, l) : Nat}:
1592	  %Equal.sym(Nat, 1n+Nat.add(p, l), Nat.add(p, 1n+l), add_succ_r(p, l)) : {s == 1n+_ : Nat}
1593	  e
1594	
1595	# s == n + 0  ->  s == n
1596	def zero_out(-s: Nat, +n: Nat, e: {s == Nat.add(n, 0n) : Nat}) -> {s == n : Nat}:
1597	  %add_zero_r(n) : {s == _ : Nat}
1598	  e
1599	
1600	# reversing onto acc adds the lengths
1601	def rev_len(+xs: List<&2, F32>, +acc: List<&2, F32>) -> {List.length(&2, F32, List.reverse.go(&2, F32, xs, acc)) == Nat.add(List.length(&2, F32, xs), List.length(&2, F32, acc)) : Nat}:
1602	  match xs:
1603	    case Nil{}:
1604	      {==}
1605	    case Con{h, t}:
1606	      succ_out(List.length(&2, F32, List.reverse.go(&2, F32, t, h <> acc)), List.length(&2, F32, t), List.length(&2, F32, acc), rev_len(t, h <> acc))
1607	
1608	# appending adds the lengths
1609	def app_len(+xs: List<&2, F32>, +ys: List<&2, F32>) -> {List.length(&2, F32, List.append(&2, F32, xs, ys)) == Nat.add(List.length(&2, F32, xs), List.length(&2, F32, ys)) : Nat}:
1610	  match xs:
1611	    case Nil{}:
1612	      {==}
1613	    case Con{h, t}:
1614	      %app_len(t, ys) : {1n+List.length(&2, F32, List.append(&2, F32, t, ys)) == 1n+_ : Nat}
1615	      {==}
1616	
1617	# a leaf's list: read_l reads mleft + 1 numbers onto acc and reverses them. By induction on
1618	# mleft: with 0 left it reads one number (rev_len counts the reversed list); with 1 + p left it
1619	# reads one and recurses, and len_step moves the 1 out of the sum.
1620	law leaf_len:
1621	  for -r: Nat
1622	  for -c: Nat
1623	  for w: Array<F32>
1624	  for +mleft: Nat
1625	  for +idx: U32
1626	  for +acc: List<&2, F32>
1627	  for rc: Array<F32> & F32
1628	  {bl_len(r, c, bv_g2(r, c, w, read_l(mleft, idx, acc, rc))) == Nat.add(1n+mleft, List.length(&2, F32, acc)) : Nat}
1629	
1630	def leaf_len(r, c, w, mleft, idx, acc, rc):
1631	  match mleft rc:
1632	    case 0n Tuple{d, +v}:
1633	      zero_out(List.length(&2, F32, List.reverse.go(&2, F32, v <> acc, Nil{})), 1n+List.length(&2, F32, acc), rev_len(v <> acc, Nil{}))
1634	    case 1n+p Tuple{d, +v}:
1635	      len_step(bl_len(r, c, bv_g2(r, c, w, read_l(p, (idx + 1 : U32), v <> acc, Array.get(F32, d, (idx + 1 : U32))))), p, List.length(&2, F32, acc), leaf_len(r, c, w, p, (idx + 1 : U32), v <> acc, Array.get(F32, d, (idx + 1 : U32))))
1636	
1637	# whatever the gemm wrote, reading 1 + p rows back gives 1 + p numbers (G has one constructor, so
1638	# matching g opens it, and leaf_len counts the list)
1639	def g_len(-r: Nat, -c: Nat, +p: Nat, g: G) -> {bl_len(r, c, bv_g(r, c, p, g)) == 1n+p : Nat}:
1640	  match g:
1641	    case G{a, b, y}:
1642	      zero_out(bl_len(r, c, bv_g2(r, c, b, read_l(p, 0, Nil{}, Array.get(F32, y, 0)))), 1n+p, leaf_len(r, c, b, p, 0, Nil{}, Array.get(F32, y, 0)))
1643	
1644	# LAW (band_len): the product of one band with `rows` rows gives exactly `rows` numbers.
1645	# 0 rows give the empty list; 1 + p rows go through the gemm, and g_len counts what is read back.
1646	law band_len:
1647	  for -r: Nat
1648	  for +c: Nat
1649	  for +rows: Nat
1650	  for x: Array<F32>
1651	  for w: Array<F32>
1652	  {bl_len(r, c, bv_gemm(r, c, rows, rows, x, w)) == rows : Nat}
1653	
1654	def band_len(r, c, rows, x, w):
1655	  match rows:
1656	    case 0n:
1657	      {==}
1658	    case 1n+p:
1659	      g_len(r, c, p, gemm(1n, 1n+p, c, U32.from_nat(c), 1, 1, U32.from_nat(c), x, w, Array.new(F32, cap_depth(1n+p), 0.0)))
1660	
1661	# joining two halves of h1 and h2 numbers gives h1 + h2 numbers. (The statement for the whole tree,
1662	# "Bands.matvec_l gives r numbers", would apply this to the results of the two recursive calls and
1663	# to the induction hypotheses about the same calls; that uses each band's Array twice, which Bend's
1664	# affine rules refuse. It stays open, covered by tests: NOTES.md, v3.1.)
1665	def join_len(-r: Nat, -c: Nat, -h1: Nat, -h2: Nat, p: BL<half(r), c>, q: BL<Nat.sub(r, half(r)), c>, hp: {bl_len(half(r), c, p) == h1 : Nat}, hq: {bl_len(Nat.sub(r, half(r)), c, q) == h2 : Nat}) -> {bl_len(r, c, bv_join(r, c, p, q)) == Nat.add(h1, h2) : Nat}:
1666	  match p q:
1667	    case BL{bx, +yx} BL{by, +yy}:
1668	      %hp : {List.length(&2, F32, List.append(&2, F32, yx, yy)) == Nat.add(_, h2) : Nat}
1669	      %hq : {List.length(&2, F32, List.append(&2, F32, yx, yy)) == Nat.add(List.length(&2, F32, yx), _) : Nat}
1670	      app_len(yx, yy)
1671