rustymath

rustymath gives you fast versions of the elementary maths functions, such as exp, log, sin and pow, for 16-, 32- and 64-bit floats. You can call them from Rust, where the compiler builds them into your own loops, or from Python, where they work on whole NumPy arrays. It also has APL’s tolerant comparison, number theory and polynomials.

rustymath gives up the last bits of accuracy for speed. Every 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:

pip install rustymath

For Rust, add rustymath to your Cargo.toml:

[dependencies]
rustymath = "0.1"

Using rustymath from Python

Each function is a NumPy ufunc. It takes arrays of any shape, broadcasts them as NumPy does, and gives a result of the same float width. 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 and arcsin are Rust’s ln, pow and asin.

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 compares this way. Its 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 and gcd 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 within the same tolerance. unique and isin take the arguments of numpy.unique and numpy.isin. 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)

(array([0.3, 1. ]), np.int64(0))

Using rustymath from Rust

Each function is generic over f32, f64 and half::f16. It’s an ordinary scalar function with no branches on its usual path. The compiler inlines it into the loop that calls it and vectorizes that loop, as it would for + or *:

use rustymath::{exp, ln};

let xs = [0.5f32, 1.0, 2.0];
let ys: Vec<f32> = xs.iter().map(|&x| exp(x)).collect();
assert!((ln(ys[2]) - 2.0).abs() < 1e-6);
assert_eq!(ln(1.0f64), 0.0);

sin, cos and tan reduce their argument by multiples of π. Up to a magnitude of 2^23 at f64 or 2^20 at f32, a few multiply-adds do it. Beyond that, an exact reduction in integer arithmetic keeps the same accuracy for every finite argument. That reduction needs a branch, and a branch stops the loop vectorizing. checked_sin, checked_cos and checked_tan have no branch. Beyond those magnitudes they give None:

use rustymath::checked_sin;

let sines: Option<Vec<f64>> = [0.5, 1e6].iter().map(|&x| checked_sin(x)).collect();
assert!(sines.is_some());

Tolerant comparison

The tolerant module has the tolerant comparison that the Python functions use. A Tolerance gives tolerant equality, order, floor, ceiling, residue, whole numbers and greatest common divisors, and the range of floats equal to a number. With a tolerance of 0, with_tolerance! picks Exactly instead. Its methods compile to plain IEEE comparisons.

use rustymath::tolerant::{Tolerance, Tolerant};

let t = Tolerance::DEFAULT;
assert!(t.equal(0.1 + 0.2, 0.3));
assert_eq!(t.floor(2.9999999999999996), 3.0);

Number theory

The number_theory module tests primality, factorises, finds the next and previous primes, computes Euler’s totient, and enumerates and counts primes. Each function on one integer takes a u64, a u128 or a BigUint. It computes in the narrowest of those that holds the value. A function that can run for a long time also has a _with version. Its extra argument is a check that can stop it, such as a time limit.

use rustymath::number_theory::{factor, prime_count};

assert_eq!(factor(&360u64), [(2, 3), (3, 2), (5, 1)]);
assert_eq!(prime_count(1_000_000), 78_498);

Polynomials

The polynomial module works on polynomials of one variable. A polynomial is a slice of its coefficients, in ascending order. evaluate replaces each item of a slice with the polynomial’s value there. It applies Horner’s rule to a block of items at a time, and that loop vectorizes. roots gives the eigenvalues of the companion matrix, from_roots builds a polynomial from its roots, derivative differentiates, and fit fits points by least squares. roots and fit compute in f64 through faer.

use rustymath::polynomial::{derivative, evaluate};

let c = [2.0, -3.0, 1.0];
let mut x = [0.0, 1.5, 3.0];
evaluate(&c, &mut x);
assert_eq!(x, [2.0, -0.25, 2.0]);
assert_eq!(derivative(&c), [-3.0, 2.0]);

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. The rows for sin, cos and tan call checked_sin, checked_cos and checked_tan. To print the table on your own machine, run RUSTYMATH_TIMES=1 cargo test --test speed -- --nocapture.

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.

pow computes e^(y ln |x|) with 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. root divides ln |y| by n the same way, without rounding 1/n first. At f32 and f16, pow, root and hypot compute in f64 and round once.

Licence

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