rustymath

rustymath gives you maths in plain Rust, and in Python with NumPy, from fast functions on 16-, 32- and 64-bit floats to exact rationals. It covers the elementary functions, such as exp, log, sin and pow, at 16, 32 and 64 bits. Most of them also take complex numbers. It also has special functions such as gamma and erf, probability distributions, reductions and scans, searching and sorting, linear algebra, polynomials, number theory, and APL-inspired tolerant comparison. In both languages, it also computes exactly on rationals: matrix solves and determinants, polynomials, gcds and binomial coefficients.

In Rust, you call each function on single values inside your own loops, or through an object that runs it over a whole slice. In Python, each maths function is a NumPy ufunc. The linear algebra functions run as generalized ufuncs, over stacks of matrices.

Its float functions give up the last bits of accuracy for speed. Every elementary and special function is within 4 units in the last place (ulps) of the correctly rounded result. At 32 bits, most functions run two to six times as fast as the standard library’s versions.

Install

Install the Python package from PyPI. For Rust, see Using rustymath from Rust.

pip install rustymath

Using rustymath from Python

The maths functions are NumPy ufuncs. Each takes arrays of any shape, broadcasts them as NumPy does, and gives a result of the same width. A complex array works wherever the function has complex values. An integer array works wherever NumPy’s version takes one. add, multiply, maximum and minimum have reduce and accumulate methods that run rustymath’s vectorized sums, products, maximums and minimums. Import rustymath beside NumPy:

import numpy as np, rustymath as rm

Where NumPy has a function, rustymath uses the same name and arguments. To switch, you change only the import. Here exp keeps the 32-bit width of its argument:

x = np.linspace(0, 2, 5, dtype=np.float32)
rm.exp(x)

array([1.       , 1.6487211, 2.718282 , 4.4816885, 7.389056 ],
      dtype=float32)

To check the accuracy yourself, compare with NumPy. np.spacing gives the size of one ulp at each value. Dividing the difference by it gives the distance between the two results in ulps:

xs = np.linspace(-10, 10, 100_001)
y = np.exp(xs)
(np.abs(rm.exp(xs) - y) / np.spacing(y)).max()

np.float64(2.0)

Functions that NumPy lacks keep their own names. root(n, y) gives the nth root of y:

rm.root(3, np.array([8., 27., 1000.]))

array([ 2.,  3., 10.])

The Rust names follow Rust’s standard library instead. Python’s log, power, arcsin and arctan2 are Rust’s ln, pow, asin and atan2.

Complex numbers

A complex array gives complex results, at each function’s principal value. For exp, log, power, sin, cos, sinh and cosh, rustymath computes a complex array through its vectorized real loops:

rm.sqrt(np.array([-4+0j, 3+4j])), rm.log(np.array([-1+0j]))

(array([0.+2.j, 2.+1.j]), array([0.+3.14159265j]))

Special functions and distributions

The special functions take SciPy’s names, such as gamma, gammaln, erf, erfinv, betainc and lambertw:

rm.gamma(np.array([0.5, 5, 10])), rm.erfinv(0.5)

(array([1.77245385e+00, 2.40000000e+01, 3.62880000e+05]),
 np.float64(0.4769362762044699))

rm.distribution holds sixteen distributions. Each constructor takes the distribution’s parameters, or uses standard ones where the distribution has them. The result has pdf, cdf, quantile and sample. The same seed gives the same draws:

n = rm.distribution.normal()
n.cdf([-1.96, 0, 1.96]), n.quantile(0.975), n.sample(3, seed=0)

(array([0.0249979, 0.5      , 0.9750021]),
 1.9599639845400538,
 array([ 0.11623959, -0.78285485, -2.33463531]))

Polynomials

rm.polynomial evaluates polynomials, finds roots, constructs coefficients from roots, differentiates and fits data, under the names of numpy.polynomial.polynomial. Coefficients start from the constant term. Evaluation, construction and derivatives take real or complex arrays. polyvalfromroots(x, roots, m=1) evaluates directly from roots without forming coefficients. The roots of x^3 - 6x^2 + 11x - 6 are 1, 2 and 3:

rm.polynomial.polyroots([-6, 11, -6, 1])

array([1.+0.j, 2.-0.j, 3.+0.j])

Linear algebra

rm.linalg provides real and complex matrix solves and decompositions, including solve, det, lu, cholesky, qr, svd, eigh, eig, lstsq, pinv and matrix_rank. They take stacks of matrices and broadcast the leading axes. solve and pinv also handle singular or non-square matrices through least-norm least-squares solutions. permanent takes square matrices; echelon reduces one 2-D matrix. LU factors store packed data and row positions, with .l and .u properties for the triangular arrays. Reuse the factors for solves, inverses and determinants:

a = np.array([[4., 1], [1, 3]])
f = rm.linalg.lu(a)
f.solve([1., 2]), f.determinant(), rm.linalg.eigh(a).values, rm.linalg.solve([[1., 2], [2, 4]], [1., 2])

(array([0.09090909, 0.63636364]),
 np.float64(11.0),
 array([2.38196601, 4.61803399]),
 array([0.2, 0.4]))

rm.matmul multiplies matrices as np.matmul does, with integer matrices in int64. rm.linalg.semiring(f, g) gives the matrix product whose sums use f and whose products use g, as APL’s inner product f.g does. With add and equal, it counts the matching items. Squaring a matrix of edge lengths with minimum and add gives the shortest paths of at most two edges:

d = np.array([[0, 3, np.inf], [np.inf, 0, 1], [2, np.inf, 0]])
minplus = rm.linalg.semiring(rm.minimum, rm.add)
minplus(d, d)

array([[0., 3., 4.],
       [3., 0., 1.],
       [2., 5., 0.]])

Optimization

rm.optimize takes SciPy’s names. minimize runs Brent’s method, Nelder–Mead, Powell’s method, L-BFGS or differential evolution. Without a method, it picks one from the number of variables, whether jac gives a gradient and whether there are bounds. root_scalar finds a zero of a function of one variable by Brent’s method, within a bracket or by stepping outward from a start. Here Nelder–Mead finds the minimum of Rosenbrock’s function, at (1, 1), and Brent’s method finds where cos x = x:

rosen = lambda x: (1-x[0])**2 + 100*(x[1]-x[0]**2)**2
rm.optimize.minimize(rosen, [-1.2, 1]).x, rm.optimize.root_scalar(lambda x: np.cos(x) - x, [0, 1]).root

(array([1., 1.]), 0.7390851350713885)

least_squares fits by Levenberg–Marquardt. cg, minres, gmres and lsmr solve linear systems through the matrix’s products alone, so the matrix can be an array or a SciPy LinearOperator. Here least_squares recovers a = 2 and b = 0.3 from five points on y = a e^(bx), and lsmr gives the least-squares line through five points that no line fits:

t = np.linspace(0, 2, 5)
fit = rm.optimize.least_squares(lambda p: p[0]*np.exp(p[1]*t) - 2*np.exp(0.3*t), [1, 1])
line = rm.optimize.lsmr(np.column_stack([np.ones(5), np.arange(5.)]), [1, 2.9, 5.2, 6.8, 9.1])
fit.x, line.x

Number theory

rm.number_theory tests primality, factorises, finds the next or previous prime, gives Euler’s totient and exact factorials, and lists and counts primes, on Python integers of any size. Ctrl-C stops a long call:

nt = rm.number_theory
nt.factor(2**64 + 1), nt.prime_count(10**8), nt.next_prime(10**20)

([(274177, 1), (67280421310721, 1)], 5761455, 100000000000000000039)

Exact arithmetic

rm.exact takes integers and fractions.Fraction. It provides rational gcd, lcm, remainder and binomial coefficients, matrix solves, pseudoinverses, determinants, permanents and row reduction, and polynomial construction, evaluation and derivatives. Elementwise calls broadcast and return Fractions or NumPy object arrays. Matrix calls take one 2-D matrix. Ordinary rm.linalg calls stay floating-point.

from fractions import Fraction

rm.exact.solve([[2, 1], [1, 2]], [Fraction(1, 2), 0]), rm.exact.binom(Fraction(3, 2), 2)

(array([Fraction(1, 3), Fraction(-1, 6)], dtype=object), Fraction(3, 8))

Tolerant comparison

APL compares floats within a relative tolerance. Two floats are equal when their difference is less than the tolerance times the larger magnitude. With the usual tolerance, 0.1 + 0.2 equals 0.3. rm.equal, not_equal, less, less_equal, greater and greater_equal compare this way. Their tolerance ct defaults to rm.CT, 1e-14. Pass ct=0 to compare exactly:

