NumRS2 0.4.1 Released — O(1) Copy-on-Write Arrays, Real Distributed Collectives, and a Production-Hardening Pass

The N-dimensional array core just got harder to break — O(1) clones, a real distributed transport, and a long list of bugs that were quietly returning wrong numbers, all fixed and independently verified against the working tree.

Today we released NumRS2 0.4.1 — a production-hardening pass across the whole crate: Array<T> is now Arc-backed copy-on-write, a shared dispatch layer backs the hot paths, the distributed feature’s collectives run over a real network transport instead of returning fabricated local data, lapack is on by default, and a long list of correctness bugs — some of them silently wrong numbers, not crashes — got found and fixed.

No C. No Fortran. No system BLAS/LAPACK. No hand-rolled, per-call-site SIMD dispatch to get subtly wrong under maintenance. Just clean, blazing-fast N-dimensional arrays — broadcasting, fancy indexing, SVD, FFT, autodiff — that compile to a single static binary (or WASM) and run everywhere, backed by a compute core that’s now been through a real correctness audit.

Why NumRS2 0.4.1 is a game changer

Fast array libraries are common. Array libraries that are honest about where they used to be wrong are not. NumPy is the bedrock of scientific Python, but its C/Fortran core makes this kind of audit nearly impossible for an outside contributor to perform — you can’t easily prove an upstream BLAS kernel silently mishandles a shape mismatch. NumRS2 0.4.1 spent this cycle doing exactly that kind of audit on its own codebase, in Rust, with regression tests as the receipts:

• Clone got 1000x cheaper, safely. Array<T> is now Arc-backed copy-on-write — Clone is an Arc bump (O(1)), not a deep copy, and the first mutation after a clone unshares automatically. A CI policy check pins the crate to exactly one Arc::make_mut / Arc::try_unwrap call site, so the unsafe-adjacent sharing logic can’t quietly spread.
• The distributed collectives stopped lying. broadcast/reduce/allreduce/gather/allgather/scatter/allscatter/barrier used to return the caller’s own local data with no network transport at all. They now run for real, over a new TCP Endpoint transport with LZ4 compression — and a bidirectional deadlock in the channel layer was eliminated by construction, not papered over.
• Variance was silently wrong above length 64. Statistics::var/std flipped from population to sample variance (n vs. n-1) at array length 64 — an upstream SIMD-path artifact. Fixed to population semantics everywhere, matching NumPy’s default. This changes the numeric value var/std return on arrays ≥ 64 elements versus the previous (buggy) build.
• A shape mismatch was silently truncated, not errored. The x86_64 SIMD binary ops (e.g. vectorized_add_arrays_f64) truncated a shape mismatch to min(len_a, len_b) instead of raising an error — invisible in local development because the usual dev machine here is Apple Silicon, where those kernels never even type-check. Now an explicit ShapeMismatch error.
• lapack moved to the default feature set. det/inv/svd/eig/qr/cholesky and the rest of core linear algebra are reachable on a plain cargo add numrs2 — no more opt-in feature flag required for what most users expect to just be there.
• Matmul got faster, not just leaner. The default f32/f64 dispatch tier now routes through pure-Rust matrixmultiply with an M-only row-split parallelization — measured up to 18.82x faster at 512³ and 12.70x at 256³ against the prior loop, with no regressed shape.

Technical Deep Dive: From Shared Kernels to a Real Wire Transport

Layer 1 — The kernels dispatch layer. src/kernels/ centralizes what used to be scattered, hand-rolled SIMD/parallel thresholds across the hot paths into one dtype-dispatched module: elementwise (unary/binary dispatch), gemm (2D matmul dispatch table), reduce (sum/mean/var/min/max with a deterministic accumulation order), cast (sound TypeId-guarded reinterpretation, replacing ad hoc raw-pointer type-punning), and borrow (contiguous-vs-owned operand bridging). Every ad hoc mem::transmute_copy and x as *const T as *const f64 cast in the hot paths now goes through kernels::cast, which proves soundness via TypeId::of::<T>() before reinterpreting.

Layer 2 — Expression fusion. IntoExpr::expr() builds a lazy ExprNode tree from an Array/scalar/expression; .eval() walks it in one fused pass instead of materializing every intermediate array. It’s honest about scope — eager syntax (&a + &b) doesn’t fuse on its own, only the explicit .expr()...eval() builder does (a compile-time fused! macro is deferred to 0.6.0) — and it’s honest about precision too: fused and eager evaluation are guaranteed bit-for-bit equal for every finite, infinite, and signed-zero element, but not for a NaN’s payload/sign bits when a single operation combines two distinct NaNs. That’s measured, not assumed — there’s a regression test that pins the exact case where the one-pass fused loop and the two-pass eager spelling disagree under rustc -O.

Layer 3 — A real distributed transport. src/distributed/net/ implements Endpoint, a real point-to-point TCP transport with a fixed 56-byte frame header and LZ4 payload compression, and src/distributed/collective.rs now runs all eight collectives (plus v-variants and reduce_scatter) over it for real. TSQR (Tall-Skinny QR) and block-cyclic-column Cholesky land alongside it in src/distributed/linalg/, and a LocalCluster test harness drives real loopback-TCP multi-rank tests without needing a process launcher.

Layer 4 — NumPy-parity breadth. Ufunc reduce/accumulate/outer/reduceat/.at() (in-place scatter), fftn/ifftn/rfftn/irfftn with NumPy’s exact s=/axes=/norm= conventions, all 9 NumPy≥1.22 quantile/percentile methods, histogramdd with density=, a full masked-array completion pass (std/var/prod/median/argmin/argmax/sort/cumsum/dot/concatenate), the Chebyshev/Legendre/Hermite/HermiteE/Laguerre polynomial classes, and new SeedSequence + Generator::spawn, Philox4x64, and SFC64 random generators.

Layer 5 — The correctness audit itself. Beyond variance: min/max returning a wrong finite value for certain NaN placements (a live upstream scirs2_core bug, now bypassed and pinned as a tripwire test), min_along_axis/max_along_axis panicking on every axis reduction, tile()’s N-D shape/ordering, moveaxis’s permutation for arbitrary axis pairs, einsum panicking on a HashMap miss, two confirmed scirs2_fft::fftn normalization bugs worked around at the wrapper level, and more — each with its own regression test and its own line in CHANGELOG.md pointing at the fix.

Getting Started

cargo add numrs2

use numrs2::prelude::*;

fn main() -> Result<()> {
    // N-dimensional arrays with NumPy-style broadcasting
    let a = Array::from_vec(vec![1.0, 2.0, 3.0, 4.0]).reshape(&[2, 2]);
    let b = Array::from_vec(vec![5.0, 6.0, 7.0, 8.0]).reshape(&[2, 2]);

    // Clone is O(1) now — an Arc bump, not a deep copy
    let a_clone = a.clone();
    let e = a.matmul(&b)?;             // matrix multiply, now up to 18.82x faster at 512^3
    println!("a @ b = {}", e);

    // Expression fusion: builds a lazy tree, evaluates in one fused pass
    let fused = a.expr() * 2.0 + b.expr();
    println!("fused = {}", fused.eval()?);

    // Linear algebra — lapack is a default feature now, no opt-in needed
    let (u, s, vt) = a.svd_compute()?;
    println!("singular values = {}", s);
    let _ = (u, vt, a_clone);

    Ok(())
}

The array, matmul, fusion, and linear-algebra calls above are the safe, copy-pasteable core. The new distributed collectives, the masked-array completion, and the polynomial classes are reached through their respective modules — see the Tips below and CHANGELOG.md for the full per-module list.

What’s New in 0.4.1

Added

• kernels dispatch module (elementwise, gemm, reduce, cast, borrow), exercised by three new dispatch benchmarks.
• Expression fusion (IntoExpr::expr() + .eval()).
• Ufunc reduce/accumulate/outer/reduceat/.at(); fftn/ifftn/rfftn/irfftn; 9 NumPy≥1.22 quantile methods; histogramdd with density=; multi_dot/tensorsolve/tensorinv; masked-array completion; Chebyshev/Legendre/Hermite/HermiteE/Laguerre polynomial classes; SeedSequence/Generator::spawn/Philox4x64/SFC64.
• Real TCP Endpoint transport, 8 distributed collectives + v-variants, TSQR, block-cyclic Cholesky, LocalCluster test harness.
• GPU norm_l1, N-D transpose/broadcast/slicing, Conv2D; Python N-D neural-network + random/fft bindings; WASM dlmalloc allocator.

Changed

• Array<T> is Arc-backed copy-on-write (O(1) Clone). Disclosure: now requires T: Send + Sync, which can newly constrain generic call sites.
• Default matmul dispatch moved to pure-Rust matrixmultiply — up to 18.82x faster at 512³.
• lapack is now a default feature.
• Random-distribution sampling migrated to scirs2_core::random. Breaking for direct callers of NonCentralChiSquared::sample/NonCentralF::sample/VonMises::sample/Maxwell::sample/Wald::sample, which now take an explicit &mut StdRng.
• scirs2-* 0.5.0 → 0.6.5; oxiarc-* 0.3.2 → 0.4.1; oxicode 0.2.4 → 0.2.6; wgpu 29 → 30; pyo3 0.28 → 0.29.

Fixed

• var/std sample/population flip at length ≥ 64 (numeric-value disclosure above).
• min/max/min_along_axis/max_along_axis NaN and panic bugs; tile/moveaxis/as_strided/einsum correctness; lstsq on non-square input; two scirs2_fft::fftn normalization bugs worked around.
• SIMD dispatch functions in src/simd_optimize/ now return Result<Array<T>> instead of a bare Array<T> (53 signatures). Breaking for direct callers — closes a silent shape-mismatch truncation and a fallback-path panic.
• A bidirectional deadlock in the distributed CommunicationChannel, eliminated by construction (split read/write halves instead of one shared lock).
• macOS cargo build --features python linking, via a new build.rs calling pyo3_build_config::add_extension_module_link_args().

Tips

• Reach for .expr()...eval() on multi-op chains, not eager &a + &b. Fusion only happens through the explicit builder — let e = a.expr() * 2.0 + b.expr(); e.eval()? — and it’s worth up to 1.69x at n=1,000,000 in the module’s own benchmark table. Eager syntax stays exactly as fast (and as correct) as before; it just doesn’t fuse.
• If you call the SIMD functions in numrs2::simd_optimize directly, check the 0.4.1 signature. They now return Result<Array<T>> instead of a bare Array<T> — the migration is usually just adding a ?, and it buys you a real ShapeMismatch error instead of a silently truncated result.
• Re-check any stored var/std result computed on an array of 64+ elements against a pre-0.4.1 build. The population/sample fix changes the returned number, not just an edge case — see the Changed section above before diffing against cached results.
• det/inv/svd/eig/qr/cholesky no longer need an opt-in feature. lapack is a default feature now — drop the explicit features = ["lapack"] from your Cargo.toml if you had it.
• Use Array::is_unique() before a hot mutation loop if you want to know whether you’re about to pay the unshare cost. Since Clone is now O(1) and mutation unshares lazily, a loop over many cloned handles to the same buffer pays Arc::make_mut’s deep copy on the first write, not the clone.
• The distributed collectives are real now — test them with LocalCluster, not mocks. src/distributed/testing.rs’s harness drives world_size copies of an async closure over real loopback TCP, which is what NumRS2’s own regression suite uses to catch things like the bidirectional-deadlock fix above.

This is the foundation

NumRS2 is the NumPy-class N-dimensional array core at the base of the COOLJAPAN scientific stack, and 0.4.1 is the release where the compute paths underneath every higher layer got audited and hardened rather than extended. It sits directly beneath SciRS2 (now at v0.6.5 across the full scirs2-* family), with linear algebra on OxiBLAS, serialization on OxiCode, and compression on OxiArc — and it’s a real, in-production dependency: QuantRS2 and VoiRS both pull numrs2 directly today. Around it, OptiRS and PandRS sit alongside for optimization and dataframes; the ML and applied stack spans ToRSh, sklears, and trustformers; and the acceleration tier reaches from OxiCUDA for GPU compute to OxiMedia and OxiGDAL for media and geospatial.

Repository: https://github.com/cool-japan/numrs

Star the repo if you want a NumPy-class array core that’s not just fast, but auditable — O(1) clones, a real distributed transport, and a changelog that tells you exactly which numbers used to be wrong and why.

The era of trusting an opaque C/Fortran core is over. Pure Rust numerical computing — hardened, not just extended — is here.

— KitaSan at COOLJAPAN OÜ August 29, 2026