rustymath

Rustymath provides fast elementary and special functions, statistics, searching and sorting, random sampling, linear algebra, optimization, polynomials and exact arithmetic in Rust and Python. It is the numerical core of BasedPL.

The Python package needs no NumPy. Its eager operations use the same broadcasting, selection and tolerant-comparison rules as the shared Rust core. For richer array notation, operators and fusion, use BasedPL.

Float functions target four units in the last place (ulps). Rustymath prioritizes speed over last-bit agreement and reproducibility. These functions document wider bounds:

• lgamma near its zeros at 1 and 2 is within 4 ulps of 2, not of the result. digamma near its zeros is within 4 ulps of 1, and below -20 within 4 ulps of ln(1 - x).
• binom and binom_dn of a non-whole n with a whole k up to 64 are within 2k ulps.
• binom_dk has no relative bound near n = 2k, where its two ψ terms cancel.
• Complex functions count ulps of the result’s magnitude, so a much smaller part can carry most of the error. Complex powers, Γ and binomial coefficients allow ulps of their magnitude times one more than the exponent’s magnitude. exp gives an exact direction at the float nearest a whole number of quarter turns, so e^(iπ/2) is exactly i.

Install

pip install rustymath

Rust users can add the crate to Cargo.toml. See Using rustymath from Rust for scalar and slice forms, build settings, accuracy and speed.

import rustymath as rm

Arrays and indexing

Functions take Python numbers, nested lists, rustymath arrays or buffer exporters such as NumPy arrays. Broadcasting aligns leading axes. A two-item vector adds one value to each row of this matrix:

rm.add([[1, 2, 3], [4, 5, 6]], [10, 20])

Array([[11, 12, 13], [24, 25, 26]], dtype=int64)

Arrays and indexing covers construction, dtypes, selection, reshaping, axis names, keys and interoperability.

Maths and statistics

Real and complex functions operate on whole arrays. Float comparisons use a relative tolerance by default:

rm.exp([0., 1, 2]), rm.equal(0.1 + 0.2, 0.3)

(Array([1.0, 2.7182818, 7.3890561], dtype=float64), True)

Maths and statistics covers elementary and special functions, reductions, scans, statistics and tolerant comparison.

Searching, sorting and selection

unique retains first-occurrence order and reports each value’s position, class number and count:

rm.unique([3, 1, 3, 2])

Unique(values=Array([3, 1, 2], dtype=int64), firsts=Array([0, 1, 3], dtype=int64), numbers=Array([0, 1, 0, 2], dtype=int64), counts=Array([2, 1, 1], dtype=int64))

Searching, sorting and selection distinguishes tolerant matching from exact ordering.

Random numbers and distributions

Distributions provide density, cumulative probability, quantiles and seeded sampling:

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

(Array([0.024997895, 0.5, 0.9750021], dtype=float64),
 Array([0.11623959, -0.78285485, -2.3346353], dtype=float64))

Random numbers and distributions covers sampling and probability functions.

Linear algebra

Solve real or complex systems, including least-norm solutions for singular or non-square matrices:

rm.linalg.solve([[4., 1], [1, 3]], [1., 2])

Array([0.090909091, 0.63636364], dtype=float64)

Linear algebra covers matrix products, reusable factorizations, decompositions and least squares.

Optimization and iterative solves

Minimize functions, fit residuals or find roots. Here Brent’s method finds where cos x = x:

rm.optimize.root_scalar(lambda x: rm.cos(x) - x, [0, 1]).root

0.7390851350713886

Optimization and iterative solves also covers matrix-free linear systems.

Polynomials, exact arithmetic and number theory

Polynomial coefficients start from the constant term. Factorization accepts Python integers of any size:

rm.polynomial.polyroots([-6, 11, -6, 1]), rm.number_theory.factor(2**64 + 1)

(Array([(1+0j), (2-0j), (3+0j)], dtype=complex128),
 [(274177, 1), (67280421310721, 1)])

Polynomials, exact arithmetic and number theory includes rational matrix solves and calculations with fractions.Fraction.

Accuracy and speed

The Rust page’s Accuracy and speed section compares each function’s speed with Rust’s standard library. It also explains how the 4-ulp bound holds for f16, with and without fused multiply-add, and for pow and root.

Licence

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