From 90f1128d9b67921b7294ecd968dbdca7f39fe976 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Wed, 23 Sep 2026 20:21:02 +0000
Subject: [PATCH 1/8] =?UTF-8?q?statistics:=20moments=5Fu32=20+=20Cascade::?=
 =?UTF-8?q?observe=5Fbatch=20=E2=80=94=20batch=20Welford?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`moments_u32(&[u32]) -> MomentsU32 { n, sum, sum_sq }`: exact integer
moments, eight lanes at a time through U64x8. Squares are split into
32-bit halves so every lane add stays below 2^32; lanes drain into u128
totals every 2^28 chunks, before an 8-lane reduce_sum could overflow u64.
`MomentsU32::merge` is plain integer addition — exact, associative and
commutative — so shard-parallel reductions combine to the same result as
one pass. `mean()`/`variance()` form n*sum_sq - sum^2 exactly in u128 when
it fits (always, at Hamming scale) and only round in the final division.

`Cascade::observe_batch(&[u32])` folds a whole batch into the rolling
floor with the Chan-Golub-LeVeque parallel-variance merge: same mu /
sigma / observations as per-element `observe` up to f64 rounding, and
shards may be folded in any order. Drift is judged once per batch (μ moved
more than 2σ of the pre-batch state), the test `observe` applies per step.

Tests: exact vs a u128 reference across chunk boundaries and full-range
values, u32::MAX worst case (100 003 squares of (2^32-1)^2), merge
exactness at every split in both orders, mean/variance vs two-pass f64;
observe_batch vs sequential observe (cold and warm), shard order
independence, and the alert firing on a shifted batch but not a stable one.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm
---
 src/hpc/cascade.rs    | 120 +++++++++++++++++++++++++
 src/hpc/statistics.rs | 200 ++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 320 insertions(+)

diff --git a/src/hpc/cascade.rs b/src/hpc/cascade.rs
index 4673c10d..4fb1d99b 100644
--- a/src/hpc/cascade.rs
+++ b/src/hpc/cascade.rs
@@ -208,6 +208,56 @@ impl Cascade {
         }
     }
 
+    /// Fold a whole batch of distances into the rolling floor at once.
+    ///
+    /// The batch is reduced to exact integer moments with
+    /// [`moments_u32`](crate::hpc::statistics::moments_u32) and merged into
+    /// the running `(n, μ, σ)` with the parallel-variance merge
+    /// (Chan, Golub & LeVeque), so the resulting `mu`/`sigma`/`observations`
+    /// match calling [`observe`](Self::observe) once per distance up to f64
+    /// rounding. Shards computed on different threads can each be folded in
+    /// this way, in any order.
+    ///
+    /// Drift is judged once per batch, not per element: an alert fires when
+    /// the batch moves μ by more than 2σ of the pre-batch state (and the state
+    /// had > 10 observations and σ > 0), the same test `observe` applies to a
+    /// single step.
+    pub fn observe_batch(&mut self, distances: &[u32]) -> Option<ShiftAlert> {
+        let b = crate::hpc::statistics::moments_u32(distances);
+        if b.n == 0 {
+            return None;
+        }
+        let old_mu = self.mu;
+        let old_sigma = self.sigma;
+        let old_n = self.observations;
+        let n_a = old_n as f64;
+        let n_b = b.n as f64;
+        let n = n_a + n_b;
+        let (mean_b, m2_b) = (b.mean(), b.variance() * n_b);
+        if old_n == 0 {
+            self.mu = mean_b;
+            self.sigma = (m2_b / n_b).sqrt();
+        } else {
+            let delta = mean_b - old_mu;
+            self.mu = old_mu + delta * n_b / n;
+            let m2 = old_sigma * old_sigma * n_a + m2_b + delta * delta * n_a * n_b / n;
+            self.sigma = (m2 / n).sqrt();
+        }
+        self.observations = old_n + b.n as usize;
+
+        if old_n > 10 && old_sigma > 0.0 && (self.mu - old_mu).abs() > 2.0 * old_sigma {
+            Some(ShiftAlert {
+                old_mu,
+                new_mu: self.mu,
+                old_sigma,
+                new_sigma: self.sigma,
+                observations: self.observations,
+            })
+        } else {
+            None
+        }
+    }
+
     pub fn recalibrate(&mut self, alert: &ShiftAlert) {
         self.mu = alert.new_mu;
         self.sigma = alert.new_sigma;
@@ -768,6 +818,76 @@ mod tests {
         assert_eq!(got, exact_hits(&expected, threshold));
     }
 
+    fn close(a: f64, b: f64) -> bool {
+        (a - b).abs() <= 1e-9 * a.abs().max(b.abs()).max(1.0)
+    }
+
+    fn noisy(n: usize, base: u32, spread: u32, mut s: u64) -> Vec<u32> {
+        (0..n)
+            .map(|_| {
+                s ^= s << 13;
+                s ^= s >> 7;
+                s ^= s << 17;
+                base + (s as u32) % spread
+            })
+            .collect()
+    }
+
+    /// One `observe_batch` lands on the same rolling floor as `observe`
+    /// called per distance, from an empty state and from a warm one.
+    #[test]
+    fn observe_batch_matches_sequential_observe() {
+        let warm = noisy(300, 8000, 400, 1);
+        let batch = noisy(2000, 8100, 500, 2);
+        for start in [&[][..], &warm[..]] {
+            let mut seq = Cascade::from_threshold(8000, 2048);
+            let mut bat = Cascade::from_threshold(8000, 2048);
+            for &d in start {
+                seq.observe(d);
+                bat.observe(d);
+            }
+            for &d in &batch {
+                seq.observe(d);
+            }
+            bat.observe_batch(&batch);
+            assert_eq!(seq.observations(), bat.observations());
+            assert!(close(seq.mu(), bat.mu()), "mu {} vs {}", seq.mu(), bat.mu());
+            assert!(close(seq.sigma(), bat.sigma()), "sigma {} vs {}", seq.sigma(), bat.sigma());
+        }
+    }
+
+    /// Shard-parallel use: folding shards in any order gives the same floor
+    /// as folding the whole batch.
+    #[test]
+    fn observe_batch_shards_merge_in_any_order() {
+        let x = noisy(3001, 8000, 700, 3);
+        let mut whole = Cascade::from_threshold(8000, 2048);
+        whole.observe_batch(&x);
+        let shards: Vec<&[u32]> = x.chunks(640).collect();
+        for order in [[0usize, 1, 2, 3, 4], [4, 2, 0, 3, 1]] {
+            let mut c = Cascade::from_threshold(8000, 2048);
+            for i in order {
+                c.observe_batch(shards[i]);
+            }
+            assert!(close(c.mu(), whole.mu()));
+            assert!(close(c.sigma(), whole.sigma()));
+            assert_eq!(c.observations(), whole.observations());
+        }
+    }
+
+    /// The drift alert can fire (a batch from a distribution shifted far
+    /// past 2σ) and stays silent on a batch from the same distribution.
+    #[test]
+    fn observe_batch_alerts_on_a_shift_only() {
+        let mut c = Cascade::from_threshold(8000, 2048);
+        assert!(c.observe_batch(&noisy(500, 8000, 100, 4)).is_none(), "first batch has no prior");
+        assert!(c.observe_batch(&noisy(500, 8000, 100, 5)).is_none(), "same distribution");
+        let alert = c
+            .observe_batch(&noisy(5000, 9000, 100, 6))
+            .expect("shifted batch must alert");
+        assert!(alert.new_mu > alert.old_mu + 2.0 * alert.old_sigma);
+    }
+
     #[test]
     fn packed_database_roundtrip() {
         let vec_bytes = 256;
diff --git a/src/hpc/statistics.rs b/src/hpc/statistics.rs
index 14d97010..b08b0124 100644
--- a/src/hpc/statistics.rs
+++ b/src/hpc/statistics.rs
@@ -360,3 +360,203 @@ mod tests {
         }
     }
 }