rm.equal(0.1 + 0.2, 0.3), rm.equal(0.1 + 0.2, 0.3, ct=0)

(np.True_, np.False_)

floor, ceil, remainder, gcd and lcm take ct too. A value within tolerance of a whole number counts as that whole number:

rm.floor(np.array([2.9999999999999996, 2.5]))

array([3., 2.])

remainder takes its arguments in NumPy’s order. rm.remainder(x, y) is the Rust residue(y, x).

unique, index and isin search floats within the same tolerance. They match Booleans and integers exactly. unique gives the values in order of first occurrence, with each one’s first position, each item’s class number and each class’s size. index(x, y) gives where in x each item of y first matches, or the size of x where nothing does:

rm.unique([0.3, 1.0, 0.1 + 0.2]), rm.index([0.3, 1.0], 0.1 + 0.2)

(Unique(values=array([0.3, 1. ]), firsts=array([0, 1]), numbers=array([0, 1, 0]), counts=array([2, 1])),
 np.int64(0))

rm.search has stable argsort and lexsort, searchsorted, integer bincount, class counts and group_indices. These functions are also available directly on rm. Sorting and insertion searches compare exactly and put NaNs last. rm.integer.encode and decode convert mixed-radix integers with checked int64 arithmetic; rm.polynomial.decode and rm.exact.decode accept field-valued bases and digits.

rm.minmax(x, axis=-1) computes both bounds in one pass. Boolean logical_and, logical_or, logical_xor and logical_xnor have reduction and accumulation methods. rm.logical_scan(x, table, axis=-1) scans any Boolean dyad given by its 2×2 truth table. rm.set_threads(n) sets the calling thread’s limit and returns its previous setting. Use 0 for the default or 1 for serial calls.

Using rustymath from Rust

Add rustymath to your Cargo.toml:

[dependencies]
rustymath = "0.1"

On x86-64, build for a recent CPU, with -C target-cpu=x86-64-v3 or -C target-cpu=native. You can set it once in .cargo/config.toml:

[target.'cfg(target_arch = "x86_64")']
rustflags = ["-C", "target-cpu=x86-64-v3"]

Rust’s default x86-64 target has no fused multiply-add and no vector instruction for floor or trunc. Without the setting, rustymath multiplies and adds separately. Some of its loops also run one item at a time.

Each function comes in two forms. The free function, such as exp, is generic over f32, f64 and half::f16. The compiler inlines it into the loop that calls it and vectorizes that loop:

use half::f16;
use rustymath::exp;

let e64 = exp(1.0f64); // 2.718281828459045
let e32 = exp(1.0f32); // 2.718282
let e16 = exp(f16::from_f32(1.0)); // 2.719

The function object, such as Exp or Pow, has a scalar form one and a slice form many. Both run code that rustymath compiles itself. Your build settings don’t change their speed. Each argument of many is one item, which every result pairs with, or a slice with one item for each result. many appends the results to a vector and returns whether every item succeeded. A NaN result from arguments that aren’t NaN counts as a failure. An object implements Unary, Binary or Ternary at each item type it takes: f16, f32 and f64, complex numbers where the function has complex values, and integers for some, such as Gcd at i64. Code that is generic over the float type states the bound it needs, such as where Exp: Unary<F, Out = F>. The free functions need no bound.

use rustymath::{Arg, Binary, Pow};

let mut squares = Vec::new();
Pow.many(Arg::Many(&[1.5, 2.0, 3.0]), Arg::One(2.0), 3, &mut squares); // squares: [2.25, 4.0, 9.0]

Use the slice forms of Sin, Cos and Tan for arrays. They reduce each item by multiples of π with a few multiply-adds, then give the rare items beyond 2^23 at f64, or 2^20 at f32, an exact reduction in integer arithmetic. The free sin takes a branch for those items. That branch stops a loop that calls it from vectorizing.

atan2(y, x) and its object Atan2 give the angle of the point (x, y), from -π to π. y comes first, as in f64::atan2.

