Rustymath drops precision that numerical code doesn't need, for speed and a smaller implementation. Rustymath doesn't aim for last-bit agreement with other libraries, or for reproducible results. The last bits of a result can differ between thread counts, platforms, builds and runs. Only spend speed, memory or code to make them agree if there's a real common practical use case that needs it. Tests never compare results across thread counts.
Keep the accuracy contracts that callers use. Check a change that alters arithmetic against those contracts, not against bit-identical results. Don't add sign canonicalization or compensated arithmetic without an agreed need. Each shortcut that drops precision is a hypothesis that real uses don't need that precision. When a real use needs digits that a shortcut loses, return that operation to the conservative method, or create an optional parameter.
cargo test
RUSTYMATH_COSTS=1 cargo test --test speed
RUSTYMATH_COSTS=1 pytest -q tests/test_speed.py
cargo develop && pytest -q
nbdev-test
nbdev-previewRun plain cargo test while you work. It runs the Rust integration tests, speed test and doctests. cargo develop builds the same test profile and installs the Python extension into the editable package. Rebuild before testing Python changes that depend on Rust. An already-running Python process retains its loaded extension.
When another agent's unfinished edits stop the tree compiling, run cargo wait with the Bash tool's run_in_background, then build after it exits. cargo wait comes from fastws-cli, as cargo develop does. It runs cargo check on every target each time a source file or manifest of a local package changes, and exits once the tree compiles, or after 30 minutes with the remaining errors. Use it rather than writing a wait loop.
On macOS, dev builds compile rustymath's own linear algebra beside Accelerate, through the default native feature. RUSTYMATH_NATIVE=1 cargo test runs the tests on rustymath's own kernels, which Linux and wasm use.
Build nothing else without explicit discussion: no scratch crates, alternate target directories, --release or -p builds. The exception for inspecting x86-64 machine code is described in basedpl's DEV.md, under “X86-64 machine code”. Use that procedure rather than --emit asm, which changes the codegen-unit arrangement you are inspecting. Packaging checks use ship-rs-build; pass --profile dist when checking a distributed artifact.
The root package is the Python-free crate published to crates.io. py/ is the unpublished PyO3 extension. pyproject.toml points maturin at its manifest. Editable installation through uv sync or pip install -e uses fastws.build_backend and compiles no Rust. Run cargo develop after installation. Don't use maturin develop: its build configuration differs from the shared local build.
The test profile inherits release: optimization level 2, incremental compilation, LTO off, no debug information, overflow checks or debug assertions. The configs explicitly enable LLVM's SLP vectorizer. On x86-64 they also target x86-64-v3 for AVX2 and FMA. Without those target features, some conversions and rounding operations run one item at a time. dist disables incremental compilation, enables Cargo's default within-crate LTO and strips symbols.
Keep these profiles aligned with basedpl's. Cargo uses the consuming workspace's profiles for dependencies. Concrete function forms compile in rustymath, but their build still uses that workspace's profile and target flags. LTO-off incremental builds limit how much code an edit recompiles; debug line tables and inlining across units increase that work.
tests/speed.txt holds two budgets for each function, in copies of the same number of items: one for macOS on Apple silicon and one for Linux on CI's x86 runner. The speed test checks the one that timing::NATIVE_COLUMN picks. A copy is the time to write that many items of the same type, each negated, into a second buffer. The speed test times each copy immediately before its operation. The standard-library comparisons are informational. A budget failure often identifies a loop that stopped vectorizing. RUSTYMATH_COSTS prints every cost and the comparison table in src/lib.md's "Accuracy and speed".
tests/test_speed.py times rustymath's Python API the same way, in the fresh process that pytest runs. Its cases are the Python expressions in tests/speed_python.txt, each with two budgets in copies of 65,536 items of its section's type, as in tests/speed.txt. RUSTYMATH_COSTS prints every case's cost in the speed test's format. Another agent's cargo develop can replace the installed extension at any time. To time a change against the build before it, save the Python package before your first cargo develop, then time the saved copy and your build with basedpl.bpltests.pytest_costs:
shutil.copytree(rms/'python/rustymath', '/tmp/<prefix>_rm/rustymath', ignore=shutil.ignore_patterns('__pycache__'))
pytest_costs(rms, saved='/tmp/<prefix>_rm'), pytest_costs(rms)saved puts the copy first on PYTHONPATH. set_budgets takes pytest_costs(rms) and path=rms/'tests/speed_python.txt'.
The manually triggered Speed workflow measures the Ubicloud runner used by Test and restores its cache. The Test job also sets RUSTYMATH_COSTS on every push. logged_costs reads its log as it reads a Speed job's. After a run, review it and decide which budgets should change, following basedpl's DEV.md on Speed runs. Then use basedpl.bpltests to set those budgets under the rule in tests/speed.txt:
set_budgets(native=[cargo_costs(rms), logged_costs(log)], path=rms/'tests/speed.txt', cases=cases)Here rms is this repo's Path, log is the job log and cases lists the codes of the cases you decided to change. cargo_costs(rms) on this Mac sets the macOS budgets. logged_costs(log) sets the Linux budgets. Budget measurements use parallel::set_threads(1). Measure multicore scaling separately.
On macOS, a native timing measures Accelerate wherever the Linear algebra section says Accelerate replaces rustymath's own code. To time rustymath's own versions natively on the Mac, pass native=True to basedpl's time_cases, speed_costs, cargo_costs or pytest_costs. That sets RUSTYMATH_NATIVE=1 for the process they time, which dev builds honour. basedpl's wasm build, timed with speed_costs(wasm=DIR), also runs rustymath's own versions. CI's Speed run times them natively on x86-64.
Both speed tests and basedpl's timing tools take the machine's timing lock exclusively. Rust code takes it through rustymath::timing, whose docs give its waits and messages, and Python code through rustymath.timing(), a with block around the same hold. Both speed tests time each case over its budget again after each pause in rustymath::timing::RETRY_PAUSES, and fail only for a case that's over every time. fastws.cargo.timing_lock takes the same lock in Python code that doesn't import rustymath. basedpl's DEV.md lists the builds and tests that take it shared, and the rules for agents.
Most of rustymath is generic, so most of its machine code compiles in its callers, chiefly basedpl. Measure a change's effect on build time and code size there, with the tools that basedpl's DEV.md describes under "Codegen and rebuilds":
- Code size: save a copy of basedpl's installed
bplbefore your firstcargo develop, as for timings. After the change,BASEDPL_SIZES_DIFF=/tmp/<prefix>_bpl:<new bpl> cargo test --test units machine_code_growth -- --nocapture, run in basedpl, compares the two builds' machine code by generic function, rustymath's included.BASEDPL_SIZES=patternlists the instances that match a pattern, with their sizes. - Rebuild time: wrap one value in the edited function in
std::hint::black_box, which changes the generated code but not the behaviour. Time basedpl'scargo build --lib --bins --profile test --timings, then undo the edit and time the rebuild again. The report in basedpl'starget/cargo-timings/splits each crate's time into frontend and codegen. - rustymath's own build:
cargo build --lib --profile test --timingsin this repo times the code that compiles in rustymath itself, such as concrete function forms.
A generic function compiles once for each combination of types that callers give it, a conversion closure's type included. A loop that is generic over a conversion as well as its item types multiplies its instances. gather's walks show the cost. When their fixed-width row loops took the source type, the result type and the conversion, one basedpl codegen unit held 603 instances of them, about 940 KB of bpl, and an edit to them rebuilt basedpl in 16.8 s. With those loops compiled once for each size of item, and conversions in a separate pass (see "Search and gathering"), basedpl's machine code fell by 1.1 MB on 20 October 2026, and the same kind of edit rebuilds basedpl in about 7 s. The fixed-width loops that convert Booleans and integers into f64 as they write (Item::cells_from) add about 17 KB, and leave that rebuild at about 7 s.
| Files | Responsibility |
|---|---|
lib.rs, float.rs |
Public functions and objects; Float, Native, fused/unfused arithmetic and sealed implementation traits. |
exp.rs, ln.rs, trig.rs, arc.rs, hyperbolic.rs, pow.rs, parts.rs, sqrts.rs |
Elementary scalar functions and specialized forms. |
table.rs, dd.rs |
Table-based f64 logarithms/exponentials and double-word arithmetic. |
special.rs, lambert.rs, complex.rs |
Special functions and complex forms of the same objects. |
vforce.rs |
Sines and cosines over slices through vForce, as SinCos on macOS. |
slice.rs, parallel.rs |
Scalar/slice object interfaces, output storage and elementwise scheduling. |
reduce.rs, stats.rs |
Folds, scans, axis traversal and statistical summaries. |
create.rs, structure.rs |
Typed construction, sequences and structural selection plans. |
arith.rs, integer.rs, convert.rs |
Wrapped/checked arithmetic, integer functions and vectorizable conversions. |
distribution.rs, random.rs |
Checked distributions, function objects, sampling, rolls and deals. |
number_theory.rs |
Primality, factorization, prime enumeration/counting and factorial. |
field.rs, exact.rs, linalg/exact.rs, polynomial.rs |
Generic field algorithms, rational functions, exact matrices and polynomials. |
search.rs, tolerant.rs, gather.rs |
Exact/tolerant search and classification, stable sorting, partial selection, comparisons and copying through conversions. |
error.rs, layout.rs, array.rs |
Errors that the functions on arrays return, with their kinds and messages; array shapes, axis keys and names, with their size limits; and the array type, a layout paired with the buffer that fills it. |
optimize.rs |
Minimization by five methods behind one Objective trait, and solving behind one Mapping trait: Brent's root finder and Levenberg–Marquardt for nonlinear equations, and CG, MINRES and GMRES for linear ones. |
linalg.rs, linalg/scalar.rs, linalg/factor/ |
Packed products, scalar operations and decompositions. |
signal.rs |
Correlations and convolutions, which slide a kernel along leading axes through linalg's mapped products, and overlap-add. |
py/src/, python/rustymath/ |
The Python extension: the Array class, the converter of call arguments and each module's functions, with Python modules that name them. |
buffer/ (rustymath-buffer) |
The Python buffer protocol for compact buffers, and the Python exception for each error kind. Basedpl's extension shares it. |
Free functions are generic over Float: f16 computes in f32 and rounds its result, while f32 and f64 use their own working type. Native names f32 and f64. functions! declares the free wrappers at the crate root and in fused/unfused; an as entry also declares the function object. FUSED selects hardware FMA where available. Don't introduce software FMA into targets without it.
Objects implement Unary<T>, Binary<T> or Ternary<T> at concrete item types. They can hold settings. one computes a scalar, checked also reports failure, and many appends results to a vector and returns whether every item succeeded. Arguments are Arg::One or Arg::Many. Generic code calling an object states its bound, such as Exp: Unary<F, Out = F>; free functions need no object bound.
forms! writes the sealed Forms<E, N> implementations. at real writes f32/f64 forms and an f16 adapter through f32. Explicit type lists cover complex, integer and Boolean forms. => T names an output type different from the input. Attributes before the object, such as #[inline], go on its scalar forms, scalar and checked_scalar. The integer arithmetic forms use #[inline], so that basedpl's loops over items inline them, because LTO is off. The item body returns one of three outcomes:
- A plain float/complex value fails on NaN from arguments that aren't NaN. Plain integer and Boolean values have no NaN failure.
Optionfails onNone.Totalnever fails, including where NaN is the intended result.
Arithmetic objects deliberately accept IEEE NaNs. An integer result that doesn't fit its type fails: the integer forms of Add, Subtract, Multiply and Abs return Option, as the other integer functions do. A : total entry in functions! wraps the body in Total. Test an object's actual verdict rather than applying the ordinary NaN rule to every object.
Each arity trait's TOTAL says whether an object never fails. It is true for an object whose body gives a Total output, or a plain integer or Boolean output. forms! derives it from the body's type. Basedpl's fusion runs only total objects, because a deferred call mustn't fail.
arities! gives the public traits one implementation through Forms. Keep objects concrete: their forms compile in rustymath instead of being instantiated in every caller. Forms that differ from the free functions, including Sin, Cos, Tan, Pow and Root, live beside their implementations.
kernels! selects an object's function pointers at each dtype. Fold, Scan, RightFold and LeftScan share AccumulationKernel and .run, with the probe alone choosing the operation. Its typed kernel reads an optional right tail or scalar start in the total or computing dtype.
Forms::slice writes into provided slots. many and pairs reserve spare capacity through appended, without zeroing it. written views initialized slots for later passes. Keep that unsafe storage handling in these helpers. many_into and Slots::pairs run the same loops into a parallel::Slots, which extended and Slots::append hand out. A caller such as basedpl's fused chains then writes its last step straight into its result. Each writer fills every slot it gets, and a writer that doesn't panics before the vector's length changes.
each runs the ordinary item loop. It checks outputs first and inspects arguments only for NaN outputs. two_pass computes a branch-free fast path, then recomputes its rare misses. split separates complex parts into real buffers, calls real slice forms and redoes exceptional items with the scalar function. CHUNK governs the cache-sized pieces used by these forms and their adapters.
Keep these performance rules when changing loops:
- Pass item bodies through
each!,two_pass!orsplit!. They mark closures#[inline(always)]; leaving a large body out of line can stop vectorization. - Write separate closures for functions selected at runtime. A shared closure can leave the selector inside every iteration.
trig.rsshows the pattern. - Keep
pairsat#[inline(never)]. Inlining a function whose loop runs behindparallel::split's&dynonly repeats its setup at every call site. - Check code in basedpl as well as in this crate when changing a generic function that it instantiates. The small test crate is not the same compilation context.
- Fill a large per-call buffer through
appended's spare capacity, or withVec::with_capacityandresizeaslinalg'sfilleddoes. On macOS,vec![0; n]allocates through calloc. Calloc'smadvisecosts more than writing the items. Its cost grows with the number of threads that touch the buffer. - Write each hot loop as a plain loop over slices, with each value's representation chosen before the loop, or as a call to one of our own
#[inline]helpers, nested helpers included. Local builds of thetestprofile are incremental, with LTO off, in up to 256 codegen units. LLVM inlines nothing across units. A closure or an#[inline]function gets a copy in every unit that calls it. In an incremental build, a generic function without#[inline]gets one copy per crate, in one unit. std's iterator plumbing, such asMap::foldandVec'sspec_extend, has no#[inline]. When its copy compiles outside the caller's unit, the loop runs as calls, sometimes one for each item:filterinstructure::scaledcost about a nanosecond per item. Which unit holds a copy changes with unrelated edits. Thedistprofile's thin-local LTO inlines across units. CI builds thetestprofile withCARGO_INCREMENTAL=0. Its builds split code into units differently from local builds. The out-of-line calls above were measured in local incremental builds. A loop throughSteps::offsetsalso runsget's variant match and reloads its fields for every item, because the compiler can't rule out that the output buffer aliases theSteps. In basedpl,BASEDPL_CALLS=pattern cargo test --test units -- --nocapturelists each call from a matching function into another unit, in basedpl's build and in the rustymath rlib it links.- Keep
Steps::offsetsaRangemapped throughget. That iterator isTrustedLen, andVec::extendandzipuse their fast loops only forTrustedLeniterators, which stable Rust doesn't let a new type implement. A custom iterator type in its place took a write such asx.[i]←wfrom about 9 copies to 21. A hot loop over a table loops over its slice directly. - Keep 64-bit integer multiplies out of loops that should vectorize. NEON has no 64-bit vector multiply. Multiplying each table entry by a stride as it was read stopped
Steps::combined's loop vectorizing, and a write such asx.[i]←wcost about a copy more. A table holds finished offsets, whichTable::scaledmultiplies once. - Loop over an axis's positions through
Steps::for_each, which matches the kind of step once, before its loop. LLVM can keepgetout of line despite its#[inline], even within the caller's unit. A loop over rows throughgetoroffsetsthen makes a call for each row. Beforefor_each, each short row cost about 3.5 ns inscatter_intoandwalk_into:[∞;0 1]⌷mon 65,536 rows of 4 floats cost 15 copies, against 2.2 with it. - Don't build a std
zipfor each short row. Its constructor,ZipImpl::new, has no#[inline], so each row pays a call to its one copy. An index loop over two slices of the same known length needs no bounds checks, asrowingather.rsshows.
- Keep
Elementary kernels aim for branch-free vectorizable loops. Incomplete gamma/beta and their inverses also use data-dependent series, continued fractions and iterations, with precision-appropriate stopping. Table gathers help AVX2 but are unavailable in NEON/wasm. Short dependency chains and vectorizable conversions matter as much as instruction count.
Pow/Root use table-based f64 arithmetic and narrow-width vector paths. Their scalar object forms have separate fast paths; f32 scalar calls use tables that the vector paths avoid. A whole exponent from 2 to 8 squares and multiplies, in the scalar form and in a slice form with a constant exponent. A degree of 2 or 3 takes sqrt or cbrt. A table-free f64 Pow measured about 30% slower: with ln and exp in one loop, their polynomial constants don't fit in registers. The loop then rebuilds them for every item. cbrt computes at each width's own precision. Trig slice forms redo large arguments after their fast reduction. f64 inverse-error slice forms group inputs by approximation region. Erfc chooses a single formula for blocks wholly on one side of its split. Normal/log-normal CDF and quantile objects reuse those optimized forms through composed.
parallel::ranges schedules independent ranges. split applies it to disjoint output slices; many, pairs and extended use that interface. A long call times its first unit, then cuts the rest into pieces of at least PIECE_TIME of work, at most four for each thread. pieces runs fewer pieces than threads on a pool of the largest power of two up to their number: memory-bound calls of a million items ran faster on small pools than on all 18 threads of an M5 Max. Chunked forms keep their temporary buffers in cache. Work and result initialization finish for every piece even if another piece reports failure. parallel::extended lets a caller such as basedpl's fusion fill a vector's spare capacity in pieces, through Slots.
parallel::stack_with runs a function on each cell of a stack of matrices, for basedpl's f⍤2 and the binding's matrix functions. It times its first cells on the calling thread, then splits the rest across threads in rounds of about ROUND_TIME for each thread, with the caller's check between rounds and one scratch value, such as an LuWork, for each piece. A cell on another thread doesn't split again.
The default pool follows RAYON_NUM_THREADS, otherwise Rayon's available-core default. parallel::set_threads(n) changes calls made on the current thread and returns the previous setting. Use 0 to restore the default or 1 for serial calls. Calls already on a Rayon worker do not split again. Wasm runs serially. Pools are created on demand and rebuilt in a child after fork.
This scheduler covers elementwise forms, folds and matrix products. Scans don't split. reduce::along folds a lane, or a group of inner lanes, that fits in one piece of reduce::fold_pieces whole, and runs of them split across threads. Larger arrays fold in one piece per thread, each of at least 2^18 items and 256 rows. A lane that crosses pieces merges its parts' folds. Splitting a single matrix's columns across threads made +⌿ twice as slow, because each thread read part of every row. A product splits its rows into bands of about BAND_WORK multiply-adds, each with its own packed block of a. It splits once for each panel of b of up to PANEL items, or SERIAL_PANEL bytes when it runs on one thread, across at most √(work / THREAD_WORK) threads. Each thread it wakes costs about as much as THREAD_WORK multiply-adds. Decompositions split in the products of their blocked updates. The Hermitian and bidiagonal reductions also cut each pass's rows into one piece per thread, each of at least ROW_PIECE rows, through eigen::summed_rows. Divide and conquer splits each merge's secular roots, corrected z and vector rows. Panels and the recursion's two halves run on one thread. Callers that submit only small converted chunks also limit the work available to it. The shared non-generic scheduler keeps Rayon's machinery out of each kernel instantiation. Item types and forms are Send + Sync through the sealed traits.
Reduction defines the unit, scalar step, checked step, contiguous fold, merge and scan for Sum, Product, Max and Min. Integer sums/products wrap and report whether they fit. f16 sums and products have an f32 Total on every fold path, column folds and fold_pieces merges included. The caller rounds it once. f16 extrema total in f16. Float sums/products can regroup algebraic operations. Beside NaN, float extrema follow each target's instructions. AArch64's fmax and fmin give NaN. Every other step passes the running result as maximum's or minimum's y, which a NaN item leaves unchanged.
slice::running chooses the partial-total count for the target and item size. Explicit vector extrema are used where they outperform ordinary loops. On AArch64 and in wasm, f32 running sums scan four-lane vectors through sum_scan, which makes +\ about 1.3 times faster than grouped's groups of four items. Wasm uses max/min, which keep NaN. They take about eight instructions each on x86-64 hosts, but run faster than pmax/pmin in a browser on an Arm host, such as a Mac. Integer extrema scans skip blocks that cannot change their running result.
Logical represents any Boolean dyad as a truth table and scans by composing its two-state maps. Xor/Xnor use that path too. Lanes of at most 64 consecutive items scan in one call of mapped, which restarts each lane inside its word. SumAs<C> counts Booleans or sums integers in the caller's chosen total type. The caller must choose a width that holds the totals.
Alternating implements right folds by subtraction or division. Float/complex pairs sum differences or multiply ratios. Integer subtraction keeps separate totals for even and odd positions. A floating quotient reports failure when its intermediate products or result leave the normal range; basedpl then folds item by item.
along, scan_along and fold_into traverse [outer, len, inner] arrays. Contiguous lanes of at least SHORT_LANE items use slice reductions, and shorter ones fold item by item. Along another axis, rows no wider than slice::running's count of running totals fold a block of whole rows at a time, into one running total for each item of the block. Each column's totals then merge. Wider rows fold in strips of that many columns, each down every row, where there are at least 4 rows and the totals are narrower than complex numbers. Other wide rows update a row of totals. fold_into accepts blocks from callers such as fusion. extremes and extremes_of compute both bounds in one pass. along_mapped folds an array that a mapping lays over storage, without building it. Where the mapping reads a block of storage, as Mapping::block finds, Block::fold gives the dimensions to fold the block along in its own order, and the tables that read the totals in the array's order. Other mappings are read in pieces.
sliding folds each window of width rows along the middle axis of an [outer, len, inner] array. Windows of up to DIRECT rows fold directly, one pass for each row. Wider windows take van Herk-Gil-Werman's method: each window folds one block's suffix with the next block's prefix, at a cost for each item that doesn't grow with the width. For lanes of single items, one const parameter compiles each loop over a row to scalar code and keeps the running fold in a register. A running sum, recomputed every 1,024 windows, ran float sums two to three times faster on the Mac. Its error grows with every item that passes through its run, not with the window's own items. A spike of 1e20 among ones then zeroes the windows after it until the run ends. Sums take van Herk-Gil-Werman too. sliding_picks gives the row of each window's first maximum or minimum, through the same code over pairs of an item and its row. view_sliding_picks runs it at a view's own type, for basedpl's ⌈ and ⌊ rules and the Python sliding_argmax and sliding_argmin. kernels!(Sliding of r) runs sliding at each dtype, for the Python sliding.
Complex forms attach to the same objects as real forms. Abs and Angle return real values; real Cis returns complex values. On macOS, vForce implements SinCos, in vforce.rs. Real Cis calls its vvcosisin, and the complex slice forms that split into real parts take their sines and cosines from one vvsincos call. Most expensive complex slice forms split into vectorized real operations. Abs, Angle and Sqrt keep direct loops. Exceptional magnitudes and branch-cut cases take scalar paths. Inverse circle/hyperbolic functions use Kahan's formulas; complex ln Gamma uses Stirling, recurrence and reflection with SciPy's branch convention.
Complex errors are measured against a magnitude-based scale, with function-specific scales for ln Gamma, powers, Gamma and binomial. Quarter-turn directions are deliberately exact. See complex.rs for that contract rather than applying a separate ulp bound to each component.
Distribution validates parameters and holds a Law. Pdf, Cdf and Quantile select the law outside their item loops. Formulas use f64 working arithmetic and narrow their outputs. Sampling, rolls and deals take a caller-owned RNG. Keep RNG ownership separate from the numerical algorithms.
Overflowing gives wrapped integer operations and fit flags without branching in the hot loops. The integer forms of Add, Subtract, Multiply, Abs and Pow give that pair as their outcome. forms! derives their WRAPS constant from it. Kernels::wraps reads that constant for the dtype an argument computes in. Every reduction wraps too. A result that doesn't fit fails, with the wrapped result in its slot. The Python API keeps that result, as NumPy does. BPL's exact integers widen instead. Width::float and DType::inexact take the widest float that an integer may become: f64 for NumPy's rule, or a session's narrower default float. DType::common gives an integer beside a float or complex dtype that dtype, as JAX does, where NumPy widens the float. convert::int/float use bit-based conversions where target instructions do not vectorize. Check NATIVE_64 and SATURATING before replacing them with casts. Divisor turns repeated division by one integer into multiplication and shifts.
integer provides checked gcd, lcm, binomial and power functions. residue returns an integer directly, including the dividend for a zero divisor. Objects also have narrow and unsigned forms. encode/decode handle mixed radix; polynomial decoding handles fields.
number_theory::Natural narrows BigUint through u128 to u64 before computing. The fixed-width tiers use Montgomery arithmetic; BigUint uses ordinary residues. Pass cancellation checks by value through narrowing calls. Passing &mut check recursively creates new closure types and unbounded monomorphization. Primes is a segmented wheel sieve.
exact operates on rational numerators/denominators. linalg::exact supplies row reduction, solve, pseudoinverse and determinant for fields; least-norm solves use a rank factorization. Permanent uses Ryser's Gray-code formula. Basedpl uses rational matrix solves/determinants and native floating decompositions for approximate inputs. Generic permanent and polynomial helpers also run on f64 and complex floats. Long-running functions expose _with cancellation variants.
polynomial provides generic value, factored, from_roots and derivative, vectorized Horner evaluation, companion-matrix roots and least-squares fitting. Keep Horner's blocks fixed-size: variable slices can introduce copies and iterator setup in downstream codegen. Measure changes in basedpl too.
search owns exact search, classification, counts, interval index, stable counting/radix sorting and grouped positions. Positions are u32; slices must fit that range. Whole-slice algorithms stay out of line to avoid duplicating them in callers. Each buffer of a Grouped holds whole classes. Two classes take a buffer each, which into_classes hands over without a copy. More classes share one buffer, because a buffer for each class costs an allocation that callers reading only the positions don't need.
Sorting, partitioning, rank selection and quantiles share exact ordering keys. Float ordering equates signed zeros and puts NaNs last ascending. Partial selection leaves unrequested ranks unordered. Median and quantiles select their required ranks without a full sort.
tolerant owns relative comparisons, rounding and tolerant searches. Search can index whole-number keys, but validates returned matches unless the haystack is exactly whole and falls back to buckets when needed. Classification intentionally permits tolerant chaining and whole-key groups; it is not the same relation as matching every member against a representative. Near and runs implement the non-whole paths. Tolerance-sensitive tests need independent inside/outside cases and exact controls.
gather copies through a caller's conversion for cycle, repeat, compress, indexed take and scatter. Mask compression skips empty blocks and copies full blocks. In a mixed block, compress uses a branch-free loop that reads every item, which suits items loaded from a slice. positions, and the split into two classes in search, write only the selected items through the 8-bit SET table, which suits computed positions. The table loop measured 13% slower for items loaded from a slice. walk copies the items that one Steps per axis reaches: strides, clipped ranges, rotations and offset tables. Build a clipped range with Steps::clipped and a rotation with Steps::rotated. Each gives a Stride where nothing clips or wraps. Every Clipped entry reads at least one fill. Every Rotated entry wraps. Positions outside the source take the caller's fill. It merges neighbouring runs, copies long runs whole, transposes through linalg::transpose_into's tiles and writes short windows through windows. walk and scatter_walk take a Convert. With Same, they move items straight between source and result. Any other conversion runs as a separate pass over blocks of about CHUNK items: a walk gathers a block of source items, then converts it, and a scatter converts a block of items, then writes it. Their walks therefore compile once for each item type (Item), and their fixed-width loops, windows and windows_into, once for each Unit: an unsigned integer, or an array of two, with an item type's size and alignment. Two slice casts, units and written_as, move items as units. Both casts are sound by the promises of two unsafe traits. An Item type has no padding and has its Unit's size and alignment, which units! checks at compile time. Every bit pattern of a Unit is a value of it, so items read as units are valid units. written_as writes only units that units read from items of the same type, so each slot it fills holds a valid item. That matters for bool and char, which have bit patterns that aren't values. A conversion's own pass costs about 0.25 ns per item on this Mac, measured on 20 October 2026. A conversion inside the walk cost nothing extra, because the walk waits on memory. So scatter_walk writes short cells along a stride from Booleans and integers into f64 through fixed-width loops that convert as they write (Item::cells_from). units! gives those five source types alone such loops, which add about 17 KB of basedpl's machine code. Steps::through composes two walks into one only where the result stays regular, or where either walk is already an offset table. Elsewhere callers walk the inner copy into a buffer first. That costs one pass, where building and reading a table of every offset costs two. Basedpl reuses these paths for structural operations. Basedpl's reverse mode transposes structural copies, without scattering, through the reads that Steps::transposed gives. Measure search/gather changes there as well.
agreement pairs two layouts by BPL's broadcasting rule and runs kernels through the pairing. block and blocks split the frame into ranges through parallel::extended, so broadcasts run in parallel as contiguous calls do, and their leaf writes each range's results into its slots. Within a range, blocks hands contiguous arguments, and broadcast blocks of at least SHORT_REPEATED or SHORT_CYCLED items, to the leaf whole, cutting a block where the range ends. It expands other runs through expanded, CHUNK items at a time, or the whole rows that Mapping::line_block asks for where an argument's innermost axis skips items, as a transpose's does. map1 and map2 run a function object's many_into through block and blocks where the arguments' items already have the computing type. Otherwise each range reads runs of the same length through read_mapped, which expands and converts them. The leaf runs on every chunk even after a failure, because extended checks that every slot is written. Agreement is generic over its layout type, which basedpl uses to carry its renderer. blocks, map1, map2 and expanded are #[inline] and generic, and compile in each caller's kernel units. An edit to blocks rebuilds basedpl in about 12 s.
optimize has four fronts. minimize takes an Objective, and solve, solve_linear and solve_affine take a Mapping. An objective gives a value, a gradient where has_gradient says it has one, and a batch of values in one call through values where batch says it evaluates batches. Without batch, Search calls value on each point of a batch. A mapping writes its values into a slice, gives vector–Jacobian products where has_vjp says it can, and gives its Jacobian through jacobian. By default jacobian builds the rows from vector–Jacobian products when the mapping has no more values than variables, and otherwise leaves Levenberg–Marquardt to forward differences. Where has_preconditioner says it has one, its precondition gives solve_linear a preconditioner. A solve without one uses each residual as it is, with no copy. Each optional method has a flag beside it and an unreachable default. An override then holds only the case it adds. A closure that returns a Result is an objective or a mapping, and Differentiable pairs a value closure with its derivative. The function's error type is the front's error. A caller's interrupt or failed call therefore stops the search without a separate check, which matters because basedpl's and Python's functions can fail on any evaluation.
Every method of minimize and solve runs on a private Search. It counts evaluations against the limit, moves each point inside the bounds, in the method's own buffer, before evaluating it, reads NaN as +∞ and keeps the best point. A method returns whether its stopping test passed. Running out of evaluations unwinds through Halt::Spent, and passed then gives converged false for the best point. solve wraps its mapping in Norm, an objective whose value is the norm of the residuals, and which keeps the residuals for the method to read. solve and solve_affine take first, the mapping's values at the start, which give the number of residuals. The caller computes them, because only an evaluation tells a binding how many there are. Search::seed counts that evaluation and records it, and the methods start from it, so it isn't repeated. Bounds, the limit, NaN and the result therefore live in one place for minimizing and solving, and a new method should run on a Search too.
One line search serves Brent's method and Powell's. least_along brackets a minimum with golden-ratio steps, which stop at the bounds' limits along the line, and brent searches the bracket. Brent's method is one such line search along its variable. Every Nelder–Mead move puts the worst vertex at a point on the line from the centroid of the others through it, with Gao and Han's coefficients for the number of variables. L-BFGS uses linalg's axpy and dot. Differential evolution evaluates each generation through one values call, which lets basedpl evaluate a generation as one array, and then polishes its best point with Powell's method. Its generator is random::SplitMix64, seeded from Options::seed, which the truncated SVD's test matrix also uses.
root finds a zero of one equation in one variable, from finite bounds as its bracket or by stepping outward from the start. Levenberg–Marquardt solves each damped step through linalg::lstsq on the stacked matrix [J; √λ D], where D holds each variable's largest Jacobian column norm so far. lstsq keeps the step accurate where JᵀJ would square J's condition number.
solve_linear runs CG, MINRES and restarted GMRES through Products. It counts products against the limit, and for least squares it applies AᵀA as one product with A and one with Aᵀ. Task::least_squares decides when a system is least squares. Calls of Mapping::precondition count no product. Each method's doc comment gives its stopping test. solve_linear reports the true ‖b − A x‖ from one more product. MINRES and GMRES build their plane rotations through linalg's givens, which the tridiagonal and bidiagonal iterations also use. solve_affine passes solve_linear an Affine mapping, which computes each product from the values of an affine residual, as (f(t v) − f(0)) ÷ t. Its doc comment says how it picks t.
basedpl's ⇄ and the Python binding share Task, Method::choose's default policy, Options::check's rules for options that a front rejects, and Method::NAMES, the methods' names. Every tolerance is relative, with an absolute floor near the tolerance squared. On Rosenbrock's function from (-1.2, 1), Nelder–Mead took 364 evaluations, Powell's method 423, differential evolution 1767 including its polish, and L-BFGS 53, each ending within 10⁻⁸ of the minimum.
Inner<T> objects compile matrix products in rustymath. On macOS 15+, Accelerate computes float and complex plus-times products. Elsewhere the native product blocks for caches and calls register micro-kernels, with explicit vectors on supported AArch64, AVX2/FMA and wasm SIMD targets. tiles owns target tile sizes. Integer products select a narrower exact kernel when their bounds permit it; Inner::checked rejects possible overflow. Boolean products use packed bits, and equality-count products use kernels appropriate to their item width. Inner::inner_mapped takes operands whose items an agreement::Mapping reads, as a pending copy in basedpl holds them. Its default reads each operand whole. In blocked, float and complex PlusTimes use an operand where it lies when its mapping reads it at positive strides with one stride 1, as a strided or transposed matrix. They read any other operand a block of stored rows at a time, and multiply each block as it's read. A product with a vector, whose other operand's rows are adjacent in storage in runs of at least 64, as the windows of a vector or a single-channel image are, goes to by_columns. It adds each column of that operand, read where it lies and scaled by its item of the vector, into the result. With 4 or more columns, or runs of 3 or 16 rows, gemm over gathered blocks measured faster. Python's matmul, and inner with add and multiply, call it for a float or complex view where Mapped::in_place says gemm reads both operands where they lie. They copy any other view first, such as reversed rows, because one copy and one gemm took less time than gemm over blocks. signal::correlation gives the reads of a correlation or convolution: the array with the fills its mode adds, its windows, and the kernel, reversed for a convolution. signal::correlate runs them through inner_mapped, which builds no windows matrix for float and complex PlusTimes. Python's correlate and convolve use it for sums of products of floats and complex numbers. Other pairs and dtypes gather the windows matrix and call the same product as inner.
One private Matrix view and Tiled's multiply, gemm, lower-triangle gram and in-place triangular serve public products and decompositions. Views carry leading dimensions and transpose/conjugate operations. Tiled also holds each item type's dense drivers, svd, eigh and eig, behind fronts in factor.rs that handle nonfinite and empty matrices. On macOS, accelerate.rs serves all of these except complex eig from one table of BLAS (GEMM, SYRK/HERK, TRSM) and LAPACK (GESDD, SYEVD/HEEVD) symbols, plus a real-only GEEV table, through the ILP64 interface, which takes any size. Complex eig uses the native solver on every target, because Accelerate's ZGEEV and CGEEV took 2.4 and 2.9 times as long. Real eig takes GEEV from accelerate::GEEV_FROM (256) rows, and the native solver below that, which was 1.2 to 2.7 times as fast up to 128. build.rs sets the accelerate cfg on macOS, and the native cfg wherever rustymath's own kernels compile: on other targets, and on macOS with the native feature. native is a default feature. Dev builds keep it. Mac release builds turn it off with --no-default-features. Dependent crates, such as rustymath-py, depend with default-features = false and make native their own default. The native plus-times products, triangular solves, SVD and Hermitian solver compile only where native holds. Where both compile, dispatch! uses Accelerate unless RUSTYMATH_NATIVE=1 selects rustymath's own kernels for the whole process. accelerated() says which one serves it. LAPACK reads a row-major matrix as its column-major transpose. GESDD factors a tall matrix through its adjoint, because Accelerate's GESDD took nine times as long on a 4096×16 matrix that it saw as wide. HEEVD reads the lower triangle as its upper one. GEEV takes a column-major copy. Every LAPACK working copy takes its leading dimension from stride. Unpadded, a 512×512 matrix took twice as long in GEEV as one of 511 or 513, 27% longer in GESDD and 8% longer in SYEVD. Each thread's BLAS threading follows rustymath's thread setting: single-threaded where parallel::serial() holds, as on a batch's worker threads and after set_threads(1), and Accelerate's own threads otherwise. It changes only when that choice does. NumPy's calls into Accelerate on the same thread use the same setting. Rustymath needs macOS 15 at run time, for Accelerate's BLASSetThreading. Its pyproject.toml gives the wheels that minimum. Cargo never reads a dependency's .cargo/config.toml. A crate that ships a binary therefore sets MACOSX_DEPLOYMENT_TARGET=15.0 in its own, as basedpl does. Without it, the binary still links but records Rust's default minimum, 11.0.
Scalar and Factor cover f32/f64 and complex numbers with either parts' width. Scalar::Real is the type of singular values and Hermitian eigenvalues. Scalar method names avoid clashes with float traits. Complex division uses scaled Smith arithmetic; public divide keeps its faster, narrower-range formula. Native complex products form each part from two real products of the parts, as scalar multiplication does. Accelerate uses complex GEMM directly. Nonfinite complex products need not match scalar arithmetic. Singular/eigenvector signs and phases are not canonicalized.
linalg/factor.rs owns LU, Cholesky, QR and their solves. LU factors a column-major copy of the whole matrix in place. Lu keeps the factors in that layout. lu_panel halves the columns until LEAF or fewer remain. It factors the left half, applies its row swaps to the right half across threads through swap_rows, updates the right half through substitution and a packed product, then factors the right half. LU's solves go through triangular. Its inverse solves from the permuted identity where Accelerate serves the process, and from the identity on rustymath's own kernels, which skip the zeros above the diagonal of L⁻¹. Cholesky works in blocks of cholesky_block columns. Trailing updates and blocked triangular solves use packed products. QR stores Householder blocks in WY form; one Block::apply serves QR updates, Q construction and applying Q or its adjoint. QR factors each block of qr_width columns through panel_qr, which halves the block and updates the right half and T through products. Block::new builds T for the eigen and SVD blocks through the same halving. LU, Cholesky and QR work in copies whose rows, or LU's columns, lie stride(n) items apart. stride adds 8 items to a length that's a multiple of 128, which keeps a power-of-two size from mapping every row to the same cache sets. LU and Cholesky close their copies up to n items apart before they return. Cholesky reads only the lower triangle. Thin qr and lstsq of a tall matrix first try tall_qr, CholeskyQR2 with a partial Gram matrix for each thread, above the cutoffs in tall_for_qr and tall_for_lstsq. The cutoffs admit the sizes where it measured faster on one thread, and larger sizes where its threaded passes are faster on the default thread count. Where Accelerate runs the products, the larger sizes are admitted on one thread too. There, 1000000×10 is about 3 times as fast as Householder. Elsewhere they're admitted only where the call can split across threads. Wasm and set_threads(1) calls use Householder for them. It falls back to Householder where a Cholesky factor fails or Q₁ᴴQ₁ is far from I. For tall or square matrices, lstsq reduces b to Qᴴ b without storing Q, then solves through R's SVD. Wide matrices solve through their SVD directly.
solve and pinv use LU when the square matrix's pivots pass the scale-dependent eligibility check. Otherwise they use the least-norm SVD path. The SVD cutoff is shared with rank; LU applies no singular-value cutoff. Consequently rank can be smaller than the dimensions used by a successful LU solve. Keep this distinction rather than forcing both routes to agree at a rank boundary. Lu supplies solves, inverse, factors, determinant and log-determinant. Lu and Svd are generic over their storage: LuRef and SvdRef are the borrowed forms, which new builds from slices and as_ref gives from an owned factorization, with the same methods. qr_complete forms full Q.
The native SVD and Hermitian solvers compile only where the native cfg holds. The general solver compiles everywhere: complex matrices use it at every size, and real ones use it on macOS below GEEV_FROM rows. factor/eigen.rs reduces Hermitian matrices to real tridiagonal form. It reads only the lower triangle. Small vector problems use implicit QL; larger ones use divide and conquer with deflation and secular roots. Reconstructing the secular update from the roots maintains orthogonality. Values-only calls omit rotations of eigenvectors. factor/general.rs reduces to Hessenberg/Schur form and back-substitutes eigenvectors: real matrices use Francis double shifts, complex matrices use single shifts. Padding avoids cache-set conflicts in column reads. Eigen routines return None on failure. GEEV balances and permutes first, so on macOS a real triangular matrix of GEEV_FROM or more rows needn't give its eigenvalues in diagonal order.
factor/svd.rs reduces to a real bidiagonal problem. Small vector problems use Golub-Kahan QR; larger ones use divide and conquer. Values-only calls use QR at every size without constructing or rotating singular vectors. Tall matrices first reduce to R, through tall_qr where qr or lstsq would use it. Power-of-two normalization protects the shift/rotation arithmetic. factor/truncated.rs holds truncated_svd, which runs on every target: it approximates the leading singular values and vectors with a randomized range finder from a fixed seed, through products, thin QR factorizations and one small SVD. Nonfinite input or nonconvergence gives NaN outputs. Pseudoinverse forms the scaled vector product directly. Source constants own block widths and algorithm cutoffs; test their boundaries through the hidden exports rather than duplicating values in tests. Tests compare a scaled matrix's singular values, Hermitian eigenvalues and least-squares solutions with the unscaled ones to rounding, because LAPACK's internal scaling isn't a power of two. Scaling can change the signs, phases and order of singular vectors and eigenvectors. Tests check the scaled matrix's vectors through their own residuals and orthogonality or unit length. They check its general eigenvalues through their sum, which must equal the trace.
A change that adds a user-facing public Rust function also adds its Python binding, following the existing API's names, argument order and dtypes. Every binding follows the conventions of the Rust core that BasedPL shares. Related bindings agree with each other. A function with NumPy's name needn't behave as NumPy's does. remainder and floor_divide are tolerant, set functions keep first-occurrence order, orderings put NaN last, and diagonal and trace take the last two axes, as the matrix functions do.
Python operations are eager. python/rustymath/create.py gives fill constructors their names/defaults over one native fill primitive, and names the dtype facts that Rust returns. Sequence generation stays in create. Structural functions and bracket indexing in py/src/structure.rs delegate selector decoding to buffer/src/selection.rs, shared with BasedPL, for the selection and writeback paths. Brackets use structure::select with broadcasting, as BasedPL's .[] does: index arrays pair by broadcasting, by arithmetic's rule for names. A selector with one axis name selects along the axis with that name (structure::addressed). Basic affine selections share storage. Assignment detaches shared storage, preserving independent values. rm.at combines through the same selection with Array::combine_at and its function's Scatter kernels, sharing Array::assign's checks. A small view retains its whole base until compact() copies its logical items, including when it is a contiguous subrange. Methods call the same binding functions as package-level operations. python/rustymath/structure.py composes conveniences such as trace and repeated differences. Differences run subtraction kernels on two mappings of the input, not gathered copies. python/rustymath/random.py subclasses the native generator for distribution conveniences and multivariate-normal sampling. py/src/random.rs owns its RNG. Sampling, rolls and deals stay in the shared core.
py/src/array.rs holds rustymath.Array, which wraps a rustymath array::Array. numpy.asarray receives a read-only snapshot through rustymath_buffer::export. Each export pins its exact backing allocation and signed-stride geometry until __releasebuffer__, independently of later assignment to the Python owner. rm.array(x) shares existing Rustymath storage with independent values and copies foreign inputs. py/src/convert.rs retains the call's ownership, scalar classification and NEP-50 promotion policy. Its arguments expose core ArrayViews. Dense callers prepare once through core Items. Elementwise callers keep their mappings. Arrays of different dtypes promote through DType::common. Nested lists, tuples, range, dict and str follow the same converter. give returns a rank-0 result as a Python number when every data argument was a Python scalar. rustymath_buffer::read_array reads sibling snapshots with keys and names. The shared buffer reader normalizes signed bytes, wider unsigned integers and foreign byte order, borrows contiguous aligned items, and uses CPython's PyBuffer_ToContiguous to copy other layouts into C order. It doesn't use PyO3's PyUntypedBuffer, which rejects rank-0 exports. The crate has no Rust tests, because PyO3's extension-module feature unifies across the workspace, and a test that embeds Python can't link. tests/test_array.py tests it through the class.
py/src/functions.rs holds the elementwise functions, reductions and scans, one entry each in the functions! and reductions! macros. A function converts its data arguments through convert, finds its object's kernels with kernels!, and runs them at the dtype that Kernels::computes gives. convert promotes the arrays' dtypes with DType::common. It then applies DType::common_weak for the Python values, as NEP 50 does for weak scalars. BPL shares both rules. An argument whose dtype has no kernel computes at the first dtype that has one: int64, then float64, then complex128 for Booleans and integers, and a wider float, then a complex dtype at least as wide, for floats. Kernels casts the argument to that dtype first. The object's slice form then reads the items in place and splits independent ranges across threads, as it does for Rust callers. checked turns a kernel's result into the call's: the missing error for a dtype that the function doesn't take, the wrapped integer result where Kernels::wraps says that the kernel wrote one, as NumPy does, and OverflowError for any other integer result that doesn't fit its dtype. The tolerant functions take one scalar ct.
With no axis, a reduction reduces the whole array and gives a Python number, and a scan runs over the flattened array. Reductions take axis positions, names or tuples of them. accumulate(f, x, axis=None, start=None) scans any native binary function from the left, with an optional scalar seed. The same functions! entries supply its callbacks and inner's. Known reductions retain their optimized scans. Comparisons keep the computing dtype, reading later results as 0 or 1. A sort, argsort or lexsort, takes one axis and works along the last axis by default. Grading shares native-width key construction and stable sorting with BPL, including heterogeneous lexsort keys. Folds and scans run through reduce::along, scan_along and left_scan. A fold reads the items that ArrayView::fold_items gives: a view that reads a block of its storage, such as a transpose or a reversal, lends the block, which folds in place, and any other view gathers its items first. minmax over the whole array takes reduce::extremes' single pass for Booleans and integers, and two folds otherwise. logical_scan takes a 2×2 truth table.
Statistical bindings in py/src/functions.rs share conversion and axis preparation through stats::Axes. Numerical loops live in stats, including two-pass variance, weighted averages, covariance and histogram counting. Float16 statistics accumulate and return float32. Complex variance and norms return real values.
py/src/linalg.rs runs rustymath's matrix functions on stacks of matrices. Stack converts a call's data arguments through convert and chooses the dtype that the function computes in. Each float or complex width computes at its own width. float16 computes in float32 and gives float16 results. Booleans and integers compute in the float that DType::inexact gives them: float32 for Booleans, uint8 and int16, and float64 for int32 and int64. rustymath::agreement::Frames pairs the arguments' frames, the axes before each argument's cell, by rustymath's broadcasting rule. It gives each argument's cell for each cell of the result. run runs a function on each cell of the result through parallel::stack_with, with scratch such as an LuWork for each thread, and checks for Ctrl-C between rounds. A right-hand side of one axis is one vector.
Singular values and Hermitian eigenvalues take the real dtype of the matrices' width. General eigenvalues and eigenvectors take the complex dtype. Row positions and ranks are int64. A matrix that a decomposition can't factor gives NaN results, with no report. lu returns the packed factors in rustymath's own order, where row j holds column j. lu_solve, lu_inverse, lu_det and lu_slogdet read them in place. lu_lower and lu_upper build the separate triangular matrices that Lu.l and Lu.u give. product shares frame and keyed-axis preparation for matmul and inner(f, g). Recognized pairs keep their specialized kernels. Other native pairs call linalg::inner_with per cell through run, then join the results. The existing functions! entries generate their bulk binary and right-fold callbacks; no Python callbacks run inside a contraction. Empty contractions use existing reduction units or raise when none is known. Integer products raise OverflowError for a result that doesn't fit int64. python/rustymath/linalg.py names the functions and builds the Lu, Svd, Eigen and Echelon results, whose methods call the Rust functions. Preserve the intentional API differences: svd returns V, not NumPy's V-adjoint, and complete QR is selected with complete=True.
py/src/exact.rs reads integers, Fractions, nested lists of them and integer buffers through the shared walk in convert.rs, computes with BigRational, and returns Fractions, in nested lists for arrays. Its elementwise functions broadcast through Agreement, and its matrix functions check their own shapes, so python/rustymath/exact.py only names the functions. Exact matrix calls take one matrix, and solve takes a vector or a matrix as its right-hand side. The ordinary linalg module also exposes permanent and 2-D row reduction. Both row-reduction APIs return Echelon(rows, pivots, determinant).
The remaining binding files match their Python modules. The polynomial functions compute in float64 for Booleans and integers, and in float32 for float16, rounding the result back. Integer mixed-radix conversion checks int64 overflow. Number theory takes Python integers through BigUint and polls for Ctrl-C. The search functions sort, find insertion positions, classify and group through the search module's functions on views, which basedpl's searches also call. search.py adds the Unique tuple and bins, and the package root reexports the search functions. Sorts compare exactly, equate signed zeros and put NaNs last. Distribution constructors live in a submodule, so their names don't collide with the elementwise functions'.
py/src/optimize.rs runs rustymath::optimize on Python functions, and python/rustymath/optimize.py gives it SciPy's names: minimize, root_scalar, least_squares, and the linear solvers cg, minres, gmres, cgls and lsmr. Each returns a named tuple: OptimizeResult, RootResults, LeastSquaresResult or LinearResult. The module stays out of the package root, where root is the nth-root ufunc. A float32 or float16 start, or b for the linear solvers, computes in float32, and anything else in float64. Each evaluation checks for Ctrl-C, then calls the Python function with a rustymath array. read takes the items of what it returns through the shared converter, cast to the search's float type, in row-major order. A function can therefore return any shape, such as a 2-D Jacobian. An exception from the function stops the search. The wrapper turns SciPy's bounds pairs and its maxfev, popsize and seed options into the binding's arguments. The binding's first call of least_squares's fun, at x0, gives the number of residuals and becomes the search's first evaluation. least_squares's jac gives Vectors::jacobian the whole Jacobian in place of vector–Jacobian products. The linear solvers take a matrix, whose products Product computes in Rust through PlusTimes.inner_ops, or anything with a shape and a matvec, as SciPy's LinearOperator has. cgls and lsmr use the transpose, through rmatvec or the matrix, because the core solves least squares for any matrix with a transpose. M takes the same forms and becomes the core's preconditioner. Our gmres applies M on the right, where SciPy's applies it on the left. SciPy's lsmr takes no M. maxiter keeps SciPy's meaning, iterations or GMRES's restart cycles, and linear_at converts it to the core's limit on products. Options::check and Method's FromStr give the argument errors, so the binding repeats none of the core's rules. With vectorized, the function gets points as the columns of a 2-D array, as SciPy's differential_evolution passes them.
lib.rs adds each module's functions, and set_threads, thread_setting and threads, to the extension module.
Rust tests use the public API in tests/accuracy.rs, slices.rs, linalg.rs, numeric.rs and speed.rs. The library sets test = false to avoid a separate cfg(test) build. Because no build runs #[cfg(test)] code, the library holds no inline tests. Any internal behaviour worth a test is worth knowing for callers too. These tests check it through the public API, as walks_compose checks Steps::through. Shared support::random and seeded use SplitMix64 without additional rand features. Other helpers stay local. Combine related checks into readable tests rather than splitting every assertion into a test.
Python tests are grouped by topic in tests/test_ufuncs.py, test_linalg.py, test_search.py, test_exact.py and test_optimize.py. A test that checks threaded paths takes the threads fixture from tests/conftest.py, which runs it on one thread and on four.
Elementary real accuracy checks use a four-ulp bound with an additional allowance for f64 libm references. Gamma and incomplete functions use scaled errors and function-specific bounds; f64 iterative checks have a wider allowance. tests/accuracy.rs defines those bounds. Real checks cover fused/unfused arithmetic, f16 inputs, sampled wider inputs and reference files. Complex checks exercise public scalar/slice forms at the build's FMA mode, including branch cuts. testdata/references.py generates the mpmath references; preserve their notices and THIRD-PARTY.md.
nbs/index.ipynb generates the README and home page, with short examples linking to feature notebooks. nbs/sidebar.yml lists those pages. src/lib.md supplies the Rust documentation and runnable doctests; nbs/rust.qmd symlinks it as a standalone Quarto page. Doc comments in .rs files avoid doctests. They show an example as a text block. Edit those sources, not the generated README. nbdev-test checks notebook code; nbdev-preview serves saved outputs because _quarto.yml disables execution. .github/workflows/docs.yml deploys the site.
Test runs on pushes to main and pull requests. It uses the shared test-profile build before Python/notebook tests. Speed is manual. Release runs separately on v* tags, builds dist wheels/sdist, publishes rustymath then rustymath-buffer to crates.io and the Python artifacts to PyPI, and attaches artifacts to the GitHub release. It does not wait for Test or run tests itself.
Before release approval, verify successful tests for the release commit and run cargo develop && pytest -q for the installed extension. Confirm [workspace.package].version in Cargo.toml; Python gets it through dynamic = ["version"]. After Jeremy approves, run ship-release. It pushes the current version's tag to start publication, then bumps Cargo's version. Never commit, push, tag or publish without approval.