From e12b82c2f139fab6f840081737a7da6f613d217f Mon Sep 17 00:00:00 2001
From: tannevaled <tannevaled@users.noreply.github.com>
Date: Wed, 7 Oct 2026 12:52:21 +0200
Subject: [PATCH 1/4] arm64: VFSQRT/VFADD/VFSUB/VFMUL mnemonics, NEON div and
 max/min kernels (Go 1.27)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 benchmarks/bench_test.go             |  13 ++
 internal/kernels/asmgen/arm64/gen.go | 194 +++++++++++----------------
 internal/kernels/reduce_arm64.go     |  73 +++++-----
 internal/kernels/reduce_simd_test.go |  26 ++++
 internal/kernels/sum_arm64.s         | 144 ++++++++++++++++----
 5 files changed, 275 insertions(+), 175 deletions(-)

diff --git a/benchmarks/bench_test.go b/benchmarks/bench_test.go
index 15bab9a..ed43bb5 100644
--- a/benchmarks/bench_test.go
+++ b/benchmarks/bench_test.go
@@ -170,6 +170,19 @@ func BenchmarkMulInto(b *testing.B) {
 	}
 }
 
+func BenchmarkDivInto(b *testing.B) {
+	for _, n := range elemSizes {
+		x, y := vec(n), vec(n)
+		out, _ := nd.New(n)
+		b.Run(fmt.Sprintf("n=%d", n), func(b *testing.B) {
+			b.SetBytes(int64(n * 8))
+			for i := 0; i < b.N; i++ {
+				_ = x.DivInto(out, y)
+			}
+		})
+	}
+}
+
 func BenchmarkSqrtInto(b *testing.B) {
 	for _, n := range elemSizes {
 		x := vec(n)
diff --git a/internal/kernels/asmgen/arm64/gen.go b/internal/kernels/asmgen/arm64/gen.go
index f2c925d..a9d74fa 100644
--- a/internal/kernels/asmgen/arm64/gen.go
+++ b/internal/kernels/asmgen/arm64/gen.go
@@ -43,9 +43,12 @@ func main() {
 	f.Add(sumKernel())
 	f.Add(dotKernel())
 	f.Add(sqrtKernel())
-	f.Add(addKernel())
-	f.Add(subKernel())
-	f.Add(mulKernel())
+	f.Add(binKernel("addNEON", "VFADD", "FADDD"))
+	f.Add(binKernel("subNEON", "VFSUB", "FSUBD"))
+	f.Add(binKernel("mulNEON", "VFMUL", "FMULD"))
+	f.Add(binKernel("divNEON", "VFDIV", "FDIVD"))
+	f.Add(extremumKernel("maxNEON", "VFMAX", "FMAXD"))
+	f.Add(extremumKernel("minNEON", "VFMIN", "FMIND"))
 	f.Add(gemmKernel())
 	if err := os.WriteFile("sum_arm64.s", []byte(f.String()), 0o644); err != nil {
 		fmt.Fprintln(os.Stderr, err)
@@ -65,39 +68,27 @@ func binSig() arm64.Signature {
 	)
 }
 
-// Go's arm64 assembler exposes NO plain vector double ADD/SUB/MUL/DIV — only the
-// fused multiply-accumulates VFMLA/VFMLS (the same wall the sum reduction hit).
-// So the elementwise add/sub/mul kernels reach the op through an FMA against a
-// constant vector, each of which rounds bit-identically to the plain op:
+// The elementwise kernels use the plain vector double ops VFADD/VFSUB/VFMUL/
+// VFDIV, which Go's arm64 assembler names since Go 1.27 (this module requires
+// 1.27.1). Before that only the fused VFMLA/VFMLS existed, and add/sub/mul went
+// through an FMA against an exact constant (b + a*1.0, a - b*1.0, -0.0 + a*b),
+// bit-identical but harder to read, while division stayed scalar. Each lane is
+// one correctly rounded IEEE-754 operation, the same the scalar FADDD/FSUBD/
+// FMULD/FDIVD perform, so the kernels are bit-identical to the scalar oracle
+// (TestBinSIMD). Go's operand order is VFSUB Vm, Vn, Vd: Vd = Vn - Vm.
 //
-//	add:  dst = b + a*1.0   via VFMLA (a, 1.0, acc=b)   — FMA(a,1,b) == a+b
-//	sub:  dst = a - b*1.0   via VFMLS (b, 1.0, acc=a)   — FMA(-b,1,a) == a-b
-//	mul:  dst = 0 + a*b     via VFMLA (a, b,   acc=0)   — FMA(a,b,0)  == a*b
-//
-// In each case the FMA's extra addend is exact (a 1.0 multiply, or a +0.0 add),
-// so the single IEEE rounding lands on exactly the value the plain operation
-// would, making these bit-identical to the scalar Add/Sub/Mul oracle (validated
-// per-arch in CI). Division has neither an FMA form nor a vector-double divide on
-// arm64, so Div stays on the scalar FDIVD oracle (divBin in reduce_arm64.go).
-//
-// VFMLA Va, Vb, Vacc computes Vacc += Va*Vb; VFMLS computes Vacc -= Va*Vb (the
-// Go/asmgen operand order). The loop does 8 doubles (64 bytes) per iteration
-// across four D2 lanes, then a scalar FP tail for the (n mod 8) remainder.
+// The loop does 8 doubles (64 bytes) per iteration across four D2 lanes, then a
+// scalar FP tail for the (n mod 8) remainder.
 
 // sqrtKernel builds sqrtNEON(dst, src *float64, n int): dst[i] = sqrt(src[i])
 // with the packed NEON double square root FSQRT Vd.2D (2 doubles/op, unrolled 4x
 // = 8 lanes/iter) plus a scalar FSQRTD tail.
 //
-// Go's arm64 assembler exposes only the *scalar* FSQRTD — there is no VFSQRT
-// mnemonic for the vector form — so the vector instruction is emitted directly as
-// its fixed 32-bit encoding via WORD. FSQRT Vd.2D, Vn.2D is the Advanced-SIMD
-// two-register-misc fp64 square root: 0110_1110_1110_0001_1111_10nn_nnnd_dddd,
-// i.e. base 0x6EE1F800 | (Vn<<5) | Vd. This is the same hardware NEON double
-// sqrt OpenBLAS/numpy use; it is correctly-rounded IEEE-754, identical lane-by-
-// lane to the scalar FSQRTD (and to math.Sqrt), so the kernel is bit-identical to
-// the scalar sqrt oracle including negatives->NaN, ±Inf, and signed zeros
-// (validated bit-for-bit per-arch in CI). Using raw WORD keeps this pure Go asm
-// (no cgo); it is the one instruction the assembler cannot name.
+// FSQRT Vd.2D is VFSQRT Vn.D2, Vd.D2 in Go since 1.27; it was emitted as a raw
+// WORD before (0x6EE1F800 | Vn<<5 | Vd, the same encoding, checked against Go's
+// own disassembly). It is correctly rounded IEEE-754, identical lane by lane to
+// the scalar FSQRTD and math.Sqrt, so the kernel is bit-identical to the scalar
+// oracle including negatives->NaN, ±Inf and signed zeros.
 //
 // This replaces the previous arm64 sqrt path (the compiler's scalar FSQRTD loop,
 // one lane at a time): the packed form does two lanes per instruction, which is
@@ -117,10 +108,10 @@ func sqrtKernel() *emit.Function {
 		Raw("CMP $8, R2").
 		Raw("BLT tail").
 		Raw("VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2]").
-		Raw("WORD $0x6EE1F800"). // FSQRT V0.2D, V0.2D
-		Raw("WORD $0x6EE1F821"). // FSQRT V1.2D, V1.2D
-		Raw("WORD $0x6EE1F842"). // FSQRT V2.2D, V2.2D
-		Raw("WORD $0x6EE1F863"). // FSQRT V3.2D, V3.2D
+		Raw("VFSQRT V0.D2, V0.D2").
+		Raw("VFSQRT V1.D2, V1.D2").
+		Raw("VFSQRT V2.D2, V2.D2").
+		Raw("VFSQRT V3.D2, V3.D2").
 		Raw("VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0)").
 		Raw("SUB $8, R2").
 		Raw("B block").
@@ -137,91 +128,63 @@ func sqrtKernel() *emit.Function {
 	return b.Func()
 }
 
-// onesVec emits V31 = [1.0, 1.0] (the IEEE bit pattern of 1.0 broadcast).
-func onesVec(b *arm64.Builder) {
-	b.Raw("MOVD $0x3FF0000000000000, R7")
-	b.Raw("VDUP R7, V31.D2")
-}
-
-// addKernel builds addNEON(dst, a, b *float64, n int): dst = a + b via FMA.
-func addKernel() *emit.Function {
-	b := arm64.NewFunc("addNEON", binSig(), 0)
-	b.LoadArg("dst", "R0").LoadArg("a", "R1").LoadArg("b", "R2").LoadArg("n", "R3")
-	onesVec(b)
-	b.Raw("block:").Raw("CMP $8, R3").Raw("BLT tail").
-		Raw("VLD1.P 64(R2), [V0.D2, V1.D2, V2.D2, V3.D2]"). // acc = b
-		Raw("VLD1.P 64(R1), [V4.D2, V5.D2, V6.D2, V7.D2]"). // a
-		Raw("VFMLA V4.D2, V31.D2, V0.D2").                  // acc += a*1.0 => a+b
-		Raw("VFMLA V5.D2, V31.D2, V1.D2").
-		Raw("VFMLA V6.D2, V31.D2, V2.D2").
-		Raw("VFMLA V7.D2, V31.D2, V3.D2").
-		Raw("VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0)").
-		Raw("SUB $8, R3").Raw("B block").
-		Raw("tail:").Raw("CBZ R3, done").
-		Raw("FMOVD.P 8(R1), F0").Raw("FMOVD.P 8(R2), F1").
-		Raw("FADDD F1, F0, F0").
-		Raw("FMOVD.P F0, 8(R0)").
-		Raw("SUB $1, R3").Raw("B tail").
-		Raw("done:").Ret()
-	return b.Func()
-}
-
-// subKernel builds subNEON(dst, a, b *float64, n int): dst = a - b via FMA.
-func subKernel() *emit.Function {
-	b := arm64.NewFunc("subNEON", binSig(), 0)
+// binKernel builds NAME(dst, a, b *float64, n int): dst[i] = a[i] OP b[i],
+// with vop the vector form (Go operand order VOP Vm, Vn, Vd: Vd = Vn OP Vm) and
+// sop the scalar one for the tail (SOP Fm, Fn, Fd: Fd = Fn OP Fm).
+func binKernel(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")
-	onesVec(b)
 	b.Raw("block:").Raw("CMP $8, R3").Raw("BLT tail").
-		Raw("VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2]"). // acc = a
-		Raw("VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2]"). // b
-		Raw("VFMLS V4.D2, V31.D2, V0.D2").                  // acc -= b*1.0 => a-b
-		Raw("VFMLS V5.D2, V31.D2, V1.D2").
-		Raw("VFMLS V6.D2, V31.D2, V2.D2").
-		Raw("VFMLS V7.D2, V31.D2, V3.D2").
-		Raw("VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0)").
+		Raw("VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2]"). // a
+		Raw("VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2]")  // b
+	for r := 0; r < 4; r++ {
+		b.Raw(fmt.Sprintf("%s V%d.D2, V%d.D2, V%d.D2", vop, 4+r, r, r)) // a OP b
+	}
+	b.Raw("VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0)").
 		Raw("SUB $8, R3").Raw("B block").
 		Raw("tail:").Raw("CBZ R3, done").
 		Raw("FMOVD.P 8(R1), F0").Raw("FMOVD.P 8(R2), F1").
-		Raw("FSUBD F1, F0, F0").
+		Raw(sop + " F1, F0, F0").
 		Raw("FMOVD.P F0, 8(R0)").
 		Raw("SUB $1, R3").Raw("B tail").
 		Raw("done:").Ret()
 	return b.Func()
 }
 
-// mulKernel builds mulNEON(dst, a, b *float64, n int): dst = a * b via FMA into a
-// -0.0 accumulator: FMA(a, b, -0.0) == a*b *bit-for-bit, including sign-of-zero*.
-// Adding +0.0 would NOT be safe — when the product is -0.0, (-0.0)+(+0.0) rounds
-// to +0.0 under round-to-nearest, flipping the sign that plain a*b keeps. Adding
-// -0.0 is the correct identity: -0.0+(-0.0) = -0.0 and (+0.0)+(-0.0) = +0.0, so
-// the result sign matches the product's, and for any nonzero product adding -0.0
-// changes nothing. Hence the tile is bit-identical to the scalar Mul oracle
-// (asserted in TestBinSIMD, which pairs -0.0*0.0 etc.).
-func mulKernel() *emit.Function {
-	b := arm64.NewFunc("mulNEON", binSig(), 0)
-	b.LoadArg("dst", "R0").LoadArg("a", "R1").LoadArg("b", "R2").LoadArg("n", "R3")
-	// V30 = [-0.0, -0.0] (IEEE bits 0x8000000000000000), the FMA accumulator seed.
-	b.Raw("MOVD $0x8000000000000000, R7")
-	b.Raw("VDUP R7, V30.D2")
-	b.Raw("block:").Raw("CMP $8, R3").Raw("BLT tail").
-		Raw("VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2]"). // a
-		Raw("VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2]"). // b
-		Raw("VMOV V30.B16, V8.B16").                        // acc = -0.0
-		Raw("VMOV V30.B16, V9.B16").
-		Raw("VMOV V30.B16, V10.B16").
-		Raw("VMOV V30.B16, V11.B16").
-		Raw("VFMLA V0.D2, V4.D2, V8.D2"). // acc += a*b => a*b (sign-exact)
-		Raw("VFMLA V1.D2, V5.D2, V9.D2").
-		Raw("VFMLA V2.D2, V6.D2, V10.D2").
-		Raw("VFMLA V3.D2, V7.D2, V11.D2").
-		Raw("VST1.P [V8.D2, V9.D2, V10.D2, V11.D2], 64(R0)").
-		Raw("SUB $8, R3").Raw("B block").
-		Raw("tail:").Raw("CBZ R3, done").
-		Raw("FMOVD.P 8(R1), F0").Raw("FMOVD.P 8(R2), F1").
-		Raw("FMULD F1, F0, F0").
-		Raw("FMOVD.P F0, 8(R0)").
-		Raw("SUB $1, R3").Raw("B tail").
-		Raw("done:").Ret()
+// extremumKernel builds NAME(a *float64, n int) float64, the maximum (or
+// minimum) of a, for n >= 8: four D2 accumulators seeded with the first eight
+// elements, folded together, then the two lanes, then the tail one at a time.
+// FMAX/FMIN return NaN when either operand is NaN and order -0 below +0, as
+// the scalar FMAXD/FMIND that Go's builtin max/min compile to do, and the
+// operation is associative under those rules, so the result is bit-identical
+// to the serial oracle (TestMaxMinSIMD).
+func extremumKernel(name, vop, sop string) *emit.Function {
+	b := arm64.NewFunc(name, reduceSig(), 0)
+	b.LoadArg("a", "R0").LoadArg("n", "R1").
+		Raw("VLD1.P 64(R0), [V0.D2, V1.D2, V2.D2, V3.D2]").
+		Raw("SUB $8, R1").
+		Raw("block:").Raw("CMP $8, R1").Raw("BLT fold").
+		Raw("VLD1.P 64(R0), [V4.D2, V5.D2, V6.D2, V7.D2]")
+	for r := 0; r < 4; r++ {
+		b.Raw(fmt.Sprintf("%s V%d.D2, V%d.D2, V%d.D2", vop, 4+r, r, r))
+	}
+	b.Raw("SUB $8, R1").Raw("B block").
+		Raw("fold:").
+		Raw(vop+" V1.D2, V0.D2, V0.D2").
+		Raw(vop+" V3.D2, V2.D2, V2.D2").
+		Raw(vop+" V2.D2, V0.D2, V0.D2").
+		Raw("VMOV V0.D[0], R3").
+		Raw("VMOV V0.D[1], R4").
+		Raw("FMOVD R3, F0").
+		Raw("FMOVD R4, F9").
+		Raw(sop+" F9, F0, F0").
+		Raw("tail:").Raw("CBZ R1, done").
+		Raw("FMOVD.P 8(R0), F10").
+		Raw(sop+" F10, F0, F0").
+		Raw("SUB $1, R1").Raw("B tail").
+		Raw("done:").
+		StoreRet("F0", "ret").
+		Ret()
 	return b.Func()
 }
 
@@ -290,8 +253,8 @@ func fmlaElem(acc, breg, areg, idx int) string {
 //     ~1.28x faster on the L1-resident micro-kernel (44 -> 58 GFLOP/s/core);
 //     the earlier note claiming the indexed form is "slower" was wrong (it was
 //     never actually measured against a correct encoding). See docs/perf.md.
-//  2. There is no plain vector FP add (only VFMLA/VFMLS), so the closing C += acc
-//     is done as C = C + acc*1.0 via VFMLA against V31=[1.0,1.0]. Multiplying by
+//  2. The closing C += acc is done as C = C + acc*1.0 via VFMLA against
+//     V31=[1.0,1.0], written when Go's assembler had no VFADD (it has since 1.27). Multiplying by
 //     exactly 1.0 is exact and FMA(acc,1.0,C) rounds identically to C+acc (C is
 //     the only inexact addend), so this is a true add — the tile result is
 //     bit-identical to the scalar oracle's ikj accumulation order.
@@ -347,8 +310,8 @@ func gemmKernel() *emit.Function {
 		Raw(fmlaElem(15, 19, 21, 1)).
 		Raw("SUB $1, R0").
 		Raw("CBNZ R0, gkloop").
-		// C += acc, row by row, via VFMLA against V31=[1.0,1.0] (no vector FP add
-		// available — see header note 2). Bit-identical to C+acc.
+		// C += acc, row by row, via VFMLA against V31=[1.0,1.0]: exact, so
+		// bit-identical to C+acc (written before Go 1.27 named VFADD).
 		Raw("gstore:").
 		Raw("MOVD $0x3FF0000000000000, R6").
 		Raw("VDUP R6, V31.D2").
@@ -389,9 +352,8 @@ func gemmKernel() *emit.Function {
 
 // sumKernel builds sumNEON(a *float64, n int) float64.
 //
-// Go's arm64 assembler exposes no plain vector double ADD (only VFMLA/VFMLS,
-// the fused multiply-accumulates). So the lane-parallel accumulation is done
-// with VFMLA against a vector of 1.0: acc += src * 1.0. Multiplying a double by
+// The lane-parallel accumulation is VFMLA against a vector of 1.0: acc += src *
+// 1.0 (written when Go's assembler had no VFADD; it has since Go 1.27). Multiplying a double by
 // exactly 1.0 is exact, and FMA(src, 1.0, acc) rounds identically to acc+src
 // (the addend is the only inexact step), so this is a true vector add, not an
 // approximation — it simply reaches it through the only vector-FP op available.
@@ -399,7 +361,7 @@ func gemmKernel() *emit.Function {
 // with eight D2 accumulators V0..V7 (16 doubles per iteration): two loads per
 // FMA make the loop load-bound, and eight chains keep the four FP pipes of an
 // Apple core fed through the FMA latency. The accumulators are folded with
-// FMLA against 1.0 (exact; Go's assembler has no plain vector FP add), the two
+// FMLA against 1.0 (exact; written before Go 1.27 named VFADD), the two
 // lanes added, and the tail fused in one element at a time with FMADDD, as gc
 // fuses the scalar loop.
 func dotKernel() *emit.Function {
diff --git a/internal/kernels/reduce_arm64.go b/internal/kernels/reduce_arm64.go
index 02f1317..f2460aa 100644
--- a/internal/kernels/reduce_arm64.go
+++ b/internal/kernels/reduce_arm64.go
@@ -5,14 +5,10 @@ package kernels
 //
 //go:generate sh -c "cd asmgen/arm64 && GOFLAGS=-mod=mod go run gen.go && mv sum_arm64.s ../../sum_arm64.s"
 //
-// sqrt/max/min are NOT hand-vectorized on arm64: Go's arm64 assembler exposes
-// no packed vector-double sqrt/max/min (only VFMLA/VFMLS), but it DOES intrinsify
-// the *scalar* math.Sqrt -> FSQRTD and math.Max/math.Min -> FMAXD/FMIND. So the
-// direct scalar oracle loops already compile to the hardware instructions with no
-// per-element call overhead — the win over the old func-pointer Map path (which
-// could not intrinsify math.Sqrt and paid a real call per element, ~2.7x slower)
-// plus the multicore fan-out is what beats numpy. arm64's FMAXD/FMIND propagate
-// NaN, matching the math.Max-based oracle and numpy.
+// Every kernel here is NEON: sum, dot, sqrt, add/sub/mul/div, max/min and the
+// GEMM micro-kernel. The vector sqrt, divide and max/min use mnemonics Go's
+// assembler gained in 1.27 (VFSQRT, VFDIV, VFMAX, VFMIN); before that, sqrt was
+// a WORD and divide and max/min stayed scalar.
 
 // HaveReduceSIMD reports that this build ships at least one hand-vectorized SIMD
 // kernel (the NEON sum). The per-arch CI execution job validates it against the
@@ -29,8 +25,7 @@ func sumSIMD(a []float64) float64 {
 }
 
 // sqrtSIMD writes sqrt(src[i]) into dst[i] with the packed NEON FSQRT V.2D kernel
-// (two doubles per instruction, emitted via WORD because Go's arm64 assembler has
-// no vector-sqrt mnemonic — see asmgen/arm64/gen.go). FSQRT V.2D is correctly-
+// (two doubles per instruction, VFSQRT in Go). FSQRT V.2D is correctly-
 // rounded IEEE-754, identical lane-by-lane to the scalar FSQRTD/math.Sqrt, so it
 // is bit-identical to sqrtScalar including negatives->NaN, ±Inf and signed zeros.
 // It supersedes the prior scalar-FSQRTD loop, halving the per-element sqrt count.
@@ -42,22 +37,28 @@ func sqrtSIMD(dst, src []float64) {
 	sqrtNEON(&dst[0], &src[0], len(dst))
 }
 
-// maxSIMD / minSIMD are the four-accumulator NaN-propagating reducers: the
-// builtin max/min lower to the FMAXD/FMIND hardware instructions and the four
-// independent chains hide their latency (~3.6x over one accumulator). They are
-// bit-identical to the serial Max/Min oracle (max/min is associative for the
-// NaN-propagating rule). A packed vector max is not expressible on arm64 (Go's
-// assembler exposes no VFMAX.D2, only VFMLA/VFMLS), so this register-unrolled
-// scalar form is the arm64 SIMD path.
-func maxSIMD(a []float64) float64 { return maxUnrolled(a) }
-func minSIMD(a []float64) float64 { return minUnrolled(a) }
-
-// addBin/subBin/mulBin write a[i] OP b[i] into dst[i] with the NEON D2 kernels
-// (8 doubles/iter + scalar tail). Go's arm64 assembler has no plain vector
-// double add/sub/mul (only VFMLA/VFMLS), so each reaches the op through an FMA
-// against an exact constant — add: b+a*1.0, sub: a-b*1.0, mul: 0+a*b — which
-// rounds bit-identically to the plain op (the extra addend is exact), so the
-// result is bit-identical to the scalar Add/Sub/Mul oracle. They are the serial
+// maxSIMD / minSIMD are the NaN-propagating extrema of a (non-empty): the NEON
+// FMAX/FMIN kernel with four D2 accumulators from eight elements up, the
+// four-chain scalar reducer below that. Bit-identical to the serial Max/Min
+// oracle (see extremumKernel in asmgen/arm64/gen.go).
+func maxSIMD(a []float64) float64 {
+	if len(a) < 8 {
+		return maxUnrolled(a)
+	}
+	return maxNEON(&a[0], len(a))
+}
+
+func minSIMD(a []float64) float64 {
+	if len(a) < 8 {
+		return minUnrolled(a)
+	}
+	return minNEON(&a[0], len(a))
+}
+
+// addBin/subBin/mulBin/divBin write a[i] OP b[i] into dst[i] with the NEON D2
+// kernels (VFADD/VFSUB/VFMUL/VFDIV, 8 doubles/iter + scalar tail). Each lane is
+// one correctly rounded IEEE-754 operation, so the result is bit-identical to
+// the scalar Add/Sub/Mul/Div oracle. They are the serial
 // elementwise inner loop the per-op fast path and runBinaryP use. Empty slices
 // are a no-op (no &a[0]).
 func addBin(dst, a, b []float64) {
@@ -84,10 +85,13 @@ func mulBin(dst, a, b []float64) {
 	mulNEON(&dst[0], &a[0], &b[0], len(dst))
 }
 
-// divBin stays on the scalar oracle: arm64 has neither a vector-double divide nor
-// an FMA form for division, and the scalar Div loop already lowers to the FDIVD
-// hardware instruction per element (the same throughput-bound op numpy runs).
-func divBin(dst, a, b []float64) { Div(dst, a, b) }
+func divBin(dst, a, b []float64) {
+	if len(dst) == 0 {
+		return
+	}
+	_, _ = a[len(dst)-1], b[len(dst)-1] // the kernel reads len(dst) of each
+	divNEON(&dst[0], &a[0], &b[0], len(dst))
+}
 
 //go:noescape
 func sumNEON(a *float64, n int) float64
@@ -103,3 +107,12 @@ func subNEON(dst, a, b *float64, n int)
 
 //go:noescape
 func mulNEON(dst, a, b *float64, n int)
+
+//go:noescape
+func divNEON(dst, a, b *float64, n int)
+
+//go:noescape
+func maxNEON(a *float64, n int) float64
+
+//go:noescape
+func minNEON(a *float64, n int) float64
diff --git a/internal/kernels/reduce_simd_test.go b/internal/kernels/reduce_simd_test.go
index 30973da..ca4caa9 100644
--- a/internal/kernels/reduce_simd_test.go
+++ b/internal/kernels/reduce_simd_test.go
@@ -255,3 +255,29 @@ func BenchmarkSumSIMD(b *testing.B) {
 		}
 	})
 }
+
+// TestMaxMinSIMDSpecialsInVectorLanes puts signed zeros and infinities at every
+// position of arrays long enough for the vector kernels (8 and up), where they
+// land in each accumulator lane and in the tail: -0 must order below +0 and the
+// extremes must come through, bit for bit as the oracle gives them.
+func TestMaxMinSIMDSpecialsInVectorLanes(t *testing.T) {
+	for n := 8; n <= 40; n++ {
+		for pos := 0; pos < n; pos++ {
+			for _, fill := range []float64{-0.0, 0.0, -1, 1} {
+				for _, v := range []float64{-0.0, 0.0, math.Inf(1), math.Inf(-1)} {
+					a := make([]float64, n)
+					for i := range a {
+						a[i] = fill
+					}
+					a[pos] = v
+					if g, w := maxSIMD(a), Max(a); !bitEq(g, w) || math.Signbit(g) != math.Signbit(w) {
+						t.Fatalf("n=%d fill=%v a[%d]=%v: maxSIMD=%v, Max=%v", n, fill, pos, v, g, w)
+					}
+					if g, w := minSIMD(a), Min(a); !bitEq(g, w) || math.Signbit(g) != math.Signbit(w) {
+						t.Fatalf("n=%d fill=%v a[%d]=%v: minSIMD=%v, Min=%v", n, fill, pos, v, g, w)
+					}
+				}
+			}
+		}
+	}
+}
diff --git a/internal/kernels/sum_arm64.s b/internal/kernels/sum_arm64.s
index 49c984a..a3c1791 100644
--- a/internal/kernels/sum_arm64.s
+++ b/internal/kernels/sum_arm64.s
@@ -105,10 +105,10 @@ TEXT ·sqrtNEON(SB), NOSPLIT, $0-24
 	CMP $8, R2
 	BLT tail
 	VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2]
-	WORD $0x6EE1F800
-	WORD $0x6EE1F821
-	WORD $0x6EE1F842
-	WORD $0x6EE1F863
+	VFSQRT V0.D2, V0.D2
+	VFSQRT V1.D2, V1.D2
+	VFSQRT V2.D2, V2.D2
+	VFSQRT V3.D2, V3.D2
 	VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0)
 	SUB $8, R2
 	B block
@@ -127,17 +127,15 @@ TEXT ·addNEON(SB), NOSPLIT, $0-32
 	MOVD a+8(FP), R1
 	MOVD b+16(FP), R2
 	MOVD n+24(FP), R3
-	MOVD $0x3FF0000000000000, R7
-	VDUP R7, V31.D2
 	block:
 	CMP $8, R3
 	BLT tail
-	VLD1.P 64(R2), [V0.D2, V1.D2, V2.D2, V3.D2]
-	VLD1.P 64(R1), [V4.D2, V5.D2, V6.D2, V7.D2]
-	VFMLA V4.D2, V31.D2, V0.D2
-	VFMLA V5.D2, V31.D2, V1.D2
-	VFMLA V6.D2, V31.D2, V2.D2
-	VFMLA V7.D2, V31.D2, V3.D2
+	VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2]
+	VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2]
+	VFADD V4.D2, V0.D2, V0.D2
+	VFADD V5.D2, V1.D2, V1.D2
+	VFADD V6.D2, V2.D2, V2.D2
+	VFADD V7.D2, V3.D2, V3.D2
 	VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0)
 	SUB $8, R3
 	B block
@@ -157,17 +155,15 @@ TEXT ·subNEON(SB), NOSPLIT, $0-32
 	MOVD a+8(FP), R1
 	MOVD b+16(FP), R2
 	MOVD n+24(FP), R3
-	MOVD $0x3FF0000000000000, R7
-	VDUP R7, V31.D2
 	block:
 	CMP $8, R3
 	BLT tail
 	VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2]
 	VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2]
-	VFMLS V4.D2, V31.D2, V0.D2
-	VFMLS V5.D2, V31.D2, V1.D2
-	VFMLS V6.D2, V31.D2, V2.D2
-	VFMLS V7.D2, V31.D2, V3.D2
+	VFSUB V4.D2, V0.D2, V0.D2
+	VFSUB V5.D2, V1.D2, V1.D2
+	VFSUB V6.D2, V2.D2, V2.D2
+	VFSUB V7.D2, V3.D2, V3.D2
 	VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0)
 	SUB $8, R3
 	B block
@@ -187,22 +183,16 @@ TEXT ·mulNEON(SB), NOSPLIT, $0-32
 	MOVD a+8(FP), R1
 	MOVD b+16(FP), R2
 	MOVD n+24(FP), R3
-	MOVD $0x8000000000000000, R7
-	VDUP R7, V30.D2
 	block:
 	CMP $8, R3
 	BLT tail
 	VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2]
 	VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2]
