Memory-safe, clean-room NumPy reimplementation in Rust.
Zero unsafe code. Extensive conformance coverage. Bit-exact RNG parity. Every feature family green.

NumPy is the backbone of scientific Python, but it carries 30 years of C/C++ legacy: buffer overflows in parsers, undefined behavior in edge cases, opaque stride semantics, and no formal compatibility contracts. Rewriting pieces in Cython or C++ perpetuates the same class of memory bugs.

FrankenNumPy rebuilds NumPy's semantics from scratch in safe Rust with two non-negotiable goals:

1. Absolute behavioral compatibility with legacy NumPy. Not a subset, not "inspired by." The full API, edge cases and all.
2. Rigorous architecture with formal contracts, dual-mode runtime (strict/hardened), and differential conformance against a NumPy oracle.

use fnp_ufunc::{UFuncArray, DType, BinaryOp};

// Create arrays
let data = UFuncArray::new(vec![5], vec![1.0, 2.0, 3.0, 4.0, 5.0], DType::F64)?;

// Z-score normalization: (data - mean) / std
let mean = data.reduce_mean(None, false)?;
let std = data.reduce_std(None, false, 0)?;
let z = data.elementwise_binary(&mean.broadcast_to(&[5])?, BinaryOp::Sub)?
            .elementwise_binary(&std.broadcast_to(&[5])?, BinaryOp::Div)?;
// z = [-1.414, -0.707, 0.0, 0.707, 1.414]

// Sort, cumsum, percentile chain
let sorted = data.sort(None, None)?;           // [1, 2, 3, 4, 5]
let cumsum = sorted.cumsum(None)?;             // [1, 3, 6, 10, 15]
let median = data.percentile(50.0, None)?;     // 3.0

// Linear algebra
let a = UFuncArray::new(vec![2, 2], vec![3.0, 1.0, 1.0, 2.0], DType::F64)?;
let b = UFuncArray::new(vec![2], vec![9.0, 8.0], DType::F64)?;
let x = a.solve(&b)?;                         // [2.0, 3.0]

use fnp_random::Generator;

// Bit-exact NumPy-compatible RNG
let mut rng = Generator::from_pcg64_dxsm(12345)?;
let normals = rng.standard_normal(1000);        // Identical to NumPy's output
let samples = rng.binomial(10, 0.5, 100);       // BTPE algorithm, same as NumPy
let perm = rng.permutation(&data)?;             // Fisher-Yates via random_interval

Parity debt, not feature cuts. Every behavioral gap with NumPy is tracked as debt to be closed, not as an accepted scope reduction.

Stride Calculus Engine (SCE). All shape transformations (broadcast, reshape, transpose, view aliasing) flow through a single deterministic engine. This is the non-negotiable compatibility kernel.

Dual-mode runtime. Strict mode maximizes observable NumPy compatibility. Hardened mode adds safety guards and bounded defensive recovery. Both modes log decisions to an evidence ledger.

Fail-closed by default. Unknown wire formats, unrecognized metadata, and incompatible features cause explicit errors, not silent corruption.

Oracle-verified. Every RNG distribution, every linalg decomposition, and every reduction edge case is tested against NumPy's actual output from the same seed.

Over 1,000 public functions across 9 crates covering the full NumPy API:

# Clone and build
git clone https://github.com/Dicklesworthstone/franken_numpy.git
cd franken_numpy
cargo build --workspace

# Run all tests
cargo test --workspace

# Run with all features
cargo test --workspace --all-features

Requirements: Rust nightly (pinned in rust-toolchain.toml to nightly-2026-02-20).

           ┌──────────────────────────────────────────────┐
           │              User API Layer                  │
           │    UFuncArray · MaskedArray · StringArray    │
           └───────────────────────┬──────────────────────┘
                                  │
      ┌───────────┬───────────────┼──────────────┬───────────┐
      ▼           ▼               ▼              ▼           ▼
 ┌──────────┐ ┌──────────┐ ┌──────────┐ ┌──────────┐ ┌──────────┐
 │ fnp-ufunc│ │fnp-linalg│ │fnp-random│ │  fnp-io  │ │ fnp-fft  │
 │ array ops│ │ linalg   │ │ RNG      │ │ NPY/NPZ  │ │ (in-ufunc│
 │ dispatch │ │ kernels  │ │ engine   │ │ + text   │ │  module) │
 └─────┬────┘ └──────────┘ └──────────┘ └──────────┘ └──────────┘
       │
 ┌─────┴────────────────────────────────────────────────────┐
 │                                                          │
 ▼                                                          ▼
 ┌──────────┐ ┌────────────┐ ┌──────────┐ ┌──────────────┐
 │ fnp-dtype│ │ fnp-ndarray│ │ fnp-iter │ │ fnp-runtime  │
 │ dtype    │ │ SCE core   │ │ transfer │ │ strict/      │
 │ rules    │ │ strides    │ │ semantics│ │ hardened     │
 └──────────┘ └────────────┘ └──────────┘ └──────────────┘
                    ▲
                    │
         Stride Calculus Engine (SCE)
         - shape -> strides (C/F)
         - broadcast legality
         - reshape with -1 inference
         - alias-safe view transitions

The Stride Calculus Engine is the non-regression kernel. Every shape transformation (broadcast, reshape, transpose, view creation) is validated through SCE before execution.

The type system is the foundation everything else builds on. 18 dtype variants cover the full NumPy type hierarchy:

Bool  I8  I16  I32  I64  U8  U16  U32  U64
F16  F32  F64  Complex64  Complex128
Str  DateTime64  TimeDelta64  Structured

