From c57acbedb32691307c2a6b8920b02e8ac31bb28c Mon Sep 17 00:00:00 2001
From: tannevaled <tannevaled@users.noreply.github.com>
Date: Sun, 18 Oct 2026 22:57:56 +0200
Subject: [PATCH 1/3] perf: float32 Add/Sub/Mul/Div kernels on six targets;
 *Into writes typed results in place

Against NumPy 2.5.3 on one POWER9 core, float32 Add was 7.1x slower at
2^14 (33.8 us against 4.74). A profile put 57% in binNum[float32]: a
scalar add per element with two bounds checks.

- typed.go: the elementwise loops reslice their operands to len(d), so
  the compiler drops the per-element bounds checks (checked with
  -d=ssa/check_bce): float32/int32 Add -7..-12% on POWER9.
- kernels.Bin32: 16 floats per iteration with SSE2 (amd64), NEON .S4
  (arm64), VSX xvaddsp... (ppc64le), z14 vfasb... (s390x, behind
  useVXE), LSX vfadd.s... (loong64, behind useLASX) and RVV E32
  (riscv64, behind useRVV). One IEEE operation per lane: bit-identical to
  the scalar loop. binFloat becomes binFloat32 (only float32 reaches it).
- arithInto: when out has the result dtype, write into it directly
  instead of computing a new array and copying it.
- asmgen/regen: one Go command regenerates all six .s files.
- benchmarks: BenchmarkDTypeAddInto.

POWER9: float32 Add -55% (2^14) / -68% (2^20). float32 AddInto 4.35 us
(NumPy out=: 4.23) and 294 us (NumPy: 366). The kernel tests and the
PROT_NONE fence pass on Zen 3, Apple M, POWER9, IBM Z, SpacemiT X60 and
Loongson 3C5000L; a sabotaged NEON kernel fails them.

Still open, written in docs/perf.md: the allocating Add is 1.5-3.2x
slower than NumPy's a+b, float64 included. The allocation and the GC
account for the difference, not the kernel.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 README.md                              |   7 +-
 SECURITY.md                            |  17 +++
 arith.go                               |  40 +++++-
 benchmarks/dtype_bench_test.go         |  29 +++-
 docs/perf.md                           |  56 ++++++++
 docs/plan-ndarray.md                   |   6 +-
 dtype_api_test.go                      |  18 +++
 gen_dtypes.go                          |   2 +-
 internal/kernels/asmgen/amd64/gen.go   |  40 ++++++
 internal/kernels/asmgen/arm64/gen.go   |  28 ++++
 internal/kernels/asmgen/loong64/gen.go |  31 ++++
 internal/kernels/asmgen/ppc64le/gen.go |  40 ++++++
 internal/kernels/asmgen/regen/main.go  |  40 ++++++
 internal/kernels/asmgen/riscv64/gen.go |  23 +++
 internal/kernels/asmgen/s390x/gen.go   |  38 +++++
 internal/kernels/bin32.go              |  54 +++++++
 internal/kernels/bin32_amd64.go        |  20 +++
 internal/kernels/bin32_arm64.go        |  20 +++
 internal/kernels/bin32_generic.go      |   7 +
 internal/kernels/bin32_loong64.go      |  25 ++++
 internal/kernels/bin32_ppc64le.go      |  20 +++
 internal/kernels/bin32_riscv64.go      |  24 ++++
 internal/kernels/bin32_s390x.go        |  25 ++++
 internal/kernels/bin32_test.go         |  83 +++++++++++
 internal/kernels/fence_test.go         |  16 +++
 internal/kernels/sum_amd64.s           | 176 +++++++++++++++++++++++
 internal/kernels/sum_arm64.s           | 112 +++++++++++++++
 internal/kernels/sum_loong64.s         | 172 ++++++++++++++++++++++
 internal/kernels/sum_ppc64le.s         | 188 +++++++++++++++++++++++++
 internal/kernels/sum_riscv64.s         |  76 ++++++++++
 internal/kernels/sum_s390x.s           | 140 ++++++++++++++++++
 typed.go                               |  21 ++-
 zz_dtypes_gen.go                       |   2 +-
 33 files changed, 1574 insertions(+), 22 deletions(-)
 create mode 100644 internal/kernels/asmgen/regen/main.go
 create mode 100644 internal/kernels/bin32.go
 create mode 100644 internal/kernels/bin32_amd64.go
 create mode 100644 internal/kernels/bin32_arm64.go
 create mode 100644 internal/kernels/bin32_generic.go
 create mode 100644 internal/kernels/bin32_loong64.go
 create mode 100644 internal/kernels/bin32_ppc64le.go
 create mode 100644 internal/kernels/bin32_riscv64.go
 create mode 100644 internal/kernels/bin32_s390x.go
 create mode 100644 internal/kernels/bin32_test.go

diff --git a/README.md b/README.md
index 3e524e5..60ee7e3 100644
--- a/README.md
+++ b/README.md
@@ -14,7 +14,9 @@
 - **Thirteen dtypes** — bool, int8…int64, uint8…uint64, float32/float64,
   complex64/complex128, chosen at run time with **NumPy 2 promotion** (`AsType`,
   `ZerosOf`, `FromSlice`, typed access with `Data[T]`/`Item[T]`, weak Go
-  scalars with `Scalar`). Float64 is the default and the SIMD-accelerated one
+  scalars with `Scalar`). Float64 is the default and the most SIMD-accelerated
+  one. Float32 has SIMD `Sum` and `Add`/`Sub`/`Mul`/`Div` on every 64-bit
+  target, bit-identical to its scalar loops
   ([docs/design-dtypes.md](docs/design-dtypes.md)).
 - **Creation** — `New`/`Zeros`/`Ones`/`Full`/`FromData`/`Arange`/`Linspace`/
   `Eye`/`Identity`.
@@ -150,7 +152,8 @@ and `Max` by 92–96% ([docs/perf.md](docs/perf.md), *SIMD coverage*). It is a
 > pair Go supports; v0.1.0 and v0.2.0 were also run on real amd64, arm64,
 > ppc64le, riscv64 and loong64 hardware, v0.3.0 on real amd64, arm64, ppc64le
 > and loong64. See **[docs/plan-ndarray.md](docs/plan-ndarray.md)**
-> for the roadmap (SIMD for float32, linear algebra decompositions).
+> for the roadmap (float32 SIMD for sqrt, max/min and the products; linear
+> algebra decompositions).
 
 ## Why this module?
 
diff --git a/SECURITY.md b/SECURITY.md
index e3fcd76..477f677 100644
--- a/SECURITY.md
+++ b/SECURITY.md
@@ -66,6 +66,23 @@ standard library's known vulnerabilities; `govulncheck ./...` reports none.
 
 ## Audit log
 
+**2026-10-18 (v0.18.0).** The float32 elementwise kernels (`kernels.Bin32`,
+on all six 64-bit targets) take raw pointers and a count. `Bin32` reslices `a`
+and `b` to `len(dst)` before the call, so a short operand panics in Go rather
+than being read past its end, and an empty `dst` never reaches a kernel.
+`TestBin32StaysInsideItsOperands` puts `dst`, `a` and `b` against PROT_NONE
+pages at every length from 1 to 40, and passes on real hardware for each
+kernel: Zen 3, Apple M, POWER9, IBM Z, SpacemiT X60 and Loongson 3C5000L.
+A sabotaged kernel failed the bit-identity test. The `*Into` direct path
+keeps the contiguity and shape checks of the copying one, and computes in
+place, as the float64 path does. A partial overlap between `out` and an
+operand cannot reach it. `out` must be contiguous, which a view with an
+offset (`x[1:]`, even `x[:39]`) is not, so it is refused. An operand sharing
+`out`'s storage with `out`'s shape is then `out` itself, the exact alias,
+which a lane-wise kernel handles. Any other view of it is gathered into a
+copy first. A probe of `x[1:] = x[:-1] + 1` and `x[:39] = x[1:40] + 1` got
+the contiguity error in both directions, for float64 and float32.
+
 **2026-10-18 (v0.16.0).** The float32 pairwise-sum block kernels (SSE2, NEON,
 RVV, LSX since v0.14.0/v0.15.0, VSX now) take raw pointers and a block count.
 `AddBlocks8F32` refuses a length that is not a positive multiple of 8, and