-	VMOV V30.B16, V8.B16
-	VMOV V30.B16, V9.B16
-	VMOV V30.B16, V10.B16
-	VMOV V30.B16, V11.B16
-	VFMLA V0.D2, V4.D2, V8.D2
-	VFMLA V1.D2, V5.D2, V9.D2
-	VFMLA V2.D2, V6.D2, V10.D2
-	VFMLA V3.D2, V7.D2, V11.D2
-	VST1.P [V8.D2, V9.D2, V10.D2, V11.D2], 64(R0)
+	VFMUL V4.D2, V0.D2, V0.D2
+	VFMUL V5.D2, V1.D2, V1.D2
+	VFMUL V6.D2, V2.D2, V2.D2
+	VFMUL V7.D2, V3.D2, V3.D2
+	VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0)
 	SUB $8, R3
 	B block
 	tail:
@@ -216,6 +206,102 @@ TEXT ·mulNEON(SB), NOSPLIT, $0-32
 	done:
 	RET
 
+TEXT ·divNEON(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 $8, R3
+	BLT tail
+	VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2]
+	VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2]
+	VFDIV V4.D2, V0.D2, V0.D2
+	VFDIV V5.D2, V1.D2, V1.D2
+	VFDIV V6.D2, V2.D2, V2.D2
+	VFDIV V7.D2, V3.D2, V3.D2
+	VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0)
+	SUB $8, R3
+	B block
+	tail:
+	CBZ R3, done
+	FMOVD.P 8(R1), F0
+	FMOVD.P 8(R2), F1
+	FDIVD F1, F0, F0
+	FMOVD.P F0, 8(R0)
+	SUB $1, R3
+	B tail
+	done:
+	RET
+
+TEXT ·maxNEON(SB), NOSPLIT, $0-24
+	MOVD a+0(FP), R0
+	MOVD n+8(FP), R1
+	VLD1.P 64(R0), [V0.D2, V1.D2, V2.D2, V3.D2]
+	SUB $8, R1
+	block:
+	CMP $8, R1
+	BLT fold
+	VLD1.P 64(R0), [V4.D2, V5.D2, V6.D2, V7.D2]
+	VFMAX V4.D2, V0.D2, V0.D2
+	VFMAX V5.D2, V1.D2, V1.D2
+	VFMAX V6.D2, V2.D2, V2.D2
+	VFMAX V7.D2, V3.D2, V3.D2
+	SUB $8, R1
+	B block
+	fold:
+	VFMAX V1.D2, V0.D2, V0.D2
+	VFMAX V3.D2, V2.D2, V2.D2
+	VFMAX V2.D2, V0.D2, V0.D2
+	VMOV V0.D[0], R3
+	VMOV V0.D[1], R4
+	FMOVD R3, F0
+	FMOVD R4, F9
+	FMAXD F9, F0, F0
+	tail:
+	CBZ R1, done
+	FMOVD.P 8(R0), F10
+	FMAXD F10, F0, F0
+	SUB $1, R1
+	B tail
+	done:
+	FMOVD F0, ret+16(FP)
+	RET
+
+TEXT ·minNEON(SB), NOSPLIT, $0-24
+	MOVD a+0(FP), R0
+	MOVD n+8(FP), R1
+	VLD1.P 64(R0), [V0.D2, V1.D2, V2.D2, V3.D2]
+	SUB $8, R1
+	block:
+	CMP $8, R1
+	BLT fold
+	VLD1.P 64(R0), [V4.D2, V5.D2, V6.D2, V7.D2]
+	VFMIN V4.D2, V0.D2, V0.D2
+	VFMIN V5.D2, V1.D2, V1.D2
+	VFMIN V6.D2, V2.D2, V2.D2
+	VFMIN V7.D2, V3.D2, V3.D2
+	SUB $8, R1
+	B block
+	fold:
+	VFMIN V1.D2, V0.D2, V0.D2
+	VFMIN V3.D2, V2.D2, V2.D2
+	VFMIN V2.D2, V0.D2, V0.D2
+	VMOV V0.D[0], R3
+	VMOV V0.D[1], R4
+	FMOVD R3, F0
+	FMOVD R4, F9
+	FMIND F9, F0, F0
+	tail:
+	CBZ R1, done
+	FMOVD.P 8(R0), F10
+	FMIND F10, F0, F0
+	SUB $1, R1
+	B tail
+	done:
+	FMOVD F0, ret+16(FP)
+	RET
+
 TEXT ·gemmMicro4x8(SB), NOSPLIT, $0-40
 	MOVD kc+0(FP), R0
 	MOVD pa+8(FP), R1