Each dtype maps to a type-safe ArrayStorage variant with native Rust containers. No flattened Vec<u8> reinterpretation. I64 values live in Vec<i64>, Complex128 values live in Vec<(f64, f64)>, and F16 values use the half crate's f16 type. This means integer fidelity is preserved (no silent truncation through f64 for i64 values > 2^53), f32 identity is maintained, and complex numbers are stored as native interleaved pairs.

Promotion table. The promote(lhs, rhs) function implements NumPy's exact promotion rules as a deterministic const fn:

The U64+signed→F64 promotion is the most counterintuitive rule in the table, but it's what NumPy does: there simply isn't a 128-bit integer type in NumPy's type system.

Cast policy. can_cast(src, dst, casting) supports all five NumPy casting modes: "no", "equiv", "safe", "same_kind", and "unsafe".

The SCE owns all shape transformation rules and is the correctness backbone of the entire project:

1. Shape → strides. Given a shape [d0, d1, ..., dn] and memory order (C or F), compute contiguous strides. C-order: rightmost dimension has stride = item_size, each dimension to the left multiplies by the next dimension's size. F-order: leftmost dimension has stride = item_size.

2. Broadcast legality. Two shapes are broadcast-compatible when, aligned from the right, each pair of dimensions is either equal or one is 1. The output shape takes the maximum at each position. FrankenNumPy handles mixed-rank broadcasting (different number of dimensions) by left-padding the shorter shape with 1s.

3. Reshape with -1 inference. At most one dimension can be -1, meaning "infer from the total element count and the other dimensions." The engine divides the total element count by the product of known dimensions and validates that the result is exact (no remainder).

4. View safety. as_strided() and sliding_window_view() compute the required byte span from the minimum to maximum reachable offset, then verify it fits within the base allocation. Negative strides (for reverse slicing) are fully supported.

5. Broadcast strides. When a dimension has size 1 but the broadcast output has size > 1, the stride for that dimension is set to 0. This creates a "virtual repeat" without copying data.

The runtime implements a risk-aware decision engine with posterior probability estimation that chooses between three actions for each compatibility check:

Risk scoring with posterior estimation. Each decision starts with class-specific prior odds (1% incompatibility for KnownCompatible, 99% for KnownIncompatible, 50% for Unknown). A risk score is computed as the log-likelihood ratio of the evidence against a threshold. The posterior incompatibility probability determines the action via a loss model.

The loss model encodes the asymmetric costs:

allow_if_compatible:       0.0    (correct accept: no cost)
allow_if_incompatible:   100.0    (silent corruption: catastrophic)
full_validate_if_compatible: 4.0  (unnecessary work: small cost)
full_validate_if_incompatible: 2.0 (caught by validation: acceptable)
fail_closed_if_compatible: 125.0  (false rejection: high cost)
fail_closed_if_incompatible: 1.0  (correct rejection: minimal cost)

Every decision is logged to an EvidenceLedger with timestamp, mode, class, evidence terms, and action taken. Override requests are tracked separately in OverrideAuditEvent records with audit references.

The ufunc engine is the largest crate (~30,000 lines) and implements the array operation layer:

35 binary operations — Add, Sub, Mul, Div, Power, FloorDivide, Remainder, Fmod, Minimum, Maximum, Fmax, Fmin, Copysign, Heaviside, Nextafter, Arctan2, Hypot, Logaddexp, Logaddexp2, Ldexp, FloatPower, BitwiseAnd/Or/Xor, LeftShift, RightShift, LogicalAnd/Or/Xor, and 6 comparison ops.

42 unary operations — Abs, Negative, Sqrt, Exp, Log, all trig and hyperbolic functions, Floor, Ceil, Round, Cbrt, Expm1, Log1p, Spacing, Signbit, Isnan, Isinf, Isfinite, Invert, and more.

Every binary operation handles broadcasting through the same code path: compute the output shape via SCE's broadcast_shape, map output indices to source indices using broadcast strides (0-stride for size-1 dimensions), and apply the operation elementwise.

NaN semantics are explicit. After fixing 18 NaN-handling bugs found through systematic oracle testing:

• reduce_min, reduce_max: propagate NaN (use custom nan_min/nan_max, not f64::min/f64::max)
• cummin, cummax: NaN propagates through the accumulator
• sort, argsort: NaN sorts to end via nan_last_cmp (not partial_cmp().unwrap_or(Equal))
• median, percentile, quantile: NaN early-returns before sorting
• ptp: NaN propagates in both UFuncArray and MaskedArray variants

Float error handling. A thread-local FloatErrorState tracks divide-by-zero, overflow, underflow, and invalid operation events. Each mode (Ignore, Warn, Raise, Call, Print, Log) can be configured independently. seterr(), geterr(), and errstate() match NumPy's API.

Generalized ufuncs. GufuncSignature parses signatures like "(n?,k),(k,m?)->(n?,m?)" for operations that act on sub-arrays (e.g., matrix multiply). The parser normalizes signatures, extracts core dimensions, and validates operand shapes against the pattern.

einsum. Three implementations: basic einsum() for direct evaluation, einsum_path() for contraction path optimization, and einsum_optimized() for greedy or brute-force strategy selection.

The RNG crate achieves bit-exact parity with NumPy by porting every algorithm from NumPy's C source code.

Bit generators. Three implementations:

• Pcg64DxsmRng — the default. State is 128-bit with 128-bit increment. Uses the DXSM output function and cheap multiplier for generation. Seeding goes through SeedSequence which applies Melissa O'Neill's seed_seq design with SplitMix64 mixing.
• Mt19937Rng — Mersenne Twister for legacy RandomState compatibility. 624-word state array with Matsumoto-Nishimura twist.
• DeterministicRng — simple SplitMix64-based generator for testing.

SeedSequence. Hierarchical seeding with spawn/generate_state contracts. spawn(n) creates child SeedSequences by incorporating a monotonic spawn counter into the entropy pool. generate_state(words) cycles through the pool with hash transformations to produce initialization vectors. This exactly matches NumPy's numpy.random.SeedSequence.

Bounded integers. NumPy's random_bounded_uint64() dispatch:

• Range fits in 32 bits → Lemire's method via buffered next_uint32() (each u64 is split into two u32s, low first)
• Range exceeds 32 bits → Lemire's method with 128-bit multiplication
• Range = 0 → return 0; Range = u32::MAX → raw next_uint32(); Range = u64::MAX → raw next_u64()

Shuffle/permutation. Uses random_interval() (masked bit rejection, not Lemire), matching NumPy's _shuffle_raw code path. For each index from n-1 down to 1, generate a uniform random integer in [0, i] by masking to the smallest bit-width covering i and rejecting values > i.

92 public functions organized into three tiers:

2x2 fast paths. solve_2x2, det_2x2, inv_2x2, qr_2x2, svd_2x2, eigh_2x2, cholesky_2x2 avoid the overhead of general NxN algorithms for the smallest matrix size.

NxN general algorithms:

• QR decomposition — Householder reflections with qr_nxn (reduced) and qr_mxn (rectangular)
• SVD — Golub-Kahan bidiagonalization + implicit shifted QR
• Eigenvalues — Hessenberg reduction + implicit shifted QR iteration for real Schur form
• Symmetric eigenvalues — Tridiagonal reduction + implicit QL shifts (eigvalsh_nxn, eigh_nxn)
• LU factorization — Partial pivoting with lu_factor_nxn and lu_solve
• Cholesky — Column-wise lower-triangular factorization
• Least squares — Normal equations via A^T A for lstsq_nxn

Spectral methods: expm_nxn (matrix exponential via Pade approximation), sqrtm_nxn (matrix square root), logm_nxn (matrix logarithm), funm_nxn (general matrix function via Schur decomposition), polar_nxn (polar decomposition U*P), schur_nxn (Schur triangular form).

Batch operations. 14 batch functions (batch_inv, batch_det, batch_solve, batch_svd, batch_qr, batch_cholesky, batch_eig, batch_eigvalsh, batch_eigh, batch_slogdet, batch_svd_full, batch_matrix_norm, batch_matrix_rank, batch_trace) operate on stacked matrices with leading batch dimensions, matching NumPy's broadcasting semantics for linear algebra.

Complex number support. 10 functions for complex-valued matrices: complex_solve_nxn, complex_det_nxn, complex_inv_nxn, complex_cholesky_nxn, complex_qr_mxn, complex_matmul, complex_matvec, complex_conjugate_transpose, complex_matrix_norm_frobenius, complex_trace_nxn.

NPY format. Complete implementation of NumPy's .npy binary format:

• Magic prefix \x93NUMPY + version (1.0 or 2.0) + header length + Python dict header + padding to 16-byte alignment + raw data payload
• 24 supported dtype descriptors including little/big-endian variants for all integer and float types, complex types, fixed-width byte strings, Unicode strings, and object type
• Header parsing validates all three required fields (descr, fortran_order, shape) and rejects unknown keys
• Hardened with MAX_HEADER_BYTES = 65,536 to prevent allocation bombs

NPZ format. ZIP-based container for multiple arrays:

• savez (stored) and savez_compressed (DEFLATE) via the flate2 crate
• MAX_ARCHIVE_MEMBERS = 4,096 and MAX_ARCHIVE_UNCOMPRESSED_BYTES = 2 GB limits prevent archive bombs
• Each member is a complete .npy file within the ZIP

Text I/O. loadtxt, savetxt, and genfromtxt handle delimiter-separated text files with configurable dtypes, missing value handling, and column selection (usecols).

Pickle policy. Object dtype arrays require explicit opt-in via allow_pickle=true, matching NumPy's security posture against arbitrary code execution from untrusted .npy files.

The conformance crate is the quality backbone with four layers:

1. Differential harness. Captures NumPy's output for a corpus of input fixtures, then runs the same inputs through FrankenNumPy and compares shapes, dtypes, and values (with configurable tolerance). Covers ufunc, linalg, FFT, polynomial, string, masked array, datetime, RNG, and I/O operations.

2. Metamorphic testing. Verifies algebraic identities that must hold regardless of input: a + b = b + a, a * 1 = a, sum(a) = sum(sort(a)), etc. 13 identities tested.

3. Adversarial fuzzing. Tests behavior on hostile inputs: NaN-filled arrays, extreme shapes (0-d, empty, very large), denormalized floats, integer overflow, and malformed NPY headers.

4. Witness stability. Hardcoded expected values for every RNG distribution ensure that code changes don't silently alter output sequences. When an algorithm is intentionally changed (e.g., porting Lemire's method), the witness values are regenerated from the new implementation.

FrankenNumPy's security posture covers more than memory safety:

• Zero unsafe Rust. All 9 crates declare #![forbid(unsafe_code)]. The 92,000+ lines of Rust contain zero unsafe blocks.
• Fail-closed by default. Unknown wire formats, unrecognized dtype descriptors, and metadata schema violations cause explicit errors, not silent fallbacks.
• Bounded resource consumption. NPY header parsing caps at 64 KB. NPZ archives cap at 4,096 members and 2 GB uncompressed. Memmap validation retries cap at 64. These prevent denial-of-service via crafted inputs.
• Pickle rejection. Object dtype arrays that could execute arbitrary code during deserialization require explicit opt-in, matching NumPy's allow_pickle security gate.
• Adversarial conformance. The security gate (run_security_gate) tests exploit scenarios from a versioned threat matrix mapped to specific parser/IO/shape-validation boundaries.

Each subsystem uses specific numerical algorithms. This section catalogs them.

Hybrid approach supporting arbitrary input lengths:

• Power-of-two lengths: Cooley-Tukey decimation-in-time (DIT). Recursively splits into even/odd subsequences, applies butterfly operations with twiddle factors exp(-2*pi*i*k/N).
• Non-power-of-two lengths: Bluestein's chirp-Z transform. Rewrites the DFT as a convolution, zero-pads to the next power of two, uses Cooley-Tukey for the convolution, then extracts the result.

All transforms (fft, ifft, fft2, ifft2, fftn, ifftn, rfft, irfft) build on these two primitives. Multi-dimensional transforms apply 1-D FFTs sequentially along each axis. fftfreq and rfftfreq compute the frequency bin centers; fftshift/ifftshift rearrange zero-frequency to the center.

• convolve(a, v): Direct O(n*m) convolution, full mode (output length n+m-1). Flips the kernel and slides it across the input.
• correlate(a, v): Cross-correlation implemented as convolve(a, v[::-1]).
• convolve2d / correlate2d: Full 2-D convolution with output shape (h1+h2-1, w1+w2-1).

gradient(f, *varargs) computes numerical derivatives with configurable edge handling:

Non-uniform spacing uses Lagrange polynomial coefficients (3-point at edges, 2-point for edge_order=1).

diff(a, n) computes n-th discrete difference: diff[i] = a[i+1] - a[i], applied n times. Output length decreases by n.

interp(x, xp, fp) does 1-D piecewise linear interpolation: binary search to find the enclosing interval in xp, then fp[lo] * (1-t) + fp[hi] * t where t = (x - xp[lo]) / (xp[hi] - xp[lo]). Values outside xp are clamped to the boundary values.

Five window types for spectral analysis, all returning [1.0] for M <= 1:

histogram(a, bins) supports three automatic bin-count strategies:

Bin edges are uniformly spaced: min + i * (max - min) / bins for i in 0..=bins. Element-to-bin assignment uses O(log bins) binary search per element. histogram_bin_edges returns just the edges without counting. histogramdd generalizes to N dimensions.

11 pad modes matching NumPy's np.pad:

Ten time-value-of-money functions using closed-form annuity algebra where possible:

Five complete polynomial families, each with evaluation, arithmetic, calculus, and fitting:

All five polynomial families support the full suite of evaluation, arithmetic, calculus, root-finding, basis conversion, and fitting operations.

MaskedArray wraps a UFuncArray with an optional boolean mask and fill value:

MaskedArray { data: UFuncArray, mask: Option<UFuncArray>, fill_value: f64, hard_mask: bool }

Mask convention: 1.0 = masked (excluded from computation), 0.0 = valid. This matches NumPy's numpy.ma module. Operations that reduce masked arrays (sum, mean, min, max, var, std, median, ptp, argmin, argmax, cumsum, cumprod, count) skip masked values. Comparison and arithmetic operations propagate masks through mask_or. compressed() returns a 1-D array of only the unmasked values. filled(fill_value) replaces masked values with a fill value.

DatetimeArray and TimedeltaArray support temporal arithmetic:

• Datetime - Datetime = Timedelta
• Datetime + Timedelta = Datetime
• Timedelta + Timedelta = Timedelta
• Scalar multiplication of Timedelta

Business day functions: busday_count(start, end) counts business days between dates, busday_offset(date, offset) adds business days to a date, is_busday(date) checks if a date is a business day. Weekday mask is Monday-Friday.

33 numpy.char functions operating elementwise on string arrays: upper, lower, capitalize, title, center, ljust, rjust, zfill, strip, lstrip, rstrip, replace, find, rfind, count, startswith, endswith, isnumeric, isalpha, isdigit, isdecimal, str_len, encode, decode, translate, maketrans, partition, rpartition, split, rsplit, join, expandtabs, swapcase. String add concatenates elementwise; string multiply repeats.

• packbits(axis): Packs 8 boolean elements into 1 byte, MSB first. Output length = ceil(axis_len / 8).
• unpackbits(axis, count): Unpacks bytes into boolean bits, MSB to LSB. Output length = axis_len * 8 (or count if specified).

The einsum implementation supports three contraction strategies:

einsum_path returns the contraction order without executing, for inspection. einsum_optimized applies the selected strategy.

Every conformance artifact (fixture bundles, benchmark baselines, migration manifests) is protected by erasure-coding sidecars:

Encoding. Source data is hashed (SHA-256) and encoded into source + repair symbols using RaptorQ fountain codes. The sidecar records the codec parameters (symbol size, block count, repair overhead) alongside the encoded symbols in base64.

Scrubbing. A scrub report decodes all symbols, computes the SHA-256 of the decoded payload, and verifies it matches the source hash. It then drops one symbol and verifies that recovery from the remaining symbols still produces the correct hash.

Decode proof. An explicit artifact recording which symbol was dropped, how many repair symbols were needed to recover, and whether recovery succeeded. This provides machine-checkable evidence that the artifact can survive single-symbol loss.

The G8 CI gate (run_raptorq_gate.sh) enforces that all required bundles have valid sidecars, scrub reports with status: "ok", and decode proofs with recovery_success: true.

12 threat classes are formally mapped in security_control_checks_v1.yaml, each with assigned conformance suites, compatibility gates, and override audit policies:

Every threat log entry must include: fixture_id, seed, mode, env_fingerprint, artifact_refs, reason_code.

UFuncArrayView provides NumPy-style shared-memory views with overlap detection:

UFuncArrayView {
    shape: Vec<usize>,
    buffer: Arc<RwLock<Vec<f64>>>,  // shared backing store
    offset: isize,                   // byte offset into buffer
    strides: Vec<isize>,            // per-dimension strides (can be negative)
    writable: bool,
    dtype: DType,
}

Memory overlap detection uses a two-tier approach:

1. Fast path (may_share_memory): Checks whether the byte-offset spans of two views overlap. This is O(ndim) and conservative (may report false positives for non-contiguous views).

2. Exact path (shares_memory): First checks Arc pointer equality (different backing buffers never share). If same buffer, computes the actual set of accessed byte offsets (up to 200,000) and checks for intersection. Falls back to the fast path if offset collection exceeds the limit.

This supports safe in-place operations: if two views share memory, operations that read from one and write to the other must use temporary copies to avoid data corruption.

FrankenNumPy replicates NumPy's floating-point error handling system:

┌───────────────┐   seterr(divide='raise')   ┌──────────────┐
│ Default       │ ────────────────────────── │ Custom       │
│ divide=Warn   │                            │ divide=Raise │
│ over=Warn     │   errstate(all='ignore')   │ over=Warn    │
│ under=Ignore  │ ────────────────────────── │ under=Ignore │
│ invalid=Warn  │   (RAII guard restores)    │ invalid=Warn │
└───────────────┘                            └──────────────┘

Six error modes per category: Ignore (suppress), Warn (log and continue), Raise (return error), Call (invoke user callback), Print (stderr), Log (append to event buffer).

Four error categories: divide-by-zero, overflow, underflow, invalid operation.

errstate() returns an RAII guard that automatically restores the previous error configuration when dropped, matching NumPy's context manager semantics:

let _guard = errstate(Some(FloatErrorMode::Ignore), None, None, None, None);
// all float errors ignored in this scope
// previous state restored when _guard drops

The scimath module provides 8 functions that extend real-valued math to the complex domain for inputs outside the real function's natural domain:

These match numpy.lib.scimath and are useful in signal processing and physics where negative square roots or out-of-range inverse trig values arise naturally.

Random-generation methods available on Generator, grouped by family:

Continuous: beta, chisquare, exponential, f/f_distribution, gamma, gumbel, halfnormal, laplace, levy, logistic, lognormal, lomax, maxwell, noncentral_chisquare, noncentral_f, normal, pareto, power, rayleigh, standard_cauchy, standard_exponential, standard_gamma, standard_normal, standard_t, triangular, vonmises, wald, weibull

Discrete: binomial, geometric, hypergeometric, logseries, negative_binomial, poisson, zipf

Multivariate: dirichlet, multinomial, multivariate_normal

Uniform: random (float [0,1)), uniform (float [low,high)), integers (int [low,high))

Permutation: shuffle (in-place), permutation (copy), permuted (axis-aware)

Utility: bytes, choice/choice_weighted

State: spawn, jumped, state/set_state

FrankenNumPy implements the full set of NumPy's array construction and manipulation functions.

The iterator crate models NumPy's internal data transfer system, which decides how to move data between arrays during operations:

Transfer classes determine the copy strategy:

Overlap detection decides copy direction:

FlatIter indexing supports four modes: Single(i) for scalar access, Slice{start, stop, step} for regular ranges, Fancy(Vec<usize>) for arbitrary index arrays, and BoolMask(Vec<bool>) for boolean selection. The count_true_mask optimization processes boolean masks in 8-element chunks for vectorizable counting.

Stateful NDIter wrapper. fnp-iter now exposes a first-class Nditer state machine on top of NditerPlan, including iterindex, multi_index, reset, seek-by-index/seek-by-multi-index, and external_loop chunk iteration. It also ships a NumPy-backed nditer_python* bridge for parity checks against Python-side stepping behavior, and fnp-python now exposes that iterator state machine as a Python extension-module PyNditer.

The conformance system is organized around 9 extraction packets, each covering one domain of NumPy behavior:

Each packet produces 8 artifact files: legacy_anchor_map.md, contract_table.md, fixture_manifest.json, parity_gate.yaml, risk_note.md, parity_report.json, parity_report.raptorq.json, parity_report.decode_proof.json. The packet readiness validator checks all 8 files exist and contain required fields before a packet is marked ready.

The random number generator supports full state capture and restoration for reproducibility:

// Capture state
let payload = generator.to_pickle_payload();

// Restore state
let restored = Generator::from_pickle_payload(payload)?;
// restored produces identical sequence from this point

GeneratorPicklePayload captures the bit-generator state (seed, counter, algorithm tag, schema version) and optionally the SeedSequence snapshot for spawn lineage tracking. The RandomState wrapper provides legacy compatibility with NumPy's older numpy.random.RandomState API.

Seed material accepts multiple forms: None (random), U64(seed), U32Words(vec), SeedSequence, or direct State { seed, counter } for exact state restoration. This matches NumPy's flexible seeding interface where default_rng(12345), default_rng([1, 2, 3]), and default_rng(SeedSequence(42)) all work.

Current per-crate counts move over time; see FEATURE_PARITY.md for the live inventory. Coverage focus by crate:

The oracle suite verifies bit-exact or behaviorally equivalent parity against NumPy across the major subsystems:

• RNG oracle tests: every distribution produces identical output from the same seed, verified against numpy.random.Generator(PCG64DXSM(12345))
• Ufunc edge-case tests: NaN propagation, empty arrays, Inf arithmetic, boolean dtype promotion, and sort ordering
• Linalg oracle tests: det, inv, solve, eig, svd, cholesky, qr, norm, cond, slogdet, lstsq, and rank checks
• I/O format tests: NPY roundtrip, magic bytes, header dict format, and 16-byte alignment

Every algorithm was ported from NumPy's C source code:

Every feature family is parity_green:

See FEATURE_PARITY.md for the full matrix with evidence links.

Eight ordered gates run from fast to heavy:

G1  fmt + lint           cargo fmt --check && cargo clippy -- -D warnings
G2  unit/property        cargo test --workspace --lib
G3  oracle differential  capture_numpy_oracle → run_ufunc_differential
G4  adversarial/security run_security_policy_gate.sh
G5  test/logging contract run_test_contract_gate.sh
G6  workflow e2e         run_workflow_scenario_gate.sh
G7  performance budget   run_performance_budget_gate.sh
G8  durability/decode    run_raptorq_gate.sh

Run all gates:

scripts/e2e/run_ci_gate_topology.sh

For manual direct rch workspace checks in this repo, prefer disabling Cargo incremental artifacts to avoid noisy remote artifact-retrieval races under target/debug/incremental:

CARGO_INCREMENTAL=0 rch exec -- cargo check --workspace --all-targets
CARGO_INCREMENTAL=0 rch exec -- cargo clippy --workspace --all-targets -- -D warnings

The oracle capture pipeline runs real NumPy, captures its output, and compares against our implementation:

# Recommended: require a real NumPy oracle and let the runner bootstrap
# a repo-local .venv-numpy314 with uv if needed.
FNP_REQUIRE_REAL_NUMPY_ORACLE=1 \
  cargo run -p fnp-conformance --bin capture_numpy_oracle

# Run differential comparison
cargo run -p fnp-conformance --bin run_ufunc_differential

If you want to manage the interpreter explicitly, this still works:

uv venv --python 3.14 .venv-numpy314
uv pip install --python .venv-numpy314/bin/python numpy
FNP_ORACLE_PYTHON="$(pwd)/.venv-numpy314/bin/python3" \
  cargo run -p fnp-conformance --bin capture_numpy_oracle

franken_numpy/
├── Cargo.toml                         # Workspace root (9 crates)
├── FEATURE_PARITY.md                  # Live parity matrix + evidence links
├── PROPOSED_ARCHITECTURE.md           # High-level architecture notes
├── crates/
│   ├── fnp-dtype/                     # Dtype taxonomy, promotion table, cast policy
│   ├── fnp-ndarray/                   # Stride Calculus Engine (SCE)
│   ├── fnp-iter/                      # Transfer semantics, overlap-safe iteration
│   ├── fnp-ufunc/                     # 1,000+ array operations, reductions, einsum
│   ├── fnp-linalg/                    # solve, eig, svd, qr, cholesky, lstsq, etc.
│   ├── fnp-random/                    # Statistical distributions + generator helpers, PCG64DXSM, bit-exact parity
│   ├── fnp-io/                        # NPY/NPZ read/write, text I/O, DEFLATE
│   ├── fnp-conformance/               # Oracle capture, differential harness, benchmarks
│   └── fnp-runtime/                   # Strict/hardened mode, evidence ledger
├── legacy_numpy_code/numpy/           # Behavioral oracle (upstream NumPy source)
├── artifacts/                         # Contracts, security maps, logs, proofs
└── scripts/                           # E2E gate scripts

FrankenNumPy prioritizes correctness over speed in this phase, but the architecture is designed for future optimization:

• Release profile: opt-level = 3, LTO, single codegen unit, stripped symbols.
• Contiguous reduction kernel: Axis reductions on contiguous data avoid per-element index computation. A targeted optimization pass reduced axis-reduction latency by ~56% (p50/p95/p99 deltas of ~90% on contiguous workloads).
• Broadcast index mapping: Output-to-source index mapping uses incremental odometer updates instead of full unravel/remap per element.
• 2x2 fast paths: Linear algebra has specialized 2x2 implementations that bypass general NxN overhead for the most common small matrix case.
• Horner's method: Polynomial evaluation and Stirling series use Horner's form for numerical stability and minimal multiplications.
• Ziggurat sampling: Normal and exponential random variates use the Ziggurat method (same as NumPy), which accepts ~97% of samples on the first try.

Performance budgets are enforced by the G7 gate, which measures p50/p95/p99 latencies for ufunc and reduction sentinel workloads and rejects regressions.

What works and what doesn't:

