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/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/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/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/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)
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))
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: