perf: float32 Add/Sub/Mul/Div kernels on six targets; *Into writes typed results in place

.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

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?
 

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

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

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)
+					}
+				}
+			})
+		}
+	}
+}

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`

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

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}) {

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)
+	}
 }

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("}")

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 (

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