# rustymath rustymath gives you fast maths functions for arrays of floats, in Rust and in Python with NumPy. 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, linear algebra, polynomials, number theory, and APL-inspired tolerant comparison. 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. rustymath gives 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](#using-rustymath-from-rust). ``` sh 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: ``` python 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: ``` python 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: ``` python 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 `n`th root of `y`: ``` python 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: ``` python 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`: ``` python 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: ``` python 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.959963984540054, array([ 0.11623959, -0.78285485, -2.33463531])) ### Polynomials `rm.polynomial` has the functions of `numpy.polynomial.polynomial`, with the same names and argument order. Coefficients start from the constant term. The roots of x^3 - 6x^2 + 11x - 6 are 1, 2 and 3: ``` python rm.polynomial.polyroots([-6, 11, -6, 1]) ``` array([1., 2., 3.]) ### Linear algebra `rm.linalg` has rustymath’s linear algebra for real and complex matrices, under the names of NumPy’s similar functions, such as `solve`, `inv`, `det`, `lu`, `cholesky`, `qr`, `svd`, `eigh`, `eig`, `lstsq`, `pinv` and `matrix_rank`. Each takes a stack of matrices in its last two axes and broadcasts the others, as NumPy’s generalized ufuncs do. `solve` and `pinv` give a result for every matrix: for a singular or non-square one, the least-squares solution with the least norm, and the pseudoinverse. They run on rustymath’s own decompositions, written in plain Rust: ``` python a = np.array([[4., 1], [1, 3]]) rm.linalg.solve(a, [1., 2]), rm.linalg.eigh(a).values, rm.linalg.solve([[1., 2], [2, 4]], [1., 2]) ``` `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: ``` python 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) ``` ### Number theory `rm.number_theory` tests primality, factorises, finds the next or previous prime, gives Euler’s totient, and lists and counts primes, on Python integers of any size. Ctrl-C stops a long call: ``` python 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) ### Tolerant comparison APL compares floats within a relative tolerance. Two floats are equal when their difference is at most 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: ``` python 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: ``` python 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: ``` python rm.unique([0.3, 1.0, 0.1 + 0.2]), rm.index([0.3, 1.0], 0.1 + 0.2) ``` (array([0.3, 1. ]), np.int64(0)) ## Using rustymath from Rust Add rustymath to your `Cargo.toml`: ``` 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`: ``` 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: ``` rust 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`. The free functions need no bound. ``` rust 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. 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 at most 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. ``` rust 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. - `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. - `reduce`: sums, products, maximums and minimums over a slice, and their running versions. `Reduction` gives each fold’s unit and its step for one more item, for lanes whose items aren’t consecutive. 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. `running` scans by any `Reduction`, such as Booleans by `Max`. `rows` folds each row of a matrix, with a loop of its own for each row length below 16, such as an image’s 3 channels. `Count` counts the trues among Booleans into a byte or an `i64`. - `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 product `OrAnd`. 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. - `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. - `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_TIMES=1 cargo test --test speed -- --nocapture`.
Function f64 std f64 f32 std f32
Exp 0.84 1.51 0.33 0.90
Expm1 1.07 1.78 0.45 1.64
Ln 1.35 1.52 0.58 1.30
Log1p 1.80 2.18 0.80 1.30
Log2 1.35 1.52 0.58 1.19
Log10 1.37 1.63 0.58 1.19
Sin 0.81 1.99 0.36 1.22
Cos 0.91 2.04 0.40 1.62
Tan 1.12 2.35 0.45 1.89
Asin 0.94 2.30 0.34 1.30
Acos 1.03 2.17 0.37 1.43
Atan 0.89 2.24 0.34 2.18
Sinh 1.25 2.33 0.53 1.46
Cosh 0.95 2.28 0.38 1.48
Tanh 1.24 2.17 0.52 1.47
Asinh 1.48 2.41 0.64 1.76
Acosh 1.57 2.20 0.68 1.58
Atanh 1.40 2.34 0.58 1.52
Cbrt 1.12 1.31 0.83 1.74
Pow 3.08 4.46 1.69 2.00
Pow.one 4.11 4.19 2.26 1.93
Hypot 0.55 1.08 0.26 0.87
Root 3.26 4.00 1.79 1.98
Root.one 4.23 4.03 2.17 1.95
Maximum 0.18 0.19 0.09 0.09
reduce::sum 0.06 0.46 0.03 0.46
reduce::max 0.05 0.11 0.03 0.06
reduce::cumsum 0.14 2.62 0.14 2.50
reduce::cummax 0.12 1.27 0.11 2.54
polynomial::evaluate 0.48 3.26 0.20 2.61
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.