From ee5e247e73b0e95ab71193392aff3605421b6534 Mon Sep 17 00:00:00 2001
From: tannevaled <tannevaled@users.noreply.github.com>
Date: Wed, 7 Oct 2026 12:52:32 +0200
Subject: [PATCH 2/4] bench: DivInto without zero divisors

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 benchmarks/bench_test.go | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmarks/bench_test.go b/benchmarks/bench_test.go
index ed43bb5..a29ec7f 100644
--- a/benchmarks/bench_test.go
+++ b/benchmarks/bench_test.go
@@ -172,7 +172,7 @@ func BenchmarkMulInto(b *testing.B) {
 
 func BenchmarkDivInto(b *testing.B) {
 	for _, n := range elemSizes {
-		x, y := vec(n), vec(n)
+		x, y := vec(n), vec(n).AddScalar(1) // no zero divisors
 		out, _ := nd.New(n)
 		b.Run(fmt.Sprintf("n=%d", n), func(b *testing.B) {
 			b.SetBytes(int64(n * 8))

From 8beb837382cdcc31734eee4fe09a263c4f1b801e Mon Sep 17 00:00:00 2001
From: tannevaled <tannevaled@users.noreply.github.com>
Date: Wed, 7 Oct 2026 12:57:27 +0200
Subject: [PATCH 3/4] docs: arm64 NEON div/max/min and mnemonics; parallelFor
 scheduling measurements replace the cache-affinity guess

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 README.md            |  2 +-
 docs/perf.md         | 71 +++++++++++++++++++++++++++++++++-----------
 docs/plan-ndarray.md |  2 +-
 3 files changed, 55 insertions(+), 20 deletions(-)

diff --git a/README.md b/README.md
index 0ace2f9..441fea5 100644
--- a/README.md
+++ b/README.md
@@ -63,7 +63,7 @@ elementwise ops, reductions and products run **multicore** (across
 [go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels: sum,
 add/sub/mul/div, sqrt and the dot product (SSE2 or AVX2/FMA on amd64, NEON on
 arm64, VSX on ppc64le, LASX on loong64 when the CPU has it), max/min too on
-amd64, and a **panel-packed, cache-blocked GEMM** with an SIMD-FMA
+amd64 and arm64, and a **panel-packed, cache-blocked GEMM** with an SIMD-FMA
 micro-kernel (NEON 4×8; AVX2/FMA 6×8 chosen at run time, SSE2 fallback; VSX
 8×8; LASX 8×8), the OpenBLAS/BLIS structure. The other 64-bit targets
 (riscv64, s390x) and the 32-bit ones run the same pure-Go code those kernels
diff --git a/docs/perf.md b/docs/perf.md
index 10cfb6a..a6af69a 100644
--- a/docs/perf.md
+++ b/docs/perf.md
@@ -95,8 +95,8 @@ with no result array, so they have no separate *into* row.)
   arm64/others — plus the existing multicore fan-out is what turns the loss into
   a win. The kernel is **bit-identical** to a scalar `math.Sqrt` loop (validated
   per-arch and bit-for-bit against NumPy's `np.sqrt`). arm64 now also has a
-  **packed NEON `FSQRT V.2D`** kernel (two doubles per instruction, emitted by
-  `WORD` since Go's assembler has no vector-sqrt mnemonic), replacing the prior
+  **packed NEON `FSQRT V.2D`** kernel (two doubles per instruction; emitted by
+  `WORD` until v0.5.2, `VFSQRT` since Go 1.27 names it), replacing the prior
   one-lane-at-a-time scalar `FSQRTD` loop. Small `n=1 024`: the *alloc* form loses
   to NumPy on allocation cost (as Add/Mul), but `a.SqrtInto(out)` wins 1.8×.
 - **Max: go-ndarray now WINS — 4.0× small, 1.34× at 4 M (was 0.62 / 0.31
@@ -309,7 +309,7 @@ The SIMD kernels are also validated against the scalar oracle per-arch in CI:
 | `Add`/`Sub`/`Mul`/`Div`/`Sqrt`, **alloc** form, small n (≈1 K) | NumPy wins (Go-allocator ceiling) | the result `make` (~900 ns zeroing 8 KiB) is the whole cost; Go has no unzeroed-alloc escape and pooling would break result ownership. The `*Into` no-alloc form **wins 1.8–2.4×** — the supported parity path. Large-n alloc form already wins ~2×. |
 | `Add`/`Sub`/`Mul`/`Div` | **FIXED — go wins ~2× large n, *Into* wins all n** | new SIMD kernels: packed SSE2 `ADD/SUB/MUL/DIVPD` (amd64), NEON `VFMLA`/`VFMLS` (arm64); bit-identical to the scalar oracle |
 | `Sqrt` | **FIXED — go wins 2.8× at 4 M, *Into* 1.8× small** | packed `SQRTPD` (amd64) / packed NEON `FSQRT V.2D` (arm64) / intrinsic `FSQRTD` (others), off a non-`func`-pointer seam |
-| `Max` / `Min` | **FIXED — go wins 1.4× at 4 M, ~4× small** | builtin-`max` NaN-propagating oracle + 4-accumulator reducer + amd64 `MAXPD`+NaN-scan |
+| `Max` / `Min` | **FIXED — go wins 1.4× at 4 M, ~4× small** | builtin-`max` NaN-propagating oracle + 4-accumulator reducer + amd64 `MAXPD`+NaN-scan; NEON `FMAX`/`FMIN` on arm64 since v0.5.2 |
 | `MatMul` vs tuned BLAS (OpenBLAS) | **FIXED — parity at n=1024 (~0.99×, ~203 GFLOP/s); 0.97× at n=512** | by-lane FMLA micro-kernel (`FMLA Vd.2D,Vn.2D,Vm.D[i]` via `WORD`) closed the prior 0.76× gap. Only small n=256 trails (0.67×) on per-call overhead, not throughput — see above |
 | `Exp`, `Log`, `Log10` | **FIXED in v0.2.0** — own kernels | ports of Arm's optimized-routines, 2–3.5× `math.Exp`/`math.Log`; see the Exp and Log sections below |
 | `Sin`, `Cos`, `Tan`, … (through `Map`) | ~parity with NumPy | per element on one Zen 3 core: Go `math.Sin` 10.5 ns vs NumPy 10.4, `math.Cos` 11.1 vs 10.7, `math.Tanh` 16.5 vs 13.5 (2026-10-04); `Map`'s call per element adds 15–20% |
@@ -532,10 +532,33 @@ held to the n·ε·Σ|aᵢbᵢ| bound and exact on integer data.
 | MatVec 1024², 16 cores | 165 µs | 89 µs | — |
 
 NumPy's multi-threaded dot still wins on this part (24 µs). The go-ndarray
-version stops scaling at about 90 GB/s from 4 cores on. A likely cause, **not
-verified**, is cache affinity: OpenMP keeps each thread on one core, so it
+version stops scaling at about 90 GB/s from 4 cores on. ~~A likely cause, not
+verified, is cache affinity: OpenMP keeps each thread on one core, so it
 rereads its own chunk from its own L3 slice on every repetition, while
-goroutines move between core complexes.
+goroutines move between core complexes.~~
+
+**Measured 2026-10-07 (POWER9 bare metal, cfarm29, one thread per core):**
+the cache-affinity guess does not hold. The pattern is the fork/join of
+`parallelFor`, the same as the GEMM's (see *Parallel GEMM: start once, wait on
+work*). A 2^20 dot took 212 µs on 4 cores and 553 µs on 8, close to the 688 µs
+of one core. With 8 workers, only 2.4 CPUs were busy (`perf stat`). Three
+replacements for `parallelFor` were built and measured against it
+(8–10 interleaved rounds); **none is shipped**:
+
+| scheduler | 8 workers | 4 workers |
+|---|---|---|
+| shared queue of 4w blocks, caller works | Dot 1.55×, MatVec 2.0×, SumAxis(1) 2.0× | Sum 4 Mi 0.51×, Dot 0.55× |
+| owned runs (the static chunks) + stealing, caller works | Dot 2.18×, MatVec 2.0× | Sum 4 Mi 0.63×, Dot 0.68× |
+| owned runs + stealing, caller blocks, last block wakes it | mostly 1.0–1.4×, Exp 1.87× | Sum and AddInto at 256 Ki 0.51–0.55× |
+
+The second row keeps exactly the static chunks of the original and still
+loses at 4 workers, so affinity is not the cause. `perf stat` places the loss
+in scheduling instead: with the caller working, 2.1 of 4 CPUs were busy
+against 3.4 for the original. Its goroutines wait in the run queue of the
+caller's processor until another one is woken to steal them. With the caller
+blocking, the 4-worker case recovers, but the 8-worker gains go with it. The
+code is on the `experiment/parallel-for-scheduling` branch. Not yet solved:
+a scheduler that keeps both.
 
 ## Log and Log10 (2026-10-04)
 
@@ -731,18 +754,30 @@ stay dynamic, so its slower efficiency cores are handled as before.
   with the OS saving YMM state, else a pair of SSE2 6×4 `MULPD`+`ADDPD` tiles
   over the same packing. Generated by go-asmgen and validated per-arch in CI;
   the FMA kernel on real Cascade Lake and Zen 3 hosts.
-- **arm64 (NEON)** ships hand-vectorized `sum`, **packed `sqrt`** (`FSQRT V.2D`),
-  the **elementwise `add`/`sub`/`mul`** (via `VFMLA`/`VFMLS` against a `1.0`
-  vector — `add: b+a·1`, `sub: a−b·1`, `mul: 0+a·b`, each FMA exact so
-  bit-identical to the plain op; `div` stays on scalar `FDIVD`, no vector form),
-  and the **GEMM micro-kernel** (`gemmMicro4x8`: a 4×8 tile, 16 D2 accumulators,
-  using the **by-lane FMLA** `FMLA Vd.2D, Vn.2D, Vm.D[i]`). Three instructions are
-  emitted by raw `WORD` because Go's arm64 assembler cannot name them: the vector
-  `FSQRT V.2D`, and the indexed-element `FMLA …D[i]` (both verified vs `objdump`).
-  `max`/`min` use the four-accumulator builtin-`max`/`min` (lowering to
-  `FMAXD`/`FMIND`); the C-accumulate fold still uses `VFMLA` against `1.0` (no
-  plain vector FP add exists). All beat/parity NumPy via SIMD + multicore, no
-  `func`-pointer indirection.
+- **arm64 (NEON)** ships hand-vectorized `sum`, `dot`, **packed `sqrt`**
+  (`VFSQRT`), the **elementwise `add`/`sub`/`mul`/`div`** (`VFADD`/`VFSUB`/
+  `VFMUL`/`VFDIV`), **`max`/`min`** (`VFMAX`/`VFMIN`, four D2 accumulators from
+  8 elements up) and the **GEMM micro-kernel** (`gemmMicro4x8`: a 4×8 tile, 16 D2
+  accumulators, using the **by-lane FMLA** `FMLA Vd.2D, Vn.2D, Vm.D[i]`). Only
+  the indexed-element FMLA is still a raw `WORD` (16 of them): Go's assembler
+  does not name it yet (pending as CL 764224).
+
+  ⚠ **Corrected in v0.5.2.** Up to v0.5.1 this paragraph said Go's arm64
+  assembler had no vector `FSQRT`, no plain vector FP add, and no vector form for
+  `div`, so `add`/`sub`/`mul` went through an exact FMA against `1.0`, `sqrt` was
+  a `WORD`, and `div` and `max`/`min` stayed scalar. Go 1.27 (the floor since
+  v0.5.0) names all of them; the ISA always had them. Measured on one core,
+  interleaved against v0.5.1 (8 rounds, median; `Sum`, `SqrtInto` and `AddInto`,
+  unchanged, measured at 1.00× as the control):
+
+  | one core | `DivInto` | `Max` | `MulInto` |
+  |---|---|---|---|
+  | Apple M (local) | 1.34–2.19× | 1.87–2.24× | 1.21× at 1 Ki, 1.00× above |
+  | X-Gene (cfarm185) | 1.00× | 1.18–1.41× | 1.08–1.35× |
+
+  X-Gene's vector divide is no faster than its scalar one. The sum, dot and
+  GEMM C-accumulate still fold with an exact `VFMLA` against `1.0`, written
+  before `VFADD` existed in Go; the result is the same.
 - **ppc64le (VSX)**, since v0.2.3: sum, dot, sqrt, add/sub/mul/div and an 8×8
   GEMM micro-kernel; max/min stay scalar (the ISA's `xvmaxdp` NaN rule is not
   NumPy's). See the ppc64le section above.
diff --git a/docs/plan-ndarray.md b/docs/plan-ndarray.md
index 2b4b312..be7ad05 100644
--- a/docs/plan-ndarray.md
+++ b/docs/plan-ndarray.md
@@ -114,7 +114,7 @@ suite (plus the per-arch CI jobs already wired in `.github/workflows/ci.yml`).
     the 2026-06-23 pass.)
   - **SIMD via go-asmgen — DONE for amd64 (SSE2) + arm64 (NEON).** Hand-vectorized
     `sum` (4-accumulator), `sqrt` (packed `SQRTPD` / NEON `FSQRT V.2D`),
-    `max`/`min` (`MAXPD`/`MINPD` + NaN scan / four-accumulator builtin), the
+    `max`/`min` (`MAXPD`/`MINPD` + NaN scan / NEON `FMAX`/`FMIN` since v0.5.2), the
     elementwise `add`/`sub`/`mul`/`div`, and the **GEMM micro-kernel** (NEON 4×8
     by-lane FMLA on arm64, AVX2/FMA 6×8 with an SSE2 fallback on amd64),
     generated by go-asmgen and

From 86e384e43a51318faa47df90defe61371959cbde Mon Sep 17 00:00:00 2001
From: tannevaled <tannevaled@users.noreply.github.com>
Date: Wed, 7 Oct 2026 12:58:18 +0200
Subject: [PATCH 4/4] arm64: below 8 elements max/min use the scalar oracle;
 unrolled reducers off arm64 (dead there)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 internal/kernels/reduce_arm64.go    |  6 +++---
 internal/kernels/reduce_unrolled.go | 12 +++++-------
 2 files changed, 8 insertions(+), 10 deletions(-)

diff --git a/internal/kernels/reduce_arm64.go b/internal/kernels/reduce_arm64.go
index f2460aa..c3654fd 100644
--- a/internal/kernels/reduce_arm64.go
+++ b/internal/kernels/reduce_arm64.go
@@ -39,18 +39,18 @@ func sqrtSIMD(dst, src []float64) {
 
 // maxSIMD / minSIMD are the NaN-propagating extrema of a (non-empty): the NEON
 // FMAX/FMIN kernel with four D2 accumulators from eight elements up, the
-// four-chain scalar reducer below that. Bit-identical to the serial Max/Min
+// scalar oracle below that. Bit-identical to the serial Max/Min
 // oracle (see extremumKernel in asmgen/arm64/gen.go).
 func maxSIMD(a []float64) float64 {
 	if len(a) < 8 {
-		return maxUnrolled(a)
+		return Max(a)
 	}
 	return maxNEON(&a[0], len(a))
 }
 
 func minSIMD(a []float64) float64 {
 	if len(a) < 8 {
-		return minUnrolled(a)
+		return Min(a)
 	}
 	return minNEON(&a[0], len(a))
 }
diff --git a/internal/kernels/reduce_unrolled.go b/internal/kernels/reduce_unrolled.go
index bba28f9..2f2fbd8 100644
--- a/internal/kernels/reduce_unrolled.go
+++ b/internal/kernels/reduce_unrolled.go
@@ -1,14 +1,12 @@
-//go:build !amd64
+//go:build !amd64 && !arm64
 
 package kernels
 
 // maxUnrolled / minUnrolled are the four-accumulator NaN-propagating extreme
-// reducers used by the non-amd64 SIMD dispatch (arm64 and the four scalar
-// arches). amd64 has its own packed MAXPD/MINPD kernels (reduce_amd64.go ->
-// sum_amd64.s), so these Go versions would be dead code there — hence the
-// !amd64 build tag, which keeps every build's reachable statements covered
-// (the amd64 coverage job validates the .s kernel instead; arm64/generic cover
-// these). They keep FOUR independent builtin-max/min accumulators so the
+// reducers of the arches without a vector max/min kernel. amd64 (MAXPD/MINPD)
+// and arm64 (FMAX/FMIN, since v0.5.2) have their own, so these Go versions
+// would be dead code there, hence the build tag, which keeps every build's
+// reachable statements covered. They keep FOUR independent builtin-max/min accumulators so the
 // dependency chain is broken: the builtin max/min lowers to the hardware
 // FMAXD/FMIND on arm64, and four parallel chains hide its latency (~3.6x over a
 // single accumulator). max/min is associative for the NaN-propagating rule, so