• Not a full Python package. FrankenNumPy is still primarily a Rust library. There is no pip install frankennumpy packaging story today, but the workspace now includes an early fnp-python extension module exposing PyNditer, frompyfunc, vectorize, digitize, searchsorted, interp, where, flatnonzero, argwhere, count_nonzero, isposinf, isneginf, signbit, isnan, isinf, isfinite, spacing, sign, floor, ceil, degrees, radians, sinc, copysign, nextafter, hypot, ldexp, logaddexp, logaddexp2, frexp, modf, nan_to_num, take, take_along_axis, compress, extract, select, and choose.
• No BLAS/LAPACK backend. Linear algebra uses pure-Rust implementations (Householder QR, Golub-Kahan SVD, implicit shifted QR for eigenvalues). Competitive with BLAS for small matrices; slower for large ones. Future BLAS linkage is planned.
• Complex elementwise arithmetic uses interleaved storage. Complex64/Complex128 dtypes store real/imaginary parts as interleaved floats with a trailing dimension of 2. Elementwise multiply and divide apply true complex arithmetic (a+bi)(c+di) = (ac-bd)+(ad+bc)i, but the interleaved representation adds overhead compared to native complex types.
• multivariate_normal uses Cholesky. NumPy defaults to SVD. Adding SVD would require fnp-linalg as a dependency of fnp-random (currently zero-dependency).
• multivariate_hypergeometric uses sequential draws. NumPy uses the random_mvhg_marginals algorithm.
• frompyfunc now has package-level Python exposure, but only for that surface. The Rust crate still supports closure-backed numeric frompyfunc, object-value outputs through frompyfunc_object, source-evaluated Python callables through frompyfunc_python, and imported Python callable handles through frompyfunc_python_import, while fnp-python now exposes a frompyfunc wrapper that accepts live Python callable objects and returns object-array results. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• vectorize now has package-level Python exposure, but only for that surface. fnp-python now exposes a vectorize wrapper that accepts live Python callable objects, handles broadcasting across positional arguments, supports positional excluded indices, and returns NumPy arrays or tuples of arrays with NumPy-inferred output dtypes. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• digitize now has package-level Python exposure, but only for that surface. fnp-python now exposes a digitize wrapper backed by the Rust UFuncArray::digitize_right implementation, with a reusable NumPy bridge for bool/int/uint/float arrays and regression coverage for decreasing bins, right=True, and large uint64 values. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• searchsorted now has package-level Python exposure, but only for that surface. fnp-python now exposes a searchsorted wrapper backed by the Rust UFuncArray::searchsorted implementation, including side selection and optional sorter bridging, with NumPy parity coverage for duplicate handling, probe-shape preservation, and unsorted inputs plus sorter permutations. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• interp now has package-level Python exposure, but only for that surface. fnp-python now exposes an interp wrapper backed by the Rust UFuncArray::interp_lr implementation, including default and explicit left/right fill behavior, multidimensional x shape preservation, and integer-input bridge coverage. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• where now has package-level Python exposure, but only for that surface. fnp-python now exposes both np.where(condition, x, y) and single-argument np.where(condition) behavior, backed by Rust where_select and where_nonzero, including broadcast parity, tuple-of-index-array results, and large uint64 branch preservation. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• flatnonzero and argwhere now have package-level Python exposure, but only for that surface. fnp-python now exposes flatnonzero and argwhere wrappers backed by the Rust UFuncArray::flatnonzero and UFuncArray::argwhere implementations, including multidimensional coordinate parity and the 0-D scalar shape edge cases NumPy uses for nonzero scalars. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• count_nonzero now has package-level Python exposure, but only for that surface. fnp-python now exposes a count_nonzero wrapper backed by the Rust UFuncArray::count_nonzero and count_nonzero_axes implementations, including scalar reduction, integer axis, axis-tuple handling, keepdims, and NumPy's empty-axis-tuple mask behavior. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• isposinf and isneginf now have package-level Python exposure, but only for that surface. fnp-python now exposes isposinf and isneginf wrappers backed by the Rust core helpers of the same names, including NumPy-oracle parity for scalar and multidimensional inputs containing finite values, NaNs, and signed infinities. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• signbit now has package-level Python exposure, but only for that surface. fnp-python now exposes a signbit wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for signed zero, NaNs, infinities, and multidimensional inputs. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• isnan, isinf, and isfinite now have package-level Python exposure, but only for that surface. fnp-python now exposes wrappers for the Rust unary predicates of the same names, including NumPy-oracle parity for NaNs, signed infinities, finite values, signed zero, and multidimensional inputs. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• spacing now has package-level Python exposure, but only for that surface. fnp-python now exposes a spacing wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for ordinary finite values, signed zero, NaNs, and infinities. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• sinc now has package-level Python exposure, but only for that surface. fnp-python now exposes a sinc wrapper backed by the Rust UFuncArray::sinc implementation, including NumPy-oracle parity for fractional float inputs, integer inputs, and multidimensional shape preservation. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• copysign now has package-level Python exposure, but only for that surface. fnp-python now exposes a copysign wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for broadcasted inputs, signed-zero sign transfer, and multidimensional shape preservation. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• nextafter now has package-level Python exposure, but only for that surface. fnp-python now exposes a nextafter wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for directed adjacent-float stepping, signed-zero transitions, and broadcasted multidimensional inputs. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• hypot now has package-level Python exposure, but only for that surface. fnp-python now exposes a hypot wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for elementwise Euclidean norms and broadcasted multidimensional inputs. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• ldexp now has package-level Python exposure, but only for that surface. fnp-python now exposes an ldexp wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for integer exponent scaling, signed-zero preservation, and broadcasted multidimensional inputs. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• logaddexp now has package-level Python exposure, but only for that surface. fnp-python now exposes a logaddexp wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for stable exponential summation, -inf/inf edge cases, and broadcasted multidimensional inputs. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• logaddexp2 now has package-level Python exposure, but only for that surface. fnp-python now exposes a logaddexp2 wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for base-2 exponential summation, -inf/inf edge cases, and broadcasted multidimensional inputs. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• frexp now has package-level Python exposure, but only for that surface. fnp-python now exposes a frexp wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for mantissa/exponent tuple outputs, integer exponent dtype parity, special values, and multidimensional shape preservation. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• modf now has package-level Python exposure, but only for that surface. fnp-python now exposes a modf wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for fractional/integral tuple outputs, signed zero behavior, infinities/NaNs, and multidimensional shape preservation. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• nan_to_num now has package-level Python exposure, but only for that surface. fnp-python now exposes a nan_to_num wrapper backed by the Rust core helper of the same name, including NumPy-oracle parity for default replacements, custom scalar nan/posinf/neginf substitutions, and integer-input pass-through behavior. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• sign now has package-level Python exposure, but only for that surface. fnp-python now exposes a sign wrapper backed by the Rust unary helper of the same name, including NumPy-oracle parity for signed-zero normalization, NaN propagation, integer inputs, and multidimensional shape preservation. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• floor and ceil now have package-level Python exposure, but only for that surface. fnp-python now exposes floor and ceil wrappers backed by the Rust unary helpers of the same names, including NumPy-oracle parity for floating-point special values, signed-zero behavior, integer-input dtype preservation, and multidimensional shape preservation. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• degrees and radians now have package-level Python exposure, but only for that surface. fnp-python now exposes degrees and radians wrappers backed by the Rust unary angle-conversion helpers, including NumPy-oracle parity for float and integer inputs, multidimensional shape preservation, and output-dtype matching. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• take now has package-level Python exposure, but only for that surface. fnp-python now exposes a take wrapper backed by the Rust UFuncArray::take implementation, including flat indexing, axis-aware gathers, scalar-index scalar returns, multidimensional index-shape preservation, negative-index handling, and large uint64 value preservation. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• take_along_axis now has package-level Python exposure, but only for that surface. fnp-python now exposes a take_along_axis wrapper backed by the Rust UFuncArray::take_along_axis implementation, including default axis=-1 behavior, axis=None flatten semantics, and large uint64 gather parity. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• compress and extract now have package-level Python exposure, but only for that surface. fnp-python now exposes compress and extract wrappers backed by the Rust UFuncArray::compress and UFuncArray::extract implementations, including NumPy-style short-condition truncation for compress, flat mask extraction for extract, and large uint64 selection parity. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• select now has package-level Python exposure, but only for that surface. fnp-python now exposes a select wrapper built from Rust where_select chaining, including first-match precedence, broadcast parity across condition and choice lists, integer-default behavior, and large uint64 default preservation. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• choose now has package-level Python exposure, but only for that surface. fnp-python now exposes a choose wrapper backed by the Rust UFuncArray::choose implementation, including large uint64 selection parity and full mode='raise', mode='wrap', and mode='clip' parity. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• nditer now has package-level Python exposure, but only for that surface. fnp_iter::Nditer has nditer_python / nditer_python_with_interpreter for NumPy-backed state and chunk bridging, and fnp-python now exposes a PyNditer class. Broader Python packaging and FFI coverage for the rest of FrankenNumPy remain incomplete.
• Single-threaded. All operations are single-threaded. The asupersync async runtime integration is optional and used only for conformance pipeline orchestration, not for parallel array computation.
• f64 internal representation. UFuncArray stores numeric values as Vec<f64> internally for arithmetic. For i64/u64 values > 2^53, an IntegerSidecar preserves exact integer values through storage round-trips (from_storage / to_storage). Arithmetic on large integers still uses f64 approximation.

Is this a drop-in replacement for NumPy? Not yet. It reimplements NumPy's semantics in Rust for correctness verification and eventual use as a Rust-native array library. A Python FFI bridge is planned but not implemented.

How do you verify parity with NumPy? Oracle tests: we run the same operations with the same inputs in both NumPy and FrankenNumPy, comparing outputs to floating-point tolerance. For RNG, the comparison is bit-exact.

Why Rust nightly? We use Rust Edition 2024 features. The toolchain is pinned to a specific nightly date for reproducibility.

Why zero unsafe code? Memory safety is a core value. All 9 crates declare #![forbid(unsafe_code)]. The 92,000+ lines of Rust contain zero unsafe blocks.

How fast is it? Performance is not the primary goal in this phase. Correctness and parity come first. That said, the release profile uses opt-level = 3, LTO, and single codegen unit.

Can I use just the RNG crate? Yes. fnp-random has zero dependencies and produces bit-exact NumPy-compatible random sequences from a given seed.

Please don't take this the wrong way, but I do not accept outside contributions for any of my projects. I simply don't have the mental bandwidth to review anything, and it's my name on the thing, so I'm responsible for any problems it causes; thus, the risk-reward is highly asymmetric from my perspective. I'd also have to worry about other "stakeholders," which seems unwise for tools I mostly make for myself for free. Feel free to submit issues, and even PRs if you want to illustrate a proposed fix, but know I won't merge them directly. Instead, I'll have Claude or Codex review submissions via gh and independently decide whether and how to address them. Bug reports in particular are welcome. Sorry if this offends, but I want to avoid wasted time and hurt feelings. I understand this isn't in sync with the prevailing open-source ethos that seeks community contributions, but it's the only way I can move at this velocity and keep my sanity.

MIT License (with OpenAI/Anthropic Rider). See LICENSE.