diff --git a/arith.go b/arith.go
index 6044dac..d3d8225 100644
--- a/arith.go
+++ b/arith.go
@@ -59,12 +59,9 @@ func (a *Array) arith(b *Array, op binCode, k func(dst, x, y []float64)) (*Array
 	if a.dtype == Float64 && b.dtype == Float64 {
 		return a.binOp(b, k)
 	}
-	dt := promote(a, b)
-	switch {
-	case op == opDiv && (dt == Bool || dt.IsInteger()):
-		dt = Float64 // true division
-	case op == opSub && dt == Bool:
-		return nil, fmt.Errorf("%w: subtracting bool arrays; use LogicalXor", ErrDType)
+	dt, err := arithDType(a, b, op)
+	if err != nil {
+		return nil, err
 	}
 	if dt == Float64 {
 		return a.toFloat64().binOp(b.toFloat64(), k)
@@ -78,12 +75,41 @@ func (a *Array) arith(b *Array, op binCode, k func(dst, x, y []float64)) (*Array
 	return a.result(d, append([]int(nil), shape...)), nil
 }
 
+// arithDType is the result dtype of a OP b: NumPy's promotion, with true
+// division of integers and bools in Float64, and no bool subtraction.
+func arithDType(a, b *Array, op binCode) (DType, error) {
+	dt := promote(a, b)
+	switch {
+	case op == opDiv && (dt == Bool || dt.IsInteger()):
+		dt = Float64 // true division
+	case op == opSub && dt == Bool:
+		return dt, fmt.Errorf("%w: subtracting bool arrays; use LogicalXor", ErrDType)
+	}
+	return dt, nil
+}
+
 // arithInto is arith writing into out, which must be contiguous and have the
-// broadcast shape; the result is converted to out's dtype.
+// broadcast shape; the result is converted to out's dtype. When out already
+// has the result dtype the operation writes into it directly, with no
+// intermediate array.
 func (a *Array) arithInto(out, b *Array, op binCode, k func(dst, x, y []float64)) error {
 	if a.dtype == Float64 && b.dtype == Float64 && out.dtype == Float64 {
 		return a.binOpInto(out, b, k)
 	}
+	if dt, err := arithDType(a, b, op); err == nil && dt != Float64 && dt == out.dtype {
+		shape, err := broadcastShape(a.shape, b.shape)
+		if err != nil {
+			return err
+		}
+		if !out.isContiguous() {
+			return fmt.Errorf("%w: out must be contiguous", ErrBroadcast)
+		}
+		if !sameShape(out.shape, shape) {
+			return fmt.Errorf("%w: out shape %v != result shape %v", ErrBroadcast, out.shape, shape)
+		}
+		arithStore(op, out.store(), a.operandStore(dt, shape), b.operandStore(dt, shape))
+		return nil
+	}
 	r, err := a.arith(b, op, k)
 	if err != nil {
 		return err
diff --git a/benchmarks/dtype_bench_test.go b/benchmarks/dtype_bench_test.go
index 7fdacac..c7629e0 100644
--- a/benchmarks/dtype_bench_test.go
+++ b/benchmarks/dtype_bench_test.go
@@ -7,11 +7,10 @@ import (
 	nd "github.com/go-ndarray/ndarray"
 )
 
-// The same reduction and elementwise operation on the dtypes that have SIMD
-// kernels (float64) and on those that run the pure-Go loops (float32 and the
-// integers), at a cache-resident and a memory-bound size. Run them with
-// GOMAXPROCS=1 to compare kernels rather than the multicore fan-out. They
-// measure what a float32 kernel would be worth before anyone writes one.
+// The same reduction and elementwise operation on the float dtypes, which
+// have SIMD kernels, and on the integers, which run pure-Go loops, at a
+// cache-resident and a memory-bound size. Run them with GOMAXPROCS=1 to
+// compare kernels rather than the multicore fan-out.
 
 var dtypeBenchTypes = []nd.DType{nd.Float64, nd.Float32, nd.Int64, nd.Int32}
 
@@ -42,3 +41,23 @@ func BenchmarkDTypeAdd(b *testing.B) {
 		}
 	}
 }
+
+// BenchmarkDTypeAddInto is Add without the allocation, the analogue of
+// np.add(x, y, out=out): it measures the kernel, where BenchmarkDTypeAdd also
+// measures allocating and clearing the result.
+func BenchmarkDTypeAddInto(b *testing.B) {
+	for _, dt := range dtypeBenchTypes {
+		for _, n := range []int{1 << 14, 1 << 20} {
+			x, y := vec(n).AsType(dt), vec(n).AsType(dt)
+			out := vec(n).AsType(dt)
+			b.Run(fmt.Sprintf("%v/n=%d", dt, n), func(b *testing.B) {
+				b.SetBytes(int64(n * dt.Size()))
+				for i := 0; i < b.N; i++ {
+					if err := x.AddInto(out, y); err != nil {
+						b.Fatal(err)
+					}
+				}
+			})
+		}
+	}
+}
diff --git a/docs/perf.md b/docs/perf.md
index 37daa03..7354a6d 100644
--- a/docs/perf.md
+++ b/docs/perf.md
@@ -929,6 +929,62 @@ layout, not proven.
 
 Every 64-bit target now adds the float32 block with SIMD.
 
+### float32 elementwise arithmetic (2026-10-18)
+
+Measured against NumPy 2.5.3 on one POWER9 core, float32 `Add` was
+**7.1× slower** at 2¹⁴ (33.8 µs against 4.74 µs) and 5.0× slower at 2²⁰. A
+profile put 57% of the time in the generic loop `binNum[float32]`: a scalar
+add per element, with two bounds checks, because nothing told the compiler
+that `x` and `y` were as long as `d`.
+
+Two changes:
+
+- **Bounds checks.** The elementwise loops reslice their operands to
+  `len(d)` once, so the compiler drops the per-element checks (verified with
+  `-d=ssa/check_bce`). On the POWER9, float32 and int32 `Add` gain 7–12%.
+- **Kernels.** `kernels.Bin32` adds, subtracts, multiplies or divides 16
+  floats per iteration on every 64-bit target:
+
+  | target | instructions |
+  |---|---|
+  | amd64 | SSE2 `ADDPS`… |
+  | arm64 | NEON `.S4` |
+  | ppc64le | VSX `xvaddsp`… |
+  | s390x | z14 `vfasb`… |
+  | loong64 | LSX `vfadd.s`… |
+  | riscv64 | RVV E32, strip-mined |
+
+  Each lane is one IEEE operation, so the result is the scalar loop's bit
+  for bit. The kernel tests (special values, every length to 70, in place)
+  and the PROT_NONE fence pass on Zen 3, Apple M, POWER9, IBM Z, SpacemiT X60
+  and Loongson 3C5000L. A sabotaged kernel (one `vfsub` in `add32NEON`)
+  failed them.
+
+`AddInto` and the other `*Into` methods also wrote a typed result into a new
+array and then copied it into `out`. When `out` has the result dtype they now
+write into it directly.
+
+POWER9, one core, against NumPy 2.5.3:
+
+| float32 | before | after | NumPy |
+|---|---:|---:|---:|
+| `Add`, 2¹⁴ | 33.8 µs | 15.2 µs | 4.74 µs (`a+b`) |
+| `Add`, 2²⁰ | 1.70 ms | 0.537 ms | 0.342 ms |
+| `AddInto`, 2¹⁴ | — | **4.35 µs** | 4.23 µs (`out=`) |
+| `AddInto`, 2²⁰ | — | **294 µs** | 366 µs |
+
+float64 `AddInto` against `np.add(out=)`: 5.70 against 7.75 µs (2¹⁴), and
+724 against 768 µs (2²⁰).
+
+**Open defect: the allocating operations.** `Add` is still 1.5–3.2× slower
+than NumPy's `a+b`, for float64 as well (21.9 against 8.48 µs at 2¹⁴). The
+kernel is no longer the difference, as `AddInto` shows; the allocation is.
+NumPy's allocator hands back the buffer the previous result just freed,
+still in cache. Go allocates fresh memory, zeroes it, and pays for the GC
+cycles that follow; on the POWER9 the profile of float64 `Add` is half GC.
+Until that is addressed, a hot loop should use the `*Into` methods with a
+reused `out`.
+
 ## SIMD coverage
 
 - **amd64 (SSE2)** ships hand-vectorized `sum` (4-accumulator `ADDPD`), `sqrt`
diff --git a/docs/plan-ndarray.md b/docs/plan-ndarray.md
index 13ee7d4..3aec3ba 100644
--- a/docs/plan-ndarray.md
+++ b/docs/plan-ndarray.md
@@ -154,8 +154,10 @@ suite (plus the per-arch CI jobs already wired in `.github/workflows/ci.yml`).
 - **Phase 2 — dtypes — DONE.** Thirteen dtypes (bool, signed and unsigned
   integers, float32/64, complex64/128) as a runtime property of `Array`,
   with NumPy 2.x promotion; see [design-dtypes.md](design-dtypes.md). The
-  `fft` bridge shipped in v0.10.0 (package `fft`, over go-fft). Still open:
-  SIMD kernels for `float32` (float32 runs the pure-Go code today).
+  `fft` bridge shipped in v0.10.0 (package `fft`, over go-fft). float32 has
+  SIMD kernels on all six 64-bit targets for `Sum` (v0.14.0–v0.17.0) and for
+  `Add`/`Sub`/`Mul`/`Div` (v0.18.0). Still pure Go for float32: `Sqrt`,
+  `Max`/`Min`, `Dot`/`MatMul` and the transcendental ufuncs.
 - **Phase 3 — broadcasting ufuncs — MOSTLY DONE.** Broadcasting arithmetic,
   comparisons returning Bool masks, `Maximum`/`Minimum`, `Logical*`, `Power`,
   `Where` and the per-axis reductions, in every dtype. Not done: a public
diff --git a/dtype_api_test.go b/dtype_api_test.go
index 765281d..29b90e7 100644
--- a/dtype_api_test.go
+++ b/dtype_api_test.go
@@ -267,6 +267,24 @@ func TestTypedArithmeticEdges(t *testing.T) {
 	if err := i.SqrtInto(must(Zeros(2))); !errors.Is(err, ErrDType) {
 		t.Errorf("SqrtInto int err = %v", err)
 	}
+	// Into an out of the result dtype: written in place, same checks.
+	s := must(FromSlice([]float32{1.5, -2}, 2))
+	out32 := must(ZerosOf(Float32, 2))
+	if err := s.AddInto(out32, s); err != nil || !reflect.DeepEqual(data[float32](t, out32), []float32{3, -4}) {
+		t.Errorf("float32 AddInto = %v, %v", out32, err)
+	}
+	if err := s.DivInto(s, out32); err != nil || !reflect.DeepEqual(data[float32](t, s), []float32{0.5, 0.5}) {
+		t.Errorf("float32 DivInto in place = %v, %v", s, err)
+	}
+	if err := s.AddInto(out32, must(ZerosOf(Float32, 3))); !errors.Is(err, ErrBroadcast) {
+		t.Errorf("float32 AddInto broadcast err = %v", err)
+	}
+	if err := s.AddInto(must(ZerosOf(Float32, 3)), s); !errors.Is(err, ErrBroadcast) {
+		t.Errorf("float32 AddInto shape err = %v", err)
+	}
+	if err := s.AddInto(must(must(ZerosOf(Float32, 4)).Slice(Step(2))), s); !errors.Is(err, ErrBroadcast) {
+		t.Errorf("float32 AddInto strided err = %v", err)
+	}
 	// Float64 result from an int: the float path runs.
 	f := must(FromData([]float64{0.5, 0.25}, 2))
 	if got := must(i.Mul(f)).materialize(); !reflect.DeepEqual(got, []float64{0.5, 0.5}) {
diff --git a/gen_dtypes.go b/gen_dtypes.go
index 35728b6..9324afe 100644
--- a/gen_dtypes.go
+++ b/gen_dtypes.go
@@ -266,7 +266,7 @@ func main() {
 	p("func arithStore(op binCode, d, x, y any) {")
 	p("switch d := d.(type) {")
 	each("iu", "case []$T: binInt(op, d, x.([]$T), y.([]$T))")
-	each("f-float64", "case []$T: binFloat(op, d, x.([]$T), y.([]$T))")
+	each("f-float64", "case []$T: binFloat32(op, d, x.([]$T), y.([]$T))")
 	each("c", "case []$T: binComplex(op, d, x.([]$T), y.([]$T))")
 	each("b", "case []$T: binBool(op, d, x.([]$T), y.([]$T))")
 	p("}")
diff --git a/internal/kernels/asmgen/amd64/gen.go b/internal/kernels/asmgen/amd64/gen.go
index 6a47818..c1ebc62 100644
--- a/internal/kernels/asmgen/amd64/gen.go
+++ b/internal/kernels/asmgen/amd64/gen.go
@@ -50,6 +50,10 @@ func main() {
 	f.Add(gemmSSE2HalfKernel())
 	f.Add(amd64.FeatureProbe("hasFMA", amd64.FMA))
 	f.Add(blocks8Kernel())
+	f.Add(bin32Kernel("add32SSE2", "ADDPS", "ADDSS"))
+	f.Add(bin32Kernel("sub32SSE2", "SUBPS", "SUBSS"))
+	f.Add(bin32Kernel("mul32SSE2", "MULPS", "MULSS"))
+	f.Add(bin32Kernel("div32SSE2", "DIVPS", "DIVSS"))
 
 	if err := os.WriteFile("sum_amd64.s", []byte(f.String()), 0o644); err != nil {
 		fmt.Fprintln(os.Stderr, err)
@@ -325,6 +329,42 @@ func binKernel(name, vec, sca string) *emit.Function {
 	return b.Func()
 }
 
+// bin32Kernel builds NAME(dst, a, b *float32, n int): dst[i] = a[i] OP b[i],
+// 16 floats per iteration over four XMM registers with the packed op, then a
+// scalar tail with sca. Each lane is one IEEE float32 operation, so the result
+// is the scalar loop's bit for bit.
+func bin32Kernel(name, vec, sca string) *emit.Function {
+	sig := amd64.Layout(
+		[]string{"dst", "a", "b", "n"},
+		[]amd64.Type{amd64.Ptr, amd64.Ptr, amd64.Ptr, amd64.Int64},
+		nil, nil,
+	)
+	b := amd64.NewFunc(name, sig, 0)
+	b.LoadArg("dst", "DI").LoadArg("a", "SI").LoadArg("b", "DX").LoadArg("n", "CX")
+	b.Raw("fblock:").Raw("CMPQ CX, $16").Raw("JL ftail")
+	for r := 0; r < 4; r++ {
+		b.Raw("MOVUPS %d(SI), X%d", 16*r, r)
+	}
+	for r := 0; r < 4; r++ {
+		b.Raw("MOVUPS %d(DX), X%d", 16*r, 4+r)
+	}
+	for r := 0; r < 4; r++ {
+		b.Raw("%s X%d, X%d", vec, 4+r, r) // X(r) = a OP b
+	}
+	for r := 0; r < 4; r++ {
+		b.Raw("MOVUPS X%d, %d(DI)", r, 16*r)
+	}
+	b.Raw("ADDQ $64, SI").Raw("ADDQ $64, DX").Raw("ADDQ $64, DI").Raw("SUBQ $16, CX").Raw("JMP fblock").
+		Raw("ftail:").Raw("TESTQ CX, CX").Raw("JZ fdone").
+		Raw("MOVSS (SI), X0").Raw("MOVSS (DX), X1").
+		Raw("%s X1, X0", sca).
+		Raw("MOVSS X0, (DI)").
+		Raw("ADDQ $4, SI").Raw("ADDQ $4, DX").Raw("ADDQ $4, DI").Raw("DECQ CX").Raw("JMP ftail").
+		Raw("fdone:")
+	b.Ret()
+	return b.Func()
+}
+
 // gemmMR, gemmNR are the micro-tile the packing in gemm.go produces on amd64:
 // 6 rows of packed A per k step, 8 columns of packed B.
 const (
diff --git a/internal/kernels/asmgen/arm64/gen.go b/internal/kernels/asmgen/arm64/gen.go
index 2fe444d..059d364 100644
--- a/internal/kernels/asmgen/arm64/gen.go
+++ b/internal/kernels/asmgen/arm64/gen.go
@@ -47,6 +47,10 @@ func main() {
 	f.Add(binKernel("subNEON", "VFSUB", "FSUBD"))
 	f.Add(binKernel("mulNEON", "VFMUL", "FMULD"))
 	f.Add(binKernel("divNEON", "VFDIV", "FDIVD"))
+	f.Add(bin32Kernel("add32NEON", "VFADD", "FADDS"))
+	f.Add(bin32Kernel("sub32NEON", "VFSUB", "FSUBS"))
+	f.Add(bin32Kernel("mul32NEON", "VFMUL", "FMULS"))
+	f.Add(bin32Kernel("div32NEON", "VFDIV", "FDIVS"))
 	f.Add(extremumKernel("maxNEON", "VFMAX", "FMAXD"))
 	f.Add(extremumKernel("minNEON", "VFMIN", "FMIND"))
 	f.Add(gemmKernel())
@@ -152,6 +156,30 @@ func binKernel(name, vop, sop string) *emit.Function {
 	return b.Func()
 }
 
+// bin32Kernel builds NAME(dst, a, b *float32, n int): dst[i] = a[i] OP b[i],
+// 16 floats per iteration on the .S4 arrangement (VOP Vm, Vn, Vd: Vd = Vn OP
+// Vm), then a scalar tail with sop. Each lane is one IEEE float32 operation,
+// so the result is the scalar loop's bit for bit.
+func bin32Kernel(name, vop, sop string) *emit.Function {
+	b := arm64.NewFunc(name, binSig(), 0)
+	b.LoadArg("dst", "R0").LoadArg("a", "R1").LoadArg("b", "R2").LoadArg("n", "R3")
+	b.Raw("block:").Raw("CMP $16, R3").Raw("BLT tail").
+		Raw("VLD1.P 64(R1), [V0.S4, V1.S4, V2.S4, V3.S4]"). // a
+		Raw("VLD1.P 64(R2), [V4.S4, V5.S4, V6.S4, V7.S4]")  // b
+	for r := 0; r < 4; r++ {
+		b.Raw(fmt.Sprintf("%s V%d.S4, V%d.S4, V%d.S4", vop, 4+r, r, r)) // a OP b
+	}
+	b.Raw("VST1.P [V0.S4, V1.S4, V2.S4, V3.S4], 64(R0)").
+		Raw("SUB $16, R3").Raw("B block").
+		Raw("tail:").Raw("CBZ R3, done").
+		Raw("FMOVS.P 4(R1), F0").Raw("FMOVS.P 4(R2), F1").
+		Raw(sop + " F1, F0, F0").
+		Raw("FMOVS.P F0, 4(R0)").
+		Raw("SUB $1, R3").Raw("B tail").
+		Raw("done:").Ret()
+	return b.Func()
+}
+
 // blocks8Kernel builds blocks8NEON(r, a *float32, blocks int): for each of
 // the blocks of eight float32 at a, in order, r[j] += a[8i+j]. The eight
 // accumulators live in V0 (r[0:4]) and V1 (r[4:8]); VFADD on .S4 is one IEEE
diff --git a/internal/kernels/asmgen/loong64/gen.go b/internal/kernels/asmgen/loong64/gen.go
index 6d5d157..c37003a 100644
--- a/internal/kernels/asmgen/loong64/gen.go
+++ b/internal/kernels/asmgen/loong64/gen.go
@@ -37,6 +37,8 @@ func main() {
 		{"mulLASX", "XVMULD", "MULD"}, {"divLASX", "XVDIVD", "DIVD"},
 	} {
 		f.Add(binKernel(k.name, k.vec, k.sca))
+		// "XVADDD" -> "VADDF", "ADDD" -> "ADDF": the LSX single-precision op.
+		f.Add(bin32Kernel(k.name[:3]+"32LSX", k.vec[1:len(k.vec)-1]+"F", k.sca[:len(k.sca)-1]+"F"))
 	}
 	f.Add(gemmKernel())
 	f.Add(blocks8Kernel())
@@ -168,6 +170,35 @@ func binKernel(name, vec, sca string) *emit.Function {
 	return b.Func()
 }
 
+// bin32Kernel builds name(dst, a, b *float32, n int): dst[i] = a[i] OP b[i],
+// 16 floats per iteration over four 128-bit LSX registers (four independent
+// chains, as the block kernel's lesson on one LASX register taught), a scalar
+// tail with sca. Each lane is one IEEE float32 operation, bit-identical to the
+// scalar oracle.
+func bin32Kernel(name, vec, sca string) *emit.Function {
+	sig := loong64.Layout([]string{"dst", "a", "b", "n"},
+		[]loong64.Type{loong64.Ptr, loong64.Ptr, loong64.Ptr, loong64.Int64}, nil, nil)
+	b := loong64.NewFunc(name, sig, 0)
+	b.LoadArg("dst", "R4").LoadArg("a", "R6").LoadArg("b", "R7").LoadArg("n", "R5")
+	b.Label("floop").Raw("MOVV $16, R12").Raw("BLT R5, R12, ftail")
+	for k := 0; k < 4; k++ {
+		b.Raw("VMOVQ %d(R6), V%d", 16*k, k)
+		b.Raw("VMOVQ %d(R7), V%d", 16*k, 4+k)
+	}
+	for k := 0; k < 4; k++ {
+		b.Raw("%s V%d, V%d, V%d", vec, 4+k, k, 8+k) // V(8+k) = V(k) OP V(4+k)
+	}
+	for k := 0; k < 4; k++ {
+		b.Raw("VMOVQ V%d, %d(R4)", 8+k, 16*k)
+	}
+	b.Raw("ADDV $64, R4").Raw("ADDV $64, R6").Raw("ADDV $64, R7").Raw("ADDV $-16, R5").Raw("JMP floop")
+	b.Label("ftail").Raw("BEQ R5, fdone").
+		Raw("MOVF 0(R6), F1").Raw("MOVF 0(R7), F2").Raw("%s F2, F1", sca).Raw("MOVF F1, 0(R4)").
+		Raw("ADDV $4, R4").Raw("ADDV $4, R6").Raw("ADDV $4, R7").Raw("ADDV $-1, R5").Raw("JMP ftail")
+	b.Label("fdone").Ret()
+	return b.Func()
+}
+
 // gemmKernel builds gemmMicro8x8LASX(kc int, pa, pb, c *float64, ldc int):
 //
 //	for p in [0,kc):  C[r][col] += pa[p*8+r] * pb[p*8+col]   (r<8, col<8)
diff --git a/internal/kernels/asmgen/ppc64le/gen.go b/internal/kernels/asmgen/ppc64le/gen.go
index 7c10e9b..e7b241b 100644
--- a/internal/kernels/asmgen/ppc64le/gen.go
+++ b/internal/kernels/asmgen/ppc64le/gen.go
@@ -38,6 +38,7 @@ func main() {
 		name, vec, sca string
 	}{{"addVSX", "add", "FADD"}, {"subVSX", "sub", "FSUB"}, {"mulVSX", "mul", "FMUL"}, {"divVSX", "div", "FDIV"}} {
 		f.Add(binKernel(k.name, k.vec, k.sca))
+		f.Add(bin32Kernel(k.name[:3]+"32VSX", k.vec, k.sca+"S"))
 	}
 	f.Add(gemmKernel())
 	f.Add(blocks8Kernel())
@@ -194,6 +195,45 @@ func binKernel(name, op, sca string) *emit.Function {
 	return b.Func()
 }
 
+// bin32Kernel builds name(dst, a, b *float32, n int): dst[i] = a[i] OP b[i],
+// 16 floats per iteration with the VSX single-precision op, a scalar tail with
+// sca. LXVW4X/STXVW4X put each word at the same lane in every operand and back
+// at its own index. Each lane is one IEEE float32 operation, bit-identical to
+// the scalar oracle.
+func bin32Kernel(name, op, sca string) *emit.Function {
+	sig := ppc64.Layout([]string{"dst", "a", "b", "n"},
+		[]ppc64.Type{ppc64.Ptr, ppc64.Ptr, ppc64.Ptr, ppc64.Int64}, nil, nil)
+	b := ppc64.NewFunc(name, sig, 0)
+	b.LoadArg("dst", "R3").LoadArg("a", "R5").LoadArg("b", "R6").LoadArg("n", "R4")
+	offsets(b, false)
+	b.Label("floop").Raw("CMP R4, $16").Raw("BLT ftail")
+	for k := 0; k < 4; k++ {
+		b.Raw("LXVW4X (R5)(%s), VS%d", vidx(k), 32+k)
+		b.Raw("LXVW4X (R6)(%s), VS%d", vidx(k), 36+k)
+	}
+	for k := 0; k < 4; k++ {
+		switch op {
+		case "add":
+			b.XVADDSP(40+k, 32+k, 36+k)
+		case "sub":
+			b.XVSUBSP(40+k, 32+k, 36+k)
+		case "mul":
+			b.XVMULSP(40+k, 32+k, 36+k)
+		case "div":
+			b.XVDIVSP(40+k, 32+k, 36+k)
+		}
+	}
+	for k := 0; k < 4; k++ {
+		b.Raw("STXVW4X VS%d, (R3)(%s)", 40+k, vidx(k))
+	}
+	b.Raw("ADD $64, R3").Raw("ADD $64, R5").Raw("ADD $64, R6").Raw("ADD $-16, R4").Raw("BR floop")
+	b.Label("ftail").Raw("CMP R4, $0").Raw("BEQ fdone").
+		Raw("FMOVS 0(R5), F1").Raw("FMOVS 0(R6), F2").Raw("%s F2, F1", sca).Raw("FMOVS F1, 0(R3)").
+		Raw("ADD $4, R3").Raw("ADD $4, R5").Raw("ADD $4, R6").Raw("ADD $-1, R4").Raw("BR ftail")
+	b.Label("fdone").Ret()
+	return b.Func()
+}
+
 // gemmKernel builds gemmMicro8x8VSX(kc int, pa, pb, c *float64, ldc int):
 //
 //	for p in [0,kc):  C[r][col] += pa[p*8+r] * pb[p*8+col]   (r<8, col<8)
diff --git a/internal/kernels/asmgen/regen/main.go b/internal/kernels/asmgen/regen/main.go
new file mode 100644
index 0000000..c7f5497
--- /dev/null
+++ b/internal/kernels/asmgen/regen/main.go
@@ -0,0 +1,40 @@
+// Command regen runs every architecture's generator and moves its output into
+// the kernels package: the native equivalent of the per-arch go:generate
+// lines, which `go generate` only reaches for the host GOARCH.
+//
+// Run from internal/kernels/asmgen: go run ./regen
+package main
+
+import (
+	"fmt"
+	"os"
+	"os/exec"
+	"path/filepath"
+)
+
+var arches = []string{"amd64", "arm64", "loong64", "ppc64le", "riscv64", "s390x"}
+
+func main() {
+	if err := run("."); err != nil {
+		fmt.Fprintln(os.Stderr, "regen:", err)
+		os.Exit(1)
+	}
+}
+
+// run regenerates every architecture under root (the asmgen directory).
+func run(root string) error {
+	for _, a := range arches {
+		dir := filepath.Join(root, a)
+		cmd := exec.Command("go", "run", ".")
+		cmd.Dir = dir
+		cmd.Stdout, cmd.Stderr = os.Stderr, os.Stderr
+		if err := cmd.Run(); err != nil {
+			return fmt.Errorf("%s: %w", a, err)
+		}
+		out := "sum_" + a + ".s"
+		if err := os.Rename(filepath.Join(dir, out), filepath.Join(root, "..", out)); err != nil {
+			return fmt.Errorf("%s: %w", a, err)
+		}
+	}
+	return nil
+}
diff --git a/internal/kernels/asmgen/riscv64/gen.go b/internal/kernels/asmgen/riscv64/gen.go
index 4913f50..d631abb 100644
--- a/internal/kernels/asmgen/riscv64/gen.go
+++ b/internal/kernels/asmgen/riscv64/gen.go
@@ -34,6 +34,7 @@ func main() {
 		{"mulRVV", "VFMULVV"}, {"divRVV", "VFDIVVV"},
 	} {
 		f.Add(binKernel(k.name, k.op))
+		f.Add(bin32Kernel(k.name[:3]+"32RVV", k.op))
 	}
 	f.Add(gemmKernel())
 	f.Add(blocks8Kernel())
@@ -128,6 +129,28 @@ func binKernel(name, op string) *emit.Function {
 	return b.Func()
 }
 
+// bin32Kernel builds NAME(dst, a, b *float32, n int) for n >= 1: dst = a OP b
+// on 32-bit elements, strip-mined by VSETVLI (no tail loop). Each lane is one
+// IEEE float32 operation, so the result is the scalar loop's bit for bit.
+func bin32Kernel(name, op string) *emit.Function {
+	sig := riscv64.Layout([]string{"dst", "a", "b", "n"},
+		[]riscv64.Type{riscv64.Ptr, riscv64.Ptr, riscv64.Ptr, riscv64.Int64}, nil, nil)
+	b := riscv64.NewFunc(name, sig, 0)
+	b.LoadArg("dst", "X10").LoadArg("a", "X11").LoadArg("b", "X12").LoadArg("n", "X13")
+	b.Raw("loop:").
+		Raw("VSETVLI X13, E32, M8, TA, MA, X14").
+		Raw("VLE32V (X11), V0").
+		Raw("VLE32V (X12), V8").
+		Raw("%s V8, V0, V16", op). // V16 = a OP b
+		Raw("VSE32V V16, (X10)").
+		Raw("SLLI $2, X14, X17") // vl elements of 4 bytes
+	for _, p := range []string{"X10", "X11", "X12"} {
+		b.Raw("ADD X17, %s, %s", p, p)
+	}
+	b.Raw("SUB X14, X13, X13").Raw("BNEZ X13, loop").Ret()
+	return b.Func()
+}
+
 // gemmKernel builds gemmMicro6x8RVV(kc int, pa, pb, c *float64, ldc int) for
 // kc >= 1: C[r][0:8] += sum_p pa[p*6+r] * pb[p*8:p*8+8] for r < 6, the rows of C
 // ldc float64s apart.
diff --git a/internal/kernels/asmgen/s390x/gen.go b/internal/kernels/asmgen/s390x/gen.go
index aa46b0c..fb68d03 100644
--- a/internal/kernels/asmgen/s390x/gen.go
+++ b/internal/kernels/asmgen/s390x/gen.go
@@ -35,6 +35,7 @@ func main() {
 		{"mulVX", "VFMDB", "FMUL"}, {"divVX", "VFDDB", "FDIV"},
 	} {
 		f.Add(binKernel(k.name, k.vec, k.sca))
+		f.Add(bin32Kernel(k.name[:3]+"32VXE", k.name[:3], k.sca+"S"))
 	}
 	f.Add(gemmKernel())
 	f.Add(extremumKernel("maxVXE", "VFMAXDB", "WFMAXDB"))
@@ -143,6 +144,43 @@ func binKernel(name, vec, sca string) *emit.Function {
 	return b.Func()
 }
 
+// bin32Kernel builds NAME(dst, a, b *float32, n int): dst = a OP b, 16 floats
+// per iteration over four vector registers with the vector-enhancements
+// facility 1 single-precision op (z14; go-asmgen emits it as WORD+BYTE+BYTE,
+// the cmd/asm spelling is CL 846286), then a scalar tail with sca ("sca F1,
+// F0": F0 = F0 OP F1). Each lane is one IEEE float32 operation, so the result
+// is the scalar loop's bit for bit.
+func bin32Kernel(name, op, sca string) *emit.Function {
+	sig := s390x.Layout([]string{"dst", "a", "b", "n"},
+		[]s390x.Type{s390x.Ptr, s390x.Ptr, s390x.Ptr, s390x.Int64}, nil, nil)
+	b := s390x.NewFunc(name, sig, 0)
+	b.LoadArg("dst", "R1").LoadArg("a", "R2").LoadArg("b", "R3").LoadArg("n", "R4")
+	b.Raw("loop:").Raw("CMP R4, $16").Raw("BLT tail").
+		Raw("VLM (R2), V0, V3"). // a
+		Raw("VLM (R3), V4, V7")  // b
+	for r := 0; r < 4; r++ { // V_r = a OP b
+		switch op {
+		case "add":
+			b.VFASB(r, r, 4+r)
+		case "sub":
+			b.VFSSB(r, r, 4+r)
+		case "mul":
+			b.VFMSB(r, r, 4+r)
+		case "div":
+			b.VFDSB(r, r, 4+r)
+		}
+	}
+	b.Raw("VSTM V0, V3, (R1)").
+		Raw("ADD $64, R1").Raw("ADD $64, R2").Raw("ADD $64, R3").Raw("SUB $16, R4").Raw("BR loop").
+		Raw("tail:").Raw("CMP R4, $0").Raw("BEQ done").
+		Raw("FMOVS (R2), F0").Raw("FMOVS (R3), F1").
+		Raw(sca + " F1, F0").
+		Raw("FMOVS F0, (R1)").
+		Raw("ADD $4, R1").Raw("ADD $4, R2").Raw("ADD $4, R3").Raw("SUB $1, R4").Raw("BR tail").
+		Raw("done:").Ret()
+	return b.Func()
+}
+
 // gemmKernel builds gemmMicro4x8VX(kc int, pa, pb, c *float64, ldc int) for
 // kc >= 1: C[r][0:8] += sum_p pa[p*4+r] * pb[p*8:p*8+8] for r < 4, the rows of C
 // ldc float64s apart. Sixteen two-lane accumulators V0..V15 (row r in
diff --git a/internal/kernels/bin32.go b/internal/kernels/bin32.go
new file mode 100644
index 0000000..b7d2dd8
--- /dev/null
+++ b/internal/kernels/bin32.go
@@ -0,0 +1,54 @@
+package kernels
+
+// The float32 elementwise operations, in the order of the parent package's
+// binCode: Bin32 takes these codes.
+const (
+	Add32 = iota
+	Sub32
+	Mul32
+	Div32
+)
+
+// Bin32 writes a[i] OP b[i] into dst[i] for every element, OP one of Add32,
+// Sub32, Mul32 or Div32. a and b must be at least as long as dst. A vector
+// kernel runs where the architecture has one; each lane is one IEEE float32
+// operation, so the result is the scalar loop's bit for bit.
+func Bin32(op int, dst, a, b []float32) {
+	n := len(dst)
+	if n == 0 {
+		return
+	}
+	bin32(op, dst, a[:n], b[:n]) // the kernel reads n of each
+}
+
+// bin32Go is the scalar oracle the kernels are tested against, and the
+// implementation where there is no kernel.
+func bin32Go(op int, dst, a, b []float32) {
+	a, b = a[:len(dst)], b[:len(dst)]
+	switch op {
+	case Add32:
+		for i := range dst {
+			dst[i] = a[i] + b[i]
+		}
+	case Sub32:
+		for i := range dst {
+			dst[i] = a[i] - b[i]
+		}
+	case Mul32:
+		for i := range dst {
+			dst[i] = a[i] * b[i]
+		}
+	case Div32:
+		for i := range dst {
+			dst[i] = a[i] / b[i]
+		}
+	default:
+		panic("kernels: Bin32 op out of range")
+	}
+}
+
+// bin32Func is the shape of the generated float32 kernels: n >= 1 elements.
+type bin32Func func(dst, a, b *float32, n int)
+
+// run32 calls k on dst, a and b, all of length n >= 1.
+func run32(k bin32Func, dst, a, b []float32) { k(&dst[0], &a[0], &b[0], len(dst)) }
diff --git a/internal/kernels/bin32_amd64.go b/internal/kernels/bin32_amd64.go
new file mode 100644
index 0000000..f4a75ea
--- /dev/null
+++ b/internal/kernels/bin32_amd64.go
@@ -0,0 +1,20 @@
+package kernels
+
+// bin32 runs the SSE2 kernel for op (sum_amd64.s: ADDPS, SUBPS, MULPS,
+// DIVPS, 16 floats per iteration). SSE2 is the amd64 baseline. dst, a and b
+// have the same length, at least 1.
+func bin32(op int, dst, a, b []float32) {
+	run32([...]bin32Func{add32SSE2, sub32SSE2, mul32SSE2, div32SSE2}[op], dst, a, b)
+}
+
+//go:noescape
+func add32SSE2(dst, a, b *float32, n int)
+
+//go:noescape
+func sub32SSE2(dst, a, b *float32, n int)
+
+//go:noescape
+func mul32SSE2(dst, a, b *float32, n int)
+
+//go:noescape
+func div32SSE2(dst, a, b *float32, n int)
diff --git a/internal/kernels/bin32_arm64.go b/internal/kernels/bin32_arm64.go
new file mode 100644
index 0000000..e26ac46
--- /dev/null
+++ b/internal/kernels/bin32_arm64.go
@@ -0,0 +1,20 @@
+package kernels
+
+// bin32 runs the NEON kernel for op (sum_arm64.s: VFADD, VFSUB, VFMUL, VFDIV
+// on .S4, 16 floats per iteration). NEON is the arm64 baseline. dst, a and b
+// have the same length, at least 1.
+func bin32(op int, dst, a, b []float32) {
+	run32([...]bin32Func{add32NEON, sub32NEON, mul32NEON, div32NEON}[op], dst, a, b)
+}
+
+//go:noescape
+func add32NEON(dst, a, b *float32, n int)
+
+//go:noescape
+func sub32NEON(dst, a, b *float32, n int)
+
+//go:noescape
+func mul32NEON(dst, a, b *float32, n int)
+
+//go:noescape
+func div32NEON(dst, a, b *float32, n int)
diff --git a/internal/kernels/bin32_generic.go b/internal/kernels/bin32_generic.go
new file mode 100644
index 0000000..bfb0400
--- /dev/null
+++ b/internal/kernels/bin32_generic.go
@@ -0,0 +1,7 @@
+//go:build !amd64 && !arm64 && !loong64 && !ppc64le && !riscv64 && !s390x
+
+package kernels
+
+// bin32 is the scalar loop on the architectures without float32 elementwise
+// kernels. dst, a and b have the same length, at least 1.
+func bin32(op int, dst, a, b []float32) { bin32Go(op, dst, a, b) }
diff --git a/internal/kernels/bin32_loong64.go b/internal/kernels/bin32_loong64.go
new file mode 100644
index 0000000..144dd2c
--- /dev/null
+++ b/internal/kernels/bin32_loong64.go
@@ -0,0 +1,25 @@
+package kernels
+
+// bin32 runs the LSX kernel for op (sum_loong64.s: vfadd.s, vfsub.s,
+// vfmul.s, vfdiv.s over four 128-bit registers) when the CPU has LASX, which
+// every LASX CPU's LSX comes with; the scalar loop otherwise. dst, a and b
+// have the same length, at least 1.
+func bin32(op int, dst, a, b []float32) {
+	if !useLASX {
+		bin32Go(op, dst, a, b)
+		return
+	}
+	run32([...]bin32Func{add32LSX, sub32LSX, mul32LSX, div32LSX}[op], dst, a, b)
+}
+
+//go:noescape
+func add32LSX(dst, a, b *float32, n int)
+
+//go:noescape
+func sub32LSX(dst, a, b *float32, n int)
+
+//go:noescape
+func mul32LSX(dst, a, b *float32, n int)
+
+//go:noescape
+func div32LSX(dst, a, b *float32, n int)
diff --git a/internal/kernels/bin32_ppc64le.go b/internal/kernels/bin32_ppc64le.go
new file mode 100644
index 0000000..2b080b1
--- /dev/null
+++ b/internal/kernels/bin32_ppc64le.go
@@ -0,0 +1,20 @@
+package kernels
+
+// bin32 runs the VSX kernel for op (sum_ppc64le.s: xvaddsp, xvsubsp,
+// xvmulsp, xvdivsp, 16 floats per iteration). VSX is the ppc64le baseline.
+// dst, a and b have the same length, at least 1.
+func bin32(op int, dst, a, b []float32) {
+	run32([...]bin32Func{add32VSX, sub32VSX, mul32VSX, div32VSX}[op], dst, a, b)
+}
+
+//go:noescape
+func add32VSX(dst, a, b *float32, n int)
+
+//go:noescape
+func sub32VSX(dst, a, b *float32, n int)
+
+//go:noescape
+func mul32VSX(dst, a, b *float32, n int)
+
+//go:noescape
+func div32VSX(dst, a, b *float32, n int)
diff --git a/internal/kernels/bin32_riscv64.go b/internal/kernels/bin32_riscv64.go
new file mode 100644
index 0000000..3cf60f3
--- /dev/null
+++ b/internal/kernels/bin32_riscv64.go
@@ -0,0 +1,24 @@
+package kernels
+
+// bin32 runs the RVV kernel for op (sum_riscv64.s: vfadd.vv, vfsub.vv,
+// vfmul.vv, vfdiv.vv on 32-bit elements, strip-mined) when the CPU reports V;
+// the scalar loop otherwise. dst, a and b have the same length, at least 1.
+func bin32(op int, dst, a, b []float32) {
+	if !useRVV {
+		bin32Go(op, dst, a, b)
+		return
+	}
+	run32([...]bin32Func{add32RVV, sub32RVV, mul32RVV, div32RVV}[op], dst, a, b)
+}
+
+//go:noescape
+func add32RVV(dst, a, b *float32, n int)
+
+//go:noescape
+func sub32RVV(dst, a, b *float32, n int)
+
+//go:noescape
+func mul32RVV(dst, a, b *float32, n int)
+
+//go:noescape
+func div32RVV(dst, a, b *float32, n int)
diff --git a/internal/kernels/bin32_s390x.go b/internal/kernels/bin32_s390x.go
new file mode 100644
index 0000000..42727e1
--- /dev/null
+++ b/internal/kernels/bin32_s390x.go
@@ -0,0 +1,25 @@
+package kernels
+
+// bin32 runs the vector kernel for op (sum_s390x.s: vfasb, vfssb, vfmsb,
+// vfdsb, 16 floats per iteration) on a z14 or later, whose vector-enhancements
+// facility 1 brings them (the useVXE gate of max/min); the scalar loop
+// otherwise. dst, a and b have the same length, at least 1.
+func bin32(op int, dst, a, b []float32) {
+	if !useVXE {
+		bin32Go(op, dst, a, b)
+		return
+	}
+	run32([...]bin32Func{add32VXE, sub32VXE, mul32VXE, div32VXE}[op], dst, a, b)
+}
+
+//go:noescape
+func add32VXE(dst, a, b *float32, n int)
+
+//go:noescape
+func sub32VXE(dst, a, b *float32, n int)
+
+//go:noescape
+func mul32VXE(dst, a, b *float32, n int)
+
+//go:noescape
+func div32VXE(dst, a, b *float32, n int)
diff --git a/internal/kernels/bin32_test.go b/internal/kernels/bin32_test.go
new file mode 100644
index 0000000..e777551
--- /dev/null
+++ b/internal/kernels/bin32_test.go
@@ -0,0 +1,83 @@
+package kernels
+
+import (
+	"math"
+	"testing"
+)
+
+// special32 are the values a lane-wise kernel must treat as the scalar loop
+// does: signed zeros, infinities, NaN, subnormals and the extremes.
+var special32 = []float32{
+	0, float32(math.Copysign(0, -1)), 1, -1, 0.1, -3.5,
+	float32(math.Inf(1)), float32(math.Inf(-1)), float32(math.NaN()),
+	math.SmallestNonzeroFloat32, -math.SmallestNonzeroFloat32,
+	math.MaxFloat32, -math.MaxFloat32, 1e-38, 3e38, 7,
+}
+
+func operands32(n, seed int) []float32 {
+	a := make([]float32, n)
+	for i := range a {
+		a[i] = special32[(i*7+seed)%len(special32)]
+		if (i+seed)%5 == 0 {
+			a[i] = float32(i+seed) / 3
+		}
+	}
+	return a
+}
+
+func sameBits32(x, y float32) bool {
+	if x != x && y != y {
+		return true // NaN payloads are not part of the IEEE result
+	}
+	return math.Float32bits(x) == math.Float32bits(y)
+}
+
+// TestBin32MatchesScalar checks every op against the scalar oracle, bit for
+// bit, at every length around the 16-float block and the tail.
+func TestBin32MatchesScalar(t *testing.T) {
+	for op := Add32; op <= Div32; op++ {
+		for n := 0; n <= 70; n++ {
+			a, b := operands32(n, 1), operands32(n, 4)
+			got, want := make([]float32, n), make([]float32, n)
+			Bin32(op, got, a, b)
+			bin32Go(op, want, a, b)
+			for i := range got {
+				if !sameBits32(got[i], want[i]) {
+					t.Fatalf("op %d n=%d [%d]: %v OP %v = %v, scalar %v", op, n, i, a[i], b[i], got[i], want[i])
+				}
+			}
+		}
+	}
+}
+
+// TestBin32InPlace runs dst aliasing a, the way an in-place update calls it.
+func TestBin32InPlace(t *testing.T) {
+	for op := Add32; op <= Div32; op++ {
+		a, b := operands32(53, 2), operands32(53, 3)
+		want := make([]float32, len(a))
+		bin32Go(op, want, a, b)
+		Bin32(op, a, a, b)
+		for i := range a {
+			if !sameBits32(a[i], want[i]) {
+				t.Fatalf("op %d [%d]: in place %v, scalar %v", op, i, a[i], want[i])
+			}
+		}
+	}
+}
+
+func TestBin32RefusesBadInput(t *testing.T) {
+	for _, f := range []func(){
+		func() { Bin32(Add32, make([]float32, 4), make([]float32, 3), make([]float32, 4)) },
+		func() { Bin32(Div32+1, make([]float32, 4), make([]float32, 4), make([]float32, 4)) },
+		func() { bin32Go(-1, make([]float32, 1), make([]float32, 1), make([]float32, 1)) },
+	} {
+		func() {
+			defer func() {
+				if recover() == nil {
+					t.Error("bad input was accepted")
+				}
+			}()
+			f()
+		}()
+	}
+}
diff --git a/internal/kernels/fence_test.go b/internal/kernels/fence_test.go
index 3d23895..1587dc2 100644
--- a/internal/kernels/fence_test.go
+++ b/internal/kernels/fence_test.go
@@ -80,6 +80,22 @@ func TestBlocks8StaysInsideItsOperands(t *testing.T) {
 	}
 }
 
+// TestBin32StaysInsideItsOperands fences dst, a and b of the float32
+// elementwise kernels at every length around the block and the tail.
+func TestBin32StaysInsideItsOperands(t *testing.T) {
+	for n := 1; n <= 40; n++ {
+		for _, a := range fences32(t, n) {
+			for _, b := range fences32(t, n) {
+				for _, dst := range fences32(t, n) {
+					for op := Add32; op <= Div32; op++ {
+						Bin32(op, dst, a, b)
+					}
+				}
+			}
+		}
+	}
+}
+
 func TestKernelsStayInsideTheirSlices(t *testing.T) {
 	for n := 1; n <= 70; n++ {
 		for _, a := range fences(t, n) {
diff --git a/internal/kernels/sum_amd64.s b/internal/kernels/sum_amd64.s
index 6026831..dba30ca 100644
--- a/internal/kernels/sum_amd64.s
+++ b/internal/kernels/sum_amd64.s
@@ -722,3 +722,179 @@ TEXT ·blocks8SSE2(SB), NOSPLIT, $0-24
 	MOVUPS X1, 16(DI)
 	RET
 
+TEXT ·add32SSE2(SB), NOSPLIT, $0-32
+	MOVQ dst+0(FP), DI
+	MOVQ a+8(FP), SI
+	MOVQ b+16(FP), DX
+	MOVQ n+24(FP), CX
+	fblock:
+	CMPQ CX, $16
+	JL ftail
+	MOVUPS 0(SI), X0
+	MOVUPS 16(SI), X1
+	MOVUPS 32(SI), X2
+	MOVUPS 48(SI), X3
+	MOVUPS 0(DX), X4
+	MOVUPS 16(DX), X5
+	MOVUPS 32(DX), X6
+	MOVUPS 48(DX), X7
+	ADDPS X4, X0
+	ADDPS X5, X1
+	ADDPS X6, X2
+	ADDPS X7, X3
+	MOVUPS X0, 0(DI)
+	MOVUPS X1, 16(DI)
+	MOVUPS X2, 32(DI)
+	MOVUPS X3, 48(DI)
+	ADDQ $64, SI
+	ADDQ $64, DX
+	ADDQ $64, DI
+	SUBQ $16, CX
+	JMP fblock
+	ftail:
+	TESTQ CX, CX
+	JZ fdone
+	MOVSS (SI), X0
+	MOVSS (DX), X1
+	ADDSS X1, X0
+	MOVSS X0, (DI)
+	ADDQ $4, SI
+	ADDQ $4, DX
+	ADDQ $4, DI
+	DECQ CX
+	JMP ftail
+	fdone:
+	RET
+
+TEXT ·sub32SSE2(SB), NOSPLIT, $0-32
+	MOVQ dst+0(FP), DI
+	MOVQ a+8(FP), SI
+	MOVQ b+16(FP), DX
+	MOVQ n+24(FP), CX
+	fblock:
+	CMPQ CX, $16
+	JL ftail
+	MOVUPS 0(SI), X0
+	MOVUPS 16(SI), X1
+	MOVUPS 32(SI), X2
+	MOVUPS 48(SI), X3
+	MOVUPS 0(DX), X4
+	MOVUPS 16(DX), X5
+	MOVUPS 32(DX), X6
+	MOVUPS 48(DX), X7
+	SUBPS X4, X0
+	SUBPS X5, X1
+	SUBPS X6, X2
+	SUBPS X7, X3
+	MOVUPS X0, 0(DI)
+	MOVUPS X1, 16(DI)
+	MOVUPS X2, 32(DI)
+	MOVUPS X3, 48(DI)
+	ADDQ $64, SI
+	ADDQ $64, DX
+	ADDQ $64, DI
+	SUBQ $16, CX
+	JMP fblock
+	ftail:
+	TESTQ CX, CX
+	JZ fdone
+	MOVSS (SI), X0
+	MOVSS (DX), X1
+	SUBSS X1, X0
+	MOVSS X0, (DI)
+	ADDQ $4, SI
+	ADDQ $4, DX
+	ADDQ $4, DI
+	DECQ CX
+	JMP ftail
+	fdone:
+	RET
+
+TEXT ·mul32SSE2(SB), NOSPLIT, $0-32
+	MOVQ dst+0(FP), DI
+	MOVQ a+8(FP), SI
+	MOVQ b+16(FP), DX
+	MOVQ n+24(FP), CX
+	fblock:
+	CMPQ CX, $16
+	JL ftail
+	MOVUPS 0(SI), X0
+	MOVUPS 16(SI), X1
+	MOVUPS 32(SI), X2
+	MOVUPS 48(SI), X3
+	MOVUPS 0(DX), X4
+	MOVUPS 16(DX), X5
+	MOVUPS 32(DX), X6
+	MOVUPS 48(DX), X7
+	MULPS X4, X0
+	MULPS X5, X1
+	MULPS X6, X2
+	MULPS X7, X3
+	MOVUPS X0, 0(DI)
+	MOVUPS X1, 16(DI)
+	MOVUPS X2, 32(DI)
+	MOVUPS X3, 48(DI)
+	ADDQ $64, SI
+	ADDQ $64, DX
+	ADDQ $64, DI
+	SUBQ $16, CX
+	JMP fblock
+	ftail:
+	TESTQ CX, CX
+	JZ fdone
+	MOVSS (SI), X0
+	MOVSS (DX), X1
+	MULSS X1, X0
+	MOVSS X0, (DI)
+	ADDQ $4, SI
+	ADDQ $4, DX
+	ADDQ $4, DI
+	DECQ CX
+	JMP ftail
+	fdone:
+	RET
+
+TEXT ·div32SSE2(SB), NOSPLIT, $0-32
+	MOVQ dst+0(FP), DI
+	MOVQ a+8(FP), SI
+	MOVQ b+16(FP), DX
+	MOVQ n+24(FP), CX
+	fblock:
+	CMPQ CX, $16
+	JL ftail
+	MOVUPS 0(SI), X0
+	MOVUPS 16(SI), X1
+	MOVUPS 32(SI), X2
+	MOVUPS 48(SI), X3
+	MOVUPS 0(DX), X4
+	MOVUPS 16(DX), X5
+	MOVUPS 32(DX), X6
+	MOVUPS 48(DX), X7
+	DIVPS X4, X0
+	DIVPS X5, X1
+	DIVPS X6, X2
+	DIVPS X7, X3
+	MOVUPS X0, 0(DI)
+	MOVUPS X1, 16(DI)
+	MOVUPS X2, 32(DI)
+	MOVUPS X3, 48(DI)
+	ADDQ $64, SI
+	ADDQ $64, DX
+	ADDQ $64, DI
+	SUBQ $16, CX
+	JMP fblock
+	ftail:
+	TESTQ CX, CX
+	JZ fdone
+	MOVSS (SI), X0
+	MOVSS (DX), X1
+	DIVSS X1, X0
+	MOVSS X0, (DI)
+	ADDQ $4, SI
+	ADDQ $4, DX
+	ADDQ $4, DI
+	DECQ CX
+	JMP ftail
+	fdone:
+	RET
+
diff --git a/internal/kernels/sum_arm64.s b/internal/kernels/sum_arm64.s
index dc03f0e..340c8f3 100644
--- a/internal/kernels/sum_arm64.s
+++ b/internal/kernels/sum_arm64.s
@@ -234,6 +234,118 @@ TEXT ·divNEON(SB), NOSPLIT, $0-32
 	done:
 	RET
 
+TEXT ·add32NEON(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R0
+	MOVD a+8(FP), R1
+	MOVD b+16(FP), R2
+	MOVD n+24(FP), R3
+	block:
+	CMP $16, R3
+	BLT tail
+	VLD1.P 64(R1), [V0.S4, V1.S4, V2.S4, V3.S4]
+	VLD1.P 64(R2), [V4.S4, V5.S4, V6.S4, V7.S4]
+	VFADD V4.S4, V0.S4, V0.S4
+	VFADD V5.S4, V1.S4, V1.S4
+	VFADD V6.S4, V2.S4, V2.S4
+	VFADD V7.S4, V3.S4, V3.S4
+	VST1.P [V0.S4, V1.S4, V2.S4, V3.S4], 64(R0)
+	SUB $16, R3
+	B block
+	tail:
+	CBZ R3, done
+	FMOVS.P 4(R1), F0
+	FMOVS.P 4(R2), F1
+	FADDS F1, F0, F0
+	FMOVS.P F0, 4(R0)
+	SUB $1, R3
+	B tail
+	done:
+	RET
+
+TEXT ·sub32NEON(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R0
+	MOVD a+8(FP), R1
+	MOVD b+16(FP), R2
+	MOVD n+24(FP), R3
+	block:
+	CMP $16, R3
+	BLT tail
+	VLD1.P 64(R1), [V0.S4, V1.S4, V2.S4, V3.S4]
+	VLD1.P 64(R2), [V4.S4, V5.S4, V6.S4, V7.S4]
+	VFSUB V4.S4, V0.S4, V0.S4
+	VFSUB V5.S4, V1.S4, V1.S4
+	VFSUB V6.S4, V2.S4, V2.S4
+	VFSUB V7.S4, V3.S4, V3.S4
+	VST1.P [V0.S4, V1.S4, V2.S4, V3.S4], 64(R0)
+	SUB $16, R3
+	B block
+	tail:
+	CBZ R3, done
+	FMOVS.P 4(R1), F0
+	FMOVS.P 4(R2), F1
+	FSUBS F1, F0, F0
+	FMOVS.P F0, 4(R0)
+	SUB $1, R3
+	B tail
+	done:
+	RET
+
+TEXT ·mul32NEON(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R0
+	MOVD a+8(FP), R1
+	MOVD b+16(FP), R2
+	MOVD n+24(FP), R3
+	block:
+	CMP $16, R3
+	BLT tail
+	VLD1.P 64(R1), [V0.S4, V1.S4, V2.S4, V3.S4]
+	VLD1.P 64(R2), [V4.S4, V5.S4, V6.S4, V7.S4]
+	VFMUL V4.S4, V0.S4, V0.S4
+	VFMUL V5.S4, V1.S4, V1.S4
+	VFMUL V6.S4, V2.S4, V2.S4
+	VFMUL V7.S4, V3.S4, V3.S4
+	VST1.P [V0.S4, V1.S4, V2.S4, V3.S4], 64(R0)
+	SUB $16, R3
+	B block
+	tail:
+	CBZ R3, done
+	FMOVS.P 4(R1), F0
+	FMOVS.P 4(R2), F1
+	FMULS F1, F0, F0
+	FMOVS.P F0, 4(R0)
+	SUB $1, R3
+	B tail
+	done:
+	RET
+
+TEXT ·div32NEON(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R0
+	MOVD a+8(FP), R1
+	MOVD b+16(FP), R2
+	MOVD n+24(FP), R3
+	block:
+	CMP $16, R3
+	BLT tail
+	VLD1.P 64(R1), [V0.S4, V1.S4, V2.S4, V3.S4]
+	VLD1.P 64(R2), [V4.S4, V5.S4, V6.S4, V7.S4]
+	VFDIV V4.S4, V0.S4, V0.S4
+	VFDIV V5.S4, V1.S4, V1.S4
+	VFDIV V6.S4, V2.S4, V2.S4
+	VFDIV V7.S4, V3.S4, V3.S4
+	VST1.P [V0.S4, V1.S4, V2.S4, V3.S4], 64(R0)
+	SUB $16, R3
+	B block
+	tail:
+	CBZ R3, done
+	FMOVS.P 4(R1), F0
+	FMOVS.P 4(R2), F1
+	FDIVS F1, F0, F0
+	FMOVS.P F0, 4(R0)
+	SUB $1, R3
+	B tail
+	done:
+	RET
+
 TEXT ·maxNEON(SB), NOSPLIT, $0-24
 	MOVD a+0(FP), R0
 	MOVD n+8(FP), R1
diff --git a/internal/kernels/sum_loong64.s b/internal/kernels/sum_loong64.s
index 4a88d40..9729476 100644
--- a/internal/kernels/sum_loong64.s
+++ b/internal/kernels/sum_loong64.s
@@ -184,6 +184,49 @@ addtail:
 adddone:
 	RET
 
+TEXT ·add32LSX(SB), NOSPLIT, $0-32
+	MOVV dst+0(FP), R4
+	MOVV a+8(FP), R6
+	MOVV b+16(FP), R7
+	MOVV n+24(FP), R5
+floop:
+	MOVV $16, R12
+	BLT R5, R12, ftail
+	VMOVQ 0(R6), V0
+	VMOVQ 0(R7), V4
+	VMOVQ 16(R6), V1
+	VMOVQ 16(R7), V5
+	VMOVQ 32(R6), V2
+	VMOVQ 32(R7), V6
+	VMOVQ 48(R6), V3
+	VMOVQ 48(R7), V7
+	VADDF V4, V0, V8
+	VADDF V5, V1, V9
+	VADDF V6, V2, V10
+	VADDF V7, V3, V11
+	VMOVQ V8, 0(R4)
+	VMOVQ V9, 16(R4)
+	VMOVQ V10, 32(R4)
+	VMOVQ V11, 48(R4)
+	ADDV $64, R4
+	ADDV $64, R6
+	ADDV $64, R7
+	ADDV $-16, R5
+	JMP floop
+ftail:
+	BEQ R5, fdone
+	MOVF 0(R6), F1
+	MOVF 0(R7), F2
+	ADDF F2, F1
+	MOVF F1, 0(R4)
+	ADDV $4, R4
+	ADDV $4, R6
+	ADDV $4, R7
+	ADDV $-1, R5
+	JMP ftail
+fdone:
+	RET
+
 TEXT ·subLASX(SB), NOSPLIT, $0-32
 	MOVV dst+0(FP), R4
 	MOVV a+8(FP), R6
@@ -227,6 +270,49 @@ subtail:
 subdone:
 	RET
 
+TEXT ·sub32LSX(SB), NOSPLIT, $0-32
+	MOVV dst+0(FP), R4
+	MOVV a+8(FP), R6
+	MOVV b+16(FP), R7
+	MOVV n+24(FP), R5
+floop:
+	MOVV $16, R12
+	BLT R5, R12, ftail
+	VMOVQ 0(R6), V0
+	VMOVQ 0(R7), V4
+	VMOVQ 16(R6), V1
+	VMOVQ 16(R7), V5
+	VMOVQ 32(R6), V2
+	VMOVQ 32(R7), V6
+	VMOVQ 48(R6), V3
+	VMOVQ 48(R7), V7
+	VSUBF V4, V0, V8
+	VSUBF V5, V1, V9
+	VSUBF V6, V2, V10
+	VSUBF V7, V3, V11
+	VMOVQ V8, 0(R4)
+	VMOVQ V9, 16(R4)
+	VMOVQ V10, 32(R4)
+	VMOVQ V11, 48(R4)
+	ADDV $64, R4
+	ADDV $64, R6
+	ADDV $64, R7
+	ADDV $-16, R5
+	JMP floop
+ftail:
+	BEQ R5, fdone
+	MOVF 0(R6), F1
+	MOVF 0(R7), F2
+	SUBF F2, F1
+	MOVF F1, 0(R4)
+	ADDV $4, R4
+	ADDV $4, R6
+	ADDV $4, R7
+	ADDV $-1, R5
+	JMP ftail
+fdone:
+	RET
+
 TEXT ·mulLASX(SB), NOSPLIT, $0-32
 	MOVV dst+0(FP), R4
 	MOVV a+8(FP), R6
@@ -270,6 +356,49 @@ multail:
 muldone:
 	RET
 
+TEXT ·mul32LSX(SB), NOSPLIT, $0-32
+	MOVV dst+0(FP), R4
+	MOVV a+8(FP), R6
+	MOVV b+16(FP), R7
+	MOVV n+24(FP), R5
+floop:
+	MOVV $16, R12
+	BLT R5, R12, ftail
+	VMOVQ 0(R6), V0
+	VMOVQ 0(R7), V4
+	VMOVQ 16(R6), V1
+	VMOVQ 16(R7), V5
+	VMOVQ 32(R6), V2
+	VMOVQ 32(R7), V6
+	VMOVQ 48(R6), V3
+	VMOVQ 48(R7), V7
+	VMULF V4, V0, V8
+	VMULF V5, V1, V9
+	VMULF V6, V2, V10
+	VMULF V7, V3, V11
+	VMOVQ V8, 0(R4)
+	VMOVQ V9, 16(R4)
+	VMOVQ V10, 32(R4)
+	VMOVQ V11, 48(R4)
+	ADDV $64, R4
+	ADDV $64, R6
+	ADDV $64, R7
+	ADDV $-16, R5
+	JMP floop
+ftail:
+	BEQ R5, fdone
+	MOVF 0(R6), F1
+	MOVF 0(R7), F2
+	MULF F2, F1
+	MOVF F1, 0(R4)
+	ADDV $4, R4
+	ADDV $4, R6
+	ADDV $4, R7
+	ADDV $-1, R5
+	JMP ftail
+fdone:
+	RET
+
 TEXT ·divLASX(SB), NOSPLIT, $0-32
 	MOVV dst+0(FP), R4
 	MOVV a+8(FP), R6
@@ -313,6 +442,49 @@ divtail:
 divdone:
 	RET
 
+TEXT ·div32LSX(SB), NOSPLIT, $0-32
+	MOVV dst+0(FP), R4
+	MOVV a+8(FP), R6
+	MOVV b+16(FP), R7
+	MOVV n+24(FP), R5
+floop:
+	MOVV $16, R12
+	BLT R5, R12, ftail
+	VMOVQ 0(R6), V0
+	VMOVQ 0(R7), V4
+	VMOVQ 16(R6), V1
+	VMOVQ 16(R7), V5
+	VMOVQ 32(R6), V2
+	VMOVQ 32(R7), V6
+	VMOVQ 48(R6), V3
+	VMOVQ 48(R7), V7
+	VDIVF V4, V0, V8
+	VDIVF V5, V1, V9
+	VDIVF V6, V2, V10
+	VDIVF V7, V3, V11
+	VMOVQ V8, 0(R4)
+	VMOVQ V9, 16(R4)
+	VMOVQ V10, 32(R4)
+	VMOVQ V11, 48(R4)
+	ADDV $64, R4
+	ADDV $64, R6
+	ADDV $64, R7
+	ADDV $-16, R5
+	JMP floop
+ftail:
+	BEQ R5, fdone
+	MOVF 0(R6), F1
+	MOVF 0(R7), F2
+	DIVF F2, F1
+	MOVF F1, 0(R4)
+	ADDV $4, R4
+	ADDV $4, R6
+	ADDV $4, R7
+	ADDV $-1, R5
+	JMP ftail
+fdone:
+	RET
+
 TEXT ·gemmMicro8x8LASX(SB), NOSPLIT, $0-40
 	MOVV kc+0(FP), R5
 	MOVV pa+8(FP), R6
diff --git a/internal/kernels/sum_ppc64le.s b/internal/kernels/sum_ppc64le.s
index 70ef94e..bbc8dd1 100644
--- a/internal/kernels/sum_ppc64le.s
+++ b/internal/kernels/sum_ppc64le.s
@@ -183,6 +183,53 @@ addtail:
 adddone:
 	RET
 
+TEXT ·add32VSX(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R3
+	MOVD a+8(FP), R5
+	MOVD b+16(FP), R6
+	MOVD n+24(FP), R4
+	MOVD $16, R14
+	MOVD $32, R15
+	MOVD $48, R16
+floop:
+	CMP R4, $16
+	BLT ftail
+	LXVW4X (R5)(R0), VS32
+	LXVW4X (R6)(R0), VS36
+	LXVW4X (R5)(R14), VS33
+	LXVW4X (R6)(R14), VS37
+	LXVW4X (R5)(R15), VS34
+	LXVW4X (R6)(R15), VS38
+	LXVW4X (R5)(R16), VS35
+	LXVW4X (R6)(R16), VS39
+	WORD $0xf1002207 // XVADDSP VS32, VS36, VS40
+	WORD $0xf1212a07 // XVADDSP VS33, VS37, VS41
+	WORD $0xf1423207 // XVADDSP VS34, VS38, VS42
+	WORD $0xf1633a07 // XVADDSP VS35, VS39, VS43
+	STXVW4X VS40, (R3)(R0)
+	STXVW4X VS41, (R3)(R14)
+	STXVW4X VS42, (R3)(R15)
+	STXVW4X VS43, (R3)(R16)
+	ADD $64, R3
+	ADD $64, R5
+	ADD $64, R6
+	ADD $-16, R4
+	BR floop
+ftail:
+	CMP R4, $0
+	BEQ fdone
+	FMOVS 0(R5), F1
+	FMOVS 0(R6), F2
+	FADDS F2, F1
+	FMOVS F1, 0(R3)
+	ADD $4, R3
+	ADD $4, R5
+	ADD $4, R6
+	ADD $-1, R4
+	BR ftail
+fdone:
+	RET
+
 TEXT ·subVSX(SB), NOSPLIT, $0-32
 	MOVD dst+0(FP), R3
 	MOVD a+8(FP), R5
@@ -230,6 +277,53 @@ subtail:
 subdone:
 	RET
 
+TEXT ·sub32VSX(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R3
+	MOVD a+8(FP), R5
+	MOVD b+16(FP), R6
+	MOVD n+24(FP), R4
+	MOVD $16, R14
+	MOVD $32, R15
+	MOVD $48, R16
+floop:
+	CMP R4, $16
+	BLT ftail
+	LXVW4X (R5)(R0), VS32
+	LXVW4X (R6)(R0), VS36
+	LXVW4X (R5)(R14), VS33
+	LXVW4X (R6)(R14), VS37
+	LXVW4X (R5)(R15), VS34
+	LXVW4X (R6)(R15), VS38
+	LXVW4X (R5)(R16), VS35
+	LXVW4X (R6)(R16), VS39
+	WORD $0xf1002247 // XVSUBSP VS32, VS36, VS40
+	WORD $0xf1212a47 // XVSUBSP VS33, VS37, VS41
+	WORD $0xf1423247 // XVSUBSP VS34, VS38, VS42
+	WORD $0xf1633a47 // XVSUBSP VS35, VS39, VS43
+	STXVW4X VS40, (R3)(R0)
+	STXVW4X VS41, (R3)(R14)
+	STXVW4X VS42, (R3)(R15)
+	STXVW4X VS43, (R3)(R16)
+	ADD $64, R3
+	ADD $64, R5
+	ADD $64, R6
+	ADD $-16, R4
+	BR floop
+ftail:
+	CMP R4, $0
+	BEQ fdone
+	FMOVS 0(R5), F1
+	FMOVS 0(R6), F2
+	FSUBS F2, F1
+	FMOVS F1, 0(R3)
+	ADD $4, R3
+	ADD $4, R5
+	ADD $4, R6
+	ADD $-1, R4
+	BR ftail
+fdone:
+	RET
+
 TEXT ·mulVSX(SB), NOSPLIT, $0-32
 	MOVD dst+0(FP), R3
 	MOVD a+8(FP), R5
@@ -277,6 +371,53 @@ multail:
 muldone:
 	RET
 
+TEXT ·mul32VSX(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R3
+	MOVD a+8(FP), R5
+	MOVD b+16(FP), R6
+	MOVD n+24(FP), R4
+	MOVD $16, R14
+	MOVD $32, R15
+	MOVD $48, R16
+floop:
+	CMP R4, $16
+	BLT ftail
+	LXVW4X (R5)(R0), VS32
+	LXVW4X (R6)(R0), VS36
+	LXVW4X (R5)(R14), VS33
+	LXVW4X (R6)(R14), VS37
+	LXVW4X (R5)(R15), VS34
+	LXVW4X (R6)(R15), VS38
+	LXVW4X (R5)(R16), VS35
+	LXVW4X (R6)(R16), VS39
+	WORD $0xf1002287 // XVMULSP VS32, VS36, VS40
+	WORD $0xf1212a87 // XVMULSP VS33, VS37, VS41
+	WORD $0xf1423287 // XVMULSP VS34, VS38, VS42
+	WORD $0xf1633a87 // XVMULSP VS35, VS39, VS43
+	STXVW4X VS40, (R3)(R0)
+	STXVW4X VS41, (R3)(R14)
+	STXVW4X VS42, (R3)(R15)
+	STXVW4X VS43, (R3)(R16)
+	ADD $64, R3
+	ADD $64, R5
+	ADD $64, R6
+	ADD $-16, R4
+	BR floop
+ftail:
+	CMP R4, $0
+	BEQ fdone
+	FMOVS 0(R5), F1
+	FMOVS 0(R6), F2
+	FMULS F2, F1
+	FMOVS F1, 0(R3)
+	ADD $4, R3
+	ADD $4, R5
+	ADD $4, R6
+	ADD $-1, R4
+	BR ftail
+fdone:
+	RET
+
 TEXT ·divVSX(SB), NOSPLIT, $0-32
 	MOVD dst+0(FP), R3
 	MOVD a+8(FP), R5
@@ -324,6 +465,53 @@ divtail:
 divdone:
 	RET
 
+TEXT ·div32VSX(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R3
+	MOVD a+8(FP), R5
+	MOVD b+16(FP), R6
+	MOVD n+24(FP), R4
+	MOVD $16, R14
+	MOVD $32, R15
+	MOVD $48, R16
+floop:
+	CMP R4, $16
+	BLT ftail
+	LXVW4X (R5)(R0), VS32
+	LXVW4X (R6)(R0), VS36
+	LXVW4X (R5)(R14), VS33
+	LXVW4X (R6)(R14), VS37
+	LXVW4X (R5)(R15), VS34
+	LXVW4X (R6)(R15), VS38
+	LXVW4X (R5)(R16), VS35
+	LXVW4X (R6)(R16), VS39
+	WORD $0xf10022c7 // XVDIVSP VS32, VS36, VS40
+	WORD $0xf1212ac7 // XVDIVSP VS33, VS37, VS41
+	WORD $0xf14232c7 // XVDIVSP VS34, VS38, VS42
+	WORD $0xf1633ac7 // XVDIVSP VS35, VS39, VS43
+	STXVW4X VS40, (R3)(R0)
+	STXVW4X VS41, (R3)(R14)
+	STXVW4X VS42, (R3)(R15)
+	STXVW4X VS43, (R3)(R16)
+	ADD $64, R3
+	ADD $64, R5
+	ADD $64, R6
+	ADD $-16, R4
+	BR floop
+ftail:
+	CMP R4, $0
+	BEQ fdone
+	FMOVS 0(R5), F1
+	FMOVS 0(R6), F2
+	FDIVS F2, F1
+	FMOVS F1, 0(R3)
+	ADD $4, R3
+	ADD $4, R5
+	ADD $4, R6
+	ADD $-1, R4
+	BR ftail
+fdone:
+	RET
+
 TEXT ·gemmMicro8x8VSX(SB), NOSPLIT, $0-40
 	MOVD kc+0(FP), R3
 	MOVD pa+8(FP), R4
diff --git a/internal/kernels/sum_riscv64.s b/internal/kernels/sum_riscv64.s
index 07fa8d2..461d1fd 100644
--- a/internal/kernels/sum_riscv64.s
+++ b/internal/kernels/sum_riscv64.s
@@ -84,6 +84,25 @@ TEXT ·addRVV(SB), NOSPLIT, $0-32
 	BNEZ X13, loop
 	RET
 
+TEXT ·add32RVV(SB), NOSPLIT, $0-32
+	MOV dst+0(FP), X10
+	MOV a+8(FP), X11
+	MOV b+16(FP), X12
+	MOV n+24(FP), X13
+	loop:
+	VSETVLI X13, E32, M8, TA, MA, X14
+	VLE32V (X11), V0
+	VLE32V (X12), V8
+	VFADDVV V8, V0, V16
+	VSE32V V16, (X10)
+	SLLI $2, X14, X17
+	ADD X17, X10, X10
+	ADD X17, X11, X11
+	ADD X17, X12, X12
+	SUB X14, X13, X13
+	BNEZ X13, loop
+	RET
+
 TEXT ·subRVV(SB), NOSPLIT, $0-32
 	MOV dst+0(FP), X10
 	MOV a+8(FP), X11
@@ -103,6 +122,25 @@ TEXT ·subRVV(SB), NOSPLIT, $0-32
 	BNEZ X13, loop
 	RET
 
+TEXT ·sub32RVV(SB), NOSPLIT, $0-32
+	MOV dst+0(FP), X10
+	MOV a+8(FP), X11
+	MOV b+16(FP), X12
+	MOV n+24(FP), X13
+	loop:
+	VSETVLI X13, E32, M8, TA, MA, X14
+	VLE32V (X11), V0
+	VLE32V (X12), V8
+	VFSUBVV V8, V0, V16
+	VSE32V V16, (X10)
+	SLLI $2, X14, X17
+	ADD X17, X10, X10
+	ADD X17, X11, X11
+	ADD X17, X12, X12
+	SUB X14, X13, X13
+	BNEZ X13, loop
+	RET
+
 TEXT ·mulRVV(SB), NOSPLIT, $0-32
 	MOV dst+0(FP), X10
 	MOV a+8(FP), X11
@@ -122,6 +160,25 @@ TEXT ·mulRVV(SB), NOSPLIT, $0-32
 	BNEZ X13, loop
 	RET
 
+TEXT ·mul32RVV(SB), NOSPLIT, $0-32
+	MOV dst+0(FP), X10
+	MOV a+8(FP), X11
+	MOV b+16(FP), X12
+	MOV n+24(FP), X13
+	loop:
+	VSETVLI X13, E32, M8, TA, MA, X14
+	VLE32V (X11), V0
+	VLE32V (X12), V8
+	VFMULVV V8, V0, V16
+	VSE32V V16, (X10)
+	SLLI $2, X14, X17
+	ADD X17, X10, X10
+	ADD X17, X11, X11
+	ADD X17, X12, X12
+	SUB X14, X13, X13
+	BNEZ X13, loop
+	RET
+
 TEXT ·divRVV(SB), NOSPLIT, $0-32
 	MOV dst+0(FP), X10
 	MOV a+8(FP), X11
@@ -141,6 +198,25 @@ TEXT ·divRVV(SB), NOSPLIT, $0-32
 	BNEZ X13, loop
 	RET
 
+TEXT ·div32RVV(SB), NOSPLIT, $0-32
+	MOV dst+0(FP), X10
+	MOV a+8(FP), X11
+	MOV b+16(FP), X12
+	MOV n+24(FP), X13
+	loop:
+	VSETVLI X13, E32, M8, TA, MA, X14
+	VLE32V (X11), V0
+	VLE32V (X12), V8
+	VFDIVVV V8, V0, V16
+	VSE32V V16, (X10)
+	SLLI $2, X14, X17
+	ADD X17, X10, X10
+	ADD X17, X11, X11
+	ADD X17, X12, X12
+	SUB X14, X13, X13
+	BNEZ X13, loop
+	RET
+
 TEXT ·gemmMicro6x8RVV(SB), NOSPLIT, $0-40
 	MOV kc+0(FP), X10
 	MOV pa+8(FP), X11
diff --git a/internal/kernels/sum_s390x.s b/internal/kernels/sum_s390x.s
index 04b17d5..f7af7cb 100644
--- a/internal/kernels/sum_s390x.s
+++ b/internal/kernels/sum_s390x.s
@@ -147,6 +147,41 @@ TEXT ·addVX(SB), NOSPLIT, $0-32
 	done:
 	RET
 
+TEXT ·add32VXE(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R1
+	MOVD a+8(FP), R2
+	MOVD b+16(FP), R3
+	MOVD n+24(FP), R4
+	loop:
+	CMP R4, $16
+	BLT tail
+	VLM (R2), V0, V3
+	VLM (R3), V4, V7
+	WORD $0xe7004000; BYTE $0x20; BYTE $0xe3 // VFASB V0, V4, V0
+	WORD $0xe7115000; BYTE $0x20; BYTE $0xe3 // VFASB V1, V5, V1
+	WORD $0xe7226000; BYTE $0x20; BYTE $0xe3 // VFASB V2, V6, V2
+	WORD $0xe7337000; BYTE $0x20; BYTE $0xe3 // VFASB V3, V7, V3
+	VSTM V0, V3, (R1)
+	ADD $64, R1
+	ADD $64, R2
+	ADD $64, R3
+	SUB $16, R4
+	BR loop
+	tail:
+	CMP R4, $0
+	BEQ done
+	FMOVS (R2), F0
+	FMOVS (R3), F1
+	FADDS F1, F0
+	FMOVS F0, (R1)
+	ADD $4, R1
+	ADD $4, R2
+	ADD $4, R3
+	SUB $1, R4
+	BR tail
+	done:
+	RET
+
 TEXT ·subVX(SB), NOSPLIT, $0-32
 	MOVD dst+0(FP), R1
 	MOVD a+8(FP), R2
@@ -182,6 +217,41 @@ TEXT ·subVX(SB), NOSPLIT, $0-32
 	done:
 	RET
 
+TEXT ·sub32VXE(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R1
+	MOVD a+8(FP), R2
+	MOVD b+16(FP), R3
+	MOVD n+24(FP), R4
+	loop:
+	CMP R4, $16
+	BLT tail
+	VLM (R2), V0, V3
+	VLM (R3), V4, V7
+	WORD $0xe7004000; BYTE $0x20; BYTE $0xe2 // VFSSB V4, V0, V0
+	WORD $0xe7115000; BYTE $0x20; BYTE $0xe2 // VFSSB V5, V1, V1
+	WORD $0xe7226000; BYTE $0x20; BYTE $0xe2 // VFSSB V6, V2, V2
+	WORD $0xe7337000; BYTE $0x20; BYTE $0xe2 // VFSSB V7, V3, V3
+	VSTM V0, V3, (R1)
+	ADD $64, R1
+	ADD $64, R2
+	ADD $64, R3
+	SUB $16, R4
+	BR loop
+	tail:
+	CMP R4, $0
+	BEQ done
+	FMOVS (R2), F0
+	FMOVS (R3), F1
+	FSUBS F1, F0
+	FMOVS F0, (R1)
+	ADD $4, R1
+	ADD $4, R2
+	ADD $4, R3
+	SUB $1, R4
+	BR tail
+	done:
+	RET
+
 TEXT ·mulVX(SB), NOSPLIT, $0-32
 	MOVD dst+0(FP), R1
 	MOVD a+8(FP), R2
@@ -217,6 +287,41 @@ TEXT ·mulVX(SB), NOSPLIT, $0-32
 	done:
 	RET
 
+TEXT ·mul32VXE(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R1
+	MOVD a+8(FP), R2
+	MOVD b+16(FP), R3
+	MOVD n+24(FP), R4
+	loop:
+	CMP R4, $16
+	BLT tail
+	VLM (R2), V0, V3
+	VLM (R3), V4, V7
+	WORD $0xe7004000; BYTE $0x20; BYTE $0xe7 // VFMSB V0, V4, V0
+	WORD $0xe7115000; BYTE $0x20; BYTE $0xe7 // VFMSB V1, V5, V1
+	WORD $0xe7226000; BYTE $0x20; BYTE $0xe7 // VFMSB V2, V6, V2
+	WORD $0xe7337000; BYTE $0x20; BYTE $0xe7 // VFMSB V3, V7, V3
+	VSTM V0, V3, (R1)
+	ADD $64, R1
+	ADD $64, R2
+	ADD $64, R3
+	SUB $16, R4
+	BR loop
+	tail:
+	CMP R4, $0
+	BEQ done
+	FMOVS (R2), F0
+	FMOVS (R3), F1
+	FMULS F1, F0
+	FMOVS F0, (R1)
+	ADD $4, R1
+	ADD $4, R2
+	ADD $4, R3
+	SUB $1, R4
+	BR tail
+	done:
+	RET
+
 TEXT ·divVX(SB), NOSPLIT, $0-32
 	MOVD dst+0(FP), R1
 	MOVD a+8(FP), R2
@@ -252,6 +357,41 @@ TEXT ·divVX(SB), NOSPLIT, $0-32
 	done:
 	RET
 
+TEXT ·div32VXE(SB), NOSPLIT, $0-32
+	MOVD dst+0(FP), R1
+	MOVD a+8(FP), R2
+	MOVD b+16(FP), R3
+	MOVD n+24(FP), R4
+	loop:
+	CMP R4, $16
+	BLT tail
+	VLM (R2), V0, V3
+	VLM (R3), V4, V7
+	WORD $0xe7004000; BYTE $0x20; BYTE $0xe5 // VFDSB V4, V0, V0
+	WORD $0xe7115000; BYTE $0x20; BYTE $0xe5 // VFDSB V5, V1, V1
+	WORD $0xe7226000; BYTE $0x20; BYTE $0xe5 // VFDSB V6, V2, V2
+	WORD $0xe7337000; BYTE $0x20; BYTE $0xe5 // VFDSB V7, V3, V3
+	VSTM V0, V3, (R1)
+	ADD $64, R1
+	ADD $64, R2
+	ADD $64, R3
+	SUB $16, R4
+	BR loop
+	tail:
+	CMP R4, $0
+	BEQ done
+	FMOVS (R2), F0
+	FMOVS (R3), F1
+	FDIVS F1, F0
+	FMOVS F0, (R1)
+	ADD $4, R1
+	ADD $4, R2
+	ADD $4, R3
+	SUB $1, R4
+	BR tail
+	done:
+	RET
+
 TEXT ·gemmMicro4x8VX(SB), NOSPLIT, $0-40
 	MOVD kc+0(FP), R5
 	MOVD pa+8(FP), R2
diff --git a/typed.go b/typed.go
index 65682c5..82de375 100644
--- a/typed.go
+++ b/typed.go
@@ -40,6 +40,7 @@ const (
 // binNum applies an arithmetic operation. opDiv is only used for inexact
 // types: integer division is true division, done in a float type.
 func binNum[T numT](op binCode, d, x, y []T) {
+	x, y = x[:len(d)], y[:len(d)] // so the loops below need no bounds checks
 	switch op {
 	case opAdd:
 		for i := range d {
@@ -62,6 +63,7 @@ func binNum[T numT](op binCode, d, x, y []T) {
 
 // binInt applies an arithmetic or max/min operation to integers.
 func binInt[T intT](op binCode, d, x, y []T) {
+	x, y = x[:len(d)], y[:len(d)]
 	switch op {
 	case opMax:
 		for i := range d {
@@ -76,9 +78,13 @@ func binInt[T intT](op binCode, d, x, y []T) {
 	}
 }
 
-// binFloat applies an arithmetic or max/min operation to floats. Max and min
-// propagate NaN, as NumPy's maximum and minimum do (and Go's builtins).
-func binFloat[T floatT](op binCode, d, x, y []T) {
+// binFloat32 applies an arithmetic or max/min operation to float32 (Float64
+// has its own path). Max and min propagate NaN, as NumPy's maximum and minimum
+// do (and Go's builtins). The arithmetic runs on kernels.Bin32, whose codes
+// are opAdd..opDiv: one IEEE operation per lane, the scalar loop's result bit
+// for bit.
+func binFloat32(op binCode, d, x, y []float32) {
+	x, y = x[:len(d)], y[:len(d)]
 	switch op {
 	case opMax:
 		for i := range d {
@@ -89,7 +95,7 @@ func binFloat[T floatT](op binCode, d, x, y []T) {
 			d[i] = min(x[i], y[i])
 		}
 	default:
-		binNum(op, d, x, y)
+		kernels.Bin32(int(op), d, x, y)
 	}
 }
 
@@ -97,6 +103,7 @@ func binFloat[T floatT](op binCode, d, x, y []T) {
 // Max and min order lexicographically, real part first, and propagate NaN, as
 // NumPy does.
 func binComplex[T complexT](op binCode, d, x, y []T) {
+	x, y = x[:len(d)], y[:len(d)]
 	switch op {
 	case opMax:
 		for i := range d {
@@ -136,6 +143,7 @@ func lessComplex(x, y complex128) bool {
 // multiply and minimum are "and". Subtraction is an error before this is
 // reached, and division converts to a float type first.
 func binBool(op binCode, d, x, y []bool) {
+	x, y = x[:len(d)], y[:len(d)]
 	switch op {
 	case opAdd, opMax:
 		for i := range d {
@@ -150,6 +158,7 @@ func binBool(op binCode, d, x, y []bool) {
 
 // cmpOrdered compares real values.
 func cmpOrdered[T realT](op cmpCode, d []bool, x, y []T) {
+	x, y = x[:len(d)], y[:len(d)]
 	switch op {
 	case opEq:
 		for i := range d {
@@ -181,6 +190,7 @@ func cmpOrdered[T realT](op cmpCode, d []bool, x, y []T) {
 // cmpComplex compares complex values: equality componentwise, order
 // lexicographic as NumPy defines it (false when either side has a NaN).
 func cmpComplex[T complexT](op cmpCode, d []bool, x, y []T) {
+	x, y = x[:len(d)], y[:len(d)]
 	for i := range d {
 		a, b := complex128(x[i]), complex128(y[i])
 		nan := cmplx.IsNaN(a) || cmplx.IsNaN(b)
@@ -203,6 +213,7 @@ func cmpComplex[T complexT](op cmpCode, d []bool, x, y []T) {
 
 // cmpBool compares booleans, with false < true.
 func cmpBool(op cmpCode, d []bool, x, y []bool) {
+	x, y = x[:len(d)], y[:len(d)]
 	for i := range d {
 		a, b := x[i], y[i]
 		switch op {
@@ -271,6 +282,7 @@ var complexFuncs = [...]func(complex128) complex128{
 // negation, abs, square) this is exactly the float32 result: float64 carries
 // more than twice float32's precision, so the second rounding cannot change it.
 func unaryFloat32(f func(float64) float64, d, x []float32) {
+	x = x[:len(d)]
 	for i := range d {
 		d[i] = float32(f(float64(x[i])))
 	}
@@ -279,6 +291,7 @@ func unaryFloat32(f func(float64) float64, d, x []float32) {
 // unaryComplex applies f to complex elements, computing complex64 in
 // complex128.
 func unaryComplex[T complexT](f func(complex128) complex128, d, x []T) {
+	x = x[:len(d)]
 	for i := range d {
 		d[i] = T(f(complex128(x[i])))
 	}
diff --git a/zz_dtypes_gen.go b/zz_dtypes_gen.go
index 72b5be9..b3b3f91 100644
--- a/zz_dtypes_gen.go
+++ b/zz_dtypes_gen.go
@@ -1263,7 +1263,7 @@ func arithStore(op binCode, d, x, y any) {
 	case []uint64:
 		binInt(op, d, x.([]uint64), y.([]uint64))
 	case []float32:
-		binFloat(op, d, x.([]float32), y.([]float32))
+		binFloat32(op, d, x.([]float32), y.([]float32))
 	case []complex64:
 		binComplex(op, d, x.([]complex64), y.([]complex64))
 	case []complex128:

From 41ea3848ab8b0cc03a1554192b25f873252c0418 Mon Sep 17 00:00:00 2001
From: tannevaled <tannevaled@users.noreply.github.com>
Date: Sun, 18 Oct 2026 23:01:37 +0200
Subject: [PATCH 2/3] ci: gate on uncovered blocks, not a rounded percentage

The gate read the total of `go tool cover -func`, which rounds to one
decimal: it printed 100.0% with two statements never run
(random/generator.go fill32's error return, and the complex Log2). Both
are now tested, and the gate counts uncovered blocks with fleettools'
covblocks (pinned v0.27.0), which fails on any.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .github/workflows/ci.yml | 8 ++++----
 fidelity_test.go         | 6 ++++++
 random/api_test.go       | 1 +
 3 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index f787f30..9330cce 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -67,10 +67,10 @@ jobs:
           COVERPKG=$(go list ./... | paste -sd, -)
           go test -race -coverpkg="$COVERPKG" -coverprofile=cover.out ./...
           go tool cover -func=cover.out | tail -n 40
-          cov=$(go tool cover -func=cover.out | tail -1 | awk '{print $3}' | tr -d '%')
-          echo "coverage: ${cov}%"
-          awk -v t="$cov" 'BEGIN { exit !(t+0 >= 100.0) }' || {
-            echo "::error::coverage ${cov}% < 100%"; exit 1; }
+          # Count the uncovered blocks: `go tool cover -func` rounds to one
+          # decimal, and a gate on its total passed at 99.98% with two
+          # statements never run.
+          go run github.com/go-fleettools/fleettools/cmd/covblocks@v0.27.0 cover.out
 
   # Native per-arch execution. amd64/arm64 have GitHub-hosted runners, so the
   # generated assembly runs natively. Each regenerates the committed .s via
diff --git a/fidelity_test.go b/fidelity_test.go
index c77e882..0759411 100644
--- a/fidelity_test.go
+++ b/fidelity_test.go
@@ -3,6 +3,7 @@ package ndarray
 import (
 	"errors"
 	"math"
+	"math/cmplx"
 	"testing"
 )
 
@@ -115,4 +116,9 @@ func TestLog2Accuracy(t *testing.T) {
 			t.Errorf("Log2(%v) = %v, want %v", x.At(i), got[i], want[i])
 		}
 	}
+	// Complex: log(z)/ln 2, as numpy.log2 (log2(-1) = iπ/ln 2).
+	z := data[complex128](t, mustArr(t, ok(FromSlice([]complex128{8, -1}, 2))).Log2())
+	if cmplx.Abs(z[0]-3) > 1e-15 || cmplx.Abs(z[1]-complex(0, math.Pi/math.Ln2)) > 1e-15 {
+		t.Errorf("complex Log2 = %v", z)
+	}
 }
diff --git a/random/api_test.go b/random/api_test.go
index 0122ba1..e5ce781 100644
--- a/random/api_test.go
+++ b/random/api_test.go
@@ -120,6 +120,7 @@ func TestArgumentErrors(t *testing.T) {
 	add("choice nonzero")(g.ChoiceN(3, ChoiceOptions{NoReplace: true, P: []float64{1, 0, 0}}, 2))
 	add("normal size")(g.Normal(0, 1, -1))
 	add("random32 size")(g.RandomOf(ndarray.Float32, -1))
+	add("standard_normal32 size")(g.StandardNormalOf(ndarray.Float32, -1))
 	add("poisson size")(g.Poisson(1, -1))
 	add("standard_normal size")(g.StandardNormal(-1))
 	add("integers closed zero")(g.IntegersOf(ndarray.Int8, 0, -1, true))

From 39a6c446c0211c21a6d8882e18ae150772f9f7fb Mon Sep 17 00:00:00 2001
From: tannevaled <tannevaled@users.noreply.github.com>
Date: Mon, 19 Oct 2026 11:01:11 +0200
Subject: [PATCH 3/3] linalg: test the ZLAHQR rescale below a small subdiagonal
 on every arch

The block-count gate found linalg/eig.go:718 never run on amd64. The
branch rescales h(m+2, m+1) when a QR sweep starts at m > l with rows
below it. On arm64 TestEigBranches reached it by rounding alone; on
amd64 nothing did, and the rounded percentage had hidden that.

A new case plants two consecutive 1e-9 subdiagonals in an upper
Hessenberg matrix: each is far above the deflation threshold and their
product far below the two-small-subdiagonals test, so the sweep starts
at m = 3 on any arch. Checked: the block is run on arm64 and on amd64.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 linalg/properties_test.go | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/linalg/properties_test.go b/linalg/properties_test.go
index d9241b7..4261a13 100644
--- a/linalg/properties_test.go
+++ b/linalg/properties_test.go
@@ -408,6 +408,21 @@ func TestEigBranches(t *testing.T) {
 			}
 			return 0
 		}),
+		// Two consecutive small subdiagonals in the middle of an upper
+		// Hessenberg matrix: each is far above the deflation threshold, their
+		// product far below ZLAHQR's "two small subdiagonals" test, so a QR
+		// sweep starts at m = 3 with rows below it and rescales h(m+2, m+1)
+		// to keep the subdiagonal real. Reached by rounding alone on arm64
+		// only, before this case existed.
+		"small-subdiagonal-pair": mk(8, func(i, j int) complex128 {
+			switch {
+			case i > j+1:
+				return 0
+			case (i == 3 && j == 2) || (i == 4 && j == 3):
+				return 1e-9
+			}
+			return complex(float64(1+(i*5+j*3)%7), float64((i+2*j)%3))
+		}),
 		// Badly scaled: balancing scales by large powers of two.
 		"scaled": mk(4, func(i, j int) complex128 {
 			return complex(math.Pow(1e30, float64(i-j))*float64(1+(i*3+j)%5), 0)