+
+// ── Batch moments for shard-parallel Welford ──────────────────────────────
+
+/// Exact first and second moments of a `u32` sample: count, `Σx` and `Σx²`.
+///
+/// These three integers are the sufficient statistics of a Welford rolling
+/// floor. Unlike a running `(mean, M2)` pair they merge by plain addition, so
+/// [`MomentsU32::merge`] is exact, associative and commutative: shards of a
+/// sample can be reduced in parallel, in any order, and combined into the
+/// same result as one sequential pass. Floats appear only in [`mean`] and
+/// [`variance`], at the end.
+///
+/// [`mean`]: MomentsU32::mean
+/// [`variance`]: MomentsU32::variance
+#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
+pub struct MomentsU32 {
+    /// Number of values.
+    pub n: u64,
+    /// `Σx`.
+    pub sum: u128,
+    /// `Σx²`.
+    pub sum_sq: u128,
+}
+
+impl MomentsU32 {
+    /// Moments of the union of two samples — exact integer addition.
+    #[inline]
+    #[must_use]
+    pub fn merge(self, other: Self) -> Self {
+        Self {
+            n: self.n + other.n,
+            sum: self.sum + other.sum,
+            sum_sq: self.sum_sq + other.sum_sq,
+        }
+    }
+
+    /// Sample mean; `0.0` for an empty sample.
+    pub fn mean(&self) -> f64 {
+        if self.n == 0 {
+            0.0
+        } else {
+            self.sum as f64 / self.n as f64
+        }
+    }
+
+    /// Population variance `E[(X - μ)²]`; `0.0` for an empty sample.
+    ///
+    /// Computed as `(n·Σx² − (Σx)²) / n²`. The numerator is formed exactly in
+    /// `u128` whenever it fits (always, for Hamming-scale data: n ≤ 2³², x ≤
+    /// 2¹⁷), so there is no cancellation between two large floats; only the
+    /// final division rounds. Past that range it falls back to `f64`.
+    pub fn variance(&self) -> f64 {
+        if self.n == 0 {
+            return 0.0;
+        }
+        let n = u128::from(self.n);
+        match (n.checked_mul(self.sum_sq), self.sum.checked_mul(self.sum)) {
+            (Some(a), Some(b)) => (a - b) as f64 / (self.n as f64 * self.n as f64),
+            _ => {
+                let mean = self.mean();
+                (self.sum_sq as f64 / self.n as f64 - mean * mean).max(0.0)
+            }
+        }
+    }
+}
+
+/// Exact [`MomentsU32`] of `values`, eight lanes at a time through `U64x8`.
+///
+/// Each value is widened to `u64` and its square split into 32-bit halves, so
+/// every lane add is below 2³², and the lanes are drained into `u128` totals
+/// every 2²⁸ chunks, before an 8-lane reduction could overflow `u64`. The
+/// widening multiply is plain lane-wise Rust that the compiler vectorizes
+/// (`vpmuludq` on x86); the accumulation and the final reduction go through
+/// the polyfill, so every backend runs the same code.
+///
+/// # Example
+///
+/// ```
+/// use ndarray::hpc::statistics::moments_u32;
+///
+/// let m = moments_u32(&[1, 2, 3, 4]);
+/// assert_eq!((m.n, m.sum, m.sum_sq), (4, 10, 30));
+/// assert_eq!(m.mean(), 2.5);
+/// assert_eq!(m.variance(), 1.25);
+/// ```
+pub fn moments_u32(values: &[u32]) -> MomentsU32 {
+    use crate::simd::U64x8;
+
+    /// Chunks per drain. Every lane add is below 2^32, so after 2^28 chunks
+    /// a lane is below 2^60 and the 8-lane `reduce_sum` below 2^63 — the
+    /// reduction, not the lane, is the binding limit.
+    const DRAIN: usize = 1 << 28;
+
+    let (chunks, tail) = values.as_chunks::<8>();
+    let mut out = MomentsU32 {
+        n: values.len() as u64,
+        ..MomentsU32::default()
+    };
+    for block in chunks.chunks(DRAIN) {
+        let (mut s, mut lo, mut hi) = (U64x8::splat(0), U64x8::splat(0), U64x8::splat(0));
+        for c in block {
+            let x: [u64; 8] = core::array::from_fn(|i| u64::from(c[i]));
+            let sq: [u64; 8] = core::array::from_fn(|i| x[i] * x[i]);
+            s += U64x8::from_array(x);
+            lo += U64x8::from_array(core::array::from_fn(|i| sq[i] & 0xFFFF_FFFF));
+            hi += U64x8::from_array(core::array::from_fn(|i| sq[i] >> 32));
+        }
+        out.sum += u128::from(s.reduce_sum());
+        out.sum_sq += u128::from(lo.reduce_sum()) + (u128::from(hi.reduce_sum()) << 32);
+    }
+    for &v in tail {
+        out.sum += u128::from(v);
+        out.sum_sq += u128::from(v) * u128::from(v);
+    }
+    out
+}
+
+#[cfg(test)]
+mod moments_tests {
+    use super::*;
+
+    fn xorshift(n: usize, mut s: u64, mask: u32) -> Vec<u32> {
+        (0..n)
+            .map(|_| {
+                s ^= s << 13;
+                s ^= s >> 7;
+                s ^= s << 17;
+                (s as u32) & mask
+            })
+            .collect()
+    }
+
+    fn reference(x: &[u32]) -> MomentsU32 {
+        MomentsU32 {
+            n: x.len() as u64,
+            sum: x.iter().map(|&v| u128::from(v)).sum(),
+            sum_sq: x.iter().map(|&v| u128::from(v) * u128::from(v)).sum(),
+        }
+    }
+
+    /// Exact against a u128 reference at every length across the 8-lane
+    /// chunk boundary, for small (Hamming-scale) and full-range values.
+    #[test]
+    fn moments_u32_is_exact() {
+        for mask in [0x3FFF, u32::MAX] {
+            let x = xorshift(1000, 0x9E37_79B9_7F4A_7C15, mask);
+            for n in (0..=40).chain([63, 64, 65, 999, 1000]) {
+                assert_eq!(moments_u32(&x[..n]), reference(&x[..n]), "n={n} mask={mask:#x}");
+            }
+        }
+    }
+
+    /// Worst case for the square accumulator: every square is (2^32-1)^2,
+    /// which alone nearly fills a u64. A lane that summed whole squares
+    /// would overflow on the second element.
+    #[test]
+    fn moments_u32_does_not_overflow_at_u32_max() {
+        let x = vec![u32::MAX; 100_003];
+        let m = moments_u32(&x);
+        let v = u128::from(u32::MAX);
+        assert_eq!(m.n, 100_003);
+        assert_eq!(m.sum, 100_003 * v);
+        assert_eq!(m.sum_sq, 100_003 * v * v);
+    }
+
+    /// Merging is exact integer addition: the moments of a concatenation
+    /// equal the merge of the parts' moments at any split point and in
+    /// either order — which is what makes shard-parallel statistics exact.
+    #[test]
+    fn merge_is_exact_and_order_independent() {
+        let x = xorshift(777, 42, u32::MAX);
+        let whole = moments_u32(&x);
+        for split in [0, 1, 7, 8, 9, 400, 776, 777] {
+            let (a, b) = x.split_at(split);
+            assert_eq!(moments_u32(a).merge(moments_u32(b)), whole, "split={split}");
+            assert_eq!(moments_u32(b).merge(moments_u32(a)), whole, "reversed split={split}");
+        }
+    }
+
+    /// Mean and population variance against a two-pass f64 reference.
+    #[test]
+    fn mean_and_variance_match_two_pass() {
+        for mask in [0x3FFF, u32::MAX] {
+            let x = xorshift(5000, 7, mask);
+            let n = x.len() as f64;
+            let mean = x.iter().map(|&v| f64::from(v)).sum::<f64>() / n;
+            let var = x
+                .iter()
+                .map(|&v| (f64::from(v) - mean).powi(2))
+                .sum::<f64>()
+                / n;
+            let m = moments_u32(&x);
+            assert!((m.mean() - mean).abs() <= 1e-9 * mean.abs(), "mean mask={mask:#x}");
+            assert!((m.variance() - var).abs() <= 1e-9 * var, "variance mask={mask:#x}");
+        }
+        assert_eq!(MomentsU32::default().mean(), 0.0);
+        assert_eq!(MomentsU32::default().variance(), 0.0);
+        assert_eq!(moments_u32(&[5, 5, 5]).variance(), 0.0, "constant input has zero variance");
+    }
+}

From 414aa8e92a6766400ce2a090f7531649e614f33f Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Wed, 23 Sep 2026 20:26:58 +0000
Subject: [PATCH 2/8] simd: re-export moments_u32 / MomentsU32 from the facade

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm
---
 src/simd.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/simd.rs b/src/simd.rs
index 172e5797..f7354a3f 100644
--- a/src/simd.rs
+++ b/src/simd.rs
@@ -617,6 +617,7 @@ pub use crate::hpc::fft::{wht_f32, wht_f32_new};
 pub use crate::hpc::fingerprint::{
     vector_config, Fingerprint, Fingerprint1K, Fingerprint2K, Fingerprint64K, VectorConfig, VectorWidth,
 };
+pub use crate::hpc::statistics::{moments_u32, MomentsU32};
 
 // PR-X1 — SoA carrier + const-size slice helpers, dispatched from their
 // respective `simd_{type}.rs` modules. The W1a consumer contract forbids

From 7bb45fbcfc8f6e81abb40f3c58006dd72b627043 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Wed, 23 Sep 2026 20:39:28 +0000
Subject: [PATCH 3/8] zspace: the Fisher-Z entry point and the binomial Hamming
 null
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Operator rule: anything cosine-shaped pays the Fisher-Z entry tax before
it is thresholded, averaged, banded or ranked by margin — no exceptions,
lab and calibration included; anything bitpacked goes through HDR popcount
stacking on the Hamming distance itself. ndarray had no Fisher transform at
all; this adds the one canonical entry point.

- fisher_z / fisher_z_inv / hyperbolic_depth: the helix convention
  (`helix::fisher_z::Similarity`, mirrored by `jc::stats::fisher_2z`):
  ln form, r clamped to [-1+1e-9, 1-1e-9], NaN propagates.
- fisher_z_f32_batch: 16 lanes through F32x16 + simd_ln_f32, bit-identical
  to the scalar f32 formula. The f32 rim is 1 - 2^-24 (FISHER_CLAMP_F32):
  1 - 1e-9 rounds to exactly 1.0 in f32, where ln(1 - r) is -inf.
- hamming_null_z: (d - n/2) / (sqrt(n)/2) against Binomial(n, 1/2) — the
  random-vector null needs no calibration samples.

All re-exported from ndarray::simd. Tests pin atanh(0.5) and the 1e-9 rim
value, check oddness, monotonicity, finiteness past the rim, NaN
propagation, the inverse, batch/scalar bit identity across the 16-lane
boundary, and the binomial null at 0, -1 and +3 sigma.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm
---
 src/hpc/mod.rs    |   3 +
 src/hpc/zspace.rs | 237 ++++++++++++++++++++++++++++++++++++++++++++++
 src/simd.rs       |   1 +
 3 files changed, 241 insertions(+)
 create mode 100644 src/hpc/zspace.rs

