Float32 results from the Workspace; cheaper range check in the AVX2 log

README.md

@@ -21,8 +21,9 @@
   amd64 with FMA (AVX2), ppc64le (VSX) and s390x (z14 vector), accumulating in float32 like NumPy's sgemm, and
   elsewhere on the float64 GEMM (accumulated in float64, rounded once).
   Its `Exp` and `Log` are **correctly rounded**, proven by an exhaustive test
-  over all 2³² float32 inputs, and on amd64 with AVX2 they run as fast as
-  NumPy's float32 SIMD, which is not correctly rounded
+  over all 2³² float32 inputs. On amd64 with AVX2, in a `Workspace`, `Exp`
+  is 12% faster than NumPy's float32 SIMD and `Log` 5% slower; NumPy's are
+  not correctly rounded
   ([docs/design-dtypes.md](docs/design-dtypes.md)).
 - **Creation** — `New`/`Zeros`/`Ones`/`Full`/`FromData`/`Arange`/`Linspace`/
   `Eye`/`Identity`.

SECURITY.md

@@ -66,6 +66,20 @@ standard library's known vulnerabilities; `govulncheck ./...` reports none.
 
 ## Audit log
 
+**2026-11-07 (v0.35.0).** `Workspace.alloc32` turns a slice of `(n+1)/2`
+float64 from the arena into `n` float32 with `unsafe.Slice`: the float32
+slice lies inside the float64 one (2·⌈n/2⌉·4 ≥ 4n bytes), its cap is n so an
+append reallocates instead of writing into the next result, and n = 0 takes
+the heap path so no pointer past a block is formed. The blocks hold no
+pointers, so the garbage collector keeps them alive through the interior
+pointer like any other slice. Unzeroed float32 results are requested only by
+code that writes every element; `TestWorkspace32MatchesHeap` fills the arena
+with NaN bits before each pass, so a site that read memory it had not
+written would differ from the heap result. In `log32CRAVX2`, lanes that are
+not positive normals now run through the arithmetic before being rejected;
+their table index is still `(m + 2^16) >> 17` ≤ 64, so the gathers stay in
+the tables.
+
 **2026-11-02 (v0.34.0).** The AVX2 float32 log kernel (`log32CRAVX2`) reads
 and writes at most `n` floats, a positive multiple of 8 chosen by the Go
 wrapper, and returns how many it wrote, a multiple of 8 no larger than `n`;

arith.go

@@ -72,7 +72,7 @@ func (a *Array) arith(b *Array, op binCode, k func(dst, x, y []float64)) (*Array
 	if err != nil {
 		return nil, err
 	}
-	d := makeStore(dt, prod(shape))
+	d := a.ws.store(dt, prod(shape), false) // arithStore writes every element
 	arithStore(op, d, a.operandStore(dt, shape), b.operandStore(dt, shape))
 	return a.result(d, append([]int(nil), shape...)), nil
 }
@@ -223,7 +223,7 @@ func (a *Array) mathUnary(code unaryCode, f64 func(*Array) *Array) *Array {
 		return f64(a.AsType(Float64))
 	}
 	src := a.operandStore(dt, a.shape)
-	d := makeStore(dt, a.Size())
+	d := a.ws.store(dt, a.Size(), false) // every element is written below
 	switch s := src.(type) {
 	case []float32:
 		if code == uSqrt {
@@ -302,7 +302,7 @@ func (a *Array) power(p float64) *Array {
 	src := a.contiguousStore()
 	switch s := src.(type) {
 	case []float32:
-		d := make([]float32, len(s))
+		d := a.ws.alloc32(len(s), false)
 		powFloat32(p, d, s)
 		return a.result(d, append([]int(nil), a.shape...))
 	case []complex64:
@@ -454,7 +454,7 @@ func (a *Array) clipTyped(lo, hi float64) (*Array, error) {
 	case a.dtype != Float32:
 		return a.AsType(Float64).Clip(lo, hi)
 	}
-	d := make([]float32, a.Size())
+	d := a.ws.alloc32(a.Size(), false)
 	clipOrdered(d, a.contiguousStore().([]float32), float32(lo), float32(hi))
 	return a.result(d, append([]int(nil), a.shape...)), nil
 }

benchmarks/dtype_bench_test.go

@@ -82,6 +82,33 @@ func BenchmarkDTypeExpLog(b *testing.B) {
 	}
 }
 
+// BenchmarkDTypeExpLogWS is BenchmarkDTypeExpLog for Float32 in a Workspace
+// reset every iteration: the kernels without the heap allocation and zeroing
+// of each result (4 µs at 16 Ki, 130 µs at 1 Mi on one Zen 3 core).
+func BenchmarkDTypeExpLogWS(b *testing.B) {
+	for _, n := range []int{1 << 14, 1 << 20} {
+		d := make([]float64, n)
+		for i := range d {
+			d[i] = 0.25 + float64(i%97)*0.25
+		}
+		x0, _ := nd.FromData(d, n)
+		ws := nd.NewWorkspace()
+		x := ws.Use(x0.AsType(nd.Float32))
+		for _, f := range []struct {
+			name string
+			op   func(*nd.Array) *nd.Array
+		}{{"Exp", (*nd.Array).Exp}, {"Log", (*nd.Array).Log}} {
+			b.Run(fmt.Sprintf("%s/float32/n=%d", f.name, n), func(b *testing.B) {
+				b.SetBytes(int64(n * 4))
+				for i := 0; i < b.N; i++ {
+					_ = f.op(x)
+					ws.Reset()
+				}
+			})
+		}
+	}
+}
+
 func BenchmarkDTypeMax(b *testing.B) {
 	for _, dt := range []nd.DType{nd.Float64, nd.Float32} {
 		for _, n := range []int{1 << 14, 1 << 20} {

docs/perf.md

@@ -1558,6 +1558,48 @@ The four gathers per block are not the bottleneck: replacing the `tl` gather
 by a register copy saves at most 8%. Doing better than NumPy needs a
 different algorithm, not a tweak of this one.
 
+#### Through the API, float32 results come from the Workspace (2026-11-07)
+
+The table above compares kernels. Through the API, `Log` on a Float32 array
+was still 1.25× NumPy at 16 Ki and 1.16× at 1 Mi (35.5 µs against 28.3 µs,
+2.07 ms against 1.79 ms; `BenchmarkDTypeExpLog`'s data, 0.25 + (i%97)·0.25,
+allocating on both sides). The difference is the result's allocation:
+`make` alone costs 4.1 µs at 16 Ki and 130 µs at 1 Mi on that core, because Go
+takes a large span from the heap and zeroes it, where NumPy reuses an
+unzeroed buffer. On this data the kernel hands no block back to `logCR`
+(`handed-blocks/op` in `BenchmarkLog32Kernels`).
+
+- **Float32 in a Workspace.** The Workspace used to give only float64
+  memory, so every Float32 result went to the heap even inside one.
+  `Workspace.alloc32` now carves float32 results from the same blocks (two
+  per float64 element, 64-byte aligned, cap = len). It is used for `Exp`,
+  `Log`, the other float32 ufuncs, arithmetic, `Power`, `Clip` and the
+  float32 GEMM (zeroed there, since the GEMM accumulates).
+  `TestWorkspace32MatchesHeap` holds every one of them to the heap result
+  bit for bit, also after the arena is filled with NaN bits; setting the
+  GEMM's allocation to unzeroed makes it fail.
+- **The kernel's range check** moved to the end of the block: one wrapping
+  add and one compare, ux + 0x7f800000 ≤ −0x01000001 as int32, select
+  exactly the positive normals. The separate test and branch before the
+  arithmetic are gone: −2.4% (`abbench`, 8 rounds, p ≤ 0.03 on three of four
+  sizes).
+
+| float32, one Zen 3 core, load ≤ 4.2 | 16 Ki | 1 Mi |
+|---|---:|---:|
+| `Exp`, allocating | 24.2 µs | 1.39 ms |
+| `Exp` in a Workspace | **19.4 µs** | **1.23 ms** |
+| NumPy 2.5.4, `np.exp(a)` | 22.0–22.3 µs | 1.37–1.47 ms |
+| `Log`, allocating | 34.5 µs | 2.06 ms |
+| `Log` in a Workspace | **29.9 µs** | **1.90 ms** |
+| NumPy 2.5.4, `np.log(a)` | 28.1–28.8 µs | 1.77–1.82 ms |
+
+Ours are medians of six runs in two interleaved rounds; NumPy's ranges are
+the medians of nine of every quiet round of the day (one 1 Mi `exp` round at
+2.47 ms and two 16 Ki `log` rounds at 31–34 µs are left out, which favours
+NumPy). In a Workspace, float32 `Exp` is about 12% faster than NumPy and
+`Log` about 5% slower; NumPy computes both in float32, at 2.52 and 3.83 ULP.
+Allocating, Go's zeroing of each result adds 4 µs at 16 Ki.
+
 **Open defect:** float64 `Log` is 1.23× slower than NumPy at 16 Ki on Zen 3
 (85 µs against 69 µs), while float64 `Exp` is 16% faster. The float32 kernels
 for the other architectures are still to write; there `Exp` and `Log` run

internal/kernels/asmgen/amd64/gen.go

@@ -508,14 +508,9 @@ func log32CRKernel() *emit.Function {
 	}
 	b.Raw("lloop:")
 	b.Raw("VMOVDQU (SI), Y0")
-	// Not a positive normal: v = ux - 2^23 outside [0, 0x7effffff].
-	b.Raw(fmt.Sprintf("VPSUBD %d(R8), Y0, Y4", sub))
-	b.Raw(fmt.Sprintf("VPCMPGTD %d(R8), Y4, Y5", maxv)) // v > 0x7effffff
-	b.Raw("VPXOR Y6, Y6, Y6")
-	b.Raw("VPCMPGTD Y4, Y6, Y6") // 0 > v
-	b.Raw("VPOR Y5, Y6, Y6")
-	b.Raw("VPTEST Y6, Y6")
-	b.Raw("JNZ lbail")
+	// Lanes that are not positive normals run through the arithmetic below
+	// on meaningless values (j stays at most 64, so the gathers stay inside
+	// the tables) and are rejected at the end of the block.
 	b.Raw(fmt.Sprintf("VPCMPEQD %d(R8), Y0, Y13", oneBits)) // x == 1
 	b.Raw(fmt.Sprintf("VPAND %d(R8), Y0, Y1", mmask))       // m
 	b.Raw(fmt.Sprintf("VPADDD %d(R8), Y1, Y2", jadd))
@@ -535,6 +530,11 @@ func log32CRKernel() *emit.Function {
 	b.Raw("VINSERTI128 $1, X3, Y15, Y15") // lb
 	b.Raw("VPCMPEQD Y15, Y14, Y15")       // the bounds agree
 	b.Raw("VPOR Y13, Y15, Y15")           // or x == 1
+	// Positive normals, ux in [2^23, 0x7f7fffff], are exactly the lanes with
+	// ux + 0x7f800000 (wrapping) at most -0x01000001 as int32.
+	b.Raw(fmt.Sprintf("VPADDD %d(R8), Y0, Y4", sub))
+	b.Raw(fmt.Sprintf("VPCMPGTD %d(R8), Y4, Y4", maxv)) // not a positive normal
+	b.Raw("VPANDN Y15, Y4, Y15")
 	b.Raw("VMOVMSKPS Y15, DX")
 	b.Raw("CMPL DX, $0xff")
 	b.Raw("JNE lbail")

internal/kernels/logcr_amd64.go

@@ -17,7 +17,8 @@ var logCRK = func() *logCRConsts {
 		k.expOr[i] = 1023 << 23
 	}
 	for i := range 8 {
-		k.sub[i], k.maxv[i], k.mmask[i] = 1<<23, 0x7effffff, 1<<23-1
+		// The kernel's positive-normal test: ux + sub (wrapping) <= maxv as int32.
+		k.sub[i], k.maxv[i], k.mmask[i] = 0x7f800000, 0xfeffffff, 1<<23-1
 		k.jadd[i], k.bias[i], k.oneBits[i] = 1<<16, 127, 0x3f800000
 	}
 	k.tr, k.tl = logCRTr, logCRTl

internal/kernels/logcr_test.go

@@ -1,6 +1,7 @@
 package kernels
 
 import (
+	"fmt"
 	"math"
 	"math/bits"
 	"runtime"
@@ -86,24 +87,46 @@ func logCRTakesSlow(x float32) bool {
 	return float32(r) != float32(r+0x1.f06p-33)
 }
 
-// BenchmarkLog32Kernels compares the float32 log kernels on 16 Ki inputs in
-// [0.25, 24.25), the range of benchmarks/BenchmarkDTypeExpLog.
+// BenchmarkLog32Kernels compares the float32 log kernels on two inputs in
+// [0.25, 24.25] at 16 Ki and 1 Mi: a ramp, and the multiples of 0.25 that
+// benchmarks/BenchmarkDTypeExpLog uses (0.25 + (i%97)*0.25). It reports the
+// blocks of 8 the AVX2 kernel hands back to logCR per call.
 func BenchmarkLog32Kernels(b *testing.B) {
-	src := make([]float32, 16<<10)
-	for i := range src {
-		src[i] = 0.25 + 24*float32(i)/float32(len(src))
-	}
-	dst := make([]float32, len(src))
-	for _, k := range []struct {
+	for _, data := range []struct {
 		name string
-		f    func(dst, src []float32)
-	}{{"via64", log32}, {"logCR", logCRSlice}, {"fast", log32Fast}} {
-		b.Run(k.name, func(b *testing.B) {
-			b.SetBytes(int64(4 * len(src)))
-			for b.Loop() {
-				k.f(dst, src)
+		at   func(i, n int) float32
+	}{
+		{"ramp", func(i, n int) float32 { return 0.25 + 24*float32(i)/float32(n) }},
+		{"steps", func(i, n int) float32 { return float32(0.25 + float64(i%97)*0.25) }},
+	} {
+		for _, n := range []int{16 << 10, 1 << 20} {
+			src := make([]float32, n)
+			for i := range src {
+				src[i] = data.at(i, n)
 			}
-		})
+			handed := 0
+			for i := 0; i+8 <= n; i += 8 {
+				for _, x := range src[i : i+8] {
+					if logCRTakesSlow(x) {
+						handed++
+						break
+					}
+				}
+			}
+			dst := make([]float32, n)
+			for _, k := range []struct {
+				name string
+				f    func(dst, src []float32)
+			}{{"via64", log32}, {"logCR", logCRSlice}, {"fast", log32Fast}} {
+				b.Run(fmt.Sprintf("%s/n=%d/%s", data.name, n, k.name), func(b *testing.B) {
+					b.SetBytes(int64(4 * n))
+					for b.Loop() {
+						k.f(dst, src)
+					}
+					b.ReportMetric(float64(handed), "handed-blocks/op")
+				})
+			}
+		}
 	}
 }
 

internal/kernels/sum_amd64.s

@@ -1791,13 +1791,6 @@ TEXT ·log32CRAVX2(SB), NOSPLIT, $0-40
 	XORQ AX, AX
 	lloop:
 	VMOVDQU (SI), Y0
-	VPSUBD 192(R8), Y0, Y4
-	VPCMPGTD 224(R8), Y4, Y5
-	VPXOR Y6, Y6, Y6
-	VPCMPGTD Y4, Y6, Y6
-	VPOR Y5, Y6, Y6
-	VPTEST Y6, Y6
-	JNZ lbail
 	VPCMPEQD 352(R8), Y0, Y13
 	VPAND 256(R8), Y0, Y1
 	VPADDD 288(R8), Y1, Y2
@@ -1855,6 +1848,9 @@ TEXT ·log32CRAVX2(SB), NOSPLIT, $0-40
 	VINSERTI128 $1, X3, Y15, Y15
 	VPCMPEQD Y15, Y14, Y15
 	VPOR Y13, Y15, Y15
+	VPADDD 192(R8), Y0, Y4
+	VPCMPGTD 224(R8), Y4, Y4
+	VPANDN Y15, Y4, Y15
 	VMOVMSKPS Y15, DX
 	CMPL DX, $0xff
 	JNE lbail

linalg.go

@@ -85,7 +85,7 @@ var useGemm32 = kernels.HasGemm32()
 // strides, like matView; another dtype is converted first.
 func matmulFloat32(a, b *Array) *Array {
 	m, k, n := a.shape[0], a.shape[1], b.shape[1]
-	dst := make([]float32, m*n) // MatMul32StridedP accumulates into dst
+	dst := a.ws.alloc32(m*n, true) // MatMul32StridedP accumulates into dst
 	kernels.MatMul32StridedP(dst, a.matView32(), b.matView32(), m, k, n)
 	return a.result(dst, []int{m, n})
 }

ufunc.go

@@ -88,12 +88,11 @@ func (a *Array) Log() *Array {
 	return a.unary(kernels.LogP)
 }
 
-// unary32 applies a float32 kernel to a Float32 array: the float64 kernels of
-// Exp and Log, widened and rounded once (kernels.Exp32P, Log32P), at the
-// float64 kernels' speed instead of math.Exp's and math.Log's.
+// unary32 applies a float32 kernel that writes every element (kernels.Exp32P,
+// Log32P) to a Float32 array; in a Workspace the result is not zeroed first.
 func (a *Array) unary32(f func(dst, src []float32)) *Array {
 	src := a.contiguousStore().([]float32)
-	dst := make([]float32, len(src))
+	dst := a.ws.alloc32(len(src), false)
 	f(dst, src)
 	return a.result(dst, append([]int(nil), a.shape...))
 }

workspace.go

@@ -1,6 +1,9 @@
 package ndarray
 
-import "sync"
+import (
+	"sync"
+	"unsafe"
+)
 
 // Workspace is an arena for the results of array operations, for loops that
 // compute the same shapes over and over. Bind an input with Use; every result
@@ -22,7 +25,9 @@ import "sync"
 // cycles, no zeroing of memory the kernel overwrites anyway, and the same
 // cache-warm bytes every time. This is the workspace pattern of ND4J, the
 // Java ndarray library, and the scoped counterpart of NumPy's reference-count
-// free-list, which a garbage-collected language does not have. The *Into
+// free-list, which a garbage-collected language does not have. Float64 and
+// Float32 results come from the workspace; results of the other dtypes still
+// come from the heap. The *Into
 // methods (AddInto, SqrtInto, ...) remain the way to reuse one named buffer.
 //
 // ⚠ Reset invalidates every array allocated from the workspace since the
@@ -107,6 +112,27 @@ func (w *Workspace) alloc(n int, zero bool) []float64 {
 // when it is bound to one, else from the heap (always zeroed there).
 func (a *Array) alloc(n int, zero bool) []float64 { return a.ws.alloc(n, zero) }
 
+// alloc32 is alloc for Float32 storage: n float32 carved from the same
+// blocks, two per float64 element, so Reset recycles a Float32 result like
+// any other. It keeps alloc's 64-byte start and has cap n.
+func (w *Workspace) alloc32(n int, zero bool) []float32 {
+	if w == nil || n == 0 {
+		return make([]float32, n)
+	}
+	s := w.alloc((n+1)/2, zero)
+	return unsafe.Slice((*float32)(unsafe.Pointer(unsafe.SliceData(s))), n)
+}
+
+// store returns storage for n elements of a non-Float64 dtype dt: Float32
+// from w (zeroed only if zero is set, so a caller passes false only when it
+// writes every element), the other dtypes zeroed from the heap.
+func (w *Workspace) store(dt DType, n int, zero bool) any {
+	if dt == Float32 {
+		return w.alloc32(n, zero)
+	}
+	return makeStore(dt, n)
+}
+
 // Detach returns a heap copy of a that is not bound to any workspace: the way
 // a result computed in a workspace outlives its Reset.
 func (a *Array) Detach() *Array {