The complex slice forms of Exp, Ln, Pow, Sin, Cos, Sinh and Cosh vectorize too. Each one splits a chunk of complex items into their real and imaginary parts and computes the results from the real slice forms. On macOS, the sines and cosines they need come from vForce, which computes both in one call. It then redoes the rare items with the scalar complex function: items with a tiny or non-finite part, and Pow’s special cases. Abs, Angle and Sqrt keep their loop over items, which runs faster than a split for their short formulas. Other complex slice forms, such as Lgamma’s, compute one item at a time.

divide(x, y) divides complex numbers, inline in the caller’s loop, as the free functions do. Two real arguments divide as reals, so a real zero divisor gives ±∞ or NaN with an imaginary part of 0. Its formula’s |y|² overflows beyond about 1e154 at f64 and 1e19 at f32.

Tolerant comparison

The tolerant module has APL’s tolerant comparison: x equals y when |x - y| is less than the tolerance times the larger of |x| and |y|. Its objects, such as Equal, Less, Floor, Residue and Gcd, each hold a Tolerance. With Tolerance::EXACT, a slice form runs a loop of plain IEEE comparisons, chosen once per call.

use rustymath::{tolerant::{Equal, Tolerance}, Arg, Binary};

let mut same = Vec::new();
Equal(Tolerance::DEFAULT).many(Arg::Many(&[0.3, 0.1 + 0.2, 0.31]), Arg::One(0.3), 3, &mut same); // same: [true, true, false]

Other modules

• number_theory: primality, factorisation, next and previous primes, Euler’s totient, and the enumeration and counting of primes, on a u64, a u128 or a BigUint. Each computes in the narrowest of those that holds the value. A function that can run for a long time also has a _with version, whose extra argument is a check that can stop it. factorial gives n! as a BigUint.
• optimize: minimize finds the minimum of a function of an f32 or f64 slice, by Brent’s method for one variable, Nelder–Mead, Powell’s conjugate directions, L-BFGS with a gradient, or differential evolution within bounds. solve finds a zero of a function that gives a slice, or its least-squares solution, by Brent’s method for one equation in one variable and Levenberg–Marquardt otherwise. solve_linear solves a linear system, or its least squares, through the matrix’s products with vectors, by conjugate gradients, MINRES or restarted GMRES, with an optional preconditioner. solve_affine solves f(x) = 0 for an affine function the same way, through its values alone. Method::choose picks a method where the caller doesn’t. Functions are Objectives and Mappings: closures that return a Result, whose error stops the search, or types that also give a gradient, vector–Jacobian products or batches of points.
• polynomial: evaluation over a slice by Horner’s rule, roots from the eigenvalues of the companion matrix, sorted by real part, then imaginary part, building from roots, derivatives and least-squares fitting. value and factored evaluate at one point. They, from_roots and derivative take any type with the field operations, such as exact rationals. So does decode, which gives the values of digits in a mixed radix.
• reduce: sums, products, maximums and minimums over a slice, and their running versions. Reduction gives each one’s unit, its step for one more item, its fold of a slice, its merge of the folds of two neighbouring runs, and its scan, through Sum, Product, Max and Min. Its checked is the step, and whether the result fits. Float sums and products vectorize by combining the items in any order. An integer sum or product, and each running one, wraps where it doesn’t fit, and reports whether it fits. extremes gives the least and the greatest item in one pass. Alternating gives the right fold by - through Sum, and by ÷ through Product: the sum of the pairs’ differences, or the product of their ratios. A quotient of floats reports that it doesn’t fit where a running product leaves the normal range. Booleans fold by Max and Min as or and and, by Xor and Xnor as not-equal and equal, and by SumAs as a count of trues. Logical scans Booleans by any function of two Booleans, given as its truth table. SumAs<C> also sums integers, wrapping in the type C, which the caller picks to hold every sum: a byte, i16, i32 or i64, through the traits Term and Tally. along and scan_along fold and scan along the middle axis of an array of shape [outer, len, inner] held row by row: down the columns of a matrix, for example, with [1, rows, columns]. fold_into folds the same array one block at a time, for a caller that computes the array in blocks. A row of fewer than 16 items, such as an image’s 3 channels, has a loop of its own for each length. A running maximum or minimum of integers repeats its value over each block of 256 items that leaves it unchanged.
• linalg: inner products of row-major matrices through the Inner trait: the matrix product PlusTimes, which also takes complex matrices, the max-plus and min-plus products MaxPlus and MinPlus, the counts of equal and unequal items PlusEqual and PlusNotEqual, and the Boolean products OrAnd and AndOr. Each numeric pair also takes i64 matrices, whose arithmetic wraps. Inner::checked gives None where an integer result could overflow. They block for the caches, as GotoBLAS and BLIS do. Their innermost loops run on explicit vectors on AArch64, in wasm and on x86-64 with AVX2. Its decompositions of real and complex matrices, at f32, f64 and Complex with parts of either type, are lu, cholesky, qr, qr_complete and svd, with Lu::solve, Lu::inverse, Lu::lower, Lu::upper, Lu::determinant, Lu::log_determinant, Svd::solve, Svd::pseudoinverse, singular_values and the least-squares solver lstsq. solve and pinv give the least-norm solution and the pseudoinverse of any matrix, singular or not. rank counts the singular values above lstsq’s cutoff. eigh and hermitian_eigenvalues take Hermitian matrices, and eig and eigenvalues any square matrix. transpose_into transposes a matrix through 8×8 tiles, applying a conversion to each item. Scalar is the sealed trait of the four item types. linalg::exact works in any field, such as exact rationals: echelon, solve, pinv, determinant and Ryser’s permanent, each with a _with version.
• exact: residue, gcd, lcm and binom of exact rationals, at any integer type, such as BigRational. linalg::exact solves matrices of them.
• distribution: sixteen distributions. The objects Pdf, Cdf and Quantile hold a distribution and evaluate it over a slice. Each distribution also samples from a seeded generator.
• random: rolls and deals that fill a slice from one generator.
• gather: copying a slice’s items into a buffer by position: cycle, repeat, compress by a mask, take at indices, windows, walk along a pattern of Steps with fills, and scatter at offsets, each applying the caller’s conversion to each item. walk also takes Same, which copies whole slices of one type. Steps::reach gives the source offset of one item of a walk, and Steps::run the offset of a run of items that lie in order in the source.
• convert: int converts a whole float to an integer type. It gives None outside the type’s range, and for i64 and u32 beyond a magnitude of 2^51 - 1. float converts an integer to the nearest float. Loops that call either vectorize on every target. AVX2 and wasm lack a vector instruction between f64 and 64-bit integers. On those targets a loop that converts with as runs one item at a time.

The special functions gamma, lgamma, erf, erfc, their inverses, the incomplete gamma and beta functions and Lambert W have free functions and objects at the crate root, as the elementary functions do.

Add, Subtract and Multiply are objects at the crate root, for every float, complex and integer type. Their integer forms wrap, as NumPy’s integer arrays do. Overflowing gives each integer type’s arithmetic with whether its result fits, for callers that keep integers exact.

Accuracy and speed

The table gives nanoseconds per item for a loop over 65,536 inputs on an Apple M5 Max, in October 2026, beside the standard library’s version of each function. A row named after an object times its slice form. Pow.one and Root.one time the scalar forms. To print the table on your own machine, run RUSTYMATH_COSTS=1 cargo test --test speed.

An f16 computes in f32 and rounds once. Special values, such as infinities, NaN and signed zeros, give IEEE’s results. The tests compare each function with a correctly rounded reference at every f16 input, at samples of f32 and f64 inputs, and at published hard cases.

Where the target has a fused multiply-add instruction, as on AArch64 and on x86-64 with FMA, the functions use it. Elsewhere, as in wasm, they multiply and add separately, because a call to a software fma would stop the loop vectorizing. FUSED says which the target uses. The modules fused and unfused hold every function again, with the choice made. A loop compiled with FMA for itself, such as an AVX2 version chosen when your program runs, calls fused. The 4-ulp bound holds either way.

At f64, pow and root compute ln |x| and its product with y as double words, through Arm Optimized Routines’ tables. The result stays within 4 ulps even where y ln |x| nears the overflow threshold. A power or root whose exact result is a float, such as 15^4 or ∛27, comes out exactly. Their scalar forms take one branch for special values, as Arm’s scalar pow does. Their slice forms compute every item the fast way, then redo only the rare items with special values. At f32 and f16, pow and root compute in f64 and round once. There the free pow and root have no branch. A loop that calls them vectorizes.

Licence

Apache-2.0. Some algorithms and constants come from SLEEF and from Arm Optimized Routines. THIRD-PARTY.md has their notices.