diff --git a/src/hpc/mod.rs b/src/hpc/mod.rs
index 074921b6..8067cf8f 100644
--- a/src/hpc/mod.rs
+++ b/src/hpc/mod.rs
@@ -25,6 +25,9 @@ pub mod blas_level2;
 pub mod blas_level3;
 pub mod reductions;
 pub mod statistics;
+/// z-space entry points: Fisher-Z for cosine-shaped values, the binomial
+/// null for bitpacked Hamming distances.
+pub mod zspace;
 /// Reliability & validity statistics: Pearson r, Spearman ρ, Cronbach α, ICC.
 pub mod reliability;
 /// Entropy ladder: Staunen↔Wisdom coordinate over NARS truth + Pearl-2³ SPO.
diff --git a/src/hpc/zspace.rs b/src/hpc/zspace.rs
new file mode 100644
index 00000000..3171dae1
--- /dev/null
+++ b/src/hpc/zspace.rs
@@ -0,0 +1,237 @@
+//! The z-space entry points: every similarity or distance enters the
+//! substrate as a z-score before anything thresholds, averages, bands or
+//! ranks it by margin.
+//!
+//! Two doors, no third:
+//!
+//! - **Cosine-shaped** values (cosine, Pearson r, normalized dot products —
+//!   anything bounded in `[-1, 1]`) pay the Fisher-Z entry tax:
+//!   [`fisher_z`] / [`fisher_z_f32_batch`]. No exceptions, including lab and
+//!   calibration code. A raw cosine is not a z-score: its variance shrinks
+//!   towards the rim, so equal steps do not mean equal evidence and a fixed
+//!   `cos > t` threshold means a different confidence at every `t`.
+//! - **Bitpacked** values go through HDR popcount stacking on the Hamming
+//!   distance itself, never through a cosine reinterpretation such as
+//!   `1 - 2d/N`. Their null distribution is known exactly — `Binomial(N, ½)`
+//!   for two random N-bit vectors — so [`hamming_null_z`] needs no samples.
+//!
+//! The Fisher transform here is the helix convention
+//! (`helix::fisher_z::Similarity::fisher_z`, mirrored by `jc::stats::fisher_2z`):
+//! `z = ½·(ln(1+r) − ln(1−r))` in `ln` form, with `r` clamped to
+//! `[−1+ε, 1−ε]`, `ε = `[`FISHER_CLAMP_EPS`] `= 1e-9`, so every finite input
+//! gives a finite z and NaN propagates. Hyperbolic depth `2z` is
+//! [`hyperbolic_depth`].
+
+use crate::simd::{simd_ln_f32, F32x16};
+
+/// Rim clamp for the `f64` Fisher transform — identical to
+/// `helix::fisher_z::Similarity::CLAMP_EPS`, so this and helix agree bit for
+/// bit. At `r = 1 − ε` the transform is `≈ 10.708`.
+pub const FISHER_CLAMP_EPS: f64 = 1e-9;
+
+/// Largest `f32` below `1.0` (`1 − 2⁻²⁴`): the `f32` rim clamp.
+///
+/// `1 − 1e-9` is not representable in `f32` — it rounds to exactly `1.0`,
+/// where `ln(1 − r)` is `−∞` — so the batch path clamps to the nearest
+/// representable value instead. At that rim the transform is `≈ 8.664`.
+pub const FISHER_CLAMP_F32: f32 = 1.0 - f32::EPSILON / 2.0;
+
+/// Fisher-Z transform `z = atanh(r) = ½·(ln(1+r) − ln(1−r))`.
+///
+/// `r` is clamped to `[−1+ε, 1−ε]` ([`FISHER_CLAMP_EPS`]), so the result is
+/// finite for every finite input, including exact `±1` and out-of-range
+/// values; NaN propagates. Under the usual bivariate-normal model `z` is
+/// approximately normal with variance `1/(n−3)`, which is what makes
+/// confidence bands on it mean the same thing everywhere on the scale.
+///
+/// # Example
+///
+/// ```
+/// use ndarray::hpc::zspace::fisher_z;
+///
+/// assert_eq!(fisher_z(0.0), 0.0);
+/// assert!((fisher_z(0.5) - 0.5_f64.atanh()).abs() < 1e-15);
+/// assert!(fisher_z(1.0).is_finite());
+/// ```
+#[inline]
+pub fn fisher_z(r: f64) -> f64 {
+    if r.is_nan() {
+        return f64::NAN;
+    }
+    let s = r.clamp(-1.0 + FISHER_CLAMP_EPS, 1.0 - FISHER_CLAMP_EPS);
+    0.5 * ((1.0 + s).ln() - (1.0 - s).ln())
+}
+
+/// Inverse of [`fisher_z`]: `r = tanh(z)`.
+#[inline]
+pub fn fisher_z_inv(z: f64) -> f64 {
+    z.tanh()
+}
+
+/// Hyperbolic (Poincaré-disk) depth `2·atanh(r)` — exactly `2 ×` [`fisher_z`],
+/// the "Fisher 2z" of helix and `jc`.
+#[inline]
+pub fn hyperbolic_depth(r: f64) -> f64 {
+    2.0 * fisher_z(r)
+}
+
+/// Fisher-Z over a slice of `f32` similarities, sixteen lanes at a time.
+///
+/// Each element is clamped to `[−`[`FISHER_CLAMP_F32`]`, `[`FISHER_CLAMP_F32`]`]`
+/// and transformed with the same `ln` form as [`fisher_z`], through `F32x16`
+/// and [`simd_ln_f32`]. The result is bit-identical to applying that formula
+/// to each element in scalar `f32`; NaN propagates.
+///
+/// # Panics
+///
+/// Panics if `src.len() != dst.len()`.
+///
+/// # Example
+///
+/// ```
+/// use ndarray::hpc::zspace::fisher_z_f32_batch;
+///
+/// let r = [0.0f32, 0.5, -0.5, 1.0];
+/// let mut z = [0.0f32; 4];
+/// fisher_z_f32_batch(&r, &mut z);
+/// assert_eq!(z[0], 0.0);
+/// assert_eq!(z[1], -z[2]);
+/// assert!(z[3].is_finite());
+/// ```
+pub fn fisher_z_f32_batch(src: &[f32], dst: &mut [f32]) {
+    assert_eq!(src.len(), dst.len(), "fisher_z_f32_batch: length mismatch");
+    let (cs, ts) = src.as_chunks::<16>();
+    let (cd, td) = dst.as_chunks_mut::<16>();
+    let one = F32x16::splat(1.0);
+    let half = F32x16::splat(0.5);
+    for (s, d) in cs.iter().zip(cd) {
+        let x = F32x16::from_array(core::array::from_fn(|i| clamp_f32(s[i])));
+        *d = ((simd_ln_f32(one + x) - simd_ln_f32(one - x)) * half).to_array();
+    }
+    for (s, d) in ts.iter().zip(td) {
+        *d = fisher_z_f32(*s);
+    }
+}
+
+/// Scalar `f32` Fisher-Z with the batch path's clamp — the per-element
+/// definition [`fisher_z_f32_batch`] must reproduce.
+#[inline]
+fn fisher_z_f32(r: f32) -> f32 {
+    let s = clamp_f32(r);
+    ((1.0 + s).ln() - (1.0 - s).ln()) * 0.5
+}
+
+#[inline]
+fn clamp_f32(r: f32) -> f32 {
+    // `f32::clamp` returns NaN for a NaN input, so NaN propagates.
+    r.clamp(-FISHER_CLAMP_F32, FISHER_CLAMP_F32)
+}
+
+/// z-score of a Hamming distance against the random-vector null.
+///
+/// Two independent uniformly random `n_bits`-bit vectors differ in
+/// `Binomial(n_bits, ½)` positions: `μ = n/2`, `σ = √n/2`. The result is
+/// `(d − μ)/σ`, so similar vectors score **negative** (`−3` is three sigma
+/// closer than chance) with no calibration samples at all. Returns `0.0` for
+/// `n_bits == 0`.
+///
+/// # Example
+///
+/// ```
+/// use ndarray::hpc::zspace::hamming_null_z;
+///
+/// // 16 384 bits: μ = 8192, σ = 64.
+/// assert_eq!(hamming_null_z(8192, 16_384), 0.0);
+/// assert_eq!(hamming_null_z(8192 - 3 * 64, 16_384), -3.0);
+/// ```
+#[inline]
+pub fn hamming_null_z(d: u64, n_bits: u64) -> f64 {
+    if n_bits == 0 {
+        return 0.0;
+    }
+    let n = n_bits as f64;
+    (d as f64 - n / 2.0) / (n.sqrt() / 2.0)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn fisher_z_pins_known_values() {
+        assert_eq!(fisher_z(0.0), 0.0);
+        assert!((fisher_z(0.5) - 0.549_306_144_334_054_8).abs() < 1e-15);
+        // Rim value with ε = 1e-9, computed independently in f64.
+        assert!((fisher_z(1.0) - 10.708_206_522_644_144).abs() < 1e-9);
+        assert_eq!(hyperbolic_depth(0.5), 2.0 * fisher_z(0.5));
+    }
+
+    #[test]
+    fn fisher_z_is_odd_and_monotone() {
+        let grid: Vec<f64> = (-99..=99).map(|i| f64::from(i) / 100.0).collect();
+        for w in grid.windows(2) {
+            assert!(fisher_z(w[0]) < fisher_z(w[1]), "not increasing at {}", w[0]);
+        }
+        for &r in &grid {
+            assert_eq!(fisher_z(-r), -fisher_z(r), "not odd at {r}");
+        }
+    }
+
+    #[test]
+    fn fisher_z_is_finite_everywhere_and_propagates_nan() {
+        for r in [1.0, -1.0, 2.0, -2.0, f64::INFINITY, f64::NEG_INFINITY] {
+            assert!(fisher_z(r).is_finite(), "r = {r}");
+        }
+        assert_eq!(fisher_z(1.0), fisher_z(5.0), "everything past the rim clamps to it");
+        assert!(fisher_z(f64::NAN).is_nan());
+    }
+
+    #[test]
+    fn fisher_z_inverse_round_trips() {
+        for i in -99..=99 {
+            let r = f64::from(i) / 100.0;
+            assert!((fisher_z_inv(fisher_z(r)) - r).abs() < 1e-12, "r = {r}");
+        }
+    }
+
+    /// The batch path must equal the scalar f32 definition bit for bit, at
+    /// every length across the 16-lane boundary and at the rim.
+    #[test]
+    fn batch_is_bit_identical_to_scalar_f32() {
+        let mut src: Vec<f32> = (0..200)
+            .map(|i| ((i * 37) % 199) as f32 / 99.0 - 1.0)
+            .collect();
+        src[3] = 1.0;
+        src[20] = -1.0;
+        src[33] = 7.5;
+        for n in (0..=40).chain([199, 200]) {
+            let mut dst = vec![0.0f32; n];
+            fisher_z_f32_batch(&src[..n], &mut dst);
+            for (i, (&r, &z)) in src[..n].iter().zip(&dst).enumerate() {
+                assert_eq!(z.to_bits(), fisher_z_f32(r).to_bits(), "n={n} i={i} r={r}");
+                assert!(z.is_finite(), "n={n} i={i} r={r}");
+            }
+        }
+    }
+
+    /// The f32 rim clamp is load-bearing: without it, `r = 1.0` in f32 is
+    /// `ln(0) = −∞` and the result is infinite.
+    #[test]
+    fn f32_rim_is_finite_and_close_to_f64() {
+        let mut z = [0.0f32; 2];
+        fisher_z_f32_batch(&[1.0, -1.0], &mut z);
+        assert!(z[0].is_finite() && z[1].is_finite());
+        assert!((f64::from(z[0]) - fisher_z(f64::from(FISHER_CLAMP_F32))).abs() < 1e-3);
+        let mut nan = [0.0f32; 1];
+        fisher_z_f32_batch(&[f32::NAN], &mut nan);
+        assert!(nan[0].is_nan());
+    }
+
+    #[test]
+    fn hamming_null_z_matches_the_binomial_null() {
+        assert_eq!(hamming_null_z(8192, 16_384), 0.0);
+        assert_eq!(hamming_null_z(8192 - 64, 16_384), -1.0);
+        assert_eq!(hamming_null_z(8192 + 192, 16_384), 3.0);
+        assert_eq!(hamming_null_z(5, 0), 0.0);
+    }
+}
diff --git a/src/simd.rs b/src/simd.rs
index f7354a3f..40d51667 100644
--- a/src/simd.rs
+++ b/src/simd.rs
@@ -618,6 +618,7 @@ pub use crate::hpc::fingerprint::{
     vector_config, Fingerprint, Fingerprint1K, Fingerprint2K, Fingerprint64K, VectorConfig, VectorWidth,
 };
 pub use crate::hpc::statistics::{moments_u32, MomentsU32};
+pub use crate::hpc::zspace::{fisher_z, fisher_z_f32_batch, fisher_z_inv, hamming_null_z, hyperbolic_depth};
 
 // PR-X1 — SoA carrier + const-size slice helpers, dispatched from their
 // respective `simd_{type}.rs` modules. The W1a consumer contract forbids

From b8babbd47d4728f114b7950999562cd9691e37b1 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Wed, 23 Sep 2026 20:43:30 +0000
Subject: [PATCH 4/8] reliability: FidelityReport pays the Fisher-Z entry tax

pearson / spearman / cosine stay raw for display; pearson_z / spearman_z /
cosine_z are the forms any threshold, average or comparison must use, and
the doc example now thresholds in z-space. Found by the cosine-tax census:
this was the one ndarray site consuming a raw cosine-shaped coefficient.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm
---
 src/hpc/reliability.rs | 40 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 39 insertions(+), 1 deletion(-)

diff --git a/src/hpc/reliability.rs b/src/hpc/reliability.rs
index 56adf785..a529ea32 100644
--- a/src/hpc/reliability.rs
+++ b/src/hpc/reliability.rs
@@ -215,6 +215,12 @@ pub fn icc_a1(ratings: &[&[f64]]) -> f64 {
 /// One call computes all four coefficients plus the relative-L2 error and
 /// cosine similarity, so a harness can print a row per codec/flavor without
 /// recomputing means four times.
+///
+/// `pearson`, `spearman` and `cosine` are cosine-shaped and stay raw for
+/// display; anything that thresholds, averages or compares them goes through
+/// [`pearson_z`](Self::pearson_z) / [`spearman_z`](Self::spearman_z) /
+/// [`cosine_z`](Self::cosine_z) (the Fisher-Z entry point in
+/// [`zspace`](crate::hpc::zspace)).
 #[derive(Debug, Clone, Copy, PartialEq)]
 pub struct FidelityReport {
     /// Pearson product-moment correlation (linear association).
@@ -246,7 +252,9 @@ impl FidelityReport {
     /// let truth = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
     /// let est = [0.1, 0.9, 2.1, 2.9, 4.2, 4.8];
     /// let r = FidelityReport::compute(&truth, &est);
-    /// assert!(r.pearson > 0.99 && r.spearman > 0.99);
+    /// // Cosine-shaped coefficients are consumed in Fisher-Z space.
+    /// use ndarray::hpc::zspace::fisher_z;
+    /// assert!(r.pearson_z() > fisher_z(0.99) && r.spearman_z() > fisher_z(0.99));
     /// assert!(r.rel_l2 < 0.1);
     /// // Mismatched lengths → degenerate (does NOT truncate-then-score):
     /// let bad = FidelityReport::compute(&[1.0, 2.0, 100.0], &[1.0, 2.0]);
@@ -295,12 +303,42 @@ impl FidelityReport {
             cosine,
         }
     }
+
+    /// Fisher-Z of [`pearson`](Self::pearson) — the form every threshold,
+    /// average or confidence interval over it must use.
+    pub fn pearson_z(&self) -> f64 {
+        crate::hpc::zspace::fisher_z(self.pearson)
+    }
+
+    /// Fisher-Z of [`spearman`](Self::spearman).
+    pub fn spearman_z(&self) -> f64 {
+        crate::hpc::zspace::fisher_z(self.spearman)
+    }
+
+    /// Fisher-Z of [`cosine`](Self::cosine).
+    pub fn cosine_z(&self) -> f64 {
+        crate::hpc::zspace::fisher_z(self.cosine)
+    }
 }
 
 #[cfg(test)]
 mod tests {
     use super::*;
 
+    /// The z accessors are exactly Fisher-Z of the raw coefficients, and a
+    /// perfect score stays finite (the rim clamp, not `atanh(1) = inf`).
+    #[test]
+    fn fidelity_z_accessors_are_fisher_z() {
+        use crate::hpc::zspace::fisher_z;
+        let truth = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
+        let r = FidelityReport::compute(&truth, &[0.1, 0.9, 2.1, 2.9, 4.2, 4.8]);
+        assert_eq!(r.pearson_z(), fisher_z(r.pearson));
+        assert_eq!(r.spearman_z(), fisher_z(r.spearman));
+        assert_eq!(r.cosine_z(), fisher_z(r.cosine));
+        let perfect = FidelityReport::compute(&truth, &truth);
+        assert!(perfect.pearson_z().is_finite() && perfect.pearson_z() > 10.0);
+    }
+
     #[test]
     fn pearson_perfect_and_anti() {
         let x = [1.0, 2.0, 3.0, 4.0, 5.0];

From bb14de65e1521b24fb2db0ee273989c400232a31 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Wed, 23 Sep 2026 20:51:06 +0000
Subject: [PATCH 5/8] =?UTF-8?q?zspace:=20z=20is=20never=20materialized=20?=
 =?UTF-8?q?=E2=80=94=20ZGamma=20envelope=20encodes=20cosine=20straight=20t?=
 =?UTF-8?q?o=20i8=20codes?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Remove fisher_z_f32_batch (wrote a z buffer) and fisher_z_inv (the way back
to cosine). Add ZGamma { z_min, z_range }: 8 LE bytes, the FamilyGamma
layout and code formula. fit() scans only the cosine min/max (Fisher-Z is
monotone), encode_batch() fuses clamp + transform + normalize in F32x16 and
writes only the i8 code, and threshold_code() moves a threshold into code
space once. fisher_z stays for single scalars.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm
---
 src/hpc/zspace.rs | 288 ++++++++++++++++++++++++++++++++++------------
 src/simd.rs       |   2 +-
 2 files changed, 218 insertions(+), 72 deletions(-)

diff --git a/src/hpc/zspace.rs b/src/hpc/zspace.rs
index 3171dae1..800b7a03 100644
--- a/src/hpc/zspace.rs
+++ b/src/hpc/zspace.rs
@@ -6,7 +6,7 @@
 //!
 //! - **Cosine-shaped** values (cosine, Pearson r, normalized dot products —
 //!   anything bounded in `[-1, 1]`) pay the Fisher-Z entry tax:
-//!   [`fisher_z`] / [`fisher_z_f32_batch`]. No exceptions, including lab and
+//!   [`ZGamma`] (codes) / [`fisher_z`] (one scalar). No exceptions, including lab and
 //!   calibration code. A raw cosine is not a z-score: its variance shrinks
 //!   towards the rim, so equal steps do not mean equal evidence and a fixed
 //!   `cos > t` threshold means a different confidence at every `t`.
@@ -21,6 +21,11 @@
 //! `[−1+ε, 1−ε]`, `ε = `[`FISHER_CLAMP_EPS`] `= 1e-9`, so every finite input
 //! gives a finite z and NaN propagates. Hyperbolic depth `2z` is
 //! [`hyperbolic_depth`].
+//!
+//! **z is never materialized.** Populations are encoded through a
+//! [`ZGamma`] envelope straight to `i8` codes and stay there; there is no
+//! z buffer and no inverse back to cosine. [`fisher_z`] exists for single
+//! scalars — a report statistic, a threshold — not for arrays.
 
 use crate::simd::{simd_ln_f32, F32x16};
 
@@ -62,12 +67,6 @@ pub fn fisher_z(r: f64) -> f64 {
     0.5 * ((1.0 + s).ln() - (1.0 - s).ln())
 }
 
-/// Inverse of [`fisher_z`]: `r = tanh(z)`.
-#[inline]
-pub fn fisher_z_inv(z: f64) -> f64 {
-    z.tanh()
-}
-
 /// Hyperbolic (Poincaré-disk) depth `2·atanh(r)` — exactly `2 ×` [`fisher_z`],
 /// the "Fisher 2z" of helix and `jc`.
 #[inline]
@@ -75,46 +74,133 @@ pub fn hyperbolic_depth(r: f64) -> f64 {
     2.0 * fisher_z(r)
 }
 
-/// Fisher-Z over a slice of `f32` similarities, sixteen lanes at a time.
-///
-/// Each element is clamped to `[−`[`FISHER_CLAMP_F32`]`, `[`FISHER_CLAMP_F32`]`]`
-/// and transformed with the same `ln` form as [`fisher_z`], through `F32x16`
-/// and [`simd_ln_f32`]. The result is bit-identical to applying that formula
-/// to each element in scalar `f32`; NaN propagates.
-///
-/// # Panics
-///
-/// Panics if `src.len() != dst.len()`.
+/// The z-space envelope: the metadata that makes an `i8` code mean a z-score.
 ///
-/// # Example
-///
-/// ```
-/// use ndarray::hpc::zspace::fisher_z_f32_batch;
+/// **Fisher-Z is never materialized.** There is no buffer of z values and no
+/// way back to cosine. A cosine enters, is transformed in registers, is
+/// normalized against this envelope, and leaves as an `i8` code. Everything
+/// after that — thresholds, bands, ranking — works on codes. A threshold is
+/// converted into code space once with [`ZGamma::threshold_code`]; it is never
+/// converted back.
 ///
-/// let r = [0.0f32, 0.5, -0.5, 1.0];
-/// let mut z = [0.0f32; 4];
-/// fisher_z_f32_batch(&r, &mut z);
-/// assert_eq!(z[0], 0.0);
-/// assert_eq!(z[1], -z[2]);
-/// assert!(z[3].is_finite());
-/// ```
-pub fn fisher_z_f32_batch(src: &[f32], dst: &mut [f32]) {
-    assert_eq!(src.len(), dst.len(), "fisher_z_f32_batch: length mismatch");
-    let (cs, ts) = src.as_chunks::<16>();
-    let (cd, td) = dst.as_chunks_mut::<16>();
-    let one = F32x16::splat(1.0);
-    let half = F32x16::splat(0.5);
-    for (s, d) in cs.iter().zip(cd) {
-        let x = F32x16::from_array(core::array::from_fn(|i| clamp_f32(s[i])));
-        *d = ((simd_ln_f32(one + x) - simd_ln_f32(one - x)) * half).to_array();
-    }
-    for (s, d) in ts.iter().zip(td) {
-        *d = fisher_z_f32(*s);
+/// The code is `((z − z_min) / z_range)·254 − 127`, clamped to
+/// `[−128, 127]` and truncated toward zero. `z_min`/`z_range` travel with the
+/// codes as 8 bytes of little-endian `f32` ([`ZGamma::to_le_bytes`]). The byte
+/// layout and the code formula match
+/// `bgz_tensor::fisher_z::FamilyGamma`, so an envelope written by one can be
+/// read by the other. The two clamp the cosine differently at the rim
+/// (bgz-tensor uses `0.9999`, this uses [`FISHER_CLAMP_F32`]), so codes agree
+/// bit for bit only for `|cos| ≤ 0.9999`.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct ZGamma {
+    /// z of the smallest cosine the envelope covers.
+    pub z_min: f32,
+    /// `z_max − z_min`; always `> 0`.
+    pub z_range: f32,
+}
+
+impl ZGamma {
+    /// Serialized size: two little-endian `f32`s.
+    pub const BYTES: usize = 8;
+
+    /// Fit the envelope to a population of cosines.
+    ///
+    /// Fisher-Z is monotone, so the z range is the transform of the cosine
+    /// range: only the cosine minimum and maximum are scanned, and no z value
+    /// is stored. NaN elements are ignored. An empty, all-NaN, or constant
+    /// population gets `z_range = 1.0` so that codes stay finite.
+    pub fn fit(cosines: &[f32]) -> Self {
+        let (mut lo, mut hi) = (f32::INFINITY, f32::NEG_INFINITY);
+        for &c in cosines {
+            if !c.is_nan() {
+                lo = lo.min(c);
+                hi = hi.max(c);
+            }
+        }
+        if lo > hi {
+            return Self {
+                z_min: 0.0,
+                z_range: 1.0,
+            };
+        }
+        let z_min = fisher_z_f32(lo);
+        let range = fisher_z_f32(hi) - z_min;
+        Self {
+            z_min,
+            z_range: if range > 0.0 { range } else { 1.0 },
+        }
+    }
+
+    /// Encode one cosine as a z-space code.
+    ///
+    /// NaN encodes as `0`, the code for the envelope midpoint, because Rust's
+    /// saturating `as i8` maps NaN to zero.
+    #[inline]
+    pub fn encode(&self, cosine: f32) -> i8 {
+        let n = (fisher_z_f32(cosine) - self.z_min) / self.z_range;
+        (n * 254.0 - 127.0).clamp(-128.0, 127.0) as i8
+    }
+
+    /// Encode a slice of cosines, sixteen lanes at a time.
+    ///
+    /// Clamp, Fisher transform and normalization are fused in `F32x16`
+    /// registers; the only thing written is the `i8` code. The result is
+    /// bit-identical to [`ZGamma::encode`] per element.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `cosines.len() != codes.len()`.
+    pub fn encode_batch(&self, cosines: &[f32], codes: &mut [i8]) {
+        assert_eq!(cosines.len(), codes.len(), "ZGamma::encode_batch: length mismatch");
+        let (cs, ts) = cosines.as_chunks::<16>();
+        let (cd, td) = codes.as_chunks_mut::<16>();
+        let one = F32x16::splat(1.0);
+        let half = F32x16::splat(0.5);
+        let z_min = F32x16::splat(self.z_min);
+        let z_range = F32x16::splat(self.z_range);
+        let scale = F32x16::splat(254.0);
+        let shift = F32x16::splat(127.0);
+        for (s, d) in cs.iter().zip(cd) {
+            let x = F32x16::from_array(core::array::from_fn(|i| clamp_f32(s[i])));
+            let z = (simd_ln_f32(one + x) - simd_ln_f32(one - x)) * half;
+            let v = (((z - z_min) / z_range) * scale - shift).to_array();
+            *d = core::array::from_fn(|i| v[i].clamp(-128.0, 127.0) as i8);
+        }
+        for (s, d) in ts.iter().zip(td) {
+            *d = self.encode(*s);
+        }
+    }
+
+    /// Convert a cosine threshold into code space, once.
+    ///
+    /// Because the code is monotone in the cosine, `encode(c) >= t` with
+    /// `t = threshold_code(c₀)` holds for every `c ≥ c₀`. A cosine slightly
+    /// below `c₀` can share its code, so the test is inclusive at the
+    /// resolution of one code step.
+    #[inline]
+    pub fn threshold_code(&self, cosine: f32) -> i8 {
+        self.encode(cosine)
+    }
+
+    /// The envelope as 8 little-endian bytes: `z_min` then `z_range`.
+    pub fn to_le_bytes(&self) -> [u8; Self::BYTES] {
+        let mut b = [0u8; Self::BYTES];
+        b[..4].copy_from_slice(&self.z_min.to_le_bytes());
+        b[4..].copy_from_slice(&self.z_range.to_le_bytes());
+        b
+    }
+
+    /// Read an envelope written by [`ZGamma::to_le_bytes`].
+    pub fn from_le_bytes(b: [u8; Self::BYTES]) -> Self {
+        Self {
+            z_min: f32::from_le_bytes([b[0], b[1], b[2], b[3]]),
+            z_range: f32::from_le_bytes([b[4], b[5], b[6], b[7]]),
+        }
     }
 }
 
-/// Scalar `f32` Fisher-Z with the batch path's clamp — the per-element
-/// definition [`fisher_z_f32_batch`] must reproduce.
+/// Scalar `f32` Fisher-Z with the `f32` rim clamp. Private on purpose: its
+/// result only ever lives in a register on its way to a code.
 #[inline]
 fn fisher_z_f32(r: f32) -> f32 {
     let s = clamp_f32(r);
@@ -186,45 +272,105 @@ mod tests {
         assert!(fisher_z(f64::NAN).is_nan());
     }
 
+    fn cosines() -> Vec<f32> {
+        let mut v: Vec<f32> = (0..200)
+            .map(|i| ((i * 37) % 199) as f32 / 99.0 - 1.0)
+            .collect();
+        v[3] = 1.0;
+        v[20] = -1.0;
+        v[33] = 7.5;
+        v[50] = f32::NAN;
+        v
+    }
+
+    /// The fused batch path must equal the scalar encoder bit for bit, at
+    /// every length across the 16-lane boundary, at the rim, and on NaN.
     #[test]
-    fn fisher_z_inverse_round_trips() {
-        for i in -99..=99 {
-            let r = f64::from(i) / 100.0;
-            assert!((fisher_z_inv(fisher_z(r)) - r).abs() < 1e-12, "r = {r}");
+    fn encode_batch_is_bit_identical_to_encode() {
+        let src = cosines();
+        let g = ZGamma::fit(&src[..150]);
+        for n in (0..=40).chain([199, 200]) {
+            let mut dst = vec![0i8; n];
+            g.encode_batch(&src[..n], &mut dst);
+            for (i, (&c, &k)) in src[..n].iter().zip(&dst).enumerate() {
+                assert_eq!(k, g.encode(c), "n={n} i={i} c={c}");
+            }
         }
     }
 
-    /// The batch path must equal the scalar f32 definition bit for bit, at
-    /// every length across the 16-lane boundary and at the rim.
+    /// The fitted envelope spans the full code range: the population minimum
+    /// lands at −127, the maximum at 127.
     #[test]
-    fn batch_is_bit_identical_to_scalar_f32() {
-        let mut src: Vec<f32> = (0..200)
-            .map(|i| ((i * 37) % 199) as f32 / 99.0 - 1.0)
-            .collect();
-        src[3] = 1.0;
-        src[20] = -1.0;
-        src[33] = 7.5;
-        for n in (0..=40).chain([199, 200]) {
-            let mut dst = vec![0.0f32; n];
-            fisher_z_f32_batch(&src[..n], &mut dst);
-            for (i, (&r, &z)) in src[..n].iter().zip(&dst).enumerate() {
-                assert_eq!(z.to_bits(), fisher_z_f32(r).to_bits(), "n={n} i={i} r={r}");
-                assert!(z.is_finite(), "n={n} i={i} r={r}");
+    fn fit_spans_the_code_range() {
+        let pop = [-0.2f32, 0.1, 0.4, 0.9];
+        let g = ZGamma::fit(&pop);
+        assert_eq!(g.encode(-0.2), -127);
+        assert_eq!(g.encode(0.9), 127);
+        assert!(g.encode(0.1) > -127 && g.encode(0.4) < 127);
+        // Outside the envelope saturates rather than wrapping.
+        assert_eq!(g.encode(-0.9), -128);
+        assert_eq!(g.encode(0.99), 127);
+    }
+
+    /// Codes are spaced in z, not in cosine: two cosine steps of the same
+    /// size get more codes near the rim, where the evidence is stronger.
+    #[test]
+    fn codes_are_spaced_in_z_not_in_cosine() {
+        let g = ZGamma::fit(&[0.0, 0.99]);
+        let mid = i32::from(g.encode(0.10)) - i32::from(g.encode(0.00));
+        let rim = i32::from(g.encode(0.99)) - i32::from(g.encode(0.89));
+        assert!(rim > 2 * mid, "rim step {rim} vs centre step {mid}");
+    }
+
+    #[test]
+    fn fit_is_finite_on_degenerate_populations() {
+        for pop in [&[][..], &[f32::NAN][..], &[0.5, 0.5][..], &[1.0, 1.0][..]] {
+            let g = ZGamma::fit(pop);
+            assert!(g.z_min.is_finite() && g.z_range > 0.0, "{pop:?}");
+            assert!(g.encode(0.3) >= -128);
+        }
+        assert_eq!(ZGamma::fit(&[f32::NAN]).encode(f32::NAN), 0);
+    }
+
+    #[test]
+    fn threshold_code_is_inclusive_and_monotone() {
+        let src = cosines();
+        let g = ZGamma::fit(&src);
+        let t = g.threshold_code(0.3);
+        for &c in src.iter().filter(|c| !c.is_nan()) {
+            if c >= 0.3 {
+                assert!(g.encode(c) >= t, "c={c}");
             }
         }
     }
 
-    /// The f32 rim clamp is load-bearing: without it, `r = 1.0` in f32 is
-    /// `ln(0) = −∞` and the result is infinite.
     #[test]
-    fn f32_rim_is_finite_and_close_to_f64() {
-        let mut z = [0.0f32; 2];
-        fisher_z_f32_batch(&[1.0, -1.0], &mut z);
-        assert!(z[0].is_finite() && z[1].is_finite());
-        assert!((f64::from(z[0]) - fisher_z(f64::from(FISHER_CLAMP_F32))).abs() < 1e-3);
-        let mut nan = [0.0f32; 1];
-        fisher_z_f32_batch(&[f32::NAN], &mut nan);
-        assert!(nan[0].is_nan());
+    fn envelope_round_trips_through_8_le_bytes() {
+        let g = ZGamma {
+            z_min: -0.75,
+            z_range: 3.25,
+        };
+        let b = g.to_le_bytes();
+        assert_eq!(b.len(), 8);
+        assert_eq!(&b[..4], &(-0.75f32).to_le_bytes());
+        assert_eq!(ZGamma::from_le_bytes(b), g);
+    }
+
+    /// Same formula as `bgz_tensor::fisher_z::FamilyGamma::encode`, written
+    /// out independently: codes must agree away from the rim.
+    #[test]
+    fn matches_the_family_gamma_formula() {
+        let g = ZGamma {
+            z_min: -0.5,
+            z_range: 2.0,
+        };
+        for i in -99..=99 {
+            let c = i as f32 / 100.0;
+            let z = c.clamp(-0.9999, 0.9999).atanh();
+            let want = (((z - g.z_min) / g.z_range) * 254.0 - 127.0).clamp(-128.0, 127.0) as i8;
+            let got = g.encode(c);
+            assert!((i32::from(got) - i32::from(want)).abs() <= 1, "c={c} got={got} want={want}");
+        }
     }
 
     #[test]
diff --git a/src/simd.rs b/src/simd.rs
index 40d51667..58e25f7b 100644
--- a/src/simd.rs
+++ b/src/simd.rs
@@ -618,7 +618,7 @@ pub use crate::hpc::fingerprint::{
     vector_config, Fingerprint, Fingerprint1K, Fingerprint2K, Fingerprint64K, VectorConfig, VectorWidth,
 };
 pub use crate::hpc::statistics::{moments_u32, MomentsU32};
-pub use crate::hpc::zspace::{fisher_z, fisher_z_f32_batch, fisher_z_inv, hamming_null_z, hyperbolic_depth};
+pub use crate::hpc::zspace::{fisher_z, hamming_null_z, hyperbolic_depth, ZGamma};
 
 // PR-X1 — SoA carrier + const-size slice helpers, dispatched from their
 // respective `simd_{type}.rs` modules. The W1a consumer contract forbids

From fd45c890c381822d57c60c1ab05fe4ef7536365d Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Wed, 23 Sep 2026 21:21:16 +0000
Subject: [PATCH 6/8] =?UTF-8?q?zspace:=20round,=20symmetric=20=C2=B1127,?=
 =?UTF-8?q?=20reserve=20-128=20as=20the=20NaN=20code?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The FamilyGamma code inherited from bgz-tensor truncated toward zero, which
doubles the code-0 bucket and pulls every code half a step inward (measured
+0.49 / -0.49 mean error either side of centre), used -128 as a data level
reachable only by saturating below the envelope, and mapped NaN to 0, a real
code. Now: round-ties-even, saturate to -127..=127, NaN -> ZGamma::NAN_CODE
(-128). One shared quantize() keeps scalar and batch identical. The
FamilyGamma formula-parity test is replaced by a median-bias falsifier and a
NaN-sentinel test.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm
---
 src/hpc/zspace.rs | 104 ++++++++++++++++++++++++++++++++++------------
 1 file changed, 78 insertions(+), 26 deletions(-)

diff --git a/src/hpc/zspace.rs b/src/hpc/zspace.rs
index 800b7a03..ab2a1023 100644
--- a/src/hpc/zspace.rs
+++ b/src/hpc/zspace.rs
@@ -83,14 +83,22 @@ pub fn hyperbolic_depth(r: f64) -> f64 {
 /// converted into code space once with [`ZGamma::threshold_code`]; it is never
 /// converted back.
 ///
-/// The code is `((z − z_min) / z_range)·254 − 127`, clamped to
-/// `[−128, 127]` and truncated toward zero. `z_min`/`z_range` travel with the
-/// codes as 8 bytes of little-endian `f32` ([`ZGamma::to_le_bytes`]). The byte
-/// layout and the code formula match
-/// `bgz_tensor::fisher_z::FamilyGamma`, so an envelope written by one can be
-/// read by the other. The two clamp the cosine differently at the rim
-/// (bgz-tensor uses `0.9999`, this uses [`FISHER_CLAMP_F32`]), so codes agree
-/// bit for bit only for `|cos| ≤ 0.9999`.
+/// The code is `((z − z_min) / z_range)·254 − 127`, **rounded** to the
+/// nearest integer (ties to even) and saturated to the **symmetric** range
+/// `[−127, 127]`. `z_min`/`z_range` travel with the codes as 8 bytes of
+/// little-endian `f32` ([`ZGamma::to_le_bytes`]).
+///
+/// Two deliberate departures from `bgz_tensor::fisher_z::FamilyGamma`, whose
+/// byte layout this shares but whose code does not:
+///
+/// - **Rounding, not truncation.** `as i8` truncates toward zero, which makes
+///   the code-0 bucket twice as wide as every other and pulls every code half
+///   a step toward the centre (measured: mean error `+0.49` below the centre,
+///   `−0.49` above it). Rounding removes that median bias.
+/// - **Symmetric range with a sentinel.** Two's-complement `i8` is asymmetric
+///   (`−128..=127`); using `−128` for data would give one side an extra level.
+///   Data codes stay in `−127..=127` and `−128` is reserved as
+///   [`ZGamma::NAN_CODE`], so NaN can never be mistaken for a real value.
 #[derive(Clone, Copy, Debug, PartialEq)]
 pub struct ZGamma {
     /// z of the smallest cosine the envelope covers.
@@ -103,6 +111,10 @@ impl ZGamma {
     /// Serialized size: two little-endian `f32`s.
     pub const BYTES: usize = 8;
 
+    /// The code for a NaN input. It is outside the data range `−127..=127`,
+    /// so no finite cosine ever encodes to it.
+    pub const NAN_CODE: i8 = i8::MIN;
+
     /// Fit the envelope to a population of cosines.
     ///
     /// Fisher-Z is monotone, so the z range is the transform of the cosine
@@ -133,12 +145,12 @@ impl ZGamma {
 
     /// Encode one cosine as a z-space code.
     ///
-    /// NaN encodes as `0`, the code for the envelope midpoint, because Rust's
-    /// saturating `as i8` maps NaN to zero.
+    /// NaN encodes as [`ZGamma::NAN_CODE`]; values outside the envelope
+    /// saturate to `±127`.
     #[inline]
     pub fn encode(&self, cosine: f32) -> i8 {
         let n = (fisher_z_f32(cosine) - self.z_min) / self.z_range;
-        (n * 254.0 - 127.0).clamp(-128.0, 127.0) as i8
+        quantize(n * 254.0 - 127.0)
     }
 
     /// Encode a slice of cosines, sixteen lanes at a time.
@@ -164,7 +176,7 @@ impl ZGamma {
             let x = F32x16::from_array(core::array::from_fn(|i| clamp_f32(s[i])));
             let z = (simd_ln_f32(one + x) - simd_ln_f32(one - x)) * half;
             let v = (((z - z_min) / z_range) * scale - shift).to_array();
-            *d = core::array::from_fn(|i| v[i].clamp(-128.0, 127.0) as i8);
+            *d = core::array::from_fn(|i| quantize(v[i]));
         }
         for (s, d) in ts.iter().zip(td) {
             *d = self.encode(*s);
@@ -199,6 +211,17 @@ impl ZGamma {
     }
 }
 
+/// Round to the nearest code (ties to even), saturate to the symmetric data
+/// range, and map NaN to the sentinel. Shared by the scalar and batch paths so
+/// they cannot drift.
+#[inline]
+fn quantize(v: f32) -> i8 {
+    if v.is_nan() {
+        return ZGamma::NAN_CODE;
+    }
+    v.round_ties_even().clamp(-127.0, 127.0) as i8
+}
+
 /// Scalar `f32` Fisher-Z with the `f32` rim clamp. Private on purpose: its
 /// result only ever lives in a register on its way to a code.
 #[inline]
@@ -307,8 +330,8 @@ mod tests {
         assert_eq!(g.encode(-0.2), -127);
         assert_eq!(g.encode(0.9), 127);
         assert!(g.encode(0.1) > -127 && g.encode(0.4) < 127);
-        // Outside the envelope saturates rather than wrapping.
-        assert_eq!(g.encode(-0.9), -128);
+        // Outside the envelope saturates, symmetrically, rather than wrapping.
+        assert_eq!(g.encode(-0.9), -127);
         assert_eq!(g.encode(0.99), 127);
     }
 
@@ -327,9 +350,8 @@ mod tests {
         for pop in [&[][..], &[f32::NAN][..], &[0.5, 0.5][..], &[1.0, 1.0][..]] {
             let g = ZGamma::fit(pop);
             assert!(g.z_min.is_finite() && g.z_range > 0.0, "{pop:?}");
-            assert!(g.encode(0.3) >= -128);
+            assert!(g.encode(0.3) >= -127);
         }
-        assert_eq!(ZGamma::fit(&[f32::NAN]).encode(f32::NAN), 0);
     }
 
     #[test]
@@ -356,20 +378,50 @@ mod tests {
         assert_eq!(ZGamma::from_le_bytes(b), g);
     }
 
-    /// Same formula as `bgz_tensor::fisher_z::FamilyGamma::encode`, written
-    /// out independently: codes must agree away from the rim.
+    /// No median bias: over z uniform on the envelope, code 0 holds the same
+    /// share as its neighbours and the signed error is centred on both sides.
+    /// Truncation toward zero fails this: code 0 gets twice its share and
+    /// every code is pulled half a step inward.
     #[test]
-    fn matches_the_family_gamma_formula() {
+    fn rounding_has_no_median_bias() {
         let g = ZGamma {
-            z_min: -0.5,
+            z_min: -1.0,
             z_range: 2.0,
         };
-        for i in -99..=99 {
-            let c = i as f32 / 100.0;
-            let z = c.clamp(-0.9999, 0.9999).atanh();
-            let want = (((z - g.z_min) / g.z_range) * 254.0 - 127.0).clamp(-128.0, 127.0) as i8;
-            let got = g.encode(c);
-            assert!((i32::from(got) - i32::from(want)).abs() <= 1, "c={c} got={got} want={want}");
+        let mut count = [0u32; 256];
+        let (mut err_neg, mut n_neg, mut err_pos, mut n_pos) = (0.0f64, 0u32, 0.0f64, 0u32);
+        for i in 0..=20_000 {
+            let z = -1.0 + 2.0 * i as f32 / 20_000.0;
+            let k = g.encode(z.tanh());
+            count[(i32::from(k) + 128) as usize] += 1;
+            let v = f64::from(((z - g.z_min) / g.z_range) * 254.0 - 127.0);
+            let e = f64::from(k) - v;
+            if v < -0.5 {
+                err_neg += e;
+                n_neg += 1;
+            } else if v > 0.5 {
+                err_pos += e;
+                n_pos += 1;
+            }
+        }
+        let (c0, c1) = (count[128], count[129]);
+        assert!(c0 * 10 <= c1 * 12, "code 0 over-full: {c0} vs neighbour {c1}");
+        assert!((err_neg / f64::from(n_neg)).abs() < 0.05, "bias below centre {}", err_neg / f64::from(n_neg));
+        assert!((err_pos / f64::from(n_pos)).abs() < 0.05, "bias above centre {}", err_pos / f64::from(n_pos));
+    }
+
+    /// NaN has its own code, and no finite input — however far outside the
+    /// envelope — can produce it.
+    #[test]
+    fn nan_sentinel_is_reserved() {
+        let g = ZGamma::fit(&[-0.2, 0.3]);
+        assert_eq!(g.encode(f32::NAN), ZGamma::NAN_CODE);
+        let mut out = [0i8; 1];
+        g.encode_batch(&[f32::NAN], &mut out);
+        assert_eq!(out[0], ZGamma::NAN_CODE);
+        for c in [-1.0f32, -0.99, -0.5, 0.0, 0.5, 1.0, 5.0, -5.0, f32::INFINITY, f32::NEG_INFINITY] {
+            let k = g.encode(c);
+            assert!((-127..=127).contains(&k), "c={c} -> {k}");
         }
     }
 

From f2a461ab70695e6d958461b40168e94862806f54 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Wed, 23 Sep 2026 21:34:00 +0000
Subject: [PATCH 7/8] =?UTF-8?q?zspace:=20ZGamma=20codes=20are=20bit-exact?=
 =?UTF-8?q?=20on=20every=20target=20=E2=80=94=20deterministic=20ln?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

simd_ln_f32 is per-lane libm .ln(), so the Fisher path inherited the
platform libm: measured 624 of 8539 grid values (f32) and 566 (f64) differ by
1-2 ulp between x86-64 glibc and wasm32. Codes happened to agree on the grid,
but a one-ulp difference at a rounding boundary flips a code.

ln_det: musl/FreeBSD logf from IEEE add/sub/mul/div + bit ops in fixed order,
< 1 ulp over [2^-24, 2), used by both the scalar and batch encoders. Its x86
z digest now equals the wasm (musl-derived) libm digest. Codes and z digests
are pinned as constants. fisher_z (f64 scalar) keeps libm and is documented
as not bit-exact.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm
---
 src/hpc/zspace.rs | 129 ++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 126 insertions(+), 3 deletions(-)

diff --git a/src/hpc/zspace.rs b/src/hpc/zspace.rs
index ab2a1023..df9d1d2e 100644
--- a/src/hpc/zspace.rs
+++ b/src/hpc/zspace.rs
@@ -27,7 +27,7 @@
 //! z buffer and no inverse back to cosine. [`fisher_z`] exists for single
 //! scalars — a report statistic, a threshold — not for arrays.
 
-use crate::simd::{simd_ln_f32, F32x16};
+use crate::simd::F32x16;
 
 /// Rim clamp for the `f64` Fisher transform — identical to
 /// `helix::fisher_z::Similarity::CLAMP_EPS`, so this and helix agree bit for
@@ -49,6 +49,11 @@ pub const FISHER_CLAMP_F32: f32 = 1.0 - f32::EPSILON / 2.0;
 /// approximately normal with variance `1/(n−3)`, which is what makes
 /// confidence bands on it mean the same thing everywhere on the scale.
 ///
+/// **Not bit-exact across targets.** This is the `f64` scalar for single
+/// values (a report statistic, a threshold) and uses the platform `ln`; its
+/// last ulp can differ between targets (measured: x86-64 glibc vs wasm32).
+/// The substrate path — [`ZGamma`] codes — does not use it and is bit-exact.
+///
 /// # Example
 ///
 /// ```
@@ -174,7 +179,10 @@ impl ZGamma {
         let shift = F32x16::splat(127.0);
         for (s, d) in cs.iter().zip(cd) {
             let x = F32x16::from_array(core::array::from_fn(|i| clamp_f32(s[i])));
-            let z = (simd_ln_f32(one + x) - simd_ln_f32(one - x)) * half;
+            let (p, m) = ((one + x).to_array(), (one - x).to_array());
+            let lp = F32x16::from_array(core::array::from_fn(|i| ln_det(p[i])));
+            let lm = F32x16::from_array(core::array::from_fn(|i| ln_det(m[i])));
+            let z = (lp - lm) * half;
             let v = (((z - z_min) / z_range) * scale - shift).to_array();
             *d = core::array::from_fn(|i| quantize(v[i]));
         }
@@ -224,10 +232,54 @@ fn quantize(v: f32) -> i8 {
 
 /// Scalar `f32` Fisher-Z with the `f32` rim clamp. Private on purpose: its
 /// result only ever lives in a register on its way to a code.
+///
+/// Uses [`ln_det`], not `f32::ln`: the platform libm differs in the last one
+/// or two ulp between targets (measured: 624 of 8 539 grid values differ
+/// between x86-64 glibc and wasm32), and a one-ulp difference at a rounding
+/// boundary flips a code.
 #[inline]
 fn fisher_z_f32(r: f32) -> f32 {
     let s = clamp_f32(r);
-    ((1.0 + s).ln() - (1.0 - s).ln()) * 0.5
+    (ln_det(1.0 + s) - ln_det(1.0 - s)) * 0.5
+}
+
+/// Deterministic natural log for positive normal `f32`, bit-identical on
+/// every target.
+///
+/// Built only from IEEE add, sub, mul, div and bit operations, evaluated in a
+/// fixed order (Rust never contracts into FMA), so every conforming target
+/// produces the same bits. Algorithm: musl / FreeBSD `logf` — split off the
+/// exponent `k`, reduce the mantissa to `[√½, √2)`, and evaluate a degree-8
+/// polynomial in `s = f/(2+f)`. Error is below one ulp.
+///
+/// Domain: the Fisher path only ever passes `1 ± s` with `|s| ≤ 1 − 2⁻²⁴`,
+/// i.e. `[2⁻²⁴, 2)` — positive and normal. NaN propagates. Zero, negative,
+/// subnormal and infinite inputs are outside the domain and not handled.
+#[inline]
+fn ln_det(x: f32) -> f32 {
+    const LN2_HI: f32 = f32::from_bits(0x3f31_7180); // 6.9313812256e-01
+    const LN2_LO: f32 = f32::from_bits(0x3717_f7d1); // 9.0580006145e-06
+    const LG1: f32 = f32::from_bits(0x3f2a_aaaa); // 0.66666662693
+    const LG2: f32 = f32::from_bits(0x3ecc_ce13); // 0.40000972152
+    const LG3: f32 = f32::from_bits(0x3e91_e9ee); // 0.28498786688
+    const LG4: f32 = f32::from_bits(0x3e78_9e26); // 0.24279078841
+    if x.is_nan() {
+        return x;
+    }
+    let mut ix = x.to_bits();
+    ix = ix.wrapping_add(0x3f80_0000 - 0x3f35_04f3);
+    let k = (ix >> 23) as i32 - 0x7f;
+    ix = (ix & 0x007f_ffff) + 0x3f35_04f3;
+    let f = f32::from_bits(ix) - 1.0;
+    let s = f / (2.0 + f);
+    let z = s * s;
+    let w = z * z;
+    let t1 = w * (LG2 + w * LG4);
+    let t2 = z * (LG1 + w * LG3);
+    let r = t2 + t1;
+    let hfsq = 0.5 * f * f;
+    let dk = k as f32;
+    s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI
 }
 
 #[inline]
@@ -433,3 +485,74 @@ mod tests {
         assert_eq!(hamming_null_z(5, 0), 0.0);
     }
 }
+
+#[cfg(test)]
+mod golden {
+    use super::*;
+
+    fn fnv(bytes: impl IntoIterator<Item = u8>) -> u64 {
+        let mut h = 0xcbf2_9ce4_8422_2325u64;
+        for b in bytes {
+            h = (h ^ u64::from(b)).wrapping_mul(0x0100_0000_01b3);
+        }
+        h
+    }
+
+    fn grid() -> Vec<f32> {
+        // Dense in the centre, and walking every f32 bit pattern near the rim,
+        // where ln(1 − r) is most sensitive to the last ulp.
+        let mut v: Vec<f32> = (-4096..=4096).map(|i| i as f32 / 4096.0).collect();
+        let mut r = 0.999f32;
+        while r < 1.0 {
+            v.push(r);
+            v.push(-r);
+            r = f32::from_bits(r.to_bits() + 97);
+        }
+        v
+    }
+
+    /// `ln_det` stays within one ulp of the correctly rounded result over
+    /// its whole domain `[2⁻²⁴, 2)`, sampled every 257th bit pattern.
+    #[test]
+    fn ln_det_is_within_one_ulp() {
+        let (lo, hi) = (2f32.powi(-24).to_bits(), 2f32.to_bits());
+        let mut worst = 0u32;
+        let mut b = lo;
+        while b < hi {
+            let x = f32::from_bits(b);
+            let exact = (f64::from(x).ln()) as f32;
+            let got = ln_det(x);
+            let ulp = got.to_bits().abs_diff(exact.to_bits());
+            // The sign can differ only at x = 1, where both are ±0.
+            if exact != 0.0 {
+                worst = worst.max(ulp);
+            }
+            b += 257;
+        }
+        assert!(worst <= 1, "worst error {worst} ulp");
+        assert!(ln_det(f32::NAN).is_nan());
+        assert_eq!(ln_det(1.0), 0.0);
+    }
+
+    /// Pinned digests over a fixed grid (dense centre plus every 97th f32
+    /// bit pattern near the rim). The same constants must hold on every
+    /// target — x86-64 (all realizations), aarch64 and wasm32. Measured
+    /// before `ln_det`: with libm, 624 of the 8 539 z values differed between
+    /// x86-64 glibc and wasm32.
+    pub const GOLDEN_CODES: u64 = 0xfb0b_294d_2a65_3dbb;
+    pub const GOLDEN_Z32: u64 = 0x2e92_a08e_fd1f_cb69;
+
+    #[test]
+    fn digests_are_pinned() {
+        let g = grid();
+        assert_eq!(g.len(), 8539);
+        let env = ZGamma::fit(&g);
+        let mut codes = vec![0i8; g.len()];
+        env.encode_batch(&g, &mut codes);
+        assert_eq!(fnv(codes.iter().map(|&c| c as u8)), GOLDEN_CODES, "codes drifted");
+        let z = fnv(g
+            .iter()
+            .flat_map(|&r| fisher_z_f32(r).to_bits().to_le_bytes()));
+        assert_eq!(z, GOLDEN_Z32, "z bits drifted");
+    }
+}

From aa3dfd736569d21175c38314210e377cc118789d Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Wed, 23 Sep 2026 21:37:34 +0000
Subject: [PATCH 8/8] wasm-simd-parity: pin ZGamma codes digest under node

CI's wasm job never runs cargo test, so the zspace golden digest was not
enforced on wasm32. selfcheck() now rebuilds the same 8539-value grid,
encodes it and compares against GOLDEN_CODES (rc 0x501), plus batch == scalar
per element (0x502).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019HnekoM1EidTwQLS3oFVFm
---
 crates/wasm-simd-parity/src/lib.rs | 38 ++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)

diff --git a/crates/wasm-simd-parity/src/lib.rs b/crates/wasm-simd-parity/src/lib.rs
index b84c4a6c..247e74d5 100644
--- a/crates/wasm-simd-parity/src/lib.rs
+++ b/crates/wasm-simd-parity/src/lib.rs
@@ -35,6 +35,9 @@ pub extern "C" fn selfcheck() -> u32 {
     if let Err(code) = check_i32x16_compare() {
         return code;
     }
+    if let Err(code) = check_zgamma_golden() {
+        return code;
+    }
     0
 }
 
@@ -424,3 +427,38 @@ fn check_i32x16_compare() -> Result<(), u32> {
     }
     Ok(())
 }
+
+/// `ZGamma` codes are bit-exact on every target: the same grid and pinned
+/// digest as `ndarray::hpc::zspace::golden::GOLDEN_CODES` (keep in sync).
+/// Before the deterministic `ln`, 624 of these 8 539 z values differed between
+/// x86-64 glibc and wasm32, so this is the target where drift would show.
+fn check_zgamma_golden() -> Result<(), u32> {
+    use ndarray::hpc::zspace::ZGamma;
+    const GOLDEN_CODES: u64 = 0xfb0b_294d_2a65_3dbb;
+    let mut grid: Vec<f32> = (-4096..=4096).map(|i| i as f32 / 4096.0).collect();
+    let mut r = 0.999f32;
+    while r < 1.0 {
+        grid.push(r);
+        grid.push(-r);
+        r = f32::from_bits(r.to_bits() + 97);
+    }
+    if grid.len() != 8539 {
+        return Err(0x500);
+    }
+    let env = ZGamma::fit(&grid);
+    let mut codes = vec![0i8; grid.len()];
+    env.encode_batch(&grid, &mut codes);
+    let mut h = 0xcbf2_9ce4_8422_2325u64;
+    for &c in &codes {
+        h = (h ^ u64::from(c as u8)).wrapping_mul(0x0100_0000_01b3);
+    }
+    if h != GOLDEN_CODES {
+        return Err(0x501);
+    }
+    for (&c, &k) in grid.iter().zip(&codes) {
+        if env.encode(c) != k {
+            return Err(0x502);
+        }
+    }
+    Ok(())
+}