# multicalc guide
A tour of every public module — the robotics and control layer (estimation, control, kinematics,
motion, spatial) and the calculus, autodiff, and linear-algebra core it is built on: what each does,
where to start, a snippet you can run, the errors it can return, and a link to a full demo. Read it
start to finish if you are new to the crate, or jump to a single module once you know your way
around.
Every operation is generic over the [`Numeric`](#scalars-and-automatic-differentiation) scalar
trait, which is implemented for `f32` and `f64` and defaults to `f64`. The math functions come
from `libm`, so the crate works without `std`. Methods like `f64::sin` need `std`; in a
`no_std` crate, call the `libm` version instead (`libm::sin(x)` in place of `x.sin()`). The
crate re-exports `libm` as `multicalc::libm`.
Every fallible call returns a `Result`, and the error is the module family's own enum; see
[Error handling](#error-handling).
## Contents
- [Importing](#importing)
- [Scalars and automatic differentiation](#scalars-and-automatic-differentiation)
- [Derivatives, Jacobians, and Hessians](#derivatives-jacobians-and-hessians)
- [Integration](#integration)
- [Gaussian quadrature tables](#gaussian-quadrature-tables)
- [Taylor approximation](#taylor-approximation)
- [Linear algebra](#linear-algebra)
- [Least-squares optimization](#least-squares-optimization)
- [Root finding](#root-finding)
- [Vector calculus](#vector-calculus)
- [ODE integrators](#ode-integrators)
- [Discretization](#discretization)
- [Spatial: quaternions and Lie groups](#spatial-quaternions-and-lie-groups)
- [Kinematics](#kinematics)
- [Control](#control)
- [Motion](#motion)
- [Estimation](#estimation)
- [Random](#random)
- [Error handling](#error-handling)
- [Internals](#internals)
## Importing
There is one answer: glob the prelude for the traits and one-call functions, then name the types
you need from the crate root. Every public type lives at `multicalc::Type`, so you never have to
know which file it is declared in. The examples below spell out their imports in full, but
`use multicalc::prelude::*;` covers the traits in all of them.
```rust
use multicalc::prelude::*;
use multicalc::{KalmanFilter, Matrix, Vector};
```
The one exception is a handful of free functions that stay on their own module, because their
names only make sense next to each other — `multicalc::vector_field::curl_3d` reads better than
`multicalc::curl_3d`.
### The easy path and the configurable one
Most calculus work has two ways in, and the guide uses both:
- **The one-call functions** — `derivative`, `second_derivative`, `partial`, `integral`. They need
no imported trait and no configuration, and they use exact automatic differentiation. Reach for
these first.
- **The strategy objects** — `AutoDiffSingle`, `FiniteDifferenceSingle`, `IterativeSingle`,
`GaussianSingle` and their multi-variable siblings. These are how you choose a different method,
a step size, an iteration count, or a derivative order above the second.
Both compute the same answers; the objects just expose the knobs.
## Scalars and automatic differentiation
The scalar number system that every calculus module is generic over: the `Numeric` trait, plus
the forward-mode automatic-differentiation numbers that also implement it.
- `Numeric`: the scalar trait, implemented for `f32` and `f64`.
- `Dual`, `HyperDual`, `Jet<T, N>`: autodiff scalars (dual numbers) carrying exact first,
second, and arbitrary nth-order derivatives (`Dual` is `Jet<T, 2>`).
- `ScalarFn` / `ScalarFnN` / `VectorFn`: function traits whose `eval` is generic over the
scalar, so one formula runs at `f64` or at any autodiff type.
- The `scalar_fn!` / `scalar_fn_vec!` macros build those traits from closure syntax, and `c()`
marks numeric constants inside the body (a bare `2.0 * x` cannot typecheck in a generic body).
One formula, differentiated exactly to any order:
```rust
use multicalc::AutoDiffSingle;
use multicalc::DerivatorSingleVariable;
use multicalc::scalar_fn;
let function = scalar_fn!(|x| x * x * x); // f(x) = x^3, evaluable at any Numeric
let derivator = AutoDiffSingle::default(); // forward-mode autodiff, exact
let point = 2.0;
let first = derivator.differentiate(1, &function, point).unwrap(); // 12.0
let third = derivator.differentiate(3, &function, point).unwrap(); // 6.0
```
Errors: differentiation calls return [`DiffError`](#error-handling) (for example `OrderZero`).
Credits: standard forward-mode dual numbers. Full demo:
[autodiff_scalars.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/autodiff_scalars.rs).
## Derivatives, Jacobians, and Hessians
Derivatives of any order, total and partial — exact through forward-mode autodiff, or by finite
differences for black-box functions — plus Jacobian and Hessian matrices.
- `derivative`, `second_derivative`, `partial`: one-call functions covering the common case. See
[the note on the two paths](#the-easy-path-and-the-configurable-one).
- `AutoDiffSingle` / `AutoDiffMulti`: exact derivatives, to any order.
- `FiniteDifferenceSingle` / `FiniteDifferenceMulti`: for functions you cannot author with
`scalar_fn!`.
- Both implement the `DerivatorSingleVariable` / `DerivatorMultiVariable` traits
(`differentiate`, `first_derivative`, `second_derivative`, `first_partial_derivative`).
- `Jacobian` and `Hessian` build the matrices.
For several variables, the derivative order is just the number of indices you pass:
```rust
use multicalc::AutoDiffMulti;
use multicalc::DerivatorMultiVariable;
use multicalc::scalar_fn;
// g(x, y, z) = y*sin(x) + x*cos(y) + x*y*e^z; order = number of indices passed
let g = scalar_fn!(|v: &[f64; 3]| v[1] * v[0].sin() + v[0] * v[1].cos() + v[0] * v[1] * v[2].exp());
let d = AutoDiffMulti::default();
let point = [1.0, 2.0, 3.0];
let x_index = 0;
let dx = d.first_partial_derivative(&g, x_index, &point).unwrap();
let then_by_y = [0, 1];
let mixed = d.differentiate(&g, &then_by_y, &point).unwrap(); // d(dg/dx)/dy
let twice_by_x_then_y = [0, 0, 1];
let third = d.differentiate(&g, &twice_by_x_then_y, &point).unwrap();
```
Pass a finite-difference derivator (`FiniteDifferenceSingle` / `FiniteDifferenceMulti`) instead
when the function is a black box you cannot author with `scalar_fn!`.
Write a vector-valued function with `scalar_fn_vec!` and its rows differentiate under autodiff
to give the Jacobian; a scalar field gives the Hessian:
```rust
use multicalc::Jacobian;
use multicalc::Hessian;
use multicalc::c;
use multicalc::{scalar_fn, scalar_fn_vec};
// the vector function (x*y*z, x^2 + y^2)
let jacobian: Jacobian = Jacobian::default();
let j = jacobian.evaluate(&f, &jacobian_point).unwrap(); // [[6, 3, 2], [2, 4, 0]]
// g(x, y) = y*sin(x) + 2*x*e^y
let g = scalar_fn!(|v: &[f64; 2]| v[1] * v[0].sin() + c(2.0) * v[0] * v[1].exp());
let hessian_point = [1.0, 2.0];
let hessian: Hessian = Hessian::default();
let h = hessian.evaluate(&g, &hessian_point).unwrap();
```
With the `alloc` feature, `Jacobian::evaluate_on_heap` returns a `Vec<Vec<T>>` for inputs too large
for the stack.
Errors: these calls return [`DiffError`](#error-handling): `OrderZero`, `OrderUnsupported`,
`StepSizeZero` (finite differences), or `IndexOutOfRange`. Full demos:
[differentiation.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/differentiation.rs)
and
[jacobian_hessian.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/jacobian_hessian.rs).
## Integration
Definite integration of any order: iterative Newton-Cotes rules and Gaussian quadrature, over
finite, semi-infinite, and infinite limits.
- `integral`: a one-call function covering the common case. See
[the note on the two paths](#the-easy-path-and-the-configurable-one).
- `IterativeSingle`: Boole (default), Simpson, and Trapezoidal rules; pick the rule and interval
count with `from_parameters`.
- Pairwise summation is the default; chain `.with_kahan_summation()` to opt into Kahan.
- `GaussianSingle`: Gauss-Legendre, Gauss-Hermite, and Gauss-Laguerre. Pass the **bare**
integrand; the weights already carry the weighting factor.
- Both implement the `IntegratorSingleVariable` / `IntegratorMultiVariable` traits
(`integrate`, `single_integral`, `double_integral`, …).
Iterative rules over finite and infinite limits:
```rust
use multicalc::IntegratorSingleVariable;
use multicalc::IterativeSingle;
let integrator = IterativeSingle::default(); // Boole's rule, 120 intervals
let area = integrator.single_integral(&line, &limits).unwrap(); // 4.0
// infinite / semi-infinite limits are supported for decaying integrands
let bell_curve = |x: f64| (-x * x).exp();
let real_line = [f64::NEG_INFINITY, f64::INFINITY];
let bell = integrator.single_integral(&bell_curve, &real_line).unwrap(); // sqrt(pi)
```
Choose the rule and interval count with `from_parameters`:
```rust
use multicalc::{IterativeMethod, IterativeSingle};
let interval_count = 120;
let integrator: IterativeSingle =
IterativeSingle::from_parameters(interval_count, IterativeMethod::Simpsons);
```
Each Gaussian rule integrates over a fixed domain. Pass the bare integrand `f(x)`; the weights
already carry the weighting factor:
```rust
use multicalc::IntegratorSingleVariable;
use multicalc::GaussianSingle;
use multicalc::GaussianQuadratureMethod;
// Gauss-Hermite integrates f(x) * e^(-x^2) over the whole real line.
let node_count = 5;
let hermite = GaussianSingle::from_parameters(node_count, GaussianQuadratureMethod::GaussHermite);
let val = hermite.single_integral(&square, &real_line).unwrap(); // sqrt(pi)/2
```
| Gauss-Legendre | $\int_a^b f(x)\, \mathrm{d}x$ |
| Gauss-Laguerre | $\int_0^\infty f(x)\, e^{-x}\, \mathrm{d}x$ |
| Gauss-Hermite | $\int_{-\infty}^\infty f(x)\, e^{-x^2}\, \mathrm{d}x$ |
Gaussian nodes and weights come from the [quadrature tables](#gaussian-quadrature-tables).
Errors: integration calls return [`IntegrateError`](#error-handling): `IterationsZero`,
`LimitsIllDefined`, `QuadratureOrderOutOfRange`, or `NonFinite`. Full demos:
[iterative_integration.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/iterative_integration.rs)
and
[gaussian_integration.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/gaussian_integration.rs).
## Gaussian quadrature tables
Precomputed Gauss-Legendre, Gauss-Hermite, and Gauss-Laguerre quadrature nodes and weights, up
to order `MAX_ORDER` (30). These back the Gaussian rules in [Integration](#integration); most
users reach them through `GaussianSingle` rather than directly.
- `nodes(method, order)` returns the `(weight, abscissa)` pairs as `&'static [(f64, f64)]`, or
`IntegrateError::QuadratureOrderOutOfRange` if the order is unavailable.
- Per-family data lives in the `legendre`, `hermite`, and `laguerre` submodules.
```rust
use multicalc::gaussian_tables;
use multicalc::GaussianQuadratureMethod;
let order = 5;
let pairs = gaussian_tables::nodes(GaussianQuadratureMethod::GaussHermite, order).unwrap();
for (weight, abscissa) in pairs {
// weight * f(abscissa) is one term of the 5-point Gauss-Hermite sum
}
```
Errors: an out-of-range order returns `IntegrateError::QuadratureOrderOutOfRange` (see
[Error handling](#error-handling)). Credits: generated by
`scripts/build_gaussian_integration_tables.py`.
## Taylor approximation
Local Taylor models of a function around a point (linear and quadratic) with goodness-of-fit
metrics.
- `LinearApproximator`: first-order model.
- `QuadraticApproximator`: same API, also captures curvature.
- `approximate` builds the model; `predict` evaluates it; `prediction_metrics` returns MAE, MSE,
RMSE, R², and adjusted R² against sample points.
- Metrics use pairwise summation by default; chain `.with_kahan_summation()` to opt into Kahan.
```rust
use multicalc::LinearApproximator;
use multicalc::scalar_fn;
let f = scalar_fn!(|v: &[f64; 3]| v[0] + v[1] * v[1] + v[2] * v[2] * v[2]);
let base_point = [1.0, 2.0, 3.0]; // where the model is anchored
let linear: LinearApproximator = LinearApproximator::default();
let model = linear.approximate(&f, &base_point).unwrap();
let nearby = [1.1, 2.1, 3.1];
let y = model.predict(&nearby);
// model.prediction_metrics(&samples, &f) returns RMSE, R^2, and more
```
`QuadraticApproximator` works the same way and captures curvature as well.
Errors: the underlying derivatives return [`DiffError`](#error-handling). Full demo:
[approximation.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/approximation.rs).
## Linear algebra
Fixed-size, stack-allocated `Matrix` and `Vector`. Dimensions are const generics, so a shape
mismatch is a compile error and nothing is heap-allocated. Indexing (`v[i]`, `m[(r, c)]`) is
the ergonomic path and panics on out-of-range like `Vec`. Use `get` / `get_mut` /
`try_row` / `try_column` when you want `Option` instead of a panic (panicking `row` /
`column` were removed).
- `Matrix::lu` → `Lu`: partial-pivoting Doolittle LU; `solve`, `determinant`, `inverse`.
- `Matrix::cholesky` → `Cholesky`: faster path for symmetric positive-definite matrices.
- `PivotedQr`: column-pivoted Householder QR; `solve_least_squares`.
- `Matrix::svd` → `Svd`: one-sided Jacobi SVD; `singular_values`, `condition_number`,
`pseudo_inverse`, minimum-norm `solve`.
Direct linear solves via LU and Cholesky:
```rust
use multicalc::{Matrix, Vector};
// Solve A·x = b.
let a = Matrix::<3, 3>::new([[2.0, 1.0, 1.0], [4.0, 3.0, 3.0], [8.0, 7.0, 9.0]]);
let b = Vector::new([7.0, 19.0, 49.0]);
let x = a.solve(b).unwrap(); // [1, 2, 3]
let lu = a.lu().unwrap();
let det = lu.determinant();
let inv = lu.inverse();
// A symmetric positive-definite matrix has a faster Cholesky path.
let s = Matrix::<2, 2>::new([[4.0, 2.0], [2.0, 3.0]]);
let s_inv = s.cholesky().unwrap().inverse();
```
The singular value decomposition (one-sided Jacobi) gives the pseudo-inverse, minimum-norm
least-squares solve, rank, and condition number for any shape:
```rust
use multicalc::{Matrix, Vector};
// Thin SVD of a tall matrix: A = U · diag(σ) · Vᵀ.
let a = Matrix::<3, 2>::new([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]]);
let svd = a.svd().unwrap();
let sigma = svd.singular_values(); // descending, non-negative
let cond = svd.condition_number(); // σ_max / σ_min
// Moore-Penrose pseudo-inverse: tall, square, or wide (M < N) inputs.
let a_pinv = a.pseudo_inverse().unwrap();
// Minimum-norm least-squares solve of A·x = b, without forming A⁺.
let x = svd.solve(Vector::new([1.0, 2.0, 3.0]));
```
For an overdetermined linear least-squares fit, use the column-pivoted QR directly:
```rust
use multicalc::linear_algebra::PivotedQr;
use multicalc::{Matrix, Vector};
// Least-squares fit of y = a + b*t through (0, 1), (1, 3), (2, 5): a = 1, b = 2.
let a = Matrix::<3, 2>::new([[1.0, 0.0], [1.0, 1.0], [1.0, 2.0]]);
let b = Vector::new([1.0, 3.0, 5.0]);
let x = PivotedQr::decompose(a).unwrap().solve_least_squares(b).unwrap();
```
Errors: factorizations and solves return [`LinalgError`](#error-handling): `Singular`,
`NotPositiveDefinite`, `Underdetermined` (a least-squares system with `M < N`), or `NonFinite`.
Credits: the QR factorization, damped solve, and overflow-safe norm port MINPACK's `qrfac`,
`qrsolv`, and `enorm` (Moré, Garbow, Hillstrom; public domain, netlib). LU and Cholesky follow
the standard Doolittle and Cholesky–Banachiewicz algorithms; the SVD follows Golub & Van Loan,
*Matrix Computations*, and Demmel & Veselić for high relative accuracy. Full demos:
[linear_algebra.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/linear_algebra.rs)
and
[svd.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/svd.rs).
## Least-squares optimization
Nonlinear least-squares solvers. They minimize the sum of squared residuals of a
`scalar_fn_vec!` function, differentiating it under autodiff by default.
- `LevenbergMarquardt`: the robust, damped default.
- `GaussNewton`: the faster undamped variant for well-conditioned problems.
- `minimize` returns a `MinimizationReport` whose `TerminationReason` says which convergence
test stopped the solver.
Write the residuals `model - data` with `scalar_fn_vec!` and the solver differentiates them
under autodiff:
```rust
use multicalc::LevenbergMarquardt;
use multicalc::AutoDiffMulti;
use multicalc::c;
use multicalc::scalar_fn_vec;
// Fit a*e^(b*t) to (0, 100), (1, 50), (2, 25): the minimum is a = 100, b = -ln 2.
c(-50.0) + v[0] * v[1].exp(),
c(-25.0) + v[0] * (c(2.0) * v[1]).exp(),
]);
let report = LevenbergMarquardt::<AutoDiffMulti>::default()
.minimize(&residuals, &[80.0, -0.3])
.unwrap();
// report.solution ~ [100.0, -0.693]; report.termination says which test converged
```
`GaussNewton` has the same API and suits well-conditioned problems where damping is unnecessary.
For a plain linear least-squares fit, use the QR factorization from
[Linear algebra](#linear-algebra) instead.
Errors: the solvers return [`SolveError`](#error-handling): `DidNotConverge { iters, residual }`,
`NonFinite`, or a wrapped `Linalg` / `Diff` error from a failed inner step.
Credits: the Levenberg-Marquardt driver ports MINPACK's `lmder`/`lmpar` (Moré, Garbow,
Hillstrom; public domain, netlib), following Moré (1978), "The Levenberg-Marquardt algorithm:
Implementation and theory", and Nocedal & Wright, *Numerical Optimization*, chapters 4 and 10.
Full demos:
[curve_fit.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/curve_fit.rs)
and
[optimization_solvers.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/optimization_solvers.rs).
## Root finding
Root finders for scalar equations and square systems `F(x) = 0`. Each solver takes an iteration
budget and reports why it stopped as a `RootTermination`.
- `Bisection`: brackets a scalar root and halves the interval; guaranteed to converge within
its budget.
- `Newton`: Newton's method with a derivative from any `Derivator` (exact autodiff by default,
finite differences on request); `with_backtracking(true)` adds a damped line search that
rescues far starts.
- `NewtonSystem`: Newton for square systems `F: Rⁿ → Rⁿ` with the exact Jacobian and an
optional backtracking line search on `‖F‖`.
- The scalar solvers return a `RootReport`; the system solver returns a `RootReportN`.
```rust
use multicalc::{Bisection, Newton, NewtonSystem};
use multicalc::{AutoDiffMulti, AutoDiffSingle};
use multicalc::c;
use multicalc::{scalar_fn, scalar_fn_vec};
// Bracket a scalar root: f(x) = x^2 - 2 on [0, 2].
let f = scalar_fn!(|x| c(-2.0) + x * x);
let bracketed = Bisection::default().solve(&f, 0.0, 2.0).unwrap(); // ~ sqrt(2)
// Newton with exact derivatives; damped Newton adds a backtracking line search.
let quadratic = Newton::<AutoDiffSingle>::default().solve(&f, 2.0).unwrap(); // ~ 1.41421356
let damped = Newton::<AutoDiffSingle>::default()
.with_backtracking(true)
.solve(&f, 2.0)
.unwrap();
// Square system: x^2 + y^2 = 4 and x*y = 1.
// solved.root ~ [1.9319, 0.5176]; solved.termination says which test converged
```
Errors: root finders return [`SolveError`](#error-handling): `DidNotConverge`, `InvalidBracket`
(bisection endpoints that do not enclose a sign change), `NonFinite`, or a wrapped `Linalg` /
`Diff` error.
Credits: textbook bisection and Newton–Raphson iteration; the system step reuses the crate's LU
solve and overflow-safe `enorm` from [Linear algebra](#linear-algebra). Full demo:
[root_finding.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/root_finding.rs).
## Vector calculus
Curl and divergence via autodiff, plus line and flux integrals sampled along a curve.
- `curl_2d` / `curl_3d` and `divergence_2d` / `divergence_3d` take an explicit derivator (pass
`AutoDiffMulti::default()` for exact results) and a `scalar_fn_vec!` field.
- `line_integral_2d` and `flux_integral_2d`, with their 3D and `_custom` forms, sample the field,
so they take plain closures for the field and the parametric curve.
```rust
use multicalc::AutoDiffMulti;
use multicalc::c;
use multicalc::scalar_fn_vec;
use multicalc::vector_field::{curl_2d, divergence_2d, flux_integral_2d, line_integral_2d};
// field (2xy, 3cos y)
let divergence = divergence_2d(AutoDiffMulti::default(), &field, &[1.0, 3.14]).unwrap();
// field (y, -x) along the unit circle (cos t, sin t)
let g: [&dyn Fn(&[f64; 2]) -> f64; 2] = [&(|v: &[f64; 2]| v[1]), &(|v: &[f64; 2]| -v[0])];
let curve: [&dyn Fn(f64) -> f64; 2] = [&(|t: f64| t.cos()), &(|t: f64| t.sin())];
let limit = [0.0, 2.0 * std::f64::consts::PI];
let line = line_integral_2d(&g, &curve, &limit).unwrap(); // -2*pi
let flux = flux_integral_2d(&g, &curve, &limit).unwrap(); // 0
```
The 3D curl is `(dVz/dy - dVy/dz, dVx/dz - dVz/dx, dVy/dx - dVx/dy)`.
Errors: the operators return [`DiffError`](#error-handling) from differentiation, and the
integrals return [`IntegrateError`](#error-handling). Full demo:
[vector_field.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/vector_field.rs).
## ODE integrators
Initial-value solvers for `y' = f(t, y)` systems, generic over the state dimension.
- `Rk4`: fixed-step classical Runge–Kutta. `Rk4::step` advances one step; `Rk4::integrate`
runs a fixed number of steps with a per-step callback.
- `Rk45`: adaptive Dormand–Prince 5(4) with PI step control and cubic-Hermite dense output.
`solve` integrates to a target time, `solve_on_grid` fills requested sample times via dense
output, and `for_each_step` exposes each accepted step. Tolerances are set with `with_rtol`
and `with_atol`.
```rust
use multicalc::{Rk4, Rk45};
use multicalc::Vector;
// Harmonic oscillator y'' = -y as the first-order system [position, velocity].
let timestep = 0.1;
let y1 = Rk4::step(&f, start_time, &y0, timestep); // one fixed step
// Adaptive solve over one full period returns to the start [1, 0].
let one_period = core::f64::consts::TAU;
let yf = Rk45::default().solve(&f, start_time, &y0, one_period).unwrap();
assert!((yf[0] - 1.0).abs() < 1e-6 && yf[1].abs() < 1e-6);
```
Dense output samples a whole grid in one pass, and `for_each_step` lets you track a conserved
quantity as the solver runs:
```rust
use multicalc::Rk45;
use multicalc::Vector;
let times = [0.5, 1.0, 2.0, 3.0];
let mut out = [Vector::<2, f64>::zeros(); 4];
solver.solve_on_grid(&f, start_time, &y0, ×, &mut out).unwrap();
```
Errors: the adaptive solver returns [`IntegrateError`](#error-handling): `StepSizeTooSmall`,
`DidNotConverge { steps }`, or `NonFinite`. Full demo (harmonic oscillator plus an acrobot,
a tumbling quadrotor, and an N-body model):
[ode.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/ode.rs).
## Discretization
Turn a continuous-time linear system into its discrete-time equivalent over a step `dt`.
- `zoh(a, b, dt)`: zero-order-hold discretization of `(A, B)`, returning the discrete `(F, G)`.
- `van_loan(a, qc, dt)`: Van Loan discretization of continuous process noise, returning the
discrete transition and process-noise covariance `(F, Q_d)`.
- `q_discrete_white_noise(dt, var)`: the filterpy-compatible discrete white-noise model.
Because the routines run through the matrix exponential, an autodiff scalar flows straight
through them: a single `Dual` recovers a derivative with respect to a parameter.
```rust
use multicalc::{q_discrete_white_noise, van_loan, zoh};
use multicalc::Matrix;
use multicalc::Dual;
let dt = 0.1;
// Zero-order hold of the double integrator: F = [[1, dt], [0, 1]], G = [[dt^2/2], [dt]].
let a = Matrix::<2, 2>::new([[0.0, 1.0], [0.0, 0.0]]);
let b = Matrix::<2, 1>::new([[0.0], [1.0]]);
let (f, g) = zoh::<2, 1, 3, f64>(a, b, dt).unwrap(); // f[(0, 1)] == dt, g[(1, 0)] == dt
// Van Loan process-noise discretization of continuous white noise on velocity.
let qc = Matrix::<2, 2>::new([[0.0, 0.0], [0.0, 1.0]]);
let (_f, qd) = van_loan::<2, 4, f64>(a, qc, dt).unwrap(); // qd[(1, 1)] == dt, symmetric
// Discrete white-noise model.
let q = q_discrete_white_noise::<2, f64>(dt, 2.0); // q[(1, 1)] == 2*dt^2
// d/dx expm(x·M) at x = 0 equals M, recovered by one Dual through expm.
let m = Matrix::<2, 2>::new([[0.2, 0.5], [-0.1, 0.3]]);
let ad = Matrix::<2, 2, Dual<f64>>::from_fn(|i, j| {
Dual::new(0.0, m[(i, j)])
})
.expm()
.unwrap();
// ad[(0, 1)].deriv == m[(0, 1)]
```
Errors: the matrix-exponential step returns [`LinalgError`](#error-handling) on a non-finite
input. Full demo:
[discretization.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/discretization.rs).
## Spatial: quaternions and Lie groups
Rotations, Lie groups, and rigid-body transforms for 2D and 3D. Fixed-size, stack-allocated, no
panics, and generic over the `Numeric` scalar, so `f32`, `f64`, and the autodiff duals all work.
- `Quaternion`: Hamilton quaternion, stored scalar-first `[w, x, y, z]`: the raw algebra plus
axis-angle / rotation-matrix / ZYX-Euler conversions, `slerp`, and `exp`/`ln`.
- `SO2` / `SE2`: 2D rotation and rigid-body transform.
- `SO3` / `SE3`: 3D rotation (wrapping a unit `Quaternion`, which carries the unit-rotation
invariant) and rigid-body transform.
- `Twist` / `Wrench`: typed spatial velocity and force in the linear-first `[v; ω]` /
`[force; torque]` ordering.
Every group provides `identity`, `compose` (also `*`), `inverse`, `act` on a point, `exp`/`log`,
`hat`/`vee`, `adjoint`, geodesic `interpolate`, and matrix conversions. Conventions: the tangent
ordering is `[v; ω]` (linear part first) for `SE2`/`SE3`; the retract is right-perturbation
`X · exp(ξ)`; angles are radians. `exp`/`log` Taylor-continue near θ = 0 so derivatives stay
finite at rest.
```rust
use multicalc::{SE3, SO3};
use multicalc::Vector;
// A 90° rotation about z, applied to a point.
let quarter_turn_about_z = Vector::new([0.0, 0.0, core::f64::consts::FRAC_PI_2]);
let r = SO3::<f64>::exp(quarter_turn_about_z);
let point = Vector::new([1.0, 0.0, 0.0]);
let p = r.act(point); // ≈ (0, 1, 0)
// A rigid transform: rotate, then translate.
let translation = Vector::new([1.0, 2.0, 3.0]);
let g = SE3::from_parts(r, translation);
let q = g.act(point); // ≈ (1, 3, 3)
// exp/log round trip on the tangent twist [v; ω].
let xi = g.log();
let g2 = SE3::exp(xi);
```
Because everything is generic over the scalar, a derivative with respect to a joint angle or
pose parameter flows through `act`, `compose`, and `exp`/`log` under autodiff. That is what the
inverse-kinematics showcases are built on. Full demo:
[lie_groups.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/lie_groups.rs);
worked application:
[3d_arm_ik.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/showcase/3d_arm_ik.rs).
## Kinematics
Maps between wheel motion and body motion for a differential drive, and pose integration on SE(2).
Fixed-size, no allocation, no panics, and generic over the `Numeric` scalar.
The body motion is deliberately 2-DOF, not 3. A differential drive has exactly two degrees of
freedom `(v, ω)` and exactly two wheels, so the map between them is a bijection and both round trips
are exact identities. There is no lateral term to silently drop.
- `DifferentialDrive`: the geometry, a wheel radius and a track width. Constructing it is the only fallible
operation in the module; with the geometry checked once, every map below is total.
- `WheelVelocities` / `BodyTwist`: motion per second, related by `forward` and `inverse`. A
`BodyTwist` is the se(2) twist a differential drive can realise, with the lateral term dropped.
- `WheelRotations` / `BodyArc`: motion over one tick, related by `forward_arc` and `inverse_arc`.
`WheelRotations` is what an encoder reports; a `BodyArc` is arc length and heading change, the
exponential coordinates of the relative pose.
- `integrate`: advances an `SE2` pose along the exact constant-twist arc.
- `Unicycle`: the same plant as an ODE right-hand side, for `Rk4`/`Rk45`.
- `OdometryStep`: the process model as a `VectorFn`, for autodiff Jacobians.
```rust
use multicalc::kinematics::integrate;
use multicalc::{BodyTwist, DifferentialDrive, WheelVelocities};
use multicalc::Dual;
use multicalc::SE2;
let wheel_radius = 0.036_f64; // 36 mm
let track_width = 0.235; // 235 mm between the wheels
let drive = DifferentialDrive::new(wheel_radius, track_width).unwrap();
// Wheel velocities to a body twist, and back exactly.
let wheel_speeds = WheelVelocities::new(10.0, 10.0); // rad/s on each wheel
let twist = drive.forward(wheel_speeds); // v = 0.36 m/s, ω = 0
let body_motion = BodyTwist::new(0.36, 0.0); // m/s forward, rad/s turn
let wheels = drive.inverse(body_motion); // back to (10, 10)
// The encoder path: distance travelled -> wheel rotation -> body arc -> pose.
let left_travel = 0.01; // metres rolled by each wheel
let right_travel = 0.012;
let rotations = drive.wheel_rotations_from_travel(left_travel, right_travel);
let start = SE2::identity();
let pose = integrate(start, drive.forward_arc(rotations));
// Autodiff straight through an odometry step: d(pose)/d(arc length).
let arc_length = Dual::variable(0.4); // the quantity being differentiated
let turn_rate = Dual::constant(0.3);
let duration = Dual::constant(1.0);
let step = integrate(
SE2::<Dual<f64>>::identity(),
BodyTwist::new(arc_length, turn_rate).integrate_over(duration),
);
let dx_ds = step.translation()[0].deriv;
```
Because `integrate` is built on `SE2::exp`, a straight line (ω = 0) is handled by the same code path
as an arc, with no `1/ω` to blow up: the value and its derivative stay finite at exactly zero
curvature. The arc is exact for a constant twist at any step size, so the modelling error is the
zero-order hold on the wheel velocities rather than integration error.
Full demo:
[kinematics.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/kinematics.rs).
## Control
Feedback controllers and steering laws for a mobile robot: a PID with anti-windup and a filtered
derivative, the pure-pursuit path-following law, and Follow-the-Gap reactive avoidance. Fixed-size,
no allocation, no panics, and generic over the `Numeric` scalar, so the same code runs at `f32` on a
microcontroller.
Angles are radians in the robot body frame, measured from the forward (+x) axis and positive
counter-clockwise. Every controller is configured once, with the configuration validated up front,
and every call after that is total.
- `Pid`: three gains and a fixed timestep. `with_output_limits` clamps the output and stops the
integral winding up against the clamp; `with_derivative_filter` puts a one-pole low-pass on the
derivative term, which is what makes a D gain usable on a noisy measurement.
- `OnePoleLowPass`: the filter on its own, by smoothing coefficient (`new`) or by cutoff frequency
(`from_cutoff`).
- `pure_pursuit_curvature`: the exact `κ = 2·sin(α)/L_d` steering curvature toward a lookahead
point, written in body-frame coordinates. `Curvature::to_body_twist` turns it into a command at a
chosen speed.
- `FollowTheGap`: reactive avoidance over a forward range scan. Const-generic on the beam count,
so the working buffer is stack-allocated and the beam geometry is fixed at compile time.
```rust
use multicalc::{FollowTheGap, Pid, pure_pursuit_curvature};
use multicalc::Vector;
use multicalc::SE2;
// A speed loop: PID on the forward speed, output limited, derivative filtered.
let proportional_gain = 2.0_f64;
let integral_gain = 1.0;
let derivative_gain = 0.05;
let timestep = 0.01;
let lowest_output = -1.0;
let highest_output = 1.0;
let derivative_filter_weight = 0.2;
let mut speed_loop = Pid::new(proportional_gain, integral_gain, derivative_gain, timestep)
.unwrap()
.with_output_limits(lowest_output, highest_output)
.unwrap()
.with_derivative_filter(derivative_filter_weight)
.unwrap();
let setpoint = 0.4; // m/s we want
let measurement = 0.35; // m/s we have
let command = speed_loop.update(setpoint, measurement);
// Steering toward a point 2 m ahead and 1 m to the left: a left turn, so positive curvature.
let pose = SE2::identity();
let target = Vector::new([2.0, 1.0]);
let lookahead_distance = 2.0;
let curvature = pure_pursuit_curvature(pose, target, lookahead_distance).unwrap();
let forward_speed = 0.4;
let twist = curvature.to_body_twist(forward_speed);
// Reactive avoidance over a 31-beam scan.
let field_of_view = 2.0 * core::f64::consts::PI / 3.0; // 120°
let max_range = 4.0;
let robot_radius = 0.5;
let clearance = 0.5; // a gap must beat this to count as free
let cruise_speed = 0.4;
let follower: FollowTheGap<31, f64> =
FollowTheGap::try_new(field_of_view, max_range, robot_radius, clearance, cruise_speed).unwrap();
// A clear scan drives straight ahead at cruise speed.
let goal_angle = 0.0;
let clear_scan = [4.0; 31];
let output = follower.compute(&clear_scan, goal_angle).unwrap();
assert!(output.heading().abs() < 1e-12);
// A wall all round stops, and says why.
let walled_in = [0.2; 31];
let blocked = follower.compute(&walled_in, goal_angle).unwrap();
assert!(blocked.is_blocked());
assert_eq!(blocked.body_twist().linear(), 0.0);
```
`FollowTheGap` makes two passes over the scan. First it cleans the scan: a beam that is non-finite
or non-positive counts as a dropped return and reads as free space at maximum range. Then it finds
every run of consecutive beams above the free-range threshold, throws out any run whose two bounding
returns are closer together than the chassis width, and scores the rest by
`span − goal_bias · |aim − goal_angle|`. The `span` is the run's usable arc, pulled in from each
bounded edge by the angle the robot's half-width covers at that edge's range; `aim` is the goal
angle clamped into that arc. Together they keep the robot's sides clear of the obstacles that form
the gap. It then steers toward the winning run with a yaw rate of `steering_gain · heading` and a
forward speed that scales with how far the path ahead is clear.
Measuring the gap in metres rather than in beams is what makes the width test meaningful: the same
angular gap is wide enough to pass at 4 m but too narrow at 0.4 m, and the law of cosines across the
two bounding returns settles it directly. A run that reaches either end of the field of view has no
bounding return on that side, so it counts as open: the sensor saw nothing out there, and inventing
a wall would stop the robot on no evidence.
It is a purely reactive method: with no map and no memory, it can dither in a three-sided pocket.
When no run is both clear and wide enough, it returns a stopped twist with `is_blocked()` set rather
than inventing a heading. The recovery policy — rotating in place until a gap opens, say — is left to
the caller.
Full demos:
[avoidance.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/avoidance.rs)
and
[2d_localization_obstacle_avoidance.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/showcase/2d_localization_obstacle_avoidance.rs)
(a full lap of a marked course, localizing on a map and fusing odometry, an IMU, and GPS).
## Motion
The path a controller follows. `PolylinePath` is an ordered set of waypoints joined by straight
segments, stored as a fixed array of `MAX_POINTS` with a runtime length, so it is stack-allocated
and needs no heap. It answers the two questions a path-following law asks every tick: where am I on
the path, and what point should I aim at.
- `PolylinePath<MAX_POINTS, DIMENSION, T>`: built with `try_from_points` from a slice, or `new` plus
`push` one waypoint at a time. Duplicate consecutive waypoints are accepted; every query treats a
zero-length segment as contributing no arc length.
- `total_arc_length`: the distance along the whole path, zero for fewer than two waypoints.
- `closest_point`: projects a query point onto the path, returning a `PathProjection` with the
`point` found, its `segment_index`, its `arc_length` from the start, and the `distance` from the
query point to it.
- `lookahead_point`: the point a given distance further along from a given arc length — the aim
point for `pure_pursuit_curvature`.
- `EndOfPath`: what a lookahead does when it runs off the end, `Stop` (clamp to the last waypoint,
the default) or `Loop` (wrap to the start). Set it with `with_end_of_path`.
```rust
use multicalc::{EndOfPath, PolylinePath};
use multicalc::Vector;
// An L-shaped path: three units east, then four units north.
let path: PolylinePath<3, 2, f64> = PolylinePath::try_from_points(&[
Vector::new([0.0, 0.0]),
Vector::new([3.0, 0.0]),
Vector::new([3.0, 4.0]),
])
.unwrap()
.with_end_of_path(EndOfPath::Loop);
let total = path.total_arc_length(); // 7.0
// Where is a robot sitting off to the side of the first leg?
let robot_position = Vector::new([2.0, 0.5]);
let here = path.closest_point(robot_position).unwrap();
let on_path = here.point(); // (2.0, 0.0)
let travelled = here.arc_length(); // 2.0
let cross_track = here.distance(); // 0.5
// Aim one unit further along than that.
let lookahead_distance = 1.0;
let aim = path.lookahead_point(travelled, lookahead_distance).unwrap(); // (3.0, 0.0)
```
`try_from_points` and `push` return [`MotionError::CapacityExceeded`] if there is no room for the
waypoint and [`MotionError::NonFinite`] if any coordinate is not finite. `closest_point` and
`lookahead_point` return [`MotionError::PathTooShort`] on an empty path.
Demo: `2d_localization_obstacle_avoidance` drives a lap of this kind of path under pure pursuit.
## Estimation
State estimation from noisy measurements. `KalmanFilter` is the linear filter: `predict` rolls the
state forward through a matrix model and grows the covariance by the process noise; `update` folds in
a measurement and shrinks it. Fixed-size, no allocation, and generic over the `Numeric` scalar, so a
`Dual` state differentiates the whole filter.
- `KalmanFilter<STATE_DIMENSION, MEASUREMENT_DIMENSION, T>`: built from an initial estimate and a
`KalmanModel`.
- `KalmanModel<STATE_DIMENSION, MEASUREMENT_DIMENSION, T>`: the four matrices that describe what the
filter is tracking — transition, measurement model, process noise, measurement noise. Three of
them are the same shape, so naming each field is what stops a swapped pair compiling.
- `predict` / `predict_with_control`: the time step, undriven or with a `control_model ·
control_input` term. `CONTROL_DIMENSION` lives on the method, so undriven users never meet it.
- `update`: the measurement step. The only fallible operation in the module.
- `CovarianceUpdate`: `Joseph` (the default) or `Naive`.
- `innovation` / `innovation_covariance` / `normalized_innovation_squared`: for measurement gating.
- The setters (`set_state_transition`, `set_process_noise`, …) cover the time-varying case, where a
changing timestep changes the model between steps.
```rust
use multicalc::{KalmanFilter, KalmanModel};
use multicalc::{Matrix, Vector};
// Constant velocity: position integrates velocity over a 1 s step; position is measured.
let initial_state = Vector::new([0.0, 0.0]); // [position, velocity]
let initial_covariance = Matrix::new([[1.0, 0.0], [0.0, 1.0]]);
let model = KalmanModel {
state_transition: Matrix::new([[1.0, 1.0], [0.0, 1.0]]),
measurement_model: Matrix::new([[1.0, 0.0]]), // position only
process_noise: Matrix::new([[0.01, 0.0], [0.0, 0.01]]),
measurement_noise: Matrix::new([[0.1]]),
};
let mut filter = KalmanFilter::new(initial_state, initial_covariance, model);
filter.predict();
let measurement = Vector::new([1.0]);
filter.update(measurement).unwrap();
let position = filter.state()[0];
// Gate an outlier before folding it in.
filter.predict();
let outlier = Vector::new([1.9]);
filter.update(outlier).unwrap();
let gate = filter.normalized_innovation_squared().unwrap();
```
The covariance update uses Joseph form by default — `(I − K·H)·P·(I − K·H)ᵀ + K·R·Kᵀ` — which stays
symmetric and positive definite by construction, while the naive `(I − K·H)·P` loses symmetry as
rounding builds up. Joseph form is not a guarantee at every scale: over roughly 10⁷ single-precision
updates it drifts too, and the fix there is to symmetrize and clamp the covariance.
`update` returns `EstimationError::NonFinite` for a non-finite measurement or innovation covariance,
and `EstimationError::NotPositiveDefinite` when the innovation covariance cannot be factorized — the
gain is undefined. `predict` is a cheap element-wise path and propagates non-finite values silently.
`ExtendedKalmanFilter<STATE_DIMENSION, MEASUREMENT_DIMENSION, T>` takes the process and measurement
models as functions rather than matrices — any `VectorFn` — and re-linearizes them at the current
estimate on every step. **The Jacobians come from automatic differentiation: write the model once
and its partial derivatives are exact, with no hand-derived Jacobians anywhere** — the classic source
of silent estimator bugs.
- `new` / `from_derivator`: the autodiff default, or an explicit differentiation backend (e.g.
`FiniteDifferenceMulti`).
- `predict(&process_model)` / `update(&measurement_model, measurement)`: the models are passed per
step, not stored, so the type stays `ExtendedKalmanFilter<3, 2>`. A control input or a changing
timestep lives in the model as a field the caller sets between steps — there is no
`predict_with_control`. Unlike the linear filter, `predict` here evaluates and differentiates a
model, so it returns a `Result`.
- `update_with_residual(&measurement_model, residual)`: `update` with a caller-formed residual, for
when a measurement component is an angle — plain subtraction is wrong across the ±π wrap, and only
the caller knows which components are angular.
- `CovarianceUpdate`, the accessors, and `normalized_innovation_squared` are shared with the linear
filter. `predict` and `update` also return `EstimationError::Diff` if a Jacobian step fails —
reachable only with a finite-difference backend, as the autodiff default cannot.
```rust
use multicalc::ExtendedKalmanFilter;
use multicalc::{Matrix, Vector};
use multicalc::{Numeric, VectorFn};
// Range to a landmark at (3, 4): nonlinear in the pose, so the linear filter cannot take it.
struct RangeToLandmark;
impl VectorFn<2, 1> for RangeToLandmark {
fn eval<S: Numeric>(&self, state: &[S; 2]) -> [S; 1] {
let to_landmark_x = S::from_f64(3.0) - state[0];
let to_landmark_y = S::from_f64(4.0) - state[1];
[(to_landmark_x * to_landmark_x + to_landmark_y * to_landmark_y).sqrt()]
}
}
// A stationary target: the pose carries over unchanged.
struct Stationary;
impl VectorFn<2, 2> for Stationary {
fn eval<S: Numeric>(&self, state: &[S; 2]) -> [S; 2] {
[state[0], state[1]]
}
}
let mut filter = ExtendedKalmanFilter::<2, 1>::new(
Vector::new([0.0, 0.0]), // initial pose, 5.0 from the landmark
Matrix::new([[1.0, 0.0], [0.0, 1.0]]), // initial covariance
Matrix::new([[0.01, 0.0], [0.0, 0.01]]), // process noise
Matrix::new([[0.1]]), // measurement noise
);
filter.predict(&Stationary).unwrap();
filter.update(&RangeToLandmark, Vector::new([5.5])).unwrap();
```
### Particle filter
`ParticleFilter<STATE_DIMENSION, MEASUREMENT_DIMENSION, T, R>` carries a cloud of weighted state
samples instead of a single Gaussian, so it can track a belief the Kalman filters cannot represent —
strongly nonlinear, non-Gaussian, or with several peaks at once (a robot that could be in one of two
corridors). It is the tool to reach for when a single-Gaussian belief is the thing that breaks, and
the price is running hundreds to thousands of samples every step.
- `new(particle_count, initial_mean, initial_covariance, process_noise, seed)`: samples the starting
cloud from the given Gaussian, with a seeded built-in `Pcg32`. `from_random` takes any
`RandomSource` instead. `particle_count` must be at least one.
- `predict(&process_model)`: pushes every sample through the model — any `VectorFn` — and adds a draw
of process noise. `update(&measurement_model, &likelihood, measurement)`: reweights each sample by
how well its predicted measurement matches, normalizes, and resamples if the cloud has degenerated.
- `Likelihood` scores a sample as a log-weight; `GaussianLikelihood::new(measurement_noise)` is the
batteries-included default for additive Gaussian noise. Write your own for anything else.
- `ResamplingScheme`: `Systematic` (the default), `Stratified`, `Multinomial`, or `Residual`. Set it
with `with_resampling`; tune when it fires with `with_resample_threshold`, and add post-resample
jitter with `with_roughening`.
- `mean` (the usual estimate), `maximum_a_posteriori_state` (the single heaviest sample, for when the
belief has several peaks and the mean falls between them), `effective_sample_size`, `particles`,
and `weights`.
```rust
# use multicalc::{GaussianLikelihood, ParticleFilter};
# use multicalc::{Matrix, Vector};
# use multicalc::{Numeric, VectorFn};
// A stationary 2-D point, measured directly with a little noise.
struct Stationary;
impl VectorFn<2, 2> for Stationary {
fn eval<S: Numeric>(&self, state: &[S; 2]) -> [S; 2] {
[state[0], state[1]]
}
}
let particle_count = 1000;
let initial_mean = Vector::new([0.0, 0.0]);
let initial_covariance = Matrix::new([[1.0, 0.0], [0.0, 1.0]]);
let process_noise = Matrix::new([[0.01, 0.0], [0.0, 0.01]]);
let seed = 7;
let mut filter = ParticleFilter::<2, 2>::new(
particle_count,
initial_mean,
initial_covariance,
process_noise,
seed,
)
.unwrap();
let measurement_noise = Matrix::new([[0.05, 0.0], [0.0, 0.05]]);
let sensor = GaussianLikelihood::new(measurement_noise).unwrap();
let measurement = Vector::new([1.0, 2.0]);
for _ in 0..20 {
filter.predict(&Stationary).unwrap();
filter.update(&Stationary, &sensor, measurement).unwrap();
}
assert!((filter.mean()[0] - 1.0).abs() < 0.2);
```
The particle filter is heap-backed, so it is behind the `alloc` feature and the bare-metal build does
not compile it. Its `update` returns `EstimationError::NonFinite` for a non-finite measurement and
`EstimationError::WeightsDegenerate` when no sample can explain the measurement. `GaussianLikelihood`
forms the mismatch by plain subtraction, so a measurement with an angular component needs a custom
`Likelihood` that folds the angle into a ±π band first — the same wrap the extended filter's
`update_with_residual` exists for.
## Random
A seedable generator that works without an operating system, so the stochastic parts of the library
run on bare metal. The particle filter uses it internally; it is public because process noise,
sensor models, and Monte-Carlo checks need the same thing.
- `RandomSource`: the trait a generator implements. `next_u32` is the only required method; the
trait supplies `next_u64`, `next_unit_f64` (uniform in `[0, 1)` with 53 bits of precision), and
`standard_normal` (mean 0, standard deviation 1) on top of it. Implement it to plug in a hardware
generator or your own algorithm.
- `Pcg32`: the built-in generator (PCG-XSH-RR, 32-bit output). `new(seed)` uses the default stream;
`with_stream(seed, stream)` picks another, so independent filters draw independent sequences from
the same seed. Deterministic — the same seed reproduces the same run exactly, which is what makes
a seeded simulation repeatable. Not for cryptography.
```rust
use multicalc::{Pcg32, RandomSource};
let seed = 20260722;
let mut generator = Pcg32::new(seed);
let uniform = generator.next_unit_f64(); // in [0, 1)
let noise = generator.standard_normal(); // mean 0, standard deviation 1
// The same seed replays the same sequence.
let mut replay = Pcg32::new(seed);
assert_eq!(replay.next_unit_f64(), uniform);
// A second stream from the same seed draws an independent sequence.
let stream = 1;
let mut other = Pcg32::with_stream(seed, stream);
let independent = other.standard_normal();
```
Demo: `2d_localization_obstacle_avoidance` seeds every noise source from one number, so the whole
run repeats exactly.
## Error handling
Each module family returns its own error enum. All six convert into the `CalcError` umbrella
through `From`, so a caller that spans families can hold a single type. Every enum is
`#[non_exhaustive]` and `Copy`, and implements `Display` and `core::error::Error`.
| `LinalgError` | [Linear algebra](#linear-algebra), [Discretization](#discretization) | `Singular`, `NotPositiveDefinite`, `Underdetermined`, `NonFinite` |
| `DiffError` | [Derivatives](#derivatives-jacobians-and-hessians), [Approximation](#taylor-approximation), [Vector calculus](#vector-calculus) | `OrderZero`, `OrderUnsupported`, `StepSizeZero`, `IndexOutOfRange`, `EmptyFunctionSet` |
| `IntegrateError` | [Integration](#integration), [Gaussian tables](#gaussian-quadrature-tables), [ODE](#ode-integrators) | `IterationsZero`, `LimitsIllDefined`, `QuadratureOrderOutOfRange`, `StepSizeTooSmall`, `DidNotConverge { steps }`, `NonFinite` |
| `SolveError` | [Optimization](#least-squares-optimization), [Root finding](#root-finding) | `DidNotConverge { iters, residual }`, `NonFinite`, `InvalidBracket`, `Linalg(LinalgError)`, `Diff(DiffError)` |
| `KinematicsError` | [Kinematics](#kinematics) | `NonPositiveParameter`, `NonFinite` |
| `EstimationError` | [Estimation](#estimation) | `NotPositiveDefinite`, `NonFinite`, `Diff(DiffError)`, `WeightsDegenerate` |
| `CalcError` | umbrella | `Linalg`, `Solve`, `Integrate`, `Differentiate`, `Kinematics`, `Estimation` |
`SolveError` wraps `LinalgError` and `DiffError` (a solver step can fail in either), and both
are reachable through `core::error::Error::source`. Convert up to the umbrella with `?` or
`.into()`:
```rust
use multicalc::{CalcError, Matrix, Vector};
// One return type covers a function that mixes modules: each `?` converts the module's own
// error into the umbrella on its way out.
fn solve() -> Result<(), CalcError> {
let a = Matrix::new([[2.0, 1.0, 1.0], [4.0, 3.0, 3.0], [8.0, 7.0, 9.0]]);
let b = Vector::new([7.0, 19.0, 49.0]);
let x = a.lu()?.solve(b); // LinalgError -> CalcError
assert!((a * x - b).norm() < 1e-9);
// A singular matrix returns `LinalgError::Singular` here rather than panicking.
let singular = Matrix::<3, 3>::zeros();
assert!(singular.lu().is_err());
Ok(())
}
# solve().unwrap();
```
This is the shape the converted demos use — see
[linear_algebra.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/linear_algebra.rs),
[root_finding.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/root_finding.rs),
and
[estimation.rs](https://github.com/kmolan/multicalc-rust/blob/main/demos/examples/basics/estimation.rs),
each of which returns `Result<(), CalcError>` from `main` and propagates with `?`.
## Internals
`utils` holds crate-internal numeric helpers — `pub(crate)`, not part of the public API. The main
one is the blocked pairwise summation used for long running sums, where rounding error grows like
`O(log n · eps)` instead of the naive `O(n · eps)`.