From c2afb730d64fac16fc82fd9c7de22a38fc98cc77 Mon Sep 17 00:00:00 2001
From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com>
Date: Thu, 1 Oct 2026 00:12:52 +0100
Subject: [PATCH 1/3] docs: add Signed commits section to CONTRIBUTING

Owner ruling D218. See docs/SIGNING-POLICY.adoc in hyperpolymath/standards.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WRvDivYwLSeVCJUrfjic3f
---
 .github/CONTRIBUTING.md | 17 +++++++++++++++++
 CONTRIBUTING.adoc       | 17 +++++++++++++++++
 2 files changed, 34 insertions(+)

diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md
index bb7aa60..48450a5 100644
--- a/.github/CONTRIBUTING.md
+++ b/.github/CONTRIBUTING.md
@@ -69,3 +69,20 @@ Conduct](CODE_OF_CONDUCT.md).
 
 By contributing, you agree that your contributions will be licensed
 under the same license as the project (see LICENSE).
+
+## Signed commits
+
+Every commit that reaches the default branch must be signed; a ruleset refuses
+unsigned pushes. Estate policy:
+[SIGNING-POLICY](https://github.com/hyperpolymath/standards/blob/main/docs/SIGNING-POLICY.adoc).
+
+- **People and interactive agents** sign with an SSH key registered on GitHub
+  as a *signing* key (`gpg.format=ssh`, `user.signingkey=<key>.pub`,
+  `commit.gpgsign=true`). The committer email must be verified on that account.
+- **Apps, bots and workflows** never `git push` local commits. They write
+  through the API (`createCommitOnBranch` or the estate `signed-push` action)
+  so that GitHub signs each commit.
+- Merge PRs with **squash**. The ruleset checks every commit on the PR branch,
+  not just the result, so one unsigned commit blocks the merge. Re-create such a
+  branch with signed commits (`git cherry-pick -S`) and open a new PR.
+  Rebase-merge replays commits unsigned and is disabled.
diff --git a/CONTRIBUTING.adoc b/CONTRIBUTING.adoc
index d778e8d..eeda4ed 100644
--- a/CONTRIBUTING.adoc
+++ b/CONTRIBUTING.adoc
@@ -56,3 +56,20 @@ This root document exists because the estate docs gate
 `CONTRIBUTING.md`, `CONTRIBUTING.adoc`, or `3-practice/CONTRIBUTING.adoc`
 at the repository root. Estate documentation policy: AsciiDoc by default —
 see `hyperpolymath/standards`.
+
+== Signed commits
+
+Every commit that reaches the default branch must be signed; a ruleset refuses
+unsigned pushes. Estate policy:
+https://github.com/hyperpolymath/standards/blob/main/docs/SIGNING-POLICY.adoc[SIGNING-POLICY].
+
+* **People and interactive agents** sign with an SSH key registered on GitHub
+  as a *signing* key (`gpg.format=ssh`, `user.signingkey=<key>.pub`,
+  `commit.gpgsign=true`). The committer email must be verified on that account.
+* **Apps, bots and workflows** never `git push` local commits. They write
+  through the API (`createCommitOnBranch` or the estate `signed-push` action)
+  so that GitHub signs each commit.
+* Merge PRs with **squash**. The ruleset checks every commit on the PR branch,
+  not just the result, so one unsigned commit blocks the merge. Re-create such a
+  branch with signed commits (`git cherry-pick -S`) and open a new PR.
+  Rebase-merge replays commits unsigned and is disabled.

From 529eddeac4c809907e40a962d6e4bfe892e49694 Mon Sep 17 00:00:00 2001
From: "coderabbitai[bot]"
 <136622811+coderabbitai[bot]@users.noreply.github.com>
Date: Wed, 30 Sep 2026 23:27:39 +0000
Subject: [PATCH 2/3] docs(contributing): document Git 2.34 minimum for SSH
 commit signing

---
 .github/CONTRIBUTING.md | 2 ++
 CONTRIBUTING.adoc       | 2 ++
 2 files changed, 4 insertions(+)

diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md
index 48450a5..81fca4b 100644
--- a/.github/CONTRIBUTING.md
+++ b/.github/CONTRIBUTING.md
@@ -76,6 +76,8 @@ Every commit that reaches the default branch must be signed; a ruleset refuses
 unsigned pushes. Estate policy:
 [SIGNING-POLICY](https://github.com/hyperpolymath/standards/blob/main/docs/SIGNING-POLICY.adoc).
 
+SSH commit signing requires Git 2.34 or later.
+
 - **People and interactive agents** sign with an SSH key registered on GitHub
   as a *signing* key (`gpg.format=ssh`, `user.signingkey=<key>.pub`,
   `commit.gpgsign=true`). The committer email must be verified on that account.
diff --git a/CONTRIBUTING.adoc b/CONTRIBUTING.adoc
index eeda4ed..e5ee671 100644
--- a/CONTRIBUTING.adoc
+++ b/CONTRIBUTING.adoc
@@ -63,6 +63,8 @@ Every commit that reaches the default branch must be signed; a ruleset refuses
 unsigned pushes. Estate policy:
 https://github.com/hyperpolymath/standards/blob/main/docs/SIGNING-POLICY.adoc[SIGNING-POLICY].
 
+SSH commit signing requires Git 2.34 or later.
+
 * **People and interactive agents** sign with an SSH key registered on GitHub
   as a *signing* key (`gpg.format=ssh`, `user.signingkey=<key>.pub`,
   `commit.gpgsign=true`). The committer email must be verified on that account.

From a91c5f58ffc6e85880bcbd83e758940576f8eef7 Mon Sep 17 00:00:00 2001
From: "coderabbitai[bot]"
 <136622811+coderabbitai[bot]@users.noreply.github.com>
Date: Thu, 1 Oct 2026 21:50:50 +0000
Subject: [PATCH 3/3] refactor(repro): port silicon-core investigation probes
 from Python to Julia

---
 .../B-backend-status.adoc                     |  8 +-
 .../D-numerical-parity-risks.adoc             |  6 +-
 .../E-memory-performance.adoc                 |  2 +-
 .../J-upstream-issues.adoc                    | 15 ++--
 .../2026-09-26-silicon-core-recon/README.adoc |  2 +-
 .../repro/README.adoc                         | 28 +++++--
 .../repro/activation_probe.jl                 | 40 ++++++++++
 .../repro/activation_probe.py                 | 44 -----------
 .../repro/batchnorm_probe.jl                  | 25 ++++++
 .../repro/batchnorm_probe.py                  | 22 ------
 .../repro/canary_probe.jl                     | 18 +++++
 .../repro/canary_probe.py                     | 21 -----
 .../repro/chunking_probe.jl                   | 32 ++++++++
 .../repro/chunking_probe.py                   | 29 -------
 .../repro/layout_probe.jl                     | 62 +++++++++++++++
 .../repro/layout_probe.py                     | 79 -------------------
 .../repro/matmul_cost.jl                      | 41 ++++++++++
 .../repro/matmul_cost.py                      | 34 --------
 .../repro/probe_common.jl                     | 30 +++++++
 .../repro/repro_batchnorm.zig                 |  2 +-
 20 files changed, 286 insertions(+), 254 deletions(-)
 create mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/activation_probe.jl
 delete mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/activation_probe.py
 create mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/batchnorm_probe.jl
 delete mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/batchnorm_probe.py
 create mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/canary_probe.jl
 delete mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/canary_probe.py
 create mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/chunking_probe.jl
 delete mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/chunking_probe.py
 create mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/layout_probe.jl
 delete mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/layout_probe.py
 create mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/matmul_cost.jl
 delete mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/matmul_cost.py
 create mode 100644 docs/investigation/2026-09-26-silicon-core-recon/repro/probe_common.jl

diff --git a/docs/investigation/2026-09-26-silicon-core-recon/B-backend-status.adoc b/docs/investigation/2026-09-26-silicon-core-recon/B-backend-status.adoc
index 067394e..aa03018 100644
--- a/docs/investigation/2026-09-26-silicon-core-recon/B-backend-status.adoc
+++ b/docs/investigation/2026-09-26-silicon-core-recon/B-backend-status.adoc
@@ -32,12 +32,12 @@ Built with Zig 0.15.2, `-Doptimize=ReleaseFast` (as CI). `[empirical]` unless no
 | `backend_matmul` → `axiom_matmul_checked` | yes (16×12·12×7) | correct | 9–16× slower than `axiom_matmul` (E.1). Rejects NaN/Inf with an error — different contract from the Julia reference (D.4).
 | `backend_relu` → `axiom_relu_checked` | yes (128; NaN → error) | correct | In-place `backend_relu!` uses the *unchecked* `axiom_relu_inplace` — same op, different NaN contract.
 | `backend_gelu`, `backend_sigmoid` → `axiom_gelu/sigmoid` | yes | correct | GELU uses the overflow-safe `1 − 2/(e+1)` form.
-| `backend_tanh` / `backend_tanh!` → `axiom_tanh(_inplace)` | *no* | *NaN for x ≥ 44.5 and +Inf* (Julia: 1.0); rel. error 1.3e-3 at x=1e-5, 100 % at 1e-8 | `activations.zig:103–119` uses `(e^{2x}−1)/(e^{2x}+1)`; `repro/activation_probe.py`.
+| `backend_tanh` / `backend_tanh!` → `axiom_tanh(_inplace)` | *no* | *NaN for x ≥ 44.5 and +Inf* (Julia: 1.0); rel. error 1.3e-3 at x=1e-5, 100 % at 1e-8 | `activations.zig:103–119` uses `(e^{2x}−1)/(e^{2x}+1)`; `repro/activation_probe.jl`.
 | `backend_softmax`, `backend_log_softmax` → `axiom_softmax/log_softmax` | softmax yes (8×5) | correct for `dim == ndims` | `dim` argument ignored (`zig_ffi.jl:329,351`) `[static]`. `[empirical]` B=4…17 × 2048 correct through the `.so`.
 | `backend_conv2d` → `axiom_conv2d` | yes (2×10×10×3, 3×3×3×4) | correct at that shape | Wrapper does full `permutedims` copies both ways (E.3). Stride/padding/dilation/groups other than the tested defaults: untested.
-| `backend_batchnorm` → `axiom_batchnorm` | yes (6×8) | *wrong for `num_features` 4097–8191; SIGSEGV from 8192* | Fixed `[4096]f32` stack scratch (`norm.zig:143`); `repro/batchnorm_probe.py`. Reachable from `forward(bn::BatchNorm)` (`abstract.jl:1072–1086`) with the Zig backend set. `training` flag ignored `[static]`.
-| `backend_layernorm`, `backend_rmsnorm` → `axiom_layernorm/rmsnorm` | *no* | correct through the shipped `.so` for B=4…17 × 2048/4096/65536, canary intact | `threading.zig:392,469,523` `batch_size − last_start` underflows at `batch_size = 5` (with `MAX_WORKERS = 3`, `total ≥ 8192`): panics in Debug/ReleaseSafe; the same function compiled ReleaseFast in a `zig test` binary segfaults; in the shipped `.so` LLVM happened to make it benign. Undefined behaviour either way. `repro/repro_underflow*.zig`, `repro/canary_probe.py`. 3-D inputs: wrapper treats `size(x,2)` as hidden `[static]`.
-| `backend_maxpool2d`, `backend_global_avgpool2d` → `axiom_maxpool2d`, `axiom_global_avgpool2d` (no Zig wrapper exists for `avgpool2d`) | *no* | *wrong values whenever N>1 or C>1 (and for non-square N=1,C=1)* | Wrapper passes column-major bytes to row-major NHWC kernels without `_to_row_major_vec` (`zig_ffi.jl:479–507`, `:513`). `repro/layout_probe.py`: maxpool wrong 5/8, 5/8, 49/54, 7/8 elements in four cases; global-avgpool wrong for every N>1 or C>1. Edge contracts: all-`−Inf` window → `−3.4028235e38` (Julia `−Inf`); NaN in window → `1.0` (Julia `NaN`). Not reachable from `MaxPool2d.forward` (pure Julia) — reachable from direct calls, `SmartBackend`, and benchmarks.
+| `backend_batchnorm` → `axiom_batchnorm` | yes (6×8) | *wrong for `num_features` 4097–8191; SIGSEGV from 8192* | Fixed `[4096]f32` stack scratch (`norm.zig:143`); `repro/batchnorm_probe.jl`. Reachable from `forward(bn::BatchNorm)` (`abstract.jl:1072–1086`) with the Zig backend set. `training` flag ignored `[static]`.
+| `backend_layernorm`, `backend_rmsnorm` → `axiom_layernorm/rmsnorm` | *no* | correct through the shipped `.so` for B=4…17 × 2048/4096/65536, canary intact | `threading.zig:392,469,523` `batch_size − last_start` underflows at `batch_size = 5` (with `MAX_WORKERS = 3`, `total ≥ 8192`): panics in Debug/ReleaseSafe; the same function compiled ReleaseFast in a `zig test` binary segfaults; in the shipped `.so` LLVM happened to make it benign. Undefined behaviour either way. `repro/repro_underflow*.zig`, `repro/canary_probe.jl`. 3-D inputs: wrapper treats `size(x,2)` as hidden `[static]`.
+| `backend_maxpool2d`, `backend_global_avgpool2d` → `axiom_maxpool2d`, `axiom_global_avgpool2d` (no Zig wrapper exists for `avgpool2d`) | *no* | *wrong values whenever N>1 or C>1 (and for non-square N=1,C=1)* | Wrapper passes column-major bytes to row-major NHWC kernels without `_to_row_major_vec` (`zig_ffi.jl:479–507`, `:513`). `repro/layout_probe.jl`: maxpool wrong 5/8, 5/8, 49/54, 7/8 elements in four cases; global-avgpool wrong for every N>1 or C>1. Edge contracts: all-`−Inf` window → `−3.4028235e38` (Julia `−Inf`); NaN in window → `1.0` (Julia `NaN`). Not reachable from `MaxPool2d.forward` (pure Julia) — reachable from direct calls, `SmartBackend`, and benchmarks.
 | `backend_dense` (fused) → `axiom_dense`? | no | not probed | `zig_forward(model::Dense)` (`abstract.jl:2495–2502`) drops bias and activation `[static]`.
 | attention (`axiom_scaled_dot_product_attention`, `axiom_flash_attention`) | no Julia wrapper | silent no-op for `seq_len > 64` (`attention.zig:29`), `> 4096` for flash | Output buffer left untouched, no error code.
 | 15 unwrapped exports | — | unreachable from Julia | `add(_checked)`, `mul(_checked)`, `bmm`, `fill`, `avgpool2d`, `relu`, `relu6(_checked)`, `rotary_embedding`, `scaled_dot_product_attention(_checked)`, `flash_attention(_checked)`.
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/D-numerical-parity-risks.adoc b/docs/investigation/2026-09-26-silicon-core-recon/D-numerical-parity-risks.adoc
index b256048..ebdb8bf 100644
--- a/docs/investigation/2026-09-26-silicon-core-recon/D-numerical-parity-risks.adoc
+++ b/docs/investigation/2026-09-26-silicon-core-recon/D-numerical-parity-risks.adoc
@@ -15,7 +15,7 @@ column-major `(N,H,W,C)` buffer straight to `axiom_maxpool2d` /
 Every other array wrapper in the file converts with `_to_row_major_vec`
 (`:79`); these two do not.
 
-*Measured* (`repro/layout_probe.py`, simulates exactly what the Julia wrapper
+*Measured* (`repro/layout_probe.jl`, simulates exactly what the Julia wrapper
 does, compares with correct pooling):
 
 [cols="2,1,1"]
@@ -57,7 +57,7 @@ Julia's `tanh` returns `1.0`. `gelu` in the same file already uses the
 overflow-safe `1 − 2/(e^{2x}+1)` (`:135,144`) — the fix is to use the same
 form in `tanh_inplace`/`tanh`. Near zero the formula cancels: relative error
 1.3e-3 at `x = 1e-5`, 100 % at `x = 1e-8` (absolute error negligible).
-Measured with `repro/activation_probe.py`; `axiom_tanh_inplace` behaves the
+Measured with `repro/activation_probe.jl`; `axiom_tanh_inplace` behaves the
 same. `gelu`, `sigmoid`, `silu` are fine over the probed range (sigmoid(−100)
 flushes a denormal to 0 — harmless).
 
@@ -69,7 +69,7 @@ absent from `backend_parity.jl`).
 
 `norm.zig:143` `var inv_std: [4096]f32 = undefined;` indexed by
 `num_features` with no bound check. Through the shipped ReleaseFast `.so`
-(`repro/batchnorm_probe.py`, B=2):
+(`repro/batchnorm_probe.jl`, B=2):
 
 [cols="1,2"]
 |===
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/E-memory-performance.adoc b/docs/investigation/2026-09-26-silicon-core-recon/E-memory-performance.adoc
index 9ff15aa..b4e6cf1 100644
--- a/docs/investigation/2026-09-26-silicon-core-recon/E-memory-performance.adoc
+++ b/docs/investigation/2026-09-26-silicon-core-recon/E-memory-performance.adoc
@@ -12,7 +12,7 @@ that would settle them.
 
 `backend_matmul(::ZigBackend)` always calls `axiom_matmul_checked`
 (`zig_ffi.jl:114`), which runs a scalar `O(m·n·k)` finiteness pre-pass before
-the tiled SIMD kernel. `repro/matmul_cost.py`, same `.so`, same inputs
+the tiled SIMD kernel. `repro/matmul_cost.jl`, same `.so`, same inputs
 (median of repeats):
 
 [cols="1,1,1,1,1"]
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/J-upstream-issues.adoc b/docs/investigation/2026-09-26-silicon-core-recon/J-upstream-issues.adoc
index c2e6ec0..2955481 100644
--- a/docs/investigation/2026-09-26-silicon-core-recon/J-upstream-issues.adoc
+++ b/docs/investigation/2026-09-26-silicon-core-recon/J-upstream-issues.adoc
@@ -64,8 +64,9 @@ wrong results / false claim on a documented path or major perf regression;
 == J.2 Detailed entries (evidence + reproduction)
 
 All commands assume the library was built as CI does:
-`cd zig && zig build -Doptimize=ReleaseFast` (Zig 0.15.2) and that `numpy`
-is importable. Zig test reproductions compile against the repository's
+`cd zig && zig build -Doptimize=ReleaseFast` (Zig 0.15.2) and that Julia
+is available. The probes use only Julia standard libraries; see link:repro/README.adoc[reproduction notes]
+for the port and historical measurement details. Zig test reproductions compile against the repository's
 `zig/src` tree directly. Paths are relative to `docs/investigation/2026-09-26-silicon-core-recon/repro/`.
 
 === J-01 Zig `batchnorm` above 4096 features `[E]` — Critical
@@ -73,7 +74,7 @@ is importable. Zig test reproductions compile against the repository's
 * *Where:* `zig/src/norm.zig:143` `var inv_std: [4096]f32 = undefined;` indexed by `num_features`.
 * *Reach:* `abstract.jl:1072–1086` `forward(bn::BatchNorm, x)` → `backend_batchnorm(::ZigBackend, …)` (`zig_ffi.jl:536`) whenever the current backend is not `JuliaBackend`, `!bn.training`, `bn.affine`.
 * *Observed (shipped `.so`):* 4096 correct; 4097 → 2 wrong; 4100 → 8 wrong; 8192 / 65536 / 1e6 → SIGSEGV.
-* *Repro:* `python3 batchnorm_probe.py ../../../../zig/zig-out/lib/libaxiom_zig.so`; Zig-level: `zig test --dep axiom -Mroot=repro_batchnorm.zig -Maxiom=../../../../zig/src/axiom.zig -OReleaseSafe` → index out of bounds panic.
+* *Repro:* `julia --startup-file=no batchnorm_probe.jl ../../../../zig/zig-out/lib/libaxiom_zig.so`; Zig-level: `zig test --dep axiom -Mroot=repro_batchnorm.zig -Maxiom=../../../../zig/src/axiom.zig -OReleaseSafe` → index out of bounds panic.
 * *Julia confirmation:* `set_backend!(ZigBackend(lib)); m = Sequential(Dense(4, 8192), BatchNorm(8192)); m(Tensor(randn(Float32, 2, 4)))` — expected: process terminates.
 * *Fix:* allocate `inv_std` per call (arena/heap) or compute per feature without the scratch; add `num_features ∈ {4096, 4097, 8192, 65536}` to `zig build test` and `backend_parity.jl`.
 
@@ -81,7 +82,7 @@ is importable. Zig test reproductions compile against the repository's
 
 * *Where:* `src/backends/zig_ffi.jl:479–507` (`backend_maxpool2d`), `:513` (`backend_global_avgpool2d`) pass `input` directly; compare `_to_row_major_vec` use in conv/layernorm/softmax wrappers (`:79–95`). Kernels index row-major NHWC (`zig/src/pool.zig`).
 * *Observed:* maxpool wrong 5/8 (N=1,4×4,C=2), 5/8 (N=2,C=1), 49/54 (N=2,6×6,C=3), 7/8 (N=1,5×3,C=1,k=2,s=1); global-avgpool wrong for any N>1 or C>1; only N=1,C=1,square is correct.
-* *Repro:* `python3 layout_probe.py ../../../../zig/zig-out/lib/libaxiom_zig.so`.
+* *Repro:* `julia --startup-file=no layout_probe.jl ../../../../zig/zig-out/lib/libaxiom_zig.so`.
 * *Julia confirmation:* see D.1.
 * *Reach:* direct API, `SmartBackend`, `benchmark/benchmarks.jl`; not `MaxPool2d.forward` (pure Julia). Not parity-tested.
 
@@ -89,21 +90,21 @@ is importable. Zig test reproductions compile against the repository's
 
 * *Where:* `zig/src/activations.zig:103–119` `(e^{2x}−1)/(e^{2x}+1)`; `gelu` at `:135,144` already uses `1 − 2/(e^{2x}+1)`.
 * *Observed:* `axiom_tanh(44.5) = NaN`, `(45) = NaN`, `(50) = NaN`, `(100) = NaN`, `(+Inf) = NaN`; `(44.3) = 1.0`. Near zero: rel. err 1.3e-3 at 1e-5, 100 % at 1e-8. `gelu`, `sigmoid`, `silu` fine.
-* *Repro:* `python3 activation_probe.py ../../../../zig/zig-out/lib/libaxiom_zig.so`.
+* *Repro:* `julia --startup-file=no activation_probe.jl ../../../../zig/zig-out/lib/libaxiom_zig.so`.
 * *Reach:* `backend_tanh(::ZigBackend)`, `backend_tanh!`, `SmartBackend`. `tanh` absent from `backend_parity.jl`.
 
 === J-04 Chunking underflow at `batch_size = 5` `[E]` — Medium (UB)
 
 * *Where:* `zig/src/threading.zig:392` (`parallel_layernorm`), `:469` (`parallel_rmsnorm`), `:523` (`parallel_batch`, used by batched softmax): `last_start = (num_threads−1)·chunk`, `last_count = batch_size − last_start` with `chunk = ceil(batch/num_threads)`, `num_threads = min(4, batch)`. For `batch = 5`: `chunk = 2`, `last_start = 6 > 5` → usize wrap. Reached when `total ≥ 8192` and `batch ≥ 4` (`:294–324`).
 * *Observed:* Debug/ReleaseSafe: `panic: integer overflow` at `threading.zig:392`. `zig test -OReleaseFast` binary: *segfault*. Shipped `.so` (ReleaseFast): B=5 × {2048, 4096, 65536} correct with an intact canary after the output — LLVM's treatment of the UB happened to be benign in that compilation unit. It is undefined behaviour and build-dependent.
-* *Repro:* `zig test -OReleaseSafe --dep axiom -Mroot=repro_underflow_canary.zig -Maxiom=../../../../zig/src/threading.zig` (panic); `-OReleaseFast` (crash); `python3 canary_probe.py …libaxiom_zig.so` (benign through the `.so`); `python3 chunking_probe.py …` (B = 4…17 through the `.so`).
+* *Repro:* `zig test -OReleaseSafe --dep axiom -Mroot=repro_underflow_canary.zig -Maxiom=../../../../zig/src/threading.zig` (panic); `-OReleaseFast` (crash); `julia --startup-file=no canary_probe.jl …libaxiom_zig.so` (benign through the `.so`); `julia --startup-file=no chunking_probe.jl …` (B = 4…17 through the `.so`).
 * *Fix:* `last_count = batch_size -| last_start` (saturating) or derive per-thread ranges with `@min(start + chunk, batch_size)` as the element-wise `parallel_unary` already does (`:219`). Add B ∈ {1,…,9,13} to `zig build test` (Debug would have caught it).
 
 === J-05 `axiom_matmul_checked` cost `[E]` — High (performance / claim)
 
 * *Where:* `zig_ffi.jl:114` always calls `axiom_matmul_checked`; `zig/src/axiom.zig:101` pre-pass over `m·n·k` products.
 * *Observed:* checked vs unchecked: 0.311/0.035 ms (64²), 2.82/0.18 (128²), 24.8/2.17 (256²), 284.6/21.7 (512²). Single-thread OpenBLAS 512² = 2.3 ms. `benchmark/results_2026-02-20…` reports 25.6 ms at 512² — the unchecked kernel.
-* *Repro:* `python3 matmul_cost.py ../../../../zig/zig-out/lib/libaxiom_zig.so`.
+* *Repro:* `julia --startup-file=no matmul_cost.jl ../../../../zig/zig-out/lib/libaxiom_zig.so`.
 * *Fix:* check inputs in O(m·k + k·n) (vectorised `isfinite` scans) or the output once; keep the error contract; re-run the benchmark table.
 
 === J-06 `maxpool2d` padding not passed `[S]` — Medium
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/README.adoc b/docs/investigation/2026-09-26-silicon-core-recon/README.adoc
index c226f43..d1be439 100644
--- a/docs/investigation/2026-09-26-silicon-core-recon/README.adoc
+++ b/docs/investigation/2026-09-26-silicon-core-recon/README.adoc
@@ -30,7 +30,7 @@ dependency was added. No PR is opened by this work. Everything here lives under
 | H | link:H-staged-plan.adoc[Staged plan] | Phase 0 fixes + parity tests → Phase 1 internal `Runtime` seam → Phase 2 evidence gate → Phase 3 (conditional) extraction.
 | I | link:I-not-recommended.adoc[Changes explicitly NOT recommended] | No repo split now, no `SiliconCore.jl` from stubs, no GPU deps, no replacing reference kernels, no global warning suppression, no third backend hierarchy.
 | J | link:J-upstream-issues.adoc[Upstream issues with evidence] | 40 issues, each with severity, evidence type (empirical vs static), reproduction, and bucket (fix-in-Axiom / future-standalone / not-yet-verified).
-| — | link:repro/[repro/] | Minimal reproductions (Zig test files + Python/ctypes probes against the built `.so`).
+| — | link:repro/[repro/] | Minimal reproductions (Zig test files + Julia/ccall probes against the built `.so`).
 |===
 
 == How the evidence was obtained
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/README.adoc b/docs/investigation/2026-09-26-silicon-core-recon/repro/README.adoc
index 6547a3a..ae2274f 100644
--- a/docs/investigation/2026-09-26-silicon-core-recon/repro/README.adoc
+++ b/docs/investigation/2026-09-26-silicon-core-recon/repro/README.adoc
@@ -6,22 +6,34 @@ library first exactly as CI does (Zig 0.15.2):
 
     cd zig && zig build -Doptimize=ReleaseFast     # -> zig/zig-out/lib/libaxiom_zig.so
 
-Python probes need `numpy` (`ctypes` is standard library). Run them from this
-directory; the `.so` path argument is `../../../../zig/zig-out/lib/libaxiom_zig.so`.
+Julia probes use only standard libraries (no package installation required).
+Run them from this directory, for example:
+
+    julia --startup-file=no activation_probe.jl ../../../../zig/zig-out/lib/libaxiom_zig.so
+
+The probes retain the investigated shapes, C ABI calls, comparison tolerances,
+and subprocess isolation for batchnorm/chunking crashes. Normalisation buffers
+store each C row in a Julia column; pooling deliberately passes Julia's native
+column-major NHWC layout to reproduce the wrapper defect. Random samples use
+Julia's seeded generator, so exact mismatch counts and timings can differ from
+the original Python/NumPy measurements recorded in this investigation. Those
+historical measurements are unchanged; the Julia ports reproduce the same cases.
+A successful probe exit means the probe ran, not that the backend is correct;
+inspect its mismatch, NaN, canary, and child-process crash reports.
 
 [cols="2,3,2"]
 |===
 | File | What it shows | Issue
 
-| `batchnorm_probe.py` | `axiom_batchnorm` silently wrong for 4097–8191 features, SIGSEGV from 8192 (each case in a subprocess) | J-01
+| `batchnorm_probe.jl` | `axiom_batchnorm` silently wrong for 4097–8191 features, SIGSEGV from 8192 (each case in a subprocess) | J-01
 | `repro_batchnorm.zig` | Same defect as an index-out-of-bounds panic in ReleaseSafe | J-01
-| `layout_probe.py` | Pooling wrappers' column-major buffers vs row-major kernels → wrong values; `−Inf`/NaN window contracts | J-02, J-07
-| `activation_probe.py` | `axiom_tanh` NaN for x ≥ 44.5 / +Inf; near-zero cancellation; gelu/sigmoid fine | J-03
+| `layout_probe.jl` | Pooling wrappers' column-major buffers vs row-major kernels → wrong values; `−Inf`/NaN window contracts | J-02, J-07
+| `activation_probe.jl` | `axiom_tanh` NaN for x ≥ 44.5 / +Inf; near-zero cancellation; gelu/sigmoid fine | J-03
 | `repro_underflow.zig` | `parallel_layernorm` batch=5 integer-overflow panic (ReleaseSafe) | J-04
 | `repro_underflow_canary.zig` | Same function: `-OReleaseSafe` panics, `-OReleaseFast` test binary crashes | J-04
-| `canary_probe.py` | Through the shipped `.so`, batch=5 is (accidentally) benign: rows correct, canary intact | J-04
-| `chunking_probe.py` | softmax/layernorm for batch 4…17 × 2048 through the `.so` | J-04
-| `matmul_cost.py` | `axiom_matmul_checked` vs `axiom_matmul` timing (9–16×) | J-05
+| `canary_probe.jl` | Through the shipped `.so`, batch=5 is (accidentally) benign: rows correct, canary intact | J-04
+| `chunking_probe.jl` | softmax/layernorm for batch 4…17 × 2048 through the `.so` | J-04
+| `matmul_cost.jl` | `axiom_matmul_checked` vs `axiom_matmul` timing (9–16×) | J-05
 |===
 
 Zig-level reproductions:
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/activation_probe.jl b/docs/investigation/2026-09-26-silicon-core-recon/repro/activation_probe.jl
new file mode 100644
index 0000000..7694139
--- /dev/null
+++ b/docs/investigation/2026-09-26-silicon-core-recon/repro/activation_probe.jl
@@ -0,0 +1,40 @@
+# SPDX-License-Identifier: MPL-2.0
+# Probe Zig activation overflow/cancellation against Float64 reference formulas.
+include("probe_common.jl")
+
+function activation(symbol, x)
+  y = zeros(Float32, size(x))
+  ccall(Libdl.dlsym(lib, symbol), Cvoid,
+        (Ptr{Cfloat}, Ptr{Cfloat}, Csize_t), x, y, length(x))
+  return y
+end
+
+gelu_reference(x) = Float32(0.5 * x * (1 + tanh(sqrt(2 / pi) * (x + 0.044715 * x^3))))
+xs = Float32[1e-8, 1e-5, 1e-3, 5, 9, 10, 10.5, 11, 12, 20, 44, 44.3,
+              44.5, 45, 50, 100, -50, -100, Inf, -Inf]
+zt, zg, zs = activation(:axiom_tanh, xs), activation(:axiom_gelu, xs), activation(:axiom_sigmoid, xs)
+rt = Float32.(tanh.(Float64.(xs)))
+rg = gelu_reference.(Float64.(xs))
+rs = Float32.(1 ./ (1 .+ exp.(-Float64.(xs))))
+println("         x |     zig tanh     ref tanh |     zig gelu     ref gelu |  zig sigmoid          ref")
+for i in eachindex(xs)
+  flag = ""
+  !isfinite(zt[i]) && isfinite(rt[i]) && (flag *= " TANH-NaN")
+  !isfinite(zg[i]) && isfinite(rg[i]) && (flag *= " GELU-NaN")
+  if isfinite(zt[i]) && rt[i] != 0 && abs(zt[i] - rt[i]) / abs(rt[i]) > 1e-4
+    flag *= @sprintf(" tanh-relerr=%.1e", abs(zt[i] - rt[i]) / abs(rt[i]))
+  end
+  @printf("%10.4g | %12.6g %12.6g | %12.6g %12.6g | %12.6g %12.6g%s\n",
+          xs[i], zt[i], rt[i], zg[i], rg[i], zs[i], rs[i], flag)
+end
+xt = copy(xs)
+ccall(Libdl.dlsym(lib, :axiom_tanh_inplace), Cvoid, (Ptr{Cfloat}, Csize_t), xt, length(xt))
+println("\naxiom_tanh_inplace(50) = ", xt[findfirst(==(50), xs)],
+        " axiom_tanh_inplace(11) = ", xt[findfirst(==(11), xs)])
+rng = MersenneTwister(0)
+for sigma in (1, 3, 5, 8)
+  pre = Float32.(randn(rng, 1_000_000) .* sigma)
+  g = activation(:axiom_gelu, pre)
+  @printf("gelu on N(0,%d²) x1e6: NaN count = %d, max|x|=%.1f\n",
+          sigma, count(isnan, g), maximum(abs, pre))
+end
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/activation_probe.py b/docs/investigation/2026-09-26-silicon-core-recon/repro/activation_probe.py
deleted file mode 100644
index 8f886e9..0000000
--- a/docs/investigation/2026-09-26-silicon-core-recon/repro/activation_probe.py
+++ /dev/null
@@ -1,44 +0,0 @@
-# SPDX-License-Identifier: MPL-2.0
-"""Probe Zig activation kernels (built libaxiom_zig.so) for overflow / cancellation.
-Usage: PYTHONPATH=<numpy-env> python3 activation_probe.py zig/zig-out/lib/libaxiom_zig.so
-Compares against float32 reference formulas evaluated in float64 then rounded."""
-import ctypes, sys, numpy as np
-lib = ctypes.CDLL(sys.argv[1])
-f32p = ctypes.POINTER(ctypes.c_float)
-def call(sym, x):
-    x = np.ascontiguousarray(x, dtype=np.float32)
-    y = np.zeros_like(x)
-    fn = getattr(lib, sym)
-    fn.restype = None
-    fn.argtypes = [f32p, f32p, ctypes.c_size_t]
-    fn(x.ctypes.data_as(f32p), y.ctypes.data_as(f32p), x.size)
-    return y
-def gelu_ref(x):  # same tanh-approx formula Julia's gelu uses, evaluated in float64
-    x = x.astype(np.float64)
-    c = np.sqrt(2/np.pi)
-    return (0.5*x*(1+np.tanh(c*(x+0.044715*x**3)))).astype(np.float32)
-xs = np.array([1e-8, 1e-5, 1e-3, 5, 9, 10, 10.5, 11, 12, 20, 44, 44.3, 44.5, 45, 50, 100, -50, -100, float('inf'), float('-inf')], dtype=np.float32)
-print(f"{'x':>10} | {'zig tanh':>12} {'ref tanh':>12} | {'zig gelu':>12} {'ref gelu':>12} | {'zig sigmoid':>12} {'ref':>12}")
-zt, zg, zs = call('axiom_tanh', xs), call('axiom_gelu', xs), call('axiom_sigmoid', xs)
-rt = np.tanh(xs.astype(np.float64)).astype(np.float32)
-rg = gelu_ref(xs)
-rs = (1/(1+np.exp(-xs.astype(np.float64)))).astype(np.float32)
-for i, x in enumerate(xs):
-    flag = ""
-    if not np.isfinite(zt[i]) and np.isfinite(rt[i]): flag += " TANH-NaN"
-    if not np.isfinite(zg[i]) and np.isfinite(rg[i]): flag += " GELU-NaN"
-    if np.isfinite(zt[i]) and rt[i] != 0 and abs(zt[i]-rt[i])/abs(rt[i]) > 1e-4: flag += f" tanh-relerr={abs(zt[i]-rt[i])/abs(rt[i]):.1e}"
-    print(f"{x:>10.4g} | {zt[i]:>12.6g} {rt[i]:>12.6g} | {zg[i]:>12.6g} {rg[i]:>12.6g} | {zs[i]:>12.6g} {rs[i]:>12.6g}{flag}")
-# in-place variants use the same formula
-xt = xs.copy()
-fn = lib.axiom_tanh_inplace
-fn.restype=None
-fn.argtypes=[f32p, ctypes.c_size_t]
-fn(xt.ctypes.data_as(f32p), xt.size)
-print("\naxiom_tanh_inplace(50) =", xt[list(xs).index(np.float32(50))], " axiom_tanh_inplace(11)=", xt[list(xs).index(np.float32(11))])
-# realistic batch: how often does a standard-normal*sigma pre-activation trip GELU?
-rng = np.random.default_rng(0)
-for sigma in (1, 3, 5, 8):
-    pre = (rng.standard_normal(1_000_000)*sigma).astype(np.float32)
-    g = call('axiom_gelu', pre)
-    print(f"gelu on N(0,{sigma}²) x1e6: NaN count = {np.isnan(g).sum()}, max|x|={np.abs(pre).max():.1f}")
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/batchnorm_probe.jl b/docs/investigation/2026-09-26-silicon-core-recon/repro/batchnorm_probe.jl
new file mode 100644
index 0000000..e28c9bb
--- /dev/null
+++ b/docs/investigation/2026-09-26-silicon-core-recon/repro/batchnorm_probe.jl
@@ -0,0 +1,25 @@
+# SPDX-License-Identifier: MPL-2.0
+# Probe the fixed 4096-feature batchnorm scratch, isolating stack corruption.
+include("probe_common.jl")
+
+function child(features)
+  batch = 2
+  x = Float32.(randn(MersenneTwister(0), features, batch))
+  y = zeros(Float32, size(x))
+  gamma, beta = ones(Float32, features), zeros(Float32, features)
+  mean, variance = zeros(Float32, features), ones(Float32, features)
+  ccall(Libdl.dlsym(lib, :axiom_batchnorm), Cvoid,
+        (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},
+          Csize_t, Csize_t, Cfloat),
+        x, y, gamma, beta, mean, variance, length(x), features, 1f-5)
+  reference = (x .- mean) ./ sqrt.(variance .+ 1f-5)
+  println(count(!, close_element.(y, reference; atol=1e-4)), " of ", length(x))
+end
+
+if length(ARGS) > 1 && ARGS[2] == "--child"
+  child(parse(Int, ARGS[3]))
+else
+  for features in (4096, 4097, 4100, 8192, 65536, 1_000_000)
+    @printf("num_features=%8d: %s\n", features, run_case(features))
+  end
+end
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/batchnorm_probe.py b/docs/investigation/2026-09-26-silicon-core-recon/repro/batchnorm_probe.py
deleted file mode 100644
index d716b20..0000000
--- a/docs/investigation/2026-09-26-silicon-core-recon/repro/batchnorm_probe.py
+++ /dev/null
@@ -1,22 +0,0 @@
-# SPDX-License-Identifier: MPL-2.0
-"""axiom_batchnorm uses a fixed [4096]f32 stack scratch (norm.zig:143). Probe num_features > 4096
-through the production .so, each case in a subprocess (stack corruption may crash the process).
-Usage: PYTHONPATH=<numpy-env> python3 batchnorm_probe.py zig/zig-out/lib/libaxiom_zig.so"""
-import sys, subprocess
-child = r'''
-import ctypes, sys, numpy as np
-lib = ctypes.CDLL(sys.argv[1]); F = int(sys.argv[2]); B = 2
-f32p = ctypes.POINTER(ctypes.c_float)
-rng = np.random.default_rng(0)
-x = rng.standard_normal((B, F)).astype(np.float32); y = np.zeros_like(x)
-g = np.ones(F, np.float32); b = np.zeros(F, np.float32); m = np.zeros(F, np.float32); v = np.ones(F, np.float32)
-fn = lib.axiom_batchnorm; fn.restype = None
-fn.argtypes = [f32p, f32p, f32p, f32p, f32p, f32p, ctypes.c_size_t, ctypes.c_size_t, ctypes.c_float]
-fn(x.ctypes.data_as(f32p), y.ctypes.data_as(f32p), g.ctypes.data_as(f32p), b.ctypes.data_as(f32p), m.ctypes.data_as(f32p), v.ctypes.data_as(f32p), B * F, F, 1e-5)
-ref = (x - m) / np.sqrt(v + 1e-5)
-print(int((~np.isclose(y, ref, atol=1e-4)).sum()), "of", B * F)
-'''
-for F in (4096, 4097, 4100, 8192, 65536, 1_000_000):
-    r = subprocess.run([sys.executable, "-c", child, sys.argv[1], str(F)], capture_output=True, text=True)
-    out = r.stdout.strip() if r.returncode == 0 else f"CRASH exit={r.returncode} {r.stderr.strip().splitlines()[-1][:80] if r.stderr.strip() else '(signal)'}"
-    print(f"num_features={F:>8}: {out}")
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/canary_probe.jl b/docs/investigation/2026-09-26-silicon-core-recon/repro/canary_probe.jl
new file mode 100644
index 0000000..91abedb
--- /dev/null
+++ b/docs/investigation/2026-09-26-silicon-core-recon/repro/canary_probe.jl
@@ -0,0 +1,18 @@
+# SPDX-License-Identifier: MPL-2.0
+# Probe layernorm row correctness and an output canary for batch=5.
+include("probe_common.jl")
+
+for (batch, hidden) in ((5, 2048), (5, 4096), (5, 65536))
+  x = Float32.(randn(MersenneTwister(0), hidden, batch))
+  canary = 8 * hidden
+  ybuf = fill(12345f0, batch * hidden + canary)
+  gamma, beta = ones(Float32, hidden), zeros(Float32, hidden)
+  ccall(Libdl.dlsym(lib, :axiom_layernorm), Cvoid,
+        (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t, Csize_t, Cfloat),
+        x, ybuf, gamma, beta, batch, hidden, 1f-5)
+  y = reshape(view(ybuf, 1:batch * hidden), hidden, batch)
+  reference = layernorm_reference(x)
+  bad_rows = count(any(.!close_element.(y, reference; atol=1e-3); dims=1))
+  overwritten = count(!=(12345f0), view(ybuf, batch * hidden + 1:length(ybuf)))
+  println("B=$batch H=$hidden: bad_rows=$bad_rows canary_overwritten=$overwritten/$canary")
+end
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/canary_probe.py b/docs/investigation/2026-09-26-silicon-core-recon/repro/canary_probe.py
deleted file mode 100644
index 982c964..0000000
--- a/docs/investigation/2026-09-26-silicon-core-recon/repro/canary_probe.py
+++ /dev/null
@@ -1,21 +0,0 @@
-# SPDX-License-Identifier: MPL-2.0
-"""Call the production libaxiom_zig.so axiom_layernorm with B=5,H=2048 (the chunking-underflow case)
-and check (a) row correctness and (b) a canary region after the output buffer.
-Usage: PYTHONPATH=<numpy-env> python3 canary_probe.py zig/zig-out/lib/libaxiom_zig.so"""
-import ctypes, sys, numpy as np
-lib = ctypes.CDLL(sys.argv[1])
-f32p = ctypes.POINTER(ctypes.c_float)
-for (B, H) in [(5, 2048), (5, 4096), (5, 65536)]:
-    x = np.random.default_rng(0).standard_normal((B, H)).astype(np.float32)
-    canary = 8 * H
-    ybuf = np.full(B * H + canary, 12345.0, np.float32)   # output rows followed by canary
-    g = np.ones(H, np.float32)
-    b = np.zeros(H, np.float32)
-    fn = lib.axiom_layernorm
-    fn.restype = None
-    fn.argtypes = [f32p, f32p, f32p, f32p, ctypes.c_size_t, ctypes.c_size_t, ctypes.c_float]
-    fn(x.ctypes.data_as(f32p), ybuf.ctypes.data_as(f32p), g.ctypes.data_as(f32p), b.ctypes.data_as(f32p), B, H, 1e-5)
-    y = ybuf[:B*H].reshape(B, H)
-    can = ybuf[B*H:]
-    ref = (x - x.mean(1, keepdims=True)) / np.sqrt(x.var(1, keepdims=True) + 1e-5)
-    print(f"B={B} H={H}: bad_rows={int((~np.isclose(y, ref, atol=1e-3)).any(1).sum())} canary_overwritten={int((can != 12345.0).sum())}/{canary}")
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/chunking_probe.jl b/docs/investigation/2026-09-26-silicon-core-recon/repro/chunking_probe.jl
new file mode 100644
index 0000000..c6bcfd7
--- /dev/null
+++ b/docs/investigation/2026-09-26-silicon-core-recon/repro/chunking_probe.jl
@@ -0,0 +1,32 @@
+# SPDX-License-Identifier: MPL-2.0
+# Probe softmax/layernorm chunking, isolating potential ReleaseFast crashes.
+include("probe_common.jl")
+
+function child(batch, hidden, operation)
+  x = Float32.(randn(MersenneTwister(0), hidden, batch))
+  y = fill(7f0, size(x))
+  if operation == "softmax"
+    ccall(Libdl.dlsym(lib, :axiom_softmax), Cvoid,
+          (Ptr{Cfloat}, Ptr{Cfloat}, Csize_t, Csize_t), x, y, batch, hidden)
+    reference = exp.(x .- maximum(x; dims=1))
+    reference ./= sum(reference; dims=1)
+  elseif operation == "layernorm"
+    gamma, beta = ones(Float32, hidden), zeros(Float32, hidden)
+    ccall(Libdl.dlsym(lib, :axiom_layernorm), Cvoid,
+          (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t, Csize_t, Cfloat),
+          x, y, gamma, beta, batch, hidden, 1f-5)
+    reference = layernorm_reference(x)
+  else
+    error("Unknown operation: $operation")
+  end
+  println("mismatches=", count(!, close_element.(y, reference; atol=1e-4)))
+end
+
+if length(ARGS) > 1 && ARGS[2] == "--child"
+  child(parse(Int, ARGS[3]), parse(Int, ARGS[4]), ARGS[5])
+else
+  for operation in ("softmax", "layernorm"), batch in (4, 5, 6, 7, 8, 9, 13, 17)
+    hidden = 2048
+    @printf("%9s B=%2d H=%d: %s\n", operation, batch, hidden, run_case(batch, hidden, operation))
+  end
+end
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/chunking_probe.py b/docs/investigation/2026-09-26-silicon-core-recon/repro/chunking_probe.py
deleted file mode 100644
index a2c5663..0000000
--- a/docs/investigation/2026-09-26-silicon-core-recon/repro/chunking_probe.py
+++ /dev/null
@@ -1,29 +0,0 @@
-# SPDX-License-Identifier: MPL-2.0
-"""Which batch sizes trip the chunking underflow in parallel_layernorm/rmsnorm/softmax?
-Runs each case in a subprocess because the ReleaseFast .so may segfault / corrupt memory.
-Usage: PYTHONPATH=<numpy-env> python3 chunking_probe.py zig/zig-out/lib/libaxiom_zig.so"""
-import sys, subprocess, json
-lib = sys.argv[1]
-child = r'''
-import ctypes, sys, numpy as np
-lib = ctypes.CDLL(sys.argv[1]); B, H, op = int(sys.argv[2]), int(sys.argv[3]), sys.argv[4]
-f32p = ctypes.POINTER(ctypes.c_float)
-x = np.random.default_rng(0).standard_normal((B, H)).astype(np.float32); y = np.full((B, H), 7.0, np.float32)
-if op == "softmax":
-    fn = lib.axiom_softmax; fn.argtypes = [f32p, f32p, ctypes.c_size_t, ctypes.c_size_t]; fn.restype = None
-    fn(x.ctypes.data_as(f32p), y.ctypes.data_as(f32p), B, H)
-    ref = np.exp(x - x.max(1, keepdims=True)); ref /= ref.sum(1, keepdims=True)
-elif op == "layernorm":
-    g = np.ones(H, np.float32); b = np.zeros(H, np.float32)
-    fn = lib.axiom_layernorm; fn.argtypes = [f32p, f32p, f32p, f32p, ctypes.c_size_t, ctypes.c_size_t, ctypes.c_float]; fn.restype = None
-    fn(x.ctypes.data_as(f32p), y.ctypes.data_as(f32p), g.ctypes.data_as(f32p), b.ctypes.data_as(f32p), B, H, 1e-5)
-    m = x.mean(1, keepdims=True); v = x.var(1, keepdims=True); ref = (x - m) / np.sqrt(v + 1e-5)
-bad = int((~np.isclose(y, ref, atol=1e-4)).sum())
-print(bad)
-'''
-for op in ("softmax", "layernorm"):
-    for B in (4, 5, 6, 7, 8, 9, 13, 17):
-        H = 2048
-        r = subprocess.run([sys.executable, "-c", child, lib, str(B), str(H), op], capture_output=True, text=True)
-        status = f"exit={r.returncode}" + (f" mismatches={r.stdout.strip()}" if r.returncode == 0 else f" ({r.stderr.strip().splitlines()[-1][:60] if r.stderr.strip() else 'signal'})")
-        print(f"{op:>9} B={B:>2} H={H}: {status}")
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/layout_probe.jl b/docs/investigation/2026-09-26-silicon-core-recon/repro/layout_probe.jl
new file mode 100644
index 0000000..f25312f
--- /dev/null
+++ b/docs/investigation/2026-09-26-silicon-core-recon/repro/layout_probe.jl
@@ -0,0 +1,62 @@
+# SPDX-License-Identifier: MPL-2.0
+# Pass column-major NHWC buffers directly to row-major Zig pooling kernels,
+# as the investigated Julia wrappers did. Preserve this layout mismatch.
+include("probe_common.jl")
+
+function reference_maxpool(x, kernel, stride)
+  batch, height, width, channels = size(x)
+  out_h, out_w = (height - kernel) ÷ stride + 1, (width - kernel) ÷ stride + 1
+  y = Array{Float32}(undef, batch, out_h, out_w, channels)
+  for n in 1:batch, c in 1:channels, i in 1:out_h, j in 1:out_w
+    rows = (i - 1) * stride + 1:(i - 1) * stride + kernel
+    cols = (j - 1) * stride + 1:(j - 1) * stride + kernel
+    y[n, i, j, c] = maximum(view(x, n, rows, cols, c))
+  end
+  return y
+end
+
+function julia_style_maxpool(x, kernel, stride)
+  batch, height, width, channels = size(x)
+  out_h, out_w = (height - kernel) ÷ stride + 1, (width - kernel) ÷ stride + 1
+  y = zeros(Float32, batch, out_h, out_w, channels)
+  ccall(Libdl.dlsym(lib, :axiom_maxpool2d), Cvoid,
+        (Ptr{Cfloat}, Ptr{Cfloat}, Csize_t, Csize_t, Csize_t, Csize_t,
+          Csize_t, Csize_t, Csize_t, Csize_t),
+        x, y, batch, height, width, channels, kernel, kernel, stride, stride)
+  return y
+end
+
+function julia_style_gap(x)
+  batch, height, width, channels = size(x)
+  y = zeros(Float32, batch, channels)
+  ccall(Libdl.dlsym(lib, :axiom_global_avgpool2d), Cvoid,
+        (Ptr{Cfloat}, Ptr{Cfloat}, Csize_t, Csize_t, Csize_t, Csize_t),
+        x, y, batch, height, width, channels)
+  return y
+end
+
+rng = MersenneTwister(0)
+println("case                      | max|zig-ref| | mismatched elements")
+for (batch, height, width, channels, kernel, stride) in
+    ((1, 4, 4, 1, 2, 2), (1, 4, 4, 2, 2, 2), (2, 4, 4, 1, 2, 2),
+      (2, 6, 6, 3, 2, 2), (1, 5, 3, 1, 2, 1))
+  x = Float32.(randn(rng, batch, height, width, channels))
+  reference = reference_maxpool(x, kernel, stride)
+  got = julia_style_maxpool(x, kernel, stride)
+  bad = count(!, close_element.(reference, got))
+  @printf("maxpool N=%d H=%d W=%d C=%d k=%d s=%d | %10.4f | %d/%d\n",
+          batch, height, width, channels, kernel, stride, maximum(abs.(reference .- got)), bad, length(reference))
+end
+for (batch, height, width, channels) in ((1, 3, 3, 1), (1, 3, 3, 2), (2, 3, 3, 1), (2, 3, 3, 4))
+  x = Float32.(randn(rng, batch, height, width, channels))
+  reference = dropdims(sum(x; dims=(2, 3)); dims=(2, 3)) ./ (height * width)
+  got = julia_style_gap(x)
+  bad = count(!, close_element.(reference, got))
+  @printf("gap     N=%d H=%d W=%d C=%d         | %10.4f | %d/%d\n",
+          batch, height, width, channels, maximum(abs.(reference .- got)), bad, length(reference))
+end
+println("\nEdge-case contracts (single-channel windows so layout is not a factor):")
+x = fill(-Inf32, 1, 2, 2, 1)
+println("  all -Inf window -> zig=", only(julia_style_maxpool(x, 2, 2)), " julia(maximum)=-Inf")
+x = reshape(Float32[1, NaN, 0.5, 0.25], 1, 2, 2, 1)
+println("  window with NaN -> zig=", only(julia_style_maxpool(x, 2, 2)), " julia(maximum)=NaN")
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/layout_probe.py b/docs/investigation/2026-09-26-silicon-core-recon/repro/layout_probe.py
deleted file mode 100644
index f26ea79..0000000
--- a/docs/investigation/2026-09-26-silicon-core-recon/repro/layout_probe.py
+++ /dev/null
@@ -1,79 +0,0 @@
-# SPDX-License-Identifier: MPL-2.0
-"""
-Layout probe for the Julia -> Zig pooling wrappers.
-
-src/backends/zig_ffi.jl passes Julia's column-major (N,H,W,C) buffer *directly*
-to axiom_maxpool2d / axiom_global_avgpool2d (no _to_row_major_vec), while
-zig/src/pool.zig indexes row-major NHWC.  This script reproduces exactly what
-the Julia wrapper does (column-major bytes in, column-major interpretation of
-the output) and compares with the mathematically correct pooling result.
-
-Also probes: padding silently dropped, -Inf window handling, NaN handling.
-"""
-import ctypes, sys
-import numpy as np
-
-lib = ctypes.CDLL(sys.argv[1])
-f32p = ctypes.POINTER(ctypes.c_float)
-sz = ctypes.c_size_t
-
-lib.axiom_maxpool2d.argtypes = [f32p, f32p, sz, sz, sz, sz, sz, sz, sz, sz]
-lib.axiom_maxpool2d.restype = None
-lib.axiom_global_avgpool2d.argtypes = [f32p, f32p, sz, sz, sz, sz]
-lib.axiom_global_avgpool2d.restype = None
-
-def ref_maxpool(x, k, s):
-    N, H, W, C = x.shape
-    Ho, Wo = (H - k) // s + 1, (W - k) // s + 1
-    y = np.empty((N, Ho, Wo, C), np.float32)
-    for n in range(N):
-        for c in range(C):
-            for i in range(Ho):
-                for j in range(Wo):
-                    y[n, i, j, c] = x[n, i*s:i*s+k, j*s:j*s+k, c].max()
-    return y
-
-def julia_style_maxpool(x, k, s):
-    """Mimic backend_maxpool2d(::ZigBackend): pass column-major bytes as-is,
-    read output back as column-major (N,Ho,Wo,C)."""
-    N, H, W, C = x.shape
-    Ho, Wo = (H - k) // s + 1, (W - k) // s + 1
-    xin = np.asfortranarray(x)                       # Julia memory order
-    out = np.zeros(N*Ho*Wo*C, np.float32)
-    lib.axiom_maxpool2d(xin.ctypes.data_as(f32p), out.ctypes.data_as(f32p),
-                        N, H, W, C, k, k, s, s)
-    return out.reshape((N, Ho, Wo, C), order="F")   # Julia reads back column-major
-
-def julia_style_gap(x):
-    N, H, W, C = x.shape
-    xin = np.asfortranarray(x)
-    out = np.zeros(N*C, np.float32)
-    lib.axiom_global_avgpool2d(xin.ctypes.data_as(f32p), out.ctypes.data_as(f32p), N, H, W, C)
-    return out.reshape((N, C), order="F")
-
-rng = np.random.default_rng(0)
-print("case                      | max|zig-ref| | mismatched elements")
-for (N, H, W, C, k, s) in [(1, 4, 4, 1, 2, 2), (1, 4, 4, 2, 2, 2), (2, 4, 4, 1, 2, 2), (2, 6, 6, 3, 2, 2), (1, 5, 3, 1, 2, 1)]:
-    x = rng.standard_normal((N, H, W, C)).astype(np.float32)
-    ref = ref_maxpool(x, k, s)
-    got = julia_style_maxpool(x, k, s)
-    bad = int((~np.isclose(ref, got)).sum())
-    print(f"maxpool N={N} H={H} W={W} C={C} k={k} s={s} | {np.abs(ref-got).max():10.4f} | {bad}/{ref.size}")
-
-for (N, H, W, C) in [(1, 3, 3, 1), (1, 3, 3, 2), (2, 3, 3, 1), (2, 3, 3, 4)]:
-    x = rng.standard_normal((N, H, W, C)).astype(np.float32)
-    ref = x.mean(axis=(1, 2))
-    got = julia_style_gap(x)
-    bad = int((~np.isclose(ref, got)).sum())
-    print(f"gap     N={N} H={H} W={W} C={C}         | {np.abs(ref-got).max():10.4f} | {bad}/{ref.size}")
-
-print()
-print("Edge-case contracts (row-major single-channel so layout is not a factor):")
-x = np.full((1, 2, 2, 1), -np.inf, np.float32)
-out = np.zeros(1, np.float32)
-lib.axiom_maxpool2d(np.ascontiguousarray(x).ctypes.data_as(f32p), out.ctypes.data_as(f32p), 1, 2, 2, 1, 2, 2, 2, 2)
-print(f"  all -Inf window  -> zig={out[0]!r}  julia(maximum)=-Inf   (zig initialises max at -floatmax)")
-x = np.array([[1.0, np.nan], [0.5, 0.25]], np.float32).reshape(1, 2, 2, 1)
-out = np.zeros(1, np.float32)
-lib.axiom_maxpool2d(np.ascontiguousarray(x).ctypes.data_as(f32p), out.ctypes.data_as(f32p), 1, 2, 2, 1, 2, 2, 2, 2)
-print(f"  window with NaN  -> zig={out[0]!r}  julia(maximum)=NaN    (zig `val > max` skips NaN)")
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/matmul_cost.jl b/docs/investigation/2026-09-26-silicon-core-recon/repro/matmul_cost.jl
new file mode 100644
index 0000000..950f9d8
--- /dev/null
+++ b/docs/investigation/2026-09-26-silicon-core-recon/repro/matmul_cost.jl
@@ -0,0 +1,41 @@
+# SPDX-License-Identifier: MPL-2.0
+# Measure the finiteness pre-pass in checked matmul against the SIMD kernel.
+include("probe_common.jl")
+using LinearAlgebra
+
+function matmul_unchecked(a, b, c, n)
+  ccall(Libdl.dlsym(lib, :axiom_matmul), Cvoid,
+        (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t, Csize_t, Csize_t), a, b, c, n, n, n)
+end
+
+function matmul_checked(a, b, c, n)
+  status = ccall(Libdl.dlsym(lib, :axiom_matmul_checked), UInt32,
+                (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t, Csize_t, Csize_t), a, b, c, n, n, n)
+  status == 0 || error("axiom_matmul_checked returned status $status")
+end
+
+function bench(fn, reps)
+  fn() # Warm up compilation before timing Julia closures or BLAS.
+  best = Inf
+  for _ in 1:reps
+    start = time_ns()
+    fn()
+    best = min(best, (time_ns() - start) / 1e6)
+  end
+  return best
+end
+
+rng = MersenneTwister(1)
+println("  m=k=n | unchecked ms | checked ms | ratio | Julia(BLAS) ms")
+for n in (64, 128, 256, 512)
+  a, b = Float32.(randn(rng, n, n)), Float32.(randn(rng, n, n))
+  # Explicit row-major copies for the C ABI; convert output back before comparing.
+  ar, br, cr = vec(permutedims(a)), vec(permutedims(b)), zeros(Float32, n * n)
+  reps = n >= 512 ? 5 : 20
+  tu = bench(() -> matmul_unchecked(ar, br, cr, n), reps)
+  tc = bench(() -> matmul_checked(ar, br, cr, n), reps)
+  tn = bench(() -> a * b, reps)
+  got = permutedims(reshape(cr, n, n))
+  ok = all(close_element.(got, a * b; atol=1e-3, rtol=1e-4))
+  @printf("%7d | %12.3f | %10.3f | %5.2f | %.3f (result matches BLAS: %s)\n", n, tu, tc, tc / tu, tn, ok)
+end
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/matmul_cost.py b/docs/investigation/2026-09-26-silicon-core-recon/repro/matmul_cost.py
deleted file mode 100644
index 79f3346..0000000
--- a/docs/investigation/2026-09-26-silicon-core-recon/repro/matmul_cost.py
+++ /dev/null
@@ -1,34 +0,0 @@
-# SPDX-License-Identifier: MPL-2.0
-"""Cost of axiom_matmul_checked vs axiom_matmul (the Julia wrapper always uses
-the checked one).  The checked variant runs a scalar O(m*n*k) finiteness
-pre-pass (matmulCellFinite) before the tiled SIMD kernel."""
-import ctypes, sys, time
-import numpy as np
-
-lib = ctypes.CDLL(sys.argv[1])
-f32p = ctypes.POINTER(ctypes.c_float)
-sz = ctypes.c_size_t
-lib.axiom_matmul.argtypes = [f32p, f32p, f32p, sz, sz, sz]
-lib.axiom_matmul.restype = None
-lib.axiom_matmul_checked.argtypes = [f32p, f32p, f32p, sz, sz, sz]
-lib.axiom_matmul_checked.restype = ctypes.c_uint32
-
-rng = np.random.default_rng(1)
-print(f"{'m=k=n':>7} | {'unchecked ms':>12} | {'checked ms':>10} | ratio | numpy(BLAS) ms")
-for n in (64, 128, 256, 512):
-    a = np.ascontiguousarray(rng.standard_normal((n, n)).astype(np.float32))
-    b = np.ascontiguousarray(rng.standard_normal((n, n)).astype(np.float32))
-    c = np.zeros((n, n), np.float32)
-    reps = 5 if n >= 512 else 20
-    def bench(fn):
-        best = 1e9
-        for _ in range(reps):
-            t = time.perf_counter()
-            fn()
-            best = min(best, time.perf_counter() - t)
-        return best * 1e3
-    tu = bench(lambda: lib.axiom_matmul(a.ctypes.data_as(f32p), b.ctypes.data_as(f32p), c.ctypes.data_as(f32p), n, n, n))
-    tc = bench(lambda: lib.axiom_matmul_checked(a.ctypes.data_as(f32p), b.ctypes.data_as(f32p), c.ctypes.data_as(f32p), n, n, n))
-    tn = bench(lambda: a @ b)
-    ok = np.allclose(c, a @ b, atol=1e-3, rtol=1e-4)
-    print(f"{n:>7} | {tu:12.3f} | {tc:10.3f} | {tc/tu:5.2f} | {tn:.3f}   (result matches BLAS: {ok})")
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/probe_common.jl b/docs/investigation/2026-09-26-silicon-core-recon/repro/probe_common.jl
new file mode 100644
index 0000000..654f6c9
--- /dev/null
+++ b/docs/investigation/2026-09-26-silicon-core-recon/repro/probe_common.jl
@@ -0,0 +1,30 @@
+# SPDX-License-Identifier: MPL-2.0
+# Shared standard-library support for the standalone investigation probes.
+using Libdl, Random, Printf
+
+isempty(ARGS) && error("Usage: julia <probe>.jl /path/to/libaxiom_zig.so")
+const lib = Libdl.dlopen(abspath(ARGS[1]))
+
+# Match NumPy's elementwise isclose defaults (not Julia's norm-based isapprox).
+close_element(a, b; atol=1e-8, rtol=1e-5) =
+  a == b || (isfinite(a) && isfinite(b) && abs(a - b) <= atol + rtol * abs(b))
+
+# Columns are independent rows of the C row-major (batch, hidden) buffer.
+function layernorm_reference(x)
+  centered = x .- sum(x; dims=1) ./ size(x, 1)
+  variance = sum(abs2, centered; dims=1) ./ size(x, 1)
+  return centered ./ sqrt.(variance .+ 1f-5)
+end
+
+function run_case(args...)
+  stdout_buffer, stderr_buffer = IOBuffer(), IOBuffer()
+  command = `$(Base.julia_cmd()) --startup-file=no $(abspath(PROGRAM_FILE)) $(ARGS[1]) --child $args`
+  process = run(pipeline(ignorestatus(command); stdout=stdout_buffer, stderr=stderr_buffer))
+  output = strip(String(take!(stdout_buffer)))
+  errors = strip(String(take!(stderr_buffer)))
+  if success(process)
+    return "exit=0 " * output
+  end
+  detail = isempty(errors) ? "signal" : first(last(split(errors, '\n')), 80)
+  return "CRASH exit=$(process.exitcode) signal=$(process.termsignal) ($detail)"
+end
diff --git a/docs/investigation/2026-09-26-silicon-core-recon/repro/repro_batchnorm.zig b/docs/investigation/2026-09-26-silicon-core-recon/repro/repro_batchnorm.zig
index 79d1577..5ce575c 100644
--- a/docs/investigation/2026-09-26-silicon-core-recon/repro/repro_batchnorm.zig
+++ b/docs/investigation/2026-09-26-silicon-core-recon/repro/repro_batchnorm.zig
@@ -3,7 +3,7 @@
 // indexed by `num_features` with no bound check.  num_features > 4096 writes
 // past the stack array.  In ReleaseSafe/Debug this is an index-out-of-bounds
 // panic; through the shipped ReleaseFast .so it is silently wrong (4097..8191)
-// and SIGSEGV from 8192 features (see batchnorm_probe.py).
+// and SIGSEGV from 8192 features (see batchnorm_probe.jl).
 // Run: zig test -OReleaseSafe --dep axiom -Mroot=repro_batchnorm.zig -Maxiom=../../../../zig/src/axiom.zig
 const std = @import("std");
 const axiom = @import("axiom");