precise-numpy

Interval arithmetic with guaranteed error bounds for Python, powered by Rust.

Quickstart · API · Benchmarks · Correctness · Changelog

What You Can Verify

Need guaranteed error bounds on numerical computing? Browse the capability grid below — every operation returns a rigorous (midpoint, radius) enclosure over the true mathematical value.

Stack

precise-numpy vs NumPy

Installation

pip install precise-numpy

Requires Python >= 3.10. Build from source:

git clone https://github.com/kuntal-devrat/precise-numpy.git
cd precise-numpy
maturin develop --release

Quick Start

import precise_numpy as pnp
import numpy as np

# Create interval arrays with error bounds
a = pnp.array([0.1, 0.2, 0.3], error=1e-4)
b = pnp.array([0.4, 0.5, 0.6], error=1e-4)

# Arithmetic propagates both midpoints and radii
c = a + b
print(c.radii())            # [2e-4, 2e-4, 2e-4]
print(c.max_radius())       # 2e-4

# Reductions return (midpoint, radius) tuples
mid, err = c.sum()
print(f"sum = {mid} +/- {err}")

# NumPy ufuncs work via NEP-13
r = np.exp(pnp.array([0.0, 1.0]))
print(r.values())  # [1.0, e]

# np.linalg works via NEP-18
A = pnp.array([[1.0, 2.0], [3.0, 4.0]])
print(np.linalg.norm(A))           # Frobenius norm: (mid, rad)
print(np.linalg.eigvals(A))        # eigenvalue intervals
print(np.linalg.svd(A))            # U, singular values, Vt (all interval)

API Reference

Array creation

array(values, error=0.0), asarray, zeros(shape), ones, empty, full(shape, value, error=0.0), eye(n, m=None, k=0), identity(n), diag(v), arange(start, stop, step=1.0), linspace(start, stop, num=50, endpoint=True), from_raw_parts(midpoints, radii, shape).

empty(shape) zeros memory on purpose — uninitialised radii would violate the rigorous enclosure guarantee.

to_numpy(a, dtype=None), from_numpy(x, error=scalar_or_array), astype(a, dtype), save_npy(path, arr), load_npy(path).

Properties & conversions

shape, ndim, size, dtype ("interval64"), itemsize (16), nbytes, strides, t, values(), radii(), tolist(), get(i), item(i), midpoint(i), radius(i), copy(), flatten(), ravel(), reshape(*shape), transpose(), sort(), argsort().

__array__(dtype=None) (NEP-18) — returns a plain float64 ndarray of midpoints so np.asarray(pnp_array) and np + pnp_array work naturally.

Indexing

Integer (with negatives), slices (a[1:4], a[::-1]), fancy integer lists (a[[0, 2]]), boolean masks (a[a > 2.0], a[BoolArray]), and 2D tuples (m[1, :], m[:, 0]). Assignment works for integers, slices, and lists.

Arithmetic & math

+ - * / ** @, in-place variants, scalar operands on either side, NumPy broadcasting. sin, cos, tan, exp, ln, log2, log10, sqrt, abs, floor, ceil, trunc, round, clip(a_min, a_max), sign, nan_to_num, power(x, y), maximum(x, y), minimum(x, y).

Comparisons

==, != (overlap), <, <=, >, >= (strict ordering) return BoolArray values with & | ^ ~, any(), all(), sum(), count_nonzero(), tolist().

Reductions

sum, mean, prod, var, std, min, max, argmin, argmax, all, any — each with an optional axis; plus min_val(), max_val(), cumsum(axis=None), norm(ord=None, axis=None) with numpy-compatible ord in {None/2/'fro', 1, -1, inf, -inf, 'nuc'} and int axis.

Stacking & selection

concatenate(arrays, axis=0), stack, vstack, hstack, split(a, sections_or_indices, axis=0), where(condition, x, y), nonzero(a) (tuple of index lists per axis).

Linear algebra (precise_numpy.linalg)

det, inv, solve(A, b), lstsq(A, b), pinv, eig, svd, norm, cholesky, matrix_power, matrix_rank, cond.

Rigor notes: eig_symmetric (Jacobi) — rigorous eigenvalue and eigenvector enclosures via residual bounds + Davis–Kahan. eig_general (Hessenberg QR + Schur back-substitution) — rigorous eigenvalue radii via Bauer–Fike; eigenvector radii via perturbation bound (bounded by 1.0 at defective limit; inf for fully defective matrices). svd — rigorous singular-value radii (cluster hull for degenerate clusters) + Davis–Kahan singular-vector radii. pinv — interval reciprocals of singular values, sound for all conditioning levels. solve — rigorous LU with partial pivoting + forward/back substitution.

Random (precise_numpy.random)

seed, rand(*size), random_sample(*size), random(*size), randn(*size), randint(low, high, *size), uniform(low, high, *size), normal(loc=0.0, scale=1.0, *size). Deterministic xoshiro256**; randint scalar returns a Python int.

Error-Bound Example: AI Quantization Audit

import precise_numpy as pnp

# Wrap input data with interval error bounds
X = pnp.array(floats, error=1e-4).reshape([128, 64])
W = pnp.array(weights).reshape([64, 64])

# Propagate error through a linear layer
Q = X.matmul(W)

print("Max Radius Error:", Q.max_radius())
print("Max Relative Error:", Q.max_relative_error())

Semantics

• Values: get(i), a[0], sum(), etc. return a (midpoint, radius) tuple. To get raw float arrays use a.values() / a.radii().
• == / != mean interval overlap; <, <=, >, >= use strict endpoint ordering (a.hi < b.lo, etc.).
• Division by zero raises a RuntimeWarning and produces the entire real line (radius inf) for the affected element.
• len(a) is a.shape[0], matching NumPy.
• Reductions (sum, mean, ...) support an optional axis argument; scalar results come back as (midpoint, radius) tuples, axis results as IntervalArrays.
• dtype is "interval64", itemsize is 16 (two f64s per element).
• Rounding (round) uses round-half-to-even, like Python and NumPy.
• reshape(-1, ...) infers the missing dimension; transpose/.t support 1D/2D.
• NEP-13 ufunc: np.add, np.multiply, np.negative, np.absolute, np.exp, np.log, np.sqrt, np.sin, np.cos, np.tan, np.isnan, np.isinf, np.isfinite, np.rint, np.sign, np.matmul all work on IntervalArray.
• NEP-18 function: np.linalg.norm, np.linalg.inv, np.linalg.solve, np.linalg.svd, np.linalg.eig, np.linalg.eigvals, np.linalg.pinv, np.sum, np.mean, np.std, np.var, np.max, np.min, np.all, np.any, np.concatenate, np.stack, etc.

Status

Testing

make check          # fmt + clippy + ruff + all tests
make test           # Rust + Python tests
cargo test          # Rust unit tests
python -m unittest tests/python/test_api.py   # Python integration tests

CI runs the full lint + test matrix on Linux, macOS and Windows for Python 3.10–3.13 on every push and pull request. Tagged v* releases build wheels for all platforms and publish to PyPI.

License

MIT License. See LICENSE for details.