diff --git a/README.md b/README.md
index f0a9571..004db1b 100644
--- a/README.md
+++ b/README.md
@@ -19,7 +19,7 @@ describes**:
 
 | Source | Published under |
 | --- | --- |
-| branch `main`, with `params.ndarray.version: v0.9.1` in `hugo.yaml` | `https://go-ndarray.github.io/docs/<MAJOR.MINOR>/`, here `/docs/0.9/` |
+| branch `main`, with `params.ndarray.version: v0.11.0` in `hugo.yaml` | `https://go-ndarray.github.io/docs/<MAJOR.MINOR>/`, here `/docs/0.11/` |
 | the newest version | also `https://go-ndarray.github.io/docs/latest/` |
 
 `https://go-ndarray.github.io/docs/` redirects to `latest/`. Each push to
@@ -41,22 +41,35 @@ same `versions.json` (mike's format).
 
 ## Mirrored pages
 
-Two pages are copies of files in `go-ndarray/ndarray`, kept verbatim so that
+Six pages are copies of files in `go-ndarray/ndarray`, kept verbatim so that
 re-syncing them is a copy, not an edit:
 
 | Page | Mirrors | What is ours |
 | --- | --- | --- |
+| `content/en/dtypes.md` | [`docs/design-dtypes.md`](https://github.com/go-ndarray/ndarray/blob/main/docs/design-dtypes.md), below its `#` title | the front matter |
+| `content/en/linalg.md` | [`docs/linalg.md`](https://github.com/go-ndarray/ndarray/blob/main/docs/linalg.md), below its `#` title | the front matter |
+| `content/en/random.md` | [`docs/random.md`](https://github.com/go-ndarray/ndarray/blob/main/docs/random.md), below its `#` title | the front matter |
+| `content/en/roadmap.md` | [`docs/plan-ndarray.md`](https://github.com/go-ndarray/ndarray/blob/main/docs/plan-ndarray.md), below its `#` title | the front matter |
 | `content/en/performance.md` | [`docs/perf.md`](https://github.com/go-ndarray/ndarray/blob/main/docs/perf.md), from the heading `## How a pure-Go library can beat NumPy` to the end | the front matter, and the introduction and *Where it stands* summary above that heading |
 | `content/en/security.md` | [`SECURITY.md`](https://github.com/go-ndarray/ndarray/blob/main/SECURITY.md), from `## Reporting` to the end | the front matter and the `> Mirrored from … at vX.Y.Z` line; repository-relative references are rewritten as GitHub links (`this repository` → `[the repository](https://github.com/go-ndarray/ndarray/security)`, `` `fuzz_test.go` `` → ``[`fuzz_test.go`](https://github.com/go-ndarray/ndarray/blob/main/fuzz_test.go)``, …) |
 
-To re-sync from a checkout of ndarray at tag `vX.Y.Z` (`$ND`):
+The first five are re-synced by one command, from a checkout of ndarray at tag
+`vX.Y.Z` (`$ND`); it refuses a checkout that is not at that tag:
 
 ```sh
-# performance: keep our part above the heading, replace the rest
-{ sed '/^## How a pure-Go library can beat NumPy/,$d' content/en/performance.md
-  sed -n '/^## How a pure-Go library can beat NumPy/,$p' "$ND/docs/perf.md"; } > perf.tmp
-mv perf.tmp content/en/performance.md
+scripts/sync-mirrors.sh "$ND" vX.Y.Z
+```
+
+It writes each page's `> Mirrored from … at vX.Y.Z` line, and rewrites the
+source's repository-relative links: a file this site mirrors becomes its page
+(`perf.md` → `performance.md`), a package directory (`../linalg`) becomes its
+pkg.go.dev page at the tag, any other file a GitHub link at the tag. A relative
+link it has no rule for makes it fail, naming the file and line.
 
+Security is re-synced by hand, because its references are prose (`fuzz_test.go`
+in backquotes), not links:
+
+```sh
 # security: keep our part above "## Reporting", replace the rest, then
 # rewrite the repository-relative references as GitHub links again
 { sed '/^## Reporting/,$d' content/en/security.md
diff --git a/content/en/_index.md b/content/en/_index.md
index ae15e2e..c31686f 100644
--- a/content/en/_index.md
+++ b/content/en/_index.md
@@ -4,13 +4,15 @@ linkTitle: "Home"
 type: docs
 cascade:
   type: docs
-description: "A pure-Go (no cgo) NumPy-style n-dimensional array for float64: views, broadcasting, ufuncs, reductions, linear algebra; multicore + SIMD."
+description: "A pure-Go (no cgo) NumPy-style n-dimensional array in 13 dtypes: views, broadcasting, ufuncs, reductions, numpy.linalg, numpy.fft, numpy.random; multicore + SIMD."
 ---
 
-**A pure-Go (no cgo) NumPy-style n-dimensional array** for `float64` — the
-`numpy` equivalent for Go. Creation routines, strided **views that share data**,
-broadcasting elementwise ops and ufuncs, reductions, manipulation and linear
-algebra, all with cgo disabled.
+**A pure-Go (no cgo) NumPy-style n-dimensional array** in **thirteen dtypes**
+— the `numpy` equivalent for Go. Creation routines, strided **views that share
+data**, broadcasting elementwise ops and ufuncs, reductions, manipulation,
+`numpy.linalg`, `numpy.fft` and `numpy.random` (with NumPy's numbers for the
+same seed), all with cgo disabled. `float64` is the default dtype and the one
+with SIMD kernels.
 
 Ruby has no cgo-free ndarray (`Numo::NArray`, `NMatrix` are C extensions) and
 `gonum`'s optimized assembly is amd64-only. go-ndarray pairs a portable pure-Go
@@ -34,14 +36,18 @@ WebAssembly, with nothing to install — build a pipeline step by step (creation
 
 | Area | Functions / methods |
 | --- | --- |
+| Dtypes ([page](dtypes.md)) | bool, int8…int64, uint8…uint64, float32/float64, complex64/complex128 with NumPy 2 promotion: `AsType`, `ZerosOf`, `FromSlice`, `Data[T]`/`Item[T]`, weak Go scalars with `Scalar` |
 | Creation | `New`, `Zeros`, `Ones`, `Full`, `Arange`, `Linspace`, `Eye`, `Identity`, `FromData` |
 | Views & shape | `Slice` (`All`/`R`/`Rng`/`From`/`To`/`Step`), `Reshape`, `Ravel`, `Transpose`, `Copy` |
-| Elementwise | `Add`/`Sub`/`Mul`/`Div` (+`*Scalar`, +`*Into`), `Map`, `Neg`, `Abs` |
+| Elementwise | `Add`/`Sub`/`Mul`/`Div` (+`*Scalar`, +`*Into`), `Map`, `Neg`, `Abs`, comparisons (`Greater`/`Equal`/… → Bool masks), `Maximum`/`Minimum`, `Logical*`, `Real`/`Imag`/`Conj`/`Angle` |
 | Ufuncs | `Sqrt`, `Exp`, `Log`/`Log2`/`Log10`, `Sin`/`Cos`/`Tan`, `Floor`/`Ceil`/`Round`, `Square`, `Power` |
 | Reductions | `Sum`, `Mean`, `Max`/`Min`, `Prod`, `ArgMax`/`ArgMin`, `CumSum`/`CumProd`, `Clip`, `Where` (+ per-axis) |
 | Indexing | `Slice` (basic indexing), `MaskSelect` (`a[mask]`), `Nonzero` (`flatnonzero`), `Take` (fancy indexing) |
 | Manipulation | `Flatten`, `ExpandDims`, `Squeeze`, `Concatenate`, `Stack`, `VStack`, `HStack` |
-| Linear algebra | `MatMul`, `Dot`, `Inner`, `Outer` |
+| Linear algebra | `MatMul`, `Dot`, `Inner`, `Outer`; package `linalg` ([page](linalg.md)): `Solve`, `Inv`, `Det`, `Slogdet`, `Cholesky`, `QR`, `SVD`, `Eigh`, `Eig`, `Lstsq`, `Pinv`, `MatrixRank`, `Norm`, `Cond`, `MatrixPower`, `MultiDot` |
+| Fourier transforms | package `fft`: `numpy.fft` on arrays (`FFT`/`RFFT`/`FFT2`/`FFTN`/… with `n`/`s`, `axes`, `norm`), over [go-fft](https://github.com/go-fft/fft) |
+| Random numbers | package `random` ([page](random.md)): PCG64, PCG64DXSM, MT19937, Philox, SFC64; `Generator` (`DefaultRNG`, 40 samplers) and the legacy `RandomState` (`Seed`, `Rand`) |
+| Notebooks & storage | `Repr` (prints like NumPy), `Display` (text and HTML), `MarshalBinary`/`UnmarshalBinary` (the same bytes on every architecture) |
 | Memory reuse | `NewWorkspace`, `Use`, `Reset`, `Detach` |
 
 ## Goroutines (v0.6.0)
@@ -87,13 +93,15 @@ On an AMD Zen 3 (16 cores), `x + y` on 1 024 elements goes from 2.6 µs to
 ## Performance & architectures
 
 Large elementwise ops, reductions and products run **multicore** (across
-`GOMAXPROCS`). On **amd64, arm64, ppc64le and loong64** the hot loops are
-go-asmgen 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 micro-kernel (NEON 4×8; AVX2/FMA 6×8 chosen at run time by a CPUID
-probe, SSE2 fallback; VSX 8×8; LASX 8×8). riscv64, s390x and the 32-bit
-targets run the same pure-Go code those kernels are tested against. `Exp` and `Log` are ports of Arm's
+`GOMAXPROCS`). On **amd64, arm64, ppc64le, loong64, riscv64 and s390x** the hot
+loops are go-asmgen SIMD kernels for `float64`: 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, RVV on riscv64 and the vector facility on s390x, the last three
+when the CPU has them), max/min too on amd64, arm64 and s390x (z14 and later),
+and a **panel-packed, cache-blocked GEMM** with an SIMD-FMA micro-kernel (NEON
+4×8; AVX2/FMA 6×8 chosen at run time by a CPUID probe, SSE2 fallback; VSX 8×8;
+LASX 8×8; RVV 6×8; s390x 4×8). The 32-bit targets run the same pure-Go code
+those kernels are tested against. `Exp` and `Log` are ports of Arm's
 optimized-routines (worst 0.504 and 0.508 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`
@@ -132,6 +140,9 @@ ppc64le, riscv64 and loong64 hardware.
 
 ## Where to go next
 
+- [Dtypes](dtypes.md) — the thirteen element types and NumPy 2 promotion.
+- [Linear algebra](linalg.md) — `numpy.linalg`: algorithms, accuracy, speed.
+- [Random numbers](random.md) — `numpy.random`, the same numbers for the same seed.
 - [Roadmap](roadmap.md) — the plan and what ships today.
 - [Performance](performance.md) — honest benchmarks versus NumPy / OpenBLAS.
 
diff --git a/content/en/dtypes.md b/content/en/dtypes.md
new file mode 100644
index 0000000..c9631c5
--- /dev/null
+++ b/content/en/dtypes.md
@@ -0,0 +1,183 @@
+---
+# Below the "Mirrored from" line, this page is go-ndarray/ndarray
+# docs/design-dtypes.md; scripts/sync-mirrors.sh re-syncs it.
+title: "Dtypes — thirteen element types, NumPy 2 promotion"
+linkTitle: "Dtypes"
+weight: 5
+description: "Why the dtype is a run-time property of one Array type, how arrays are stored, and where go-ndarray follows or departs from NumPy 2.x promotion."
+tags: [dtypes, numpy]
+---
+
+> Mirrored from [`docs/design-dtypes.md`](https://github.com/go-ndarray/ndarray/blob/v0.11.0/docs/design-dtypes.md) at v0.11.0.
+
+> Status: implemented (Phase 2). Supersedes decision 4 of
+> [plan-ndarray.md](roadmap.md) ("`float64` first").
+
+## Why now
+
+`Array` holds `float64` only. Three things the library is meant to do cannot be
+done on that alone:
+
+- **Spectra.** `go-fft` returns `[]complex128`. Without a complex dtype an FFT
+  result cannot be an array, so the two libraries do not compose.
+- **Masks.** Comparisons return `0/1` floats. They work for `MaskSelect` and
+  `Where`, but a mask cannot be told apart from data, and it costs 8 bytes per
+  element instead of 1.
+- **Memory-bound work.** `float32` halves the bytes moved, which is the whole
+  cost of the elementwise ops; images arrive as `uint8`; counts and indices are
+  integers.
+
+The type at the centre of the library is the one thing that gets more
+expensive to change every release, and the module is pre-1.0. That is why it
+goes first.
+
+## Decision: the dtype is a runtime property of the one `Array` type
+
+```go
+type DType uint8 // Bool, Int8…Int64, Uint8…Uint64, Float32, Float64, Complex64, Complex128
+
+a, _ := ndarray.Zeros(3, 4)          // Float64, as before
+b := a.AsType(ndarray.Complex128)
+c, _ := a.Add(b)                     // Complex128, by NumPy promotion
+m, _ := a.Greater(ndarray.Scalar(0)) // Bool
+```
+
+The alternative was a generic `Dense[T]`. It was rejected for four reasons:
+
+1. **Promotion cannot be written as a generic method.** `int32 / int32` is
+   `float64` and `int8 + float32` is `float32`: the result type depends on both
+   operands. A method on `Dense[T]` cannot return `Dense[U]` for a `U` chosen
+   from its argument, and methods cannot have their own type parameters. The
+   arithmetic would have to move to package functions, one per pair of types
+   in the general case, and NumPy's semantics would be lost or re-encoded.
+2. **The consumers are dynamic.** The notebook evaluates cells in an
+   interpreter, and the Ruby binding (`go-embedded-ruby`) is a dynamic language.
+   Both hold "an array" and ask it what it is; neither can carry a type
+   parameter.
+3. **Nothing that exists breaks.** `*Array`, its methods and the `float64` hot
+   path keep their names, signatures and speed. The Ruby binding uses method
+   expressions such as `(*nd.Array).Sqrt`; those would stop compiling on an
+   alias of an instantiated generic type.
+4. **Static typing is still available where it pays**: generic accessors
+   `Data[T]`, `FromSlice[T]` and `Item[T]` check the dtype once at the boundary
+   and hand back a `[]T`.
+
+## Storage
+
+```go
+type Array struct {
+	data    []float64 // Float64 storage — the existing hot path, unchanged
+	ext     any       // []T for every other dtype, nil for Float64
+	dtype   DType
+	shape   []int
+	strides []int
+	offset  int
+}
+```
+
+Keeping `data []float64` as a concrete field means every existing `float64`
+kernel, SIMD path and benchmark is untouched. The other dtypes are reached
+through one type assertion per operation, not per element. Shape, strides and
+views do not depend on the dtype, so slicing, reshape, transpose, broadcast and
+the manipulation functions are written once over a generic helper.
+
+Inner loops for the other dtypes are generic Go (`func add[T number](dst, x,
+y []T)`). The gc compiler stencils one copy per GC shape, and every numeric type
+has its own shape, so these are fully specialised loops. SIMD kernels come
+dtype by dtype after correctness, `float32` first.
+
+## Semantics: NumPy 2.x, checked against NumPy
+
+- **Array–array promotion** follows `numpy.result_type`. The full 13×13 table
+  is generated from NumPy 2.5.3 and committed as a test fixture
+  (`testdata/promotion.json`, embedded so the emulated CI arches can read it);
+  the implementation is tested against every cell.
+- **Go scalars are weak**, like Python scalars under NEP 50: `AddScalar(1.5)` on
+  a `float32` array stays `float32`, on an `int8` array becomes `float64`.
+  `ScalarOf(v, dtype)` makes a strong one.
+- **True division** of integers or bools gives `float64`.
+- **Reductions:** `Sum`/`Prod`/`CumSum` of bool and signed ints accumulate in
+  `int64`, of unsigned ints in `uint64`; `Mean` of integers is `float64`;
+  `Max`/`Min` keep the dtype; arg-reductions are `int64`.
+- **Float ufuncs on integers** (`Sqrt`, `Exp`, `Sin` …) give `float64` for
+  32- and 64-bit ints, `float32` for 16-bit ints.
+- **Comparisons return `Bool`.** `MaskSelect` and `Where` accept a `Bool` mask
+  and, as NumPy does, any other dtype by truthiness, so existing `0/1` masks
+  keep working.
+- **Integer overflow wraps**, as it does in both Go and NumPy.
+
+### Deliberate differences
+
+| NumPy | here | why |
+|---|---|---|
+| `float16` exists; `sqrt(int8)` is `float16` | no `float16`; `sqrt(int8)`, `sqrt(uint8)` and `sqrt(bool)` give `float32` | Go has no half-precision arithmetic; `float32` is the smallest type that holds the result exactly |
+| `round` is half-to-even | half-away-from-zero (unchanged) | documented since v0.1.0 |
+| `datetime64`, `str_`, `object`, structured dtypes | absent | out of scope for a numeric library |
+
+## Element access
+
+- `At(idx...) float64` keeps working for every real dtype, converting the
+  element. On a complex array it panics with a message that names `AtComplex`.
+  NumPy instead discards the imaginary part with a `ComplexWarning`; Go has no
+  warnings, and a silent loss is the defect that warning exists to flag.
+- `AtComplex(idx...) complex128`, and `Item[T](a, idx...) T` for exact typed
+  access.
+- `Set(v float64, idx...)` converts into the array's dtype.
+
+## Breaking changes
+
+| before | after | migration |
+|---|---|---|
+| comparisons return a `Float64` 0/1 array | they return `Bool` | `.AsType(Float64)` to get the old values; `MaskSelect` and `Where` need no change |
+| `Nonzero` returns `Float64` indices | `Int64`, as `np.flatnonzero` | `.AsType(Float64)`, or `Data[int64]` |
+| `ArgMaxAxis`/`ArgMinAxis` return `Float64` indices | `Int64`, as NumPy | same |
+| `MaskSelect` takes a 0/1 float mask | it takes any dtype (truthiness) | none |
+
+The only consumer in the fleet, `go-embedded-ruby/ruby` (`internal/vm/ndarray.go`),
+does not use comparisons; it is pinned to a June pseudo-version and will be
+moved forward in its own pull request.
+
+## What the implementation added
+
+- `Scalar(v)` makes the weak operand from any Go number; the `*Scalar`
+  methods' `float64` arguments are weak floats.
+- `SumAll`/`ProdAll` return a 0-d array in the accumulation dtype, exact for
+  every dtype; `Sum`/`Prod`/`Mean`/`Max`/`Min` keep returning `float64` and
+  panic on complex arrays (naming the method to use instead).
+- `LogicalNot`/`LogicalAnd`/`LogicalOr`/`LogicalXor`, `Any`, `All`, and
+  `Real`/`Imag`/`Conj`/`Angle`.
+- Float-to-integer casts saturate on every platform. NumPy leaves them
+  undefined: on arm64 it saturates to 64 bits and then wraps (`+Inf` to `int8`
+  is `-1`), on x86 `NaN` gives the minimum.
+- A `Workspace` serves `Float64` results only; the other dtypes allocate from
+  the heap for now.
+- Only `Float64` has SIMD kernels and multicore fan-out; the other dtypes run
+  tight generic loops (one per element type, generated by `gen_dtypes.go`).
+  `MatMul` on them is a plain triple loop.
+
+### Verification
+
+- `testdata/promotion.json`: every one of the 169 `result_type` cells, the
+  weak-scalar results, and the result dtype of the unary ufuncs and the
+  reductions, from NumPy 2.5.3. All match, with `float16` read as `float32`.
+- `testdata/values.json`: 47 operations on every dtype with edge values
+  (integer limits, NaN, ±Inf, −0) — 3,147 values. Correctly rounded
+  operations and every integer or bool result must match to the bit; the
+  others within 8 machine epsilons of NumPy's result dtype. A deliberately
+  broken integer `max` fails 48 of them.
+
+## Delivery
+
+1. **Core** — `DType`, storage, `AsType`, creation with a dtype
+   (`ZerosOf`/`OnesOf`/`FullOf`/`FromSlice[T]`), element access, and every
+   shape-only operation (reshape, transpose, slice, copy, flatten, expand,
+   squeeze, concatenate, stack, string form) for every dtype.
+2. **Arithmetic** — promotion, the four operators and the scalar forms,
+   comparisons returning `Bool`, `Where`/`MaskSelect`, reductions with their
+   accumulator types, unary ufuncs.
+3. **Complex** — `Real`/`Imag`/`Conj`/`Angle`, `Abs` of complex to float, and an
+   `fft` subpackage bridging `go-fft` on arrays (`fft.FFTN(a, axes...)`).
+4. **Speed** — `float32` SIMD kernels, then `MatMul` for `float32` and complex.
+
+Each step keeps the 100% coverage gate and the per-arch CI, and is checked
+against NumPy before it ships.
diff --git a/content/en/linalg.md b/content/en/linalg.md
new file mode 100644
index 0000000..bda0fb3
--- /dev/null
+++ b/content/en/linalg.md
@@ -0,0 +1,266 @@
+---
+# Below the "Mirrored from" line, this page is go-ndarray/ndarray
+# docs/linalg.md; scripts/sync-mirrors.sh re-syncs it.
+title: "linalg — numpy.linalg in pure Go"
+linkTitle: "Linear algebra"
+weight: 6
+description: "Solve, Inv, Det, Cholesky, QR, SVD, Eig and the rest of numpy.linalg: the algorithms, their accuracy against numpy 2.5.3, and their speed."
+tags: [linalg, numpy]
+---
+
+> Mirrored from [`docs/linalg.md`](https://github.com/go-ndarray/ndarray/blob/v0.11.0/docs/linalg.md) at v0.11.0.
+
+Package [`linalg`](https://pkg.go.dev/github.com/go-ndarray/ndarray@v0.11.0/linalg) gives ndarray arrays the functions of
+`numpy.linalg`, written in Go without cgo or assembly of its own: the
+matrix products it relies on are the ndarray package's packed SIMD GEMM.
+
+| numpy.linalg | linalg |
+|---|---|
+| `solve(a, b)` | `Solve(a, b)` |
+| `inv(a)` | `Inv(a)` |
+| `det(a)`, `slogdet(a)` | `Det(a)`, `Slogdet(a)` |
+| `cholesky(a, upper=…)` | `Cholesky(a, Lower \| Upper)` |
+| `qr(a, mode)` | `QR(a, ModeReduced \| ModeComplete \| ModeR \| ModeRaw)` |
+| `svd(a, full_matrices)`, `svdvals(a)` | `SVD(a, full)`, `SVDVals(a)` |
+| `eigh(a, UPLO)`, `eigvalsh(a, UPLO)` | `Eigh(a, uplo)`, `Eigvalsh(a, uplo)` |
+| `eig(a)`, `eigvals(a)` | `Eig(a)`, `Eigvals(a)` |
+| `lstsq(a, b, rcond)` | `Lstsq(a, b)`, `LstsqRcond(a, b, rcond)` → `LstsqResult{X, Residuals, Rank, S}` |
+| `pinv(a, rcond)` | `Pinv(a)`, `PinvRcond(a, rcond)` |
+| `matrix_rank(A, tol)` | `MatrixRank(a)`, `MatrixRankTol(a, tol)` |
+| `norm(x, ord, axis, keepdims)` | `Norm(x, ord, axes, keepdims)`; `ord` is `Ord{}` (None), `P(p)`, `Fro`, `Nuc` |
+| `cond(x, p)` | `Cond(x, p)` |
+| `matrix_power(a, n)` | `MatrixPower(a, n)` |
+| `multi_dot(arrays)` | `MultiDot(arrays...)` |
+| `LinAlgError` | `ErrLinAlg`, wrapped by `ErrSingular`, `ErrNotPosDef`, `ErrNoConvergence` |
+
+Shape problems give `ErrShape` (numpy raises `LinAlgError` for some of them
+and `ValueError` for others), bad options `ErrArgs`.
+
+## Semantics
+
+- **Stacks.** The last two axes are the matrix; leading axes are a stack,
+  processed one matrix per goroutine. `Solve` broadcasts the stack
+  dimensions of `a` and `b` as numpy's gufunc does; a 1-D `b` is one
+  vector, any other `b` a stack of matrices (numpy 2's rule).
+- **Dtypes.** As in numpy, everything is computed in double precision:
+  float64 and complex128 directly, integers and bools promoted to float64,
+  float32 and complex64 widened. Results have numpy's dtypes: single
+  precision in, single precision out (rounded from the double result, which
+  is what numpy does too — its linalg gufuncs also run in double); singular
+  values, Hermitian eigenvalues, `logabsdet` and norms are real; `Eig` is
+  always complex (numpy 2). `MatrixPower` with n ≥ 0 keeps an integer dtype.
+- **Conventions.** LAPACK's: Householder reflectors with β = −sign(α)‖x‖, so
+  R's diagonal can be negative as in numpy; ascending Hermitian eigenvalues;
+  decreasing singular values; `Eig` vectors with unit norm and their largest
+  component real, conjugate vectors for conjugate eigenvalues of a real
+  matrix (dgeev's normalization).
+- **Not unique.** Signs and phases of eigenvectors and singular vectors, and
+  any orthonormal basis of a repeated eigenvalue's space, may differ from
+  numpy's. Compare `A·v = λ·v`, `A = U·S·Vh`, not the vectors.
+
+## Algorithms, and why
+
+All are LAPACK's, re-expressed for row-major storage. Each is written once,
+generic over `float64 | complex128`; Go 1.27 stencils the two and folds the
+type switches of the helpers (`conj`, `abs`, …), so the real code pays
+nothing for the complex case (measured: same speed as a hand-specialised
+loop).
+
+| operation | algorithm | where the time goes |
+|---|---|---|
+| `Solve`, `Inv`, `Det`, `Slogdet` | LU with partial pivoting, recursive (getrf2): halves the columns, so the trailing update and the triangular solves are GEMMs | GEMM |
+| `Cholesky` | recursive potrf: `L21 = A21·L11⁻ᴴ`, `A22 −= L21·L21ᴴ` | GEMM |
+| `QR` | blocked Householder (geqrf/orgqr, panels of 32, compact WY `I − V·T·Vᴴ`) | GEMM |
+| `Eigh` | blocked tridiagonalization (latrd, panels of 32: half the flops in GEMM, half in matrix-vector products) + divide and conquer (stedc: Cuppen's split, deflation, secular equation, Gu–Eisenstat's recomputed vector) above order 25, implicit QL/QR (steqr) below | tridiagonalization |
+| `Eigvalsh` | the same reduction + QL/QR without vectors | tridiagonalization |
+| `SVD` | blocked Golub–Kahan bidiagonalization (labrd) + divide and conquer on the bidiagonal (bdsdc: Gu and Eisenstat's arrow matrix and secular equation) above order 25, implicit zero-shift/shifted QR (bdsqr, Demmel–Kahan) below; wide matrices factor Aᴴ | bidiagonalization |
+| `SVDVals`, `MatrixRank`, `Cond` (2-norm), `Norm` (2, nuc) | the same reduction + bdsqr without vectors | bidiagonalization |
+| `Eig`, `Eigvals` | balancing (gebal: permutations and powers of two), Hessenberg reduction (gehd2), Francis double-shift QR (dlahqr; single-shift zlahqr for complex input), eigenvectors by back-substitution on the Schur form with overflow scaling (trevc), back-transformed by one GEMM, gebak, dgeev's normalization | QR iteration |
+| `Lstsq` | SVD, `x = V·S⁺·Uᴴ·b` over singular values above `rcond·σ₁`; the residuals from the Householder QR of `a` (the tail of `Qᴴb`, free of the cancellation in `b − a·x`) | SVD |
+| `Pinv` | SVD, `V·S⁺·Uᴴ` | SVD |
+| `Norm`, `MatrixPower`, `MultiDot` | direct; 2-norms rescaled against overflow; repeated squaring in numpy's order; numpy's matrix-chain dynamic programme | — |
+
+**The SVD: why Golub–Kahan with QR, and divide and conquer.** The choice
+was between Golub–Kahan bidiagonalization followed by the implicit-shift QR
+iteration, and one-sided Jacobi. One-sided Jacobi is the more accurate in a
+strong sense — it computes the small singular values of a matrix whose
+columns are badly scaled to high relative accuracy (Demmel–Veselić) — but it
+needs 5–10 sweeps of O(n³) rotations, several times the work, and it is not
+what numpy computes. Golub–Kahan reduction is backward stable: the computed
+singular values are those of `A + E` with `‖E‖ ≈ ε‖A‖`, so each is within
+about `ε·σ₁` of the truth — numpy's (gesdd's) accuracy, exactly. The QR
+iteration (bdsqr, with Demmel and Kahan's zero shift) then finds the
+singular values of the bidiagonal to high relative accuracy, so it adds
+nothing to that error. It is used at the leaves of the divide and conquer
+and whenever the vectors are not wanted. With vectors, the QR iteration's
+plane rotations cost O(n³) of non-vectorized work (5.8 s at n = 1024 in the
+first version); divide and conquer, LAPACK's gesdd algorithm and the one
+numpy calls, puts that work into GEMMs (0.26 s now). Its singular vectors
+stay orthogonal because the updating vector is recomputed from the computed
+roots (Gu and Eisenstat), and the secular equation is solved with each
+distance `dᵢ − σ` formed from the nearest pole, in the squared form
+`(dᵢ − σ)(dᵢ + σ)`.
+
+**Hermitian eigenproblems** follow the same pattern: the tridiagonal QL/QR
+iteration (steqr) is used without vectors and for orders up to 25, divide
+and conquer (stedc) above, which took n = 1024 from 3.6 s to 0.13 s.
+
+**General eigenproblems** use the classical double-shift QR (dlahqr), not
+LAPACK's multishift QR with aggressive early deflation (dlaqr0–5). That is
+the largest gap to numpy, at large orders (below).
+
+## Accuracy and verification
+
+**Against numpy.** `linalg/testdata/cases.py` generates 1,040 cases with
+numpy 2.5.3 (`uv run --with numpy python linalg/testdata/cases.py >
+linalg/testdata/cases.json`; the JSON is embedded in the test binary, so the
+emulated CI architectures read it from anywhere). They cover every
+function; float64, float32, complex128, complex64, int32/int64 and bool
+inputs; 0×0, 1×1, tall, wide, stacked and broadcast shapes; Hilbert
+matrices (κ up to 10¹⁶), graded and rank-deficient matrices, exactly
+singular ones, non-symmetric matrices with complex eigenvalues, Jordan
+blocks and other defective matrices, NaN and Inf inputs; and every error
+numpy raises on them (124 cases). The test checks, for each case:
+
+- errors: an error of the same kind (`ErrSingular`, `ErrNotPosDef`,
+  `ErrNoConvergence`, `ErrLinAlg` for non-finite input to `Eig`,
+  `ErrShape`, `ErrArgs`);
+- dtype and shape of every result, exactly;
+- unique values (solutions, inverses, determinants, Cholesky factors,
+  singular values, eigenvalues, least-squares solutions and residuals,
+  pseudo-inverses, ranks, norms, condition numbers, powers, products)
+  elementwise, within `32·n·ε·κ` of the largest entry of the result, with
+  `κ` the condition number of the problem (1 for singular values, Hermitian
+  eigenvalues and norms; `cond(A)` for solves and inverses; the effective
+  condition over the kept singular values for `Lstsq` and `Pinv`, squared for
+  the residuals), plus one float32 rounding for single-precision results;
+- factors that are not unique through their identities: `‖QᴴQ − I‖`,
+  `‖A − QR‖/‖A‖`, R upper triangular; `‖UᴴU − I‖`, `‖VhVhᴴ − I‖`,
+  `‖A − U·S·Vh‖/‖A‖`; `‖VᴴV − I‖`, `‖AV − VΛ‖/‖A‖` with ascending
+  eigenvalues; for `Eig`, `‖Av − λv‖/‖A‖` per vector, unit norm, largest
+  component real, conjugate pairs exactly conjugate — all within
+  `32·n·ε` (`ε` of the result's dtype);
+- general eigenvalues matched to numpy's nearest unused one, within 100
+  times how far numpy's own eigenvalues move when each entry is perturbed by
+  10⁻¹⁴ of itself (an empirical condition number, recorded in the fixture),
+  so a Jordan block's eigenvalue, accurate only to `ε^(1/4)`, is judged by
+  what is achievable.
+
+914 cases are value-checked and 124 are errors; the 2 deliberate departures
+below are logged.
+
+**Beyond the fixture** (`properties_test.go`, `kernels_test.go`): the same
+identities at orders 100 and 131 (above the blocking and divide-and-conquer
+thresholds), real and complex, tall and wide; a stack of 16 matrices run in
+parallel equals the matrices run one by one; repeated eigenvalues and
+singular values, rank-deficient and zero matrices (the deflation paths);
+cyclic permutations, Grcar and Frank matrices, orthogonal matrices,
+high-order Jordan blocks and coupled rotation blocks (exceptional shifts,
+overflow rescaling in the back-substitution); balancing at the limits of the
+exponent range; every iteration forced to fail (`ErrNoConvergence` from
+each function); and the 2×2 kernels (laev2, las2, lasv2, lanv2, lartg,
+solve2) checked through their defining identities on the inputs that reach
+their special cases.
+
+**Negative controls.** Dropping the row interchange in the LU leaf (no
+partial pivoting) fails 64 fixture cases — 23 `solve`, 17 `inv`, 10 `det`,
+10 `slogdet`, 4 `matrix_power` — and the order-100 identity test. Loosening
+the secular-equation solver's stopping test by 10⁸ fails the four order-34
+`eigh` cases (eigenvalues off by ~10⁻¹¹ relative) and the zero-diagonal
+bidiagonal divide-and-conquer test (`‖B − USVᵀ‖` of 10⁻⁸).
+
+**Coverage** is 100% of statements under the CI gate (`go test -race
+-coverpkg=./... ./...`).
+
+## Speed against numpy
+
+Median wall time of `BenchmarkLinalg` against numpy 2.5.3 (`la.solve(a,
+ones)`, `la.inv`, `la.cholesky` of `A·Aᵀ + nI`, `la.qr` reduced,
+`la.svd(full_matrices=False)`, `la.eigh` of `A + Aᵀ`, `la.eig`), n×n
+standard normal matrices, Apple M-series (16 cores, arm64), three
+interleaved rounds. numpy's wheel uses **Apple Accelerate**. The host was
+shared and heavily loaded (load average 16–23 on 16 cores during the runs),
+so absolute numbers are noisy; the ratios held across rounds.
+
+| op | n | Go (ms) | numpy (ms) | Go / numpy |
+|---|---|---|---|---|
+| solve | 64 | 0.046 | 0.022 | 2.1 |
+| solve | 256 | 0.99 | 0.25 | 3.9 |
+| solve | 1024 | 21.2 | 10.9 | 1.9 |
+| inv | 64 | 0.11 | 0.029 | 3.7 |
+| inv | 256 | 2.19 | 0.58 | 3.8 |
+| inv | 1024 | 39.4 | 23.7 | 1.7 |
+| cholesky | 64 | 0.030 | 0.009 | 3.3 |
+| cholesky | 256 | 0.54 | 0.24 | 2.2 |
+| cholesky | 1024 | 10.9 | 6.6 | 1.7 |
+| qr | 64 | 0.15 | 0.083 | 1.8 |
+| qr | 256 | 2.32 | 1.57 | 1.5 |
+| qr | 1024 | 54.7 | 37.3 | 1.5 |
+| svd | 64 | 0.78 | 0.30 | 2.6 |
+| svd | 256 | 13.3 | 5.25 | 2.5 |
+| svd | 1024 | 263 | 123 | 2.1 |
+| eigh | 64 | 0.41 | 0.22 | 1.9 |
+| eigh | 256 | 6.54 | 3.12 | 2.1 |
+| eigh | 1024 | 130 | 78 | 1.7 |
+| eig | 64 | 1.07 | 1.06 | 1.0 |
+| eig | 256 | 55 | 19.4 | 2.8 |
+| eig | 1024 | 5060 | 712 | 7.1 |
+
+**numpy is faster on single matrices**, by 1.5–4× (7× for `eig` at 1024).
+Most of that is the matrix product underneath: on the same host the ndarray
+GEMM takes 18 ms for a 1024² product where Accelerate (which uses the
+matrix coprocessor) takes 3.8 ms, 0.37 ms against 0.10 ms at 256. The
+blocked algorithms here spend most of their time in that GEMM, so their gap
+to numpy is smaller than the GEMM's. The exceptions:
+
+- `eig` at large n: the classical double-shift QR is O(10n³) of scalar
+  work; LAPACK's multishift QR with aggressive early deflation moves it into
+  GEMMs. This is the main follow-up.
+- small orders (64): fixed costs (conversion of the input, allocation,
+  the packing of a small GEMM) weigh more than the arithmetic.
+
+**Stacks of small matrices run in parallel** — numpy's gufuncs loop over
+them on one core — and there this package is faster (median of 3, same
+host):
+
+| op | stack | Go (ms) | numpy (ms) | numpy / Go |
+|---|---|---|---|---|
+| inv | 10000 × 4×4 | 0.94 | 3.83 | 4.1 |
+| det | 10000 × 4×4 | 0.81 | 2.12 | 2.6 |
+| eigh | 10000 × 4×4 | 3.0 | 13.8 | 4.6 |
+| svd | 10000 × 4×4 | 5.3 | 29.5 | 5.6 |
+| eig | 10000 × 4×4 | 4.9 | 55.7 | 11 |
+| inv | 1000 × 32×32 | 3.4 | 8.5 | 2.5 |
+| det | 1000 × 32×32 | 1.6 | 4.0 | 2.5 |
+| eigh | 1000 × 32×32 | 21.3 | 50.4 | 2.4 |
+| svd | 1000 × 32×32 | 44.3 | 75.9 | 1.7 |
+| eig | 1000 × 32×32 | 29.7 | 232 | 7.8 |
+
+Reproduce: `go test -run XXX -bench 'Linalg|Stacked' ./linalg`, and the
+numpy side with `uv run --with numpy==2.5.3 python docs/bench_linalg_numpy.py`
+(the same sizes and distributions; the values differ, being from another
+generator).
+
+## Differences from numpy
+
+| numpy | here | why |
+|---|---|---|
+| `det` is `sign·exp(logabsdet)` | the product of the pivots, accumulated as mantissa and exponent | exact to rounding, no spurious overflow or underflow of intermediate products, and no `|logdet|·ε` error from `exp` |
+| `norm` of a vector with entries near 10²⁰⁰ overflows to `inf` | the 2-norms are rescaled (`1.414e200` for `[1e200, 1e200]`) | the true value is representable |
+| float32 norms accumulate in float32 | in float64, then rounded | more accurate; within a float32 rounding of numpy |
+| `cond(singular, 'nuc')` raises `LinAlgError` (the SVD of its NaN "inverse" fails) | `+Inf`, as for every other order | a singular matrix has infinite condition number |
+| `eigh` of a matrix with a NaN returns NaN eigenvalues | `ErrNoConvergence` | the iteration cannot proceed; numpy's `svd` and `eig` refuse NaN too |
+| `LinAlgError` for a wrong rank or a non-square matrix | `ErrShape` | a shape error, not a numerical failure |
+| `svd(…, hermitian=True)`, `pinv(…, hermitian=True)`, `matrix_rank(…, hermitian=True)`, `matrix_rank(rtol=)`, `pinv(rtol=)` | not provided (`PinvRcond(a, max(M,N)·ε)` is `rtol=None`) | follow-ups |
+| `tensorsolve`, `tensorinv`, `cross`, `outer`, `trace`, `diagonal`, `matmul`, `vecdot`, `matrix_norm`, `vector_norm`, `matrix_transpose` | not in this package (`ndarray` has `MatMul`, `Outer`, `Transpose`; `Norm` covers both norms) | follow-ups |
+| `qr` modes `'full'` and `'economic'` (deprecated) | not provided | deprecated in numpy |
+
+## What is left
+
+- Multishift QR with aggressive early deflation for `Eig` (the 7× gap at
+  n = 1024), and a blocked Hessenberg reduction.
+- A faster small-matrix path (n ≤ 64): unblocked kernels without the GEMM
+  packing, fewer copies.
+- The `hermitian` options and `rtol` arguments above, and the tensor
+  functions.
+- `Eigh` on NaN input returning NaNs like numpy, if that is preferred.
diff --git a/content/en/random.md b/content/en/random.md
new file mode 100644
index 0000000..497eb3f
--- /dev/null
+++ b/content/en/random.md
@@ -0,0 +1,230 @@
+---
+# Below the "Mirrored from" line, this page is go-ndarray/ndarray
+# docs/random.md; scripts/sync-mirrors.sh re-syncs it.
+title: "random — numpy.random, the same numbers"
+linkTitle: "Random numbers"
+weight: 7
+description: "numpy.random with NumPy's numbers for the same seed: the bit generators and samplers, how that is verified, and where a C library makes the last bit differ."
+tags: [random, numpy]
+---
+
+> Mirrored from [`docs/random.md`](https://github.com/go-ndarray/ndarray/blob/v0.11.0/docs/random.md) at v0.11.0.
+
+Package [`random`](https://pkg.go.dev/github.com/go-ndarray/ndarray@v0.11.0/random) is `numpy.random` over go-ndarray arrays. Its goal
+is migration: a notebook moved from NumPy to Go keeps its numbers for the same
+seed.
+
+```go
+rng := random.DefaultRNG(42)          // np.random.default_rng(42)
+x, _ := rng.StandardNormal(3, 4)      // rng.standard_normal((3, 4))
+k, _ := rng.Integers(0, 10, 5)        // rng.integers(0, 10, 5)
+
+random.Seed(0)                        // np.random.seed(0)
+u, _ := random.Rand(3)                // np.random.rand(3)
+```
+
+`Example`, `ExampleSeed` and `ExampleGenerator_ChoiceN` in
+[`example_test.go`](https://github.com/go-ndarray/ndarray/blob/v0.11.0/random/example_test.go) print NumPy's own numbers.
+
+## What is implemented
+
+| NumPy | here |
+|---|---|
+| `SeedSequence(entropy, spawn_key=, pool_size=)`, `.pool`, `.generate_state(n, uint32/uint64)`, `.spawn(n)` | `NewSeedSequence(v...)`, `NewSeedSequenceBig`, `NewSeedSequenceWith`, `Pool`, `GenerateState`, `GenerateState64`, `Spawn` |
+| `PCG64`, `PCG64DXSM`, `MT19937`, `Philox` (also `counter=`, `key=`), `SFC64` | `NewPCG64` … `NewSFC64`, `New…FromSeedSequence`, `NewPhiloxKey` |
+| `.jumped(n)`, `.advance(delta)`, `.random_raw(n)`, `.spawn(n)`, `.state` | `Jumped`, `Advance(*big.Int)`, `RandomRaw`, `Spawn`, `State`/`SetState` (typed structs) |
+| `default_rng(seed)`, `Generator(bitgen)`, `.spawn(n)` | `DefaultRNG(seed...)`, `NewGenerator`, `Spawn` |
+| `random(size, dtype)` | `Random(size...)`, `RandomOf(dt, size...)` |
+| `integers(low, high, size, dtype, endpoint)` | `Integers(low, high, size...)`, `IntegersOf(dt, low, high, endpoint, size...)`, `IntegersUint64` |
+| `standard_normal`, `normal`, `standard_exponential(method=)`, `exponential`, `uniform`, `standard_gamma(dtype=)`, `gamma` | `StandardNormal[Of]`, `Normal`, `StandardExponential[Of]`, `Exponential`, `Uniform`, `StandardGamma[Of]`, `Gamma` |
+| `beta`, `chisquare`, `f`, `noncentral_chisquare`, `noncentral_f`, `standard_t`, `standard_cauchy`, `lognormal`, `laplace`, `logistic`, `gumbel`, `weibull`, `pareto`, `power`, `rayleigh`, `wald`, `triangular`, `vonmises` | same names in Go case |
+| `binomial`, `negative_binomial`, `poisson`, `geometric`, `zipf`, `logseries` | same names, int64 results |
+| `choice(a, size, replace, p, axis, shuffle)` | `ChoiceN(n, opts, size...)`, `Choice(a, opts, size...)` |
+| `shuffle(x, axis)`, `permutation(n)`, `permutation(x, axis)`, `permuted(x, axis=)` | `Shuffle`, `Permutation`, `PermutationOf`, `Permuted`, `PermutedAxis` |
+| `bytes(n)`, `dirichlet(alpha, size)`, `multinomial(n, pvals, size)` | `Bytes`, `Dirichlet`, `Multinomial` |
+| `RandomState(seed)`, `.seed`, `.get_state`/`.set_state` | `NewRandomState(uint32)`, `NewRandomStateArray`, `NewRandomStateFrom(bitgen)`, `Seed`, `SeedArray`, `GetState`/`SetState` |
+| `rand`, `randn`, `randint`, `random_sample`, `normal`, `uniform`, `choice`, `shuffle`, `permutation`, `binomial`, `poisson`, `standard_normal`, `standard_exponential`, `exponential`, `standard_gamma`, `gamma`, `chisquare`, `lognormal`, `multinomial`, `bytes` | `RandomState` methods of the same names |
+| `np.random.seed`, `rand`, `randn`, `randint`, `random_sample`, `normal`, `uniform`, `choice`, `shuffle`, `permutation`, `binomial`, `poisson` | package functions over a shared `RandomState` behind a mutex |
+
+The size comes last; no size gives a 0-d array holding what NumPy returns as a
+scalar (`size=None`), which draws exactly the same values. `ChoiceOptions`'
+zero value is NumPy's defaults (`NoReplace`, `NoShuffle` are the negations).
+
+**Not implemented:** array-valued (broadcast) distribution parameters and
+integer bounds, `out=`, `hypergeometric`, `multivariate_normal`,
+`multivariate_hypergeometric`, the `RandomState` samplers not listed above
+(`beta`, `f`, `pareto` …: NumPy 1.16's versions differ from `Generator`'s and
+were not ported), `tomaxint`, `random_integers`, pickling, and `float16`.
+
+## Verification
+
+`random/testdata/cases.py` runs NumPy 2.5.3 and records:
+
+- **SeedSequence**: the pool, 9 uint32 and 5 uint64 state words, and the pools
+  and spawn keys of 4 children, for 10 entropies (0, 42, 2⁶⁴−1, `[1,2,3]`,
+  2¹²⁸+12345, a 10-word list, `[]`, pool size 8, spawn keys with a 2⁴⁰ entry,
+  a key after a 6-word entropy);
+- **bit generators**: for each of the five and six seeds (0, 1, 42, 2⁶⁴−1,
+  `[1,2,3]`, 2¹²⁸+12345), the seeded state, 100 raw outputs (700 for MT19937,
+  which crosses a regeneration), `jumped(1)`, `jumped(3)`, `advance(2¹⁰⁰+7)`,
+  `advance(-3)`, two spawned children; Philox from an explicit counter and key
+  and across a counter carry; legacy MT19937 seeding from 6 integer and array
+  seeds, and its jump;
+- **2,222 scripts, 2,492 calls, 81,806 values**: every Generator function on
+  PCG64 with every seed, every dtype and bound case of `integers` (9 dtypes,
+  open and closed, full ranges, ranges of 2³²−1, 2³², 3·2³⁰, 3·2⁶², empty
+  ranges); a subset on the four other bit generators; an interleaved script
+  per generator and seed (the buffered upper half of a 64-bit draw outlives a
+  call); `RandomState` on 6 seeds (including a 700-word array);
+- **50 long runs as SHA-256 digests** (4.49 million values): `random`,
+  `standard_normal` and `standard_exponential` (both dtypes, 200,000 each,
+  which reaches the ziggurat tails and wedges), `integers` (Lemire with
+  rejections, 8- and 16-bit buffering, bool), `binomial` (inversion and
+  BTPE), `poisson` (both methods), `choice` with and without replacement, and
+  `RandomState`'s `rand`, `randint`, `binomial`, `poisson`, `uniform`.
+
+The fixture is embedded, so CI's emulated targets check it too. It passes,
+unchanged, on darwin/arm64, darwin/amd64 (Rosetta), and linux 386, arm,
+riscv64, ppc64le and **s390x** (big-endian: `Bytes` gives NumPy's
+little-endian bytes) under qemu.
+
+**Integers, raw output, and every sampler that is plain arithmetic match NumPy
+to the bit** — including the ziggurat normal and exponential (tails included,
+over 200,000 draws), binomial, Poisson, geometric, choice, shuffle and the
+legacy `randint`/`rand`/`binomial`/`poisson`.
+
+### Negative controls
+
+Each mutation below was made on purpose and the suite run:
+
+| mutation | result |
+|---|---|
+| PCG64 multiplier, low word + 2 | 1,155 failures |
+| SeedSequence `MIX_MULT_L` + 2 | 1,998 failures |
+| SeedSequence `MULT_A` + 2 | 2,005 failures |
+| one ziggurat `ki_double` entry losing its last hex digit (an actual bug of the first table converter, which stripped a trailing `F` as if it were a C suffix) | 15 failures, all in the long runs |
+| `RandomState.randint` with Lemire instead of masked rejection | 44 failures |
+| MT19937 seeded with `pos = 624` instead of NumPy's `623` | 266 failures |
+| legacy BTPE adding the last two Stirling terms with the corrected signs | the `RandomState` binomial digests fail |
+| legacy BTPE constant 13680 → 13860 | **not detected**: the term changes by ~2·10⁻⁹ and flips no decision in 10⁶ draws at any of four (n, p) tried |
+
+## Where the numbers can differ: the C library
+
+Some samplers call `log`, `exp`, `log1p`, `expm1`, `pow`, `cos` or `acos`.
+NumPy takes them from the C library; Go's `math` package does not round them
+the same way. Measured against the libm of the x86-64 NumPy used for the
+fixture (Apple's), Go's `math` differs in 6.5% (`log`) to 71% (`pow`) of
+random arguments. Used as-is, that made 3% of `np.random.randn` values and
+10% of `lognormal` values one ULP off.
+
+So `random` rounds `log`, `log1p`, `exp`, `expm1` and `pow` **correctly**
+([`libm.go`](https://github.com/go-ndarray/ndarray/blob/v0.11.0/random/libm.go)): a fast evaluation to about 2⁻⁶⁸ with Ziv's
+rounding test, and a double-double fallback (about 2⁻¹⁰⁰) for the one argument
+in ~2¹⁵ that needs it. Checked against mpmath at 300 bits on 63,000 arguments
+(0 misrounded; the committed `libm.json` keeps 8,400), and the fast paths
+against the fallback on 400,000 more. The code converts every product
+explicitly, so Go cannot fuse it into an FMA, and the result is the same on
+every architecture.
+
+The C libraries round these functions correctly almost always — Apple's x86-64
+libm misrounds 0.06% (`log`) to 0.3% (`exp`) of arguments; `pow(x, 0.5)` is
+misrounded for 0.19% (it is not `sqrt`). What remains, on 10⁶ draws against
+x86-64 NumPy, seed 0:
+
+| sampler | values that differ (of 10⁶) | max ULP | with Go's `math` instead |
+|---|---:|---:|---:|
+| `RandomState.randn` | 294 | 3 | 29,631 |
+| `RandomState.normal(1, 3)` | 217 | 32¹ | 22,429 |
+| `RandomState.standard_gamma(3)` | 72 | 6 | 5,762 |
+| `RandomState.standard_exponential` | 664 | 1 | 68,249 |
+| `standard_gamma(0.5)` (pow) | 1,870 | 3 | 16,212 |
+| `beta(0.5, 0.5)` (pow) | 1,693 | 4 | 1,693 |
+| `lognormal(0.5, 0.25)` | 2,263 | 1 | 99,340 |
+| `laplace` / `logistic` | 690 / 601 | 1 | 68,357 / 68,518 |
+| `gumbel` | 1,297 | 2049¹ | 131,066 |
+| `standard_exponential(method="inv")` | 692 | 1 | 68,383 |
+| `pareto(3)` (expm1) | 4,645 | 1 | 65,593 |
+| `weibull(2)` (pow) | 1,777 | 1 | 1,779 |
+| `vonmises(0.5, 2)` (cos, acos) | 126,540 | 2¹⁸ ¹ | 126,540 |
+| `standard_normal`, `standard_gamma(3)`, `beta(2, 3)`, `chisquare`, `standard_t`, `binomial`, `poisson(3)`, `poisson(50)`, `geometric(0.05)` | **0** | — | 0 – 2 |
+
+¹ Large ULP counts are absolute errors of about one ULP of the inner value,
+measured relative to a result near zero (`1 + 3z` for `z ≈ −1/3`;
+`log(-log(U))` for `U` near 1/e; an angle wrapped to near 0).
+
+Each of the residual differences checked by hand (with mpmath) is a value
+NumPy's libm rounded incorrectly. `vonmises` keeps Go's `cos` and `acos`:
+Apple's `cos` itself misrounds 17% of arguments in [0, π], so rounding
+correctly would not bring it closer. No difference has ever changed an
+accept/reject decision — every integer-valued sampler matches exactly — and
+the test bounds each float sampler's ULP distance individually
+(`libmULP` in `numpy_test.go`).
+
+## NumPy does not agree with itself
+
+The same fixture generated by NumPy 2.5.3's **aarch64** macOS build differs from
+its x86-64 build: the C compiler fuses `a*b + c` into one FMA there (clang's
+default `-ffp-contract=on`), and Apple's aarch64 libm rounds differently.
+
+| | values that differ, aarch64 vs x86-64 NumPy |
+|---|---|
+| `Generator.uniform` | 20.8% of 200,000 (`low + range*U` is fused) |
+| `RandomState.randn` | 865 of 6,072 (14%) |
+| `RandomState.normal` | 1,063 of 4,048 (26%) |
+| `RandomState.standard_gamma` | 481 of 4,108 |
+| `Generator.standard_gamma`, `chisquare`, `f`, `noncentral_*`, `beta`, `gamma`, `lognormal`, `standard_t`, `vonmises` | a few % |
+| digests of `Generator.normal`, `uniform`, `standard_gamma` | differ |
+
+This package follows the C source operation by operation without fusing,
+which is what NumPy's x86-64 builds (Linux, Windows, Intel Macs) compute, and
+it gives the same numbers on every architecture. A notebook run on an Apple
+silicon Mac or an aarch64 Linux box will see NumPy's last bits move in those
+samplers; that is NumPy, not a migration error.
+
+## Deliberate differences
+
+| NumPy | here | why |
+|---|---|---|
+| `randint`'s default dtype is C `long`: int64 on Linux/macOS, **int32 on Windows** | always int64 | the same numbers on every platform; `RandintOf(ndarray.Int32, …)` for Windows' |
+| `size=None` returns a Python scalar | a 0-d array | one return type; same draws |
+| `shuffle` works on views | `Shuffle` needs an array that owns its storage (`ErrNotContiguous` for a view) | ndarray hands out a view's elements as a copy |
+| `permuted(x)` of a Fortran-ordered `x` shuffles in memory order | always C order | ndarray arrays are row-major |
+| NumPy's `(int)` of ±∞ in rejection branches is undefined C | the same branch is rejected (the `v == 0` tests) | — |
+| `zipf`'s `X > LONG_MAX` rejection | kept as a condition, unreachable short of `pow` rounding past 2⁶³ | — |
+| `exp` of an argument whose result is subnormal | rounded twice (to 53 bits, then to the subnormal) | only below 2⁻¹⁰²² |
+
+## Performance
+
+Apple M4 Max, 16 cores, shared and loaded (load average 10.7–14.4 during the
+runs). Seven interleaved rounds: the Go benchmark (`go test -bench`, 20
+iterations) then NumPy 2.5.3 aarch64 (median of 20 calls); medians of the
+rounds, in milliseconds for 10⁶ values:
+
+| | Go | NumPy | Go / NumPy | Go range | NumPy range |
+|---|---:|---:|---:|---|---|
+| `random(10**6)` | 3.51 | 3.13 | 1.12 | 3.47–3.83 | 3.09–3.41 |
+| `random(10**6, dtype=float32)` | 1.89 | 1.64 | 1.16 | 1.88–2.24 | 1.62–1.65 |
+| `integers(0, 1000, 10**6)` | 2.31 | 1.70 | 1.36 | 2.29–2.32 | 1.68–1.71 |
+| `standard_normal(10**6)` | 4.12 | 3.68 | 1.12 | 4.07–4.16 | 3.67–3.72 |
+| `RandomState(0).randn(10**6)` | 20.36 | 9.65 | 2.11 | 19.99–20.56 | 9.61–9.68 |
+
+Go is **slower everywhere**: 12–36% for the Generator, 2.1× for the legacy
+`randn`. The Generator gap is the 128-bit LCG step (3.3 ns per output, a
+dependent multiply chain) plus zeroing the result (`ZerosOf`; NumPy's
+`np.empty` does not); the default PCG64 already gets direct, inlinable calls
+for these four functions. `randn` pays for a correctly rounded `log` on every
+pair (6 ns against 3.4 ns for `math.Log`) and for MT19937 behind an
+interface. The samplers that call `pow`/`exp` per draw cost about what
+`math.Pow`/`math.Exp` would (19.5 vs 18.4 ns; 6.3 vs 3.9 ns).
+
+## Regenerating
+
+```sh
+uv run --python cpython-3.12-macos-x86_64 --with numpy==2.5.3 python random/testdata/cases.py > random/testdata/cases.json
+uv run --with mpmath python random/testdata/libm.py > random/testdata/libm.json
+uv run --with mpmath python random/testdata/libm_tables.py > random/zz_libm_tables.go
+python3 random/testdata/zig2go.py <numpy>/numpy/random/src/distributions/ziggurat_constants.h random/zz_ziggurat.go
+```
+
+`random/testdata/quant.py` and `TestQuantifyLibm` (run with
+`RANDOM_QUANT_DIR`) reproduce the 10⁶-draw table.
diff --git a/content/en/roadmap.md b/content/en/roadmap.md
index cefc8f4..07f633b 100644
--- a/content/en/roadmap.md
+++ b/content/en/roadmap.md
@@ -1,4 +1,6 @@
 ---
+# Below the "Mirrored from" line, this page is go-ndarray/ndarray
+# docs/plan-ndarray.md; scripts/sync-mirrors.sh re-syncs it.
 title: "Implementation plan — N-dimensional arrays in pure Go (go-ndarray/ndarray)"
 linkTitle: "Roadmap"
 weight: 10
@@ -6,6 +8,8 @@ description: "The plan for go-ndarray and what ships today."
 tags: [roadmap]
 ---
 
+> Mirrored from [`docs/plan-ndarray.md`](https://github.com/go-ndarray/ndarray/blob/v0.11.0/docs/plan-ndarray.md) at v0.11.0.
+
 > Goal: a **pure-Go (CGO=0)** NumPy-style N-dimensional array library, as a
 > **standalone, reusable** module with correct scalar kernels now and a clean
 > path to SIMD acceleration across all six 64-bit Go targets.
@@ -58,14 +62,7 @@ contiguous `[]float64`. These are the **reference** and the **fallback**.
 [go-asmgen](https://github.com/go-asmgen) — Plan 9 assembly emitted for all six
 64-bit Go SIMD targets: **amd64, arm64, riscv64, loong64, ppc64le (VSX),
 s390x (vector, big-endian)** — selected at runtime via build tags / CPU-feature
-detection, behind the *same* kernel signatures and the *same* tests. **Today
-(v0.2.5) amd64, arm64, ppc64le and loong64 have SIMD kernels** (loong64's LASX
-kernels run only when the CPU reports LASX); riscv64, s390x and the 32-bit
-targets run the pure-Go reference code. Go's assembler has vector-double
-arithmetic on s390x and riscv64 (V optional), so those two are work not yet
-done, not a toolchain limit; it has none on ppc64le, whose VSX kernels use the
-float64 encoders go-asmgen v0.13.0 emits as `WORD`s, and no vector FMA or
-broadcast load on loong64, which go-asmgen v0.14.0 encodes (transitionally). Because the
+detection, behind the *same* kernel signatures and the *same* tests. Because the
 public API and the test suite address the kernels through that narrow interface,
 no caller changes and correctness stays pinned by the existing 100%-coverage
 suite (plus the per-arch CI jobs already wired in `.github/workflows/ci.yml`).
@@ -114,7 +111,8 @@ suite (plus the per-arch CI jobs already wired in `.github/workflows/ci.yml`).
   Phase 4.
 - **Performance — DONE (beats NumPy on its core ops; matmul at tuned-BLAS
   parity).** Several levers behind the kernel seam, measured honestly vs NumPy
-  2.2.4 (see [docs/performance.md](performance.md)):
+  2.2.4 (see [docs/perf.md](performance.md) and the repo's
+  [BENCHMARKS.md](https://github.com/go-ndarray/ndarray/blob/v0.11.0/BENCHMARKS.md)):
   - **Multicore.** `Add/Sub/Mul/Div/Map`, `Sum/Max/Min`, the axis reductions,
     `Concatenate`/`Stack`, and the GEMM fan out across `GOMAXPROCS` above a size
     threshold (small arrays stay serial). Plus an allocation fix: the
@@ -126,35 +124,35 @@ 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
     validated per-arch in CI (bit-identical to the scalar oracle for the
     elementwise/sqrt/max ops; tight-tolerance for the sum reduction, whose
-    lane-parallel grouping is a valid reordering — the same trade-off as NumPy's
-    pairwise sum). The other targets (riscv64, s390x; ppc64le until its
-    VSX kernels in v0.2.3, loong64 until its LASX kernels in v0.2.5) keep the validated scalar oracles plus a scalar
-    4×4 GEMM micro-kernel over the packed panels, and still win via packing + cache-blocking + multicore. Split-CI: the pure-Go core + multicore
-    are held to 100% statement coverage; the generated `.s` is validated by the
-    per-arch native/qemu execution jobs.
+    lane-parallel grouping is a valid reordering — the same trade-off as
+    NumPy's pairwise sum). Since v0.2.3 ppc64le has VSX kernels too (sum,
+    dot, sqrt, elementwise, an 8×8 GEMM micro-kernel), through go-asmgen's
+    VSX encoders, and since v0.2.5 loong64 has LASX kernels (used when the CPU
+    reports LASX), and since v0.7.0 riscv64 has RVV kernels (sum, dot, sqrt,
+    elementwise, and since v0.7.1 a 6×8 GEMM micro-kernel; used when the CPU
+    reports V), and since v0.8.0 s390x has vector-facility kernels (sum, dot,
+    sqrt, elementwise, since v0.9.0 a 4×8 GEMM micro-kernel and since v0.9.1
+    max/min on z14 and later; measured on a
+    real IBM z17). The
+    32-bit targets keep the
+    validated scalar oracles plus a scalar 4×4 GEMM micro-kernel over the
+    packed panels, with packing + cache-blocking + multicore; they have not
+    been measured against NumPy. Split-CI: the pure-Go
+    core + multicore are held to 100% statement coverage; the generated `.s`
+    is validated by the per-arch native/qemu execution jobs.
   - **SIMD `Sqrt`/`Max` and the packed, cache-blocked GEMM — DONE (no longer
     TODO).** `Sqrt`/`Max`/`Min` now WIN vs single-threaded NumPy, and the
     panel-packed, cache-blocked GEMM with its SIMD-FMA micro-kernel reaches
     **parity with tuned BLAS** (single-threaded vecLib at 1024² on M4 Max,
     ≈0.99× vs multi-threaded OpenBLAS at n=1024 on an arm64 VM) and beats
-    gonum 4–10× everywhere. Remaining headroom is small-n per-call overhead.
-    (The AVX2/FMA 6×8 amd64 micro-kernel, chosen at run time by a CPUID probe,
-    shipped in v0.2.0.)
-- **Memory reuse — DONE (v0.2.0).** `Workspace`: `NewWorkspace()`, `ws.Use(x)`
-  binds an input, every result computed from a bound array is allocated from the
-  workspace arena, `ws.Reset()` recycles everything allocated since the last
-  Reset (arrays allocated from it are invalid afterwards), and `a.Detach()`
-  returns a heap copy that survives it. In loops it removes the Go GC + zeroing
-  cost of temporaries, which NumPy frees by reference count and Go cannot. On a
-  Zen 3 `x + y` on 1 024 elements goes from 2.6 µs to 0.38 µs (2.1× NumPy) and
-  `sqrt(x*y + x)` on 256 Ki from 1.18 ms to 0.39 ms (7.5× NumPy); see
-  [Performance](performance.md).
+    gonum 4–10× everywhere. Remaining headroom is small-n per-call overhead
+    and an AVX2/FMA amd64 micro-kernel.
 - **More reductions — DONE.** `ArgMax`/`ArgMin` (flat and per-axis with
   keepdims), `CumSum`/`CumProd` (per-axis scans + flat), `Clip(lo,hi)`, and the
   three-operand broadcasting `Where(cond,t,f)`. Match numpy 2.2.4.
@@ -163,8 +161,10 @@ suite (plus the per-arch CI jobs already wired in `.github/workflows/ci.yml`).
   into the flattened array, negative indices from the end). These give NumPy
   boolean indexing today using the 0/1 masks the comparison ufuncs return; a
   first-class bool dtype (so masks are real booleans) remains Phase 2.
-- **Phase 2 — dtypes.** Generalise beyond `float64` (float32, int64/int32,
-  complex128, bool) with a dtype abstraction; typed kernels.
+- **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](dtypes.md). Next:
+  SIMD kernels for `float32`, and the `fft` bridge.
 - **Phase 3 — broadcasting ufuncs.** A general ufunc framework over views
   (the axis-reduction primitives above are the first step); more elementwise
   and reduction ufuncs.
@@ -184,6 +184,8 @@ suite (plus the per-arch CI jobs already wired in `.github/workflows/ci.yml`).
 3. **Errors, not panics, for shape/broadcast problems**; panics reserved for
    programmer index errors (matching Go slice semantics). *Settled.*
 4. **`float64` first**, dtype abstraction deferred to Phase 2 but kept in mind
-   so the kernel interface generalises cleanly.
+   so the kernel interface generalises cleanly. *Superseded:* Phase 2 is now
+   designed in [design-dtypes.md](dtypes.md) — the dtype is a runtime
+   property of the one `Array` type, with NumPy 2.x promotion.
 
 BSD-3-Clause.
diff --git a/content/en/security.md b/content/en/security.md
index 95672b4..a454703 100644
--- a/content/en/security.md
+++ b/content/en/security.md
@@ -8,7 +8,7 @@ description: "How to report a vulnerability, what the library guarantees, and th
 tags: [security]
 ---
 
-> Mirrored from [`SECURITY.md`](https://github.com/go-ndarray/ndarray/blob/main/SECURITY.md) at v0.9.1.
+> Mirrored from [`SECURITY.md`](https://github.com/go-ndarray/ndarray/blob/main/SECURITY.md) at v0.11.0.
 
 ## Reporting
 
diff --git a/hugo.yaml b/hugo.yaml
index 697c906..edbc4b2 100644
--- a/hugo.yaml
+++ b/hugo.yaml
@@ -6,7 +6,7 @@
 # have used since the site was built with MkDocs and mike. The CI sets it
 # through --baseURL, from params.ndarray.version below. This default is the
 # same build.
-baseURL: "https://go-ndarray.github.io/docs/0.9/"
+baseURL: "https://go-ndarray.github.io/docs/0.11/"
 title: "go-ndarray"
 # themes/hextra is the tannevaled/hextra fork, pinned at a release tag:
 # upstream Hextra plus opt-in features (partial menu items, combos,
@@ -92,12 +92,12 @@ menu:
         icon: github
 
 params:
-  description: "A pure-Go (CGO=0) NumPy-style n-dimensional array library for float64: views, broadcasting, ufuncs, reductions, linalg; multicore + SIMD."
+  description: "A pure-Go (CGO=0) NumPy-style n-dimensional array library in 13 dtypes: views, broadcasting, ufuncs, reductions, numpy.linalg, numpy.fft, numpy.random; multicore + SIMD."
 
   # The ndarray release these pages describe. Its minor version is the
   # documentation's version label (v0.9.1 -> 0.9): the CI reads it here.
   ndarray:
-    version: "v0.9.1"
+    version: "v0.11.0"
 
   navbar:
     displayTitle: true
diff --git a/scripts/sync-mirrors.sh b/scripts/sync-mirrors.sh
new file mode 100755
index 0000000..6d53c75
--- /dev/null
+++ b/scripts/sync-mirrors.sh
@@ -0,0 +1,104 @@
+#!/bin/sh
+# Re-sync the pages that copy files of go-ndarray/ndarray.
+#
+#   scripts/sync-mirrors.sh <ndarray checkout> <tag>
+#
+# <ndarray checkout> must be at <tag> (vX.Y.Z); the script refuses otherwise,
+# since the "Mirrored from ... at <tag>" lines would then be false.
+#
+# Each page keeps what is ours -- its front matter, and for performance.md the
+# introduction above "## How a pure-Go library can beat NumPy" -- and takes the
+# rest from the source verbatim, except for the source's repository-relative
+# links, which mean nothing on this site and are rewritten:
+#
+#   a page this site mirrors     -> that page   (perf.md -> performance.md)
+#   a package directory          -> pkg.go.dev at <tag>
+#   any other file of the repo   -> github.com/go-ndarray/ndarray/blob/<tag>/...
+#
+# Then it lists any relative link left, and fails if there is one: a link the
+# rules do not know is a link that would break.
+#
+# security.md is not handled here: its source links files in ways a rule
+# cannot know (see the README), so it is re-synced by hand.
+set -eu
+
+die() { echo "sync-mirrors.sh: $*" >&2; exit 1; }
+
+[ $# -eq 2 ] || die "usage: $0 <ndarray checkout> <tag>"
+nd=$1
+tag=$2
+printf '%s\n' "$tag" | grep -Eq '^v(0|[1-9][0-9]*)\.(0|[1-9][0-9]*)\.(0|[1-9][0-9]*)$' || die "tag must be vX.Y.Z, not \"$tag\""
+head=$(git -C "$nd" rev-parse HEAD) || die "$nd is not a git checkout"
+at=$(git -C "$nd" rev-parse "$tag^{commit}") || die "$nd has no tag $tag"
+[ "$head" = "$at" ] || die "$nd is not at $tag"
+
+blob="https://github.com/go-ndarray/ndarray/blob/$tag"
+pkg="https://pkg.go.dev/github.com/go-ndarray/ndarray@$tag"
+pages=content/en
+
+# links: rewrite the source's repository-relative links (stdin -> stdout).
+# The source files live in docs/, so ../X is the repository root.
+links() {
+	sed -E \
+		-e "s#\]\(perf\.md(\#[^)]*)?\)#](performance.md\1)#g" \
+		-e "s#\]\(plan-ndarray\.md(\#[^)]*)?\)#](roadmap.md\1)#g" \
+		-e "s#\]\(design-dtypes\.md(\#[^)]*)?\)#](dtypes.md\1)#g" \
+		-e "s#\]\(\.\./(linalg|random|fft)/?\)#]($pkg/\1)#g" \
+		-e "s#\]\(\.\./([^)]+)\)#]($blob/\1)#g"
+}
+
+# front: print a page's front matter, the two --- lines included.
+front() {
+	awk 'NR == 1 && $0 == "---" { print; inside = 1; next }
+	     inside { print; if ($0 == "---") exit }' "$1"
+}
+
+# mirror <page> <source in docs/>: front matter, a "Mirrored from" line, then
+# the source without its # title (the front matter's title replaces it).
+mirror() {
+	page=$pages/$1
+	src=$nd/docs/$2
+	[ -f "$page" ] || die "$page is missing: write its front matter first"
+	[ -f "$src" ] || die "$src is missing"
+	fm=$(front "$page")
+	[ -n "$fm" ] || die "$page has no front matter"
+	{
+		printf '%s\n\n' "$fm"
+		printf '> Mirrored from [`docs/%s`](%s/docs/%s) at %s.\n' "$2" "$blob" "$2" "$tag"
+		sed 1d "$src" | links
+	} > "$page.tmp"
+	mv "$page.tmp" "$page"
+}
+
+mirror roadmap.md plan-ndarray.md
+mirror dtypes.md design-dtypes.md
+mirror linalg.md linalg.md
+mirror random.md random.md
+
+# performance: keep our part above the heading, replace the rest.
+perf=$pages/performance.md
+{
+	sed '/^## How a pure-Go library can beat NumPy/,$d' "$perf"
+	sed -n '/^## How a pure-Go library can beat NumPy/,$p' "$nd/docs/perf.md" | links
+} > "$perf.tmp"
+mv "$perf.tmp" "$perf"
+
+# Every relative link left must name a page of this site. Code is not text:
+# `Item[T](a, idx...)` is not a link, so fenced blocks and code spans are
+# blanked first (lines kept, so the numbers reported stay right).
+prose() {
+	awk '/^[[:space:]]*```/ { fence = !fence; print ""; next }
+	     fence { print ""; next }
+	     { gsub(/`[^`]*`/, ""); print }' "$1"
+}
+bad=0
+for f in roadmap.md dtypes.md linalg.md random.md performance.md; do
+	prose "$pages/$f" | grep -noE '\]\([^)h#][^)]*\)' | while IFS=: read -r line link; do
+		target=${link#](}
+		target=${target%)}
+		target=${target%%#*}
+		[ -f "$pages/$target" ] || { echo "$pages/$f:$line: $target is not a page of this site" >&2; echo x; }
+	done
+done | grep -q x && bad=1
+[ "$bad" -eq 0 ] || die "relative links left that the rules above do not rewrite"
+echo "synced from $tag; now set params.ndarray.version in hugo.yaml to $tag"