docs: correct overstated claims, cite the Go issues, state measured results

README.md

@@ -57,44 +57,67 @@ from 1.18 ms to 0.39 ms (7.5× NumPy). ⚠ An array allocated from the workspace
 is invalid after `Reset`; `Detach` is how a result leaves it. The `*Into`
 methods (`AddInto`, `SqrtInto`, …) remain the way to reuse one named buffer.
 
-The numeric inner loops are kept behind a narrow kernel API. Behind it,
-**large** elementwise ops, reductions and matmul run **multicore** (across
-`GOMAXPROCS`); the sum reduction uses a **go-asmgen SIMD kernel** (NEON on
-arm64, SSE2 on amd64); and `MatMul` is a **panel-packed, cache-blocked GEMM**
-with a go-asmgen SIMD-FMA micro-kernel (NEON 4×8 on arm64; AVX2/FMA 6×8 on
-amd64, chosen at run time, with an SSE2 fallback) —
-the OpenBLAS/BLIS structure. So go-ndarray **beats single-threaded NumPy** on
-its core ops on large arrays (Add/Mul ~2×, Sum ~2.4×), and `MatMul` reaches
-**tuned-BLAS parity at 1024²**: ≈0.99× vs multi-threaded OpenBLAS 0.3.29 on an
-arm64 VM (~203 GFLOP/s, up from an earlier ~156 GFLOP/s / 0.76× before the
-by-lane-FMLA micro-kernel fix — see **[docs/perf.md](docs/perf.md)**), and
-≈1.00× of single-threaded vecLib on an Apple M4 Max (~373 GFLOP/s), where it
-also decisively **beats the pure-Go `gonum` 4–10×** at every size — see
-**[BENCHMARKS.md](BENCHMARKS.md)**. Where a faster reference exists —
-`Sqrt`/`Max` small sizes, and tuned BLAS at small matmul sizes — it says so.
-It is a **standalone, reusable** module and the cgo-free ndarray backend behind
+The numeric inner loops are kept behind a narrow kernel API. Behind it, large
+elementwise ops, reductions and products run **multicore** (across
+`GOMAXPROCS`). On **amd64 and arm64** the hot loops are
+[go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels: sum,
+max/min, add/sub/mul/div, sqrt and the dot product (SSE2 or AVX2/FMA on amd64,
+NEON on 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), the
+OpenBLAS/BLIS structure. The other four 64-bit targets (riscv64, loong64,
+ppc64le, s390x) and the 32-bit ones run the same pure-Go code those kernels
+are tested against. `Exp` and `Log` are ports of Arm's optimized-routines
+(≤ 0.51 ULP measured, against up to 0.88 and 0.72 for Go's
+`math.Exp`/`math.Log` on arm64),
+and they are correct where amd64's `math.Exp` returns +Inf
+([golang/go#81995](https://github.com/golang/go/issues/81995)) and `math.Log`
+is wrong on subnormals
+([golang/go#56600](https://github.com/golang/go/issues/56600)).
+
+**Measured against NumPy 2.5.3 (OpenBLAS 0.3.34) on an AMD Zen 3, 16 cores**
+(full tables in **[docs/perf.md](docs/perf.md)**):
+
+- **Faster:** whole-array reductions (`Sum`/`Mean`/`Max` at 4 Mi: 2.3–4.5×),
+  row reductions (`SumAxis(1)` 2.2×, `MaxAxis(1)` 1.7×), `Exp` from 256 Ki
+  elements (2.4–5.6×), the `*Into` forms of `Add`/`Mul`/`Sqrt` (faster or at
+  parity at every size), and `MatMul` against single-threaded OpenBLAS
+  (3.7–5.9× from 512²); in a `Workspace`, `x + y` and `sqrt(x*y + x)` match or
+  beat NumPy at every size measured.
+- **At parity:** `SumAxis(0)`, `Inner` against 16-thread OpenBLAS, `Log` on
+  one core.
+- **Slower:** `MatMul` against 16-thread OpenBLAS (0.55× at 1024², 0.3× at
+  256² and below); `Dot`/mat·vec against NumPy's threaded BLAS (0.16×); `Exp`
+  below 256 Ki elements (0.82×); the *allocating* forms of elementwise ops up
+  to 256 Ki elements, `Concatenate`/`Stack` and slice copies outside a
+  `Workspace` (0.2–0.5×: Go's allocator and garbage collector, which
+  `Workspace` largely removes).
+
+On Apple silicon the GEMM reached parity with tuned BLAS at 1024² (OpenBLAS in
+an arm64 VM, single-threaded vecLib on an M4 Max), and every product beats the
+pure-Go `gonum` 4–10× (**[BENCHMARKS.md](BENCHMARKS.md)**). It is a
+**standalone, reusable** module and the cgo-free ndarray backend behind
 [go-embedded-ruby](https://github.com/go-embedded-ruby/ruby)'s `NDArray` class.
 
-> ⚠️ **Status: float64 NumPy parity for the core surface, and faster than NumPy
-> on it.** Creation, slicing/views, broadcasting elementwise + ufuncs,
-> reductions (incl. arg/cumulative/clip/where), manipulation, and linear algebra
-> are complete, **100%-covered**, six-arch CI-green, and differentially checked
-> against numpy. The hot paths are **multicore + SIMD** and beat single-threaded
-> NumPy on large arrays (**[docs/perf.md](docs/perf.md)**). See
-> **[docs/plan-ndarray.md](docs/plan-ndarray.md)** for the roadmap (more SIMD
-> kernels, more dtypes). The Ruby binding has shipped — go-embedded-ruby's
-> `NDArray` class binds this module.
+> ⚠️ **Status: float64 NumPy parity for the core surface.** Creation,
+> slicing/views, broadcasting elementwise + ufuncs, reductions (incl.
+> arg/cumulative/clip/where), manipulation, linear algebra and `Workspace` are
+> complete, **100%-covered**, and differentially checked against NumPy. CI runs
+> the suite on amd64, arm64, 386, and riscv64/loong64/ppc64le/s390x/arm under
+> qemu, on Linux, macOS and Windows, and compiles it for every `GOOS/GOARCH`
+> pair Go supports; v0.1.0 and v0.2.0 were also run on real amd64, arm64,
+> ppc64le, riscv64 and loong64 hardware. See **[docs/plan-ndarray.md](docs/plan-ndarray.md)**
+> for the roadmap (more dtypes, more SIMD targets).
 
 ## Why this module?
 
 [gonum](https://www.gonum.org/) is matrix-centric and its optimized assembly is
 **amd64-only**. More broadly, **Ruby has no cgo-free ndarray** (`Numo::NArray`,
 `NMatrix` are C extensions). A pure-Go core whose kernels are generated for every
 arch is therefore a durable foundation. The numeric loops live in
-`internal/kernels`; **Phase 1** replaces them with
-[go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels across all six
-64-bit Go targets (amd64, arm64, riscv64, loong64, ppc64le, s390x), selected at
-runtime, behind the same API and the same tests.
+`internal/kernels`, each with a pure-Go reference implementation;
+[go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels replace them
+on amd64 and arm64 today, behind the same API and held to the same tests, and
+the other 64-bit targets are next.
 
 ## Example
 

docs/perf.md

@@ -340,7 +340,7 @@ writes every element (elementwise, ufuncs, copies, broadcasts, concatenation,
 `Where`, `Clip`, scans, `Take`, `Outer`, mat·vec) skip the zeroing. Those that
 accumulate (GEMM, vec·mat, axis reductions) clear their slot.
 
-Zen 3, 16 cores pinned, NumPy 2.5.3 single-threaded, median of two runs:
+Zen 3, 16 cores pinned, NumPy 2.5.3 single-threaded, best of two runs:
 
 | op | n | heap | **Workspace** | `*Into` | NumPy | Workspace vs NumPy |
 |----|--:|--:|--:|--:|--:|:--:|
@@ -381,14 +381,17 @@ hits both alike.
 Four changes, each measured on its own before the next:
 
 1. **AVX2/FMA 6×8 micro-kernel** (the BLIS Haswell shape). In L1 it runs at
-   43 GFLOP/s on one Cascade Lake core, the FMA ports' peak at that clock; the
-   SSE2 4×4 it replaces had no FMA and half the vector width.
+   43 GFLOP/s on one Cascade Lake core (16 flops per cycle at about 2.7 GHz,
+   if the virtual machine's unknown clock is near that); the SSE2 4×4 it
+   replaces had no FMA and half the vector width.
 2. **Edge tiles through the micro-kernel.** With MR=6, every 256-row MC block
-   ends in a 4-row edge panel, and the scalar edge loop took **40%** of a 512²
-   product (profiled). A partial tile now runs the kernel into a scratch tile.
+   ends in a 4-row edge panel, and the scalar edge loop took **40%** of one
+   profile covering serial 64² and 512² products. A partial tile now runs the
+   kernel into a scratch tile.
 3. **Prefetch** of the C tile at kernel entry and of the A panel eight steps
-   ahead: serial 768² went from 27.5 to 34.7 GFLOP/s (+26%); 90% of the time was
-   already in the kernel, stalled on memory. A sweep of MC (48–252) and KC
+   ahead: serial 768² went from 27.5 to 34.7 GFLOP/s (+26%). 90% of the time was
+   already in the kernel, which ran at 28 GFLOP/s there against 43 with its
+   panels in L1, so it was waiting on memory. A sweep of MC (48–252) and KC
    (128–512) moved nothing (±3%), so the blocking is unchanged.
 4. **B packed by all workers.** MatMulP packs the shared B block once per
    (jc, pc) round, and one goroutine did it while the others waited: +30–70%
@@ -408,9 +411,11 @@ below 256²). Where the parallel loss is, measured on the Zen 3:
   kernel at 15–17 GFLOP/s per core. Per-core efficiency already drops from 40
   to ~30 GFLOP/s at 2 cores.
 - It is not scheduling overhead: 83% of the parallel profile is the
-  micro-kernel itself, which stalls on the shared packed B panel. On this
-  part the last-level cache is split between core complexes, so the panel is
-  partly remote.
+  micro-kernel itself, so the kernel runs slower when cores share work. Why is
+  **not established**. The shared packed B panel is the suspect (on this part
+  the last-level cache is split between core complexes, so the panel is partly
+  remote), but the drop already shows at 2 cores of one complex, which that
+  explanation does not cover.
 - Taller row bands (more reuse of each B micro-panel) do not help: interleaved
   over 5 runs they were worse at 1024² and equal at 2048².
 - **A 2-D tiling of C with private packing** (no barrier, nothing shared) was
@@ -420,9 +425,9 @@ below 256²). Where the parallel loss is, measured on the Zen 3:
   rows·cols·k multiply-adds, and a strided pack costs far more per element than
   an FMA, so tiles under ~256×256 lose.
 
-The next lever is therefore a shared B panel per core complex (or per cache
-domain), not a different split of C, then an AVX-512 kernel for the hosts that
-have it.
+The next step is to find where the kernel stalls at 2 cores (hardware
+counters: cache misses per FMA, serial against parallel) before choosing a
+lever; then an AVX-512 kernel for the hosts that have it.
 
 Correctness was checked on real hardware for five of the six 64-bit targets:
 amd64 (Zen 3, FMA path), arm64 (cfarm185), ppc64le (cfarm120), riscv64 (cfarm94)
@@ -482,18 +487,23 @@ time.
 - **Accuracy:** worst 0.504 ULP over 80 000 inputs across the whole finite
   range, the region near 0, the overflow edge and the subnormal range,
   measured against a 300-bit reference (`math/big`, ln 2 from its atanh
-  series), against Arm's documented 0.509. `math.Exp` measured 0.879 ULP on
-  the same inputs on arm64. The table was not trusted as copied:
+  series), against Arm's documented 0.509; the committed test samples 20 000
+  of them and finds 0.501. `math.Exp` measured 0.879 ULP on the 80 000 inputs
+  on arm64 (0.841 on the test's 20 000). The table was not trusted as copied:
   `TestExpTable` re-derives all 256 entries from 2^(k/128) at 300 bits.
-- **A Go bug it removes:** on amd64, Go 1.26's `math.Exp` returns **+Inf for
-  x ≥ 1023.5·ln 2 ≈ 709.436**, although the result is finite up to
-  ln(MaxFloat64) ≈ 709.7827 (NumPy: `exp(709.5)` = 1.3549863193146328e+308).
-  Its assembly rounds k = x·log2(e) to 1024 there, and the biased-exponent
-  check treats that as overflow although the reduced factor is below 1. Seen on
-  real Zen 3 hardware, not only Rosetta. go-ndarray's `Exp` inherited it until
-  now; `TestExpTopOfRange` pins the fix.
-- Also just above ln(2^-1075), exp now rounds up to the smallest subnormal as
-  it should, where `math.Exp` returns 0.
+- **A Go bug it removes**, reported as
+  [golang/go#81995](https://github.com/golang/go/issues/81995): on amd64,
+  `math.Exp` returns **+Inf for 709.436139303104 ≤ x < 709.782712893384**,
+  although the result is finite there (NumPy, glibc and mpmath:
+  `exp(709.5)` = 1.3549863193146328e+308). Its assembly rounds k = x·log2(e) to
+  1024 from x = 1023.5·ln 2, and the biased-exponent check treats that as
+  overflow although the reduced factor is below 1. Reproduced with Go 1.26.8 and
+  1.27.1 on AMD Zen 3, Intel Cascade Lake and Intel Haswell, with and without
+  the FMA path; the assembly is unchanged on master. go-ndarray's `Exp`
+  inherited it until v0.2.0; `TestExpTopOfRange` pins the fix.
+- Just above ln(2^-1075), exp rounds up to the smallest subnormal, as glibc
+  and NumPy do. Go's arm64 `math.Exp` returns 0 there (amd64 does not); that is
+  one unit of the last subnormal place, not a contract violation.
 
 Speed, kernel only, 16 Ki elements: M4 30 µs (`math.Exp` 65 µs); one Zen 3 core
 62 µs (`math.Exp` 215 µs, NumPy single-threaded 78 µs). Whole `Exp` on Zen 3,
@@ -517,9 +527,12 @@ held to the n·ε·Σ|aᵢbᵢ| bound and exact on integer data.
 | Dot 2^20, 16 cores | 159 µs | 99 µs | — |
 | MatVec 1024², 16 cores | 165 µs | 89 µs | — |
 
-NumPy's multi-threaded dot still wins on this part (24 µs): OpenMP pins its
-threads, so each one rereads its own chunk from its own L3 slice on every
-repetition, while goroutines move between core complexes.
+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
+rereads its own chunk from its own L3 slice on every repetition, while
+goroutines move between core complexes.
+
 ## Log and Log10 (2026-10-04)
 
 `Log` went through `Map(math.Log)`: 11.6 ns per element on a Zen 3 core
@@ -532,17 +545,20 @@ inputs within 2⁻⁴ of 1 take a degree-12 polynomial. z/c − 1 is one
 needs is not used. `Log10` is log(x)·(1/ln10), the formula of `math.Log10`, on
 this log.
 
-- **Accuracy:** worst 0.508 ULP over 10 000 inputs (whole range, near 1,
-  subnormals), against a 300-bit reference (log m = 2·atanh((m−1)/(m+1)));
-  `math.Log` measured 0.727 on arm64. `TestLogTable` checks each of the 128
+- **Accuracy:** worst 0.508 ULP over the test's 10 000 inputs (whole range,
+  near 1, subnormals), against a 300-bit reference (log m =
+  2·atanh((m−1)/(m+1))); `math.Log` measured 0.724 on the same inputs on
+  arm64. `TestLogTable` checks each of the 128
   entries against the properties `log_data.c` documents (1/invc inside its
   subinterval, logc = log(c) within the rounding of invc, 0x1.8p9 + logc
   exact).
-- **A second Go bug it removes:** on amd64, Go 1.26's `math.Log` is wrong for
-  subnormal inputs. `math.Log(5e-324)` returns −709.09 instead of −744.44, and
-  `math.Log(1e-310)` −709.09 instead of −713.80 (NumPy: −744.44, −713.80); the
-  error measured 7·10¹³ ULP. Its assembly does not normalise them. Confirmed on
-  a real Zen 3. `Log` and `Log10` inherited it; `TestLogSubnormal` pins the fix.
+- **A second Go bug it removes**, known since 2022 as
+  [golang/go#56600](https://github.com/golang/go/issues/56600) (fix pending in
+  CL 448216): on amd64, `math.Log` is wrong for every subnormal input.
+  `math.Log(5e-324)` returns −709.09 instead of −744.44, and `math.Log(1e-310)`
+  −709.09 instead of −713.80 (NumPy, glibc: −744.44, −713.80). Its assembly
+  inlines `Frexp` with bit masks that assume a normal input. Reproduced with Go
+  1.26.8 and 1.27.1 on three x86 CPUs. `Log` and `Log10` inherited it; `TestLogSubnormal` pins the fix.
   (`Log2` already normalised through `math.Frexp` and was right.)
 
 Speed, 16 Ki elements, one core: M4 37 µs (`math.Log` 72), Zen 3 73 µs
@@ -572,12 +588,18 @@ Speed, 16 Ki elements, one core: M4 37 µs (`math.Log` 72), Zen 3 73 µs
   plain vector FP add exists). All beat/parity NumPy via SIMD + multicore, no
   `func`-pointer indirection.
 - The other four 64-bit Go targets — **riscv64, loong64, ppc64le, s390x** — keep
-  the validated scalar oracles (Go's loong64/ppc64le assemblers expose no
-  vector-double arithmetic; riscv64's V extension is optional), using the same
-  four-accumulator max/min, direct sqrt loop, and a **scalar 4×4 GEMM
-  micro-kernel** over the packed panels, and still get the **packing + cache
-  blocking + multicore** win, so they also beat single-threaded NumPy on large
-  arrays. (s390x additionally exercises the big-endian path in CI.)
+  the validated scalar oracles, using the same four-accumulator max/min, direct
+  sqrt loop, and a **scalar 4×4 GEMM micro-kernel** over the packed panels, and
+  still get the **packing + cache blocking + multicore** structure (they have
+  not been measured against NumPy). What the Go assembler offers them, checked
+  on Go 1.26.4 and 1.27.1 by assembling and disassembling (2026-10-04):
+  **ppc64le** has no vector-double arithmetic (no `XVADDDP`/`XVMADDADP`);
+  **loong64** has vector-double add/sub/mul/div (`VADDD` assembles to
+  `vfadd.d`, `VMULD` to `vfmul.d`, `XVADDD` to `xvfadd.d`) but no vector FMA;
+  **s390x** has them, FMA included (`VFADB`, `VFMADB`); **riscv64** has them
+  (`VFADDVV`, `VFMACCVV`), but the V extension is optional and needs a run-time
+  check. So loong64, s390x and riscv64 kernels are work not yet done, not a
+  toolchain wall. (s390x additionally exercises the big-endian path in CI.)
 
 All six are exercised in CI (native amd64/arm64 + qemu for the rest); each
 per-arch job regenerates the committed `.s`, fails if it is stale, vets

internal/kernels/exp_test.go

@@ -144,8 +144,8 @@ func TestExpSpecial(t *testing.T) {
 		}
 	}
 	// Against the reference: the subnormal range (just above ln(2^-1075) exp
-	// rounds up to the smallest subnormal; math.Exp gives 0), and the top of
-	// the finite range, where Go 1.26's amd64 math.Exp returns +Inf from
+	// rounds up to the smallest subnormal; arm64's math.Exp gives 0), and the
+	// top of the finite range, where amd64's math.Exp (golang/go#81995) returns +Inf from
 	// x >= 1023.5*ln2 ~ 709.436 (its k = round(x*log2e) reaches 1024 and the
 	// biased exponent check calls that overflow, though fr*2^1024 is finite).
 	for _, x := range []float64{-745.13, -745.1332, -740, -720.5, -708.5, -708.39,

internal/kernels/gemm_generic.go

@@ -4,9 +4,9 @@ package kernels
 
 // Generic GEMM micro-kernel for the four targets without a vector-double SIMD
 // kernel (loong64, ppc64le, riscv64, s390x): a portable scalar 4x4 register tile
-// over the packed panels. There is no hand-vectorized .s here — Go's loong64/
-// ppc64le assemblers expose no vector-double arithmetic and riscv64's V extension
-// is optional — but the kernel still gets the two structural wins of the packed
+// over the packed panels. There is no hand-vectorized .s here yet (ppc64le's Go
+// assembler has no vector-double arithmetic; loong64, s390x and riscv64 have it,
+// see docs/perf.md) — but the kernel still gets the two structural wins of the packed
 // GEMM: the panels are contiguous unit-stride (so the loads are conflict-free,
 // the L1-collision fix), and the cache-blocking loop nest plus the multicore
 // fan-out (MatMulP) keep it well ahead of an unblocked scalar loop and beat

internal/kernels/reduce_generic.go

@@ -10,9 +10,10 @@ package kernels
 // amd64 (SSE2) and arm64 (NEON sum + intrinsic scalar sqrt/max/min) ship kernels
 // (see reduce_amd64.go / reduce_arm64.go), so this file is the
 // loong64/ppc64le/riscv64/s390x fallback. For those four no hand-vectorized
-// kernel is shipped: Go's loong64/ppc64le assemblers expose no vector-double
-// arithmetic (the same wall go-fft documented), riscv64's V extension is
-// optional, and the per-arch qemu jobs still exercise this dispatch. The scalar
+// kernel is shipped yet. Go's ppc64le assembler has no vector-double
+// arithmetic; loong64 has vector add/mul (VADDD, VMULD) but no vector FMA,
+// s390x has both, and riscv64's V extension is optional (a run-time check);
+// see docs/perf.md. The per-arch qemu jobs still exercise this dispatch. The scalar
 // math.Sqrt/math.Max/math.Min still lower to the hardware sqrt/max where the
 // target has it, and the multicore path beats single-threaded numpy on large
 // arrays regardless.