Symplex
Symbolic mathematics for Rust.
Pre-release. The API is unstable. 0.3 adds no signature-breaking changes over 0.2, but a few results are normalised differently (see Migrating from 0.2); 0.2 contained breaking changes from 0.1 (see Migrating from 0.1). Feedback welcome.
Contributing? See CONTRIBUTING.md for architecture, conventions, and how to get started. The full changelog is in CHANGELOG.md; the user guide is The Symplex Book.
What This Is
symplex is a symbolic computation library. It manipulates mathematical expressions exactly — using arbitrary-precision rational arithmetic, not floating-point — and can differentiate, integrate, sum, solve equations and systems, simplify, transform, view expressions as polynomials with symbolic coefficients, solve linear programs with exact certificates, compute integer matrix normal forms, and generate optimized Rust or C code from symbolic results.
It is designed for Rust developers working in robotics, control systems, physics simulation, signal processing, or anywhere that symbolic math feeds into numerical code.
Quick Example
use *;
cargo add symplex
When to Use This
- You need symbolic differentiation, integration, summation, or equation solving and want to stay in Rust.
- You are generating numerical code from symbolic derivations — Jacobians, transfer functions, filter coefficients, control laws — as Rust or C99.
- You need compile-time dimensional analysis for physical quantities.
- You need thread-safe symbolic computation without a GIL or global interpreter lock.
- You want exact rational arithmetic (
1/3stays as1/3, not0.33333...).
When Not to Use This
- You need a mature CAS with decades of community validation — use SymPy. It has broader coverage, more special functions, and a much larger test corpus.
- You need geometry, statistics, tensor algebra, or PDE solving — these are not available.
- You need interactive notebook-style exploration — symplex is a library, not an application. (Though see
cargo run --example replfor a basic REPL.) - You need results verified against extensive known-answer databases — symplex has ~11,000 tests including SymPy cross-validation fixtures, but SymPy has orders of magnitude more coverage.
- You need large-scale or sparse numerical optimisation — the exact simplex is dense and
O(m·n)big-rational operations per pivot (hundreds of rows, not hundreds of thousands), and thef64routines are the classic derivative-free methods, not a replacement for a dedicated optimisation library.
What You Can Do
Calculus
Differentiation handles the chain rule, product rule, all elementary functions, and the special functions (Bessel, orthogonal polynomials, digamma → polygamma). Indefinite integration uses 15+ strategies including by-parts, u-substitution, partial fractions, trig substitution, the Risch algorithm, Lazard–Rioboo–Trager log-to-real conversion, and heuristic integration. Radical coefficients (e.g., √5 from cyclotomic denominators) are handled exactly via algebraic number field arithmetic.
let ctx = new;
syms!;
expr!.diff; // 2*x*cos(x^2)
expr!.integrate; // x*exp(x) - exp(x)
expr!.limit; // 1 (Gruntz algorithm)
expr!.series; // 1 + x + x^2/2 + x^3/6 + x^4/24
// One-sided limits; the two-sided limit stays a `Limit` node when they disagree
.limit_right; // oo
.limit_left; // -oo
Definite, Improper and Numeric Integration
integrate_definite locates interior singularities, treats infinite bounds and endpoint singularities as improper integrals via one-sided limits, resolves Abs/Heaviside/DiracDelta/Piecewise integrands, and consults a table of ~30 classical improper integrals (with symbolic parameters under assumptions). Divergence is reported, never hidden.
let ctx = new;
syms!;
let = ;
x.powi.integrate_definite; // 1/3
.exp.integrate_definite; // 1
.integrate_definite; // 1/2*pi
x.ln.integrate_definite; // -1
x.abs.integrate_definite; // 13/2
// ∫₋₁¹ dx/x² diverges: F(1) − F(−1) = −2 would be wrong, so it is an error
let r = x.powi.try_integrate_definite;
assert!;
// Adaptive Gauss–Kronrod (G7/K15) quadrature when there is no closed form
let v = x.powi.exp.integrate_numeric.unwrap; // 1.4626517459…
// Residues at poles of any order, and at infinity
let z = ctx.symbol;
.residue; // 1/2
.residue_at_infinity; // 0
Summation, Products and Series
let ctx = new;
syms!;
let n = ctx.symbol_with;
let = ;
k.powi.summation; // 1/6*n^6 + 1/2*n^5 + 5/12*n^4 - 1/12*n^2
.summation; // 2^(n + 1)*(n - 1) + 2 (Gosper)
n.binomial.summation; // 2^n
k.powi.summation; // 1/6*pi^2
k.powi.summation; // zeta(3) (odd p: symbolic)
.summation; // exp(x)
.product_over; // 1/2
.is_convergent; // Some(false)
// Formal power series: lazy exact coefficients and closed-form general terms
let s = x.sin.fps_maclaurin;
s.coefficient; // -1/1551118753287382280224243016469303211063259720016986112000000000000
s.general_term; // Some(sin(1/2*k*pi)/k!)
s.reversion.unwrap.coefficients; // asin: [0, 1, 0, 1/6, 0, 3/40]
Complex Analysis and Special Functions
re, im, conjugate, arg are honest about unknown realness: with no assumption on z, z.re() is the unevaluated re(z).
let ctx = new;
let z = ctx.symbol;
let x = ctx.symbol_with;
let y = ctx.symbol_with;
let i = ctx.i_unit;
let w = &x + &i * &y;
w.conjugate; // x - y*I
w.abs_squared; // x^2 + y^2
w.exp.as_real_imag; // (cos(y)*exp(x), sin(y)*exp(x))
z.re; // re(z) — not assumed real
z.exp.re; // cos(im(z))*exp(re(z))
.eval; // zoo (complex infinity)
// New constants and special functions with exact values and evalf
ctx.int.digamma.eval; // -EulerGamma
ctx.int.zeta.eval; // 1/90*pi^4
ctx.int.polygamma.eval; // 1/6*pi^2
ctx.infinity.si.eval; // 1/2*pi
ctx.catalan.eval_decimal.unwrap; // 0.915965594177219015054603514932
Also: Gamma, log-gamma, erf/erfc, Beta, Lambert W, Bessel J/Y/I/K, Legendre/Chebyshev/Hermite/Laguerre polynomials, Si/Ci/Ei/li, Kronecker delta — all with arbitrary-precision evaluation.
Algebra and Factoring
Polynomial operations work over ℚ using arbitrary-precision rational arithmetic. Univariate factoring over ℤ uses Berlekamp–Zassenhaus (any degree); multivariate factoring uses Kronecker substitution. Gröbner bases use Buchberger's algorithm with FGLM order conversion.
let ctx = new;
syms!;
expr!.factor; // (x - 1)*(x + 1)*(x^2 + x + 1)*(x^2 + 1)*(x^2 - x + 1)*(x^4 - x^2 + 1)
expr!.factor_all; // (x + 1)*(x + y)*(x - y)
expr!.expand; // x^3 + 3*x^2 + 3*x + 1
expr!.cancel; // x + 1
// Polynomial algebra on Ex: resultant, discriminant, division, gcdex, real-root isolation, …
expr!.discriminant; // Some(4)
expr!.count_real_roots; // Some(1)
expr!.is_irreducible; // Some(true)
Polynomials as Data and Rational Normal Forms
Poly (0.3) views an expression as a sparse polynomial in an explicit list of generators. Coefficients are exact rationals or symbolic parameter expressions, terms come back in SymPy's lex-descending order, and nothing is approximated. degree/coeffs/leading_coeff on Ex accept symbolic coefficients too. ratsimp is a rational-function normal form — one cancelled fraction with integer-primitive numerator and denominator — and solve uses it for parametric linear and quadratic equations.
let ctx = new;
syms!;
// Polynomial introspection on Ex with symbolic (var-free) coefficients
let e = &a * &x.powi + &x * + 3;
e.degree; // Some(2)
e.coeffs; // Some([3, a + 1, a]) (ascending)
e.leading_coeff; // Some(a)
// Poly: sparse terms over explicit generators, exact evaluation, calculus
let p = .as_poly.unwrap;
p.terms; // [([2, 0], a), ([1, 1], 3), ([0, 1], -1), ([0, 0], 1)]
p.coeff_monomial.unwrap; // 3
p.total_degree; // Some(2)
p.eval_gen.unwrap; // Poly(5*y + 4*a + 1, y)
p.derivative.unwrap.to_ex; // 2*a*x + 3*y
// Rational normal form: nested fractions collapse to one cancelled fraction
.ratsimp; // (x + y)/(x*y + 1)
.solve; // Ok([(3*j + 1)/(5*j - 1)])
// Exact sign of a rational-coefficient polynomial on an interval (square-free part + Sturm)
.poly_is_nonnegative_on; // Some(true)
.poly_is_positive_on; // Some(false) — touches 0 at x = 1
// Linear certificates: (x + 1)² = λ₁·(x + 1) + λ₂·(x² − 1) as an exact linear system
let = ;
let goal = .powi.as_poly.unwrap;
let basis = monomial_basis.unwrap; // [[2], [1], [0]]
let m = coefficient_matrix.unwrap; // [[0, 1], [1, 0], [1, -1]]
let b = coefficient_matrix.unwrap; // [[1], [2], [1]]
linsolve_matrix; // Ok(Unique([(x1, 2), (x2, 1)]))
Also: Poly::{from_terms, all_coeffs, degree_list, eval, add/sub/mul/pow/scale, content_and_primitive, monic, to_multipoly, nroots}, and on MultiPoly a heuristic multivariate gcd/lcm, integer_content, clear_denominators.
Simplification and the Rule Engine
simplify() tries a dozen strategies and iterates to a fixpoint. The pattern-matching engine behind it is public in 0.2: build your own rules (symbols ending in _ are wildcards, rest__ absorbs the rest of a sum or product), rewrite with them, trace what fired, and interleave them with the built-in simplifier.
let ctx = new;
syms!;
let = ;
let rules = from_rules;
.rewrite; // -cos(x)^2 + 4
.rewrite; // ln(x*y)
let = .simplify_traced;
// result = 1; steps name the strategy and the rules that fired
x.powi.subs_algebraic; // y^2 (plain subs would leave x^4)
.sqrt.sqrtdenest; // sqrt(2) + sqrt(3)
Equation Solving
Polynomial equations are solved through quartic by radicals; degree ≥ 5 produces RootOf nodes with numerical evaluation. Transcendental equations use inversion peeling and Lambert W. solve never lies: identities are Err(InfiniteSolutions), contradictions and range violations are Err(NoSolution).
let ctx = new;
syms!;
expr!.solve; // Ok([3, 2])
.solve; // Ok([1/6*pi, 5/6*pi])
.solve; // Err(NoSolution)
.solve; // Err(InfiniteSolutions)
// All periodic solutions, with an integer parameter
let fam = .solve_general.unwrap;
// fam.solutions = [2*n*pi + 1/6*pi, 2*n*pi + 5/6*pi], fam.parameters = [n]
// Linear systems: unique / parametric / inconsistent, symbolic coefficients allowed
let sol = linsolve.unwrap;
// LinearSolution::Parametric { solution: [x = -z/2 + 4, y = -z/2 + 2, z = z], free: [z] }
// Polynomial systems via Gröbner bases (algebraic solutions)
solve_system_ex;
// Ok([[√2/2, √2/2], [-√2/2, -√2/2]])
// Inequalities (sign-chart method), including absolute values
expr!.solve_gt; // (-oo, -2) ∪ (2, oo)
.solve_lt; // (-1, 3)
Differential Equations and Recurrences
16 ODE classes (separable, linear, Bernoulli, Riccati, Euler–Cauchy, exact, integrating factor, Clairaut, nth-order constant-coefficient, variation of parameters, systems via matrix exponential, …), initial-value problems, and linear recurrences.
let ctx = new;
syms!;
let y = ctx.symbol;
let = ;
.solve_ode; // y = C1*sin(x) + C2*cos(x)
.solve_ode_ivp;
// Ok(sin(x))
.classify_ode; // Clairaut
// a(n+2) = a(n+1) + a(n), a(0) = 0, a(1) = 1 → Binet's formula
rsolve_linear;
Sets and Logic
SetEx and BoolEx are first-class: intervals, finite sets, unions with a normal form, three-valued queries, and boolean normal forms with a DPLL satisfiability check.
let ctx = new;
syms!;
let a = ctx.interval; // [0, 5]
let b = ctx.interval; // (3, 10]
a.intersection.simplify; // (3, 5]
a.symmetric_difference; // [0, 3] ∪ (5, 10]
a.contains; // Some(false)
a.contains; // None
a.union.measure; // Some(10)
let conds = ;
reduce_inequalities; // Ok((2, 5])
let = ;
pp.and.or.simplify; // p > 0
pp.and.not.to_nnf; // 0 >= p | 0 >= q
pp.or.is_tautology; // Some(true)
Linear Algebra
Symbolic matrices with exact decompositions. The eigen family needs no dummy variable in 0.2, structure tests are three-valued, and preconditions are Results.
use q; // exact rational literal: q(1, 2) = 1/2
let ctx = new;
syms!;
let m = matrix!;
m.det.unwrap; // 3
m.eigenvals.unwrap; // [3, 1]
m.char_poly.unwrap; // λ^2 - 4*λ + 3
m.diagonalize.unwrap; // (P, D)
m.matrix_exp_t.unwrap; // [[e^(3t)/2 + e^t/2, …], …]
m.matrix_pow_symbolic.unwrap; // [[3^n/2 + 1/2, 3^n/2 - 1/2], …]
m.matrix_sqrt.unwrap;
let spd = matrix!;
spd.cholesky.unwrap; // [[2,0,0],[6,1,0],[-8,5,3]]
spd.is_positive_definite; // Some(true)
matrix!.qr.unwrap; // exact radicals
// Irreducible characteristic polynomials give exact, evaluable RootOf eigenvalues
matrix!.eigenvals.unwrap; // [RootOf(λ^3 - λ - 1, 0), …]
// 0.3: index-list extraction, exact rationals in and out, three-valued structure tests
m.extract.unwrap; // [[1], [2]]
from_ratio.unwrap; // [[1/2, 3]]
.is_zero; // Some(true) (m is symmetric)
// 0.3.5: QMatrix / ZMatrix — plain exact matrices over ℚ / ℤ, no expression arena.
// Fraction-free (Bareiss) elimination: a 30×30 rational inverse takes 10 ms, not 470.
let h = from_fn; // Hilbert matrix
h.det.unwrap; // 1/6048000
h.inv.unwrap; // 2800 (the inverse is integral)
let = from_i64.unwrap.rref;
// r = [[1, 0, -1], [0, 1, 2]], pivots = [0, 1]
from_i64.unwrap.smith_normal_form;
Matrix::{rref, rank, nullspace, det, inv, solve}, linsolve/linsolve_matrix and the normal forms route through QMatrix/ZMatrix automatically whenever every entry is a rational literal, so existing code gets the speed-up without changes.
Also: LU, LDLᵀ, Gram–Schmidt, Jordan form, pseudo-inverse, Kronecker product, rank/nullspace/rowspace, norms, least squares, Hessian, Wronskian, quaternions, vector calculus in Cartesian/cylindrical/spherical coordinates, state-space ↔ transfer function; select_rows/select_cols/delete_row/delete_col, from_bigint/from_f64_rows, to_rational_rows/to_bigint_rows, is_integer_matrix, subs_map, nnz.
Exact Optimization and Integer Lattices
Linear programs are solved over ℚ by a two-phase simplex with Bland's rule: optima, shadow prices and Farkas infeasibility certificates are exact, never "infeasible to within tolerance". Since 0.3.5 the tableau pivots on integers with a common denominator (no rational normalisation in the inner loop), which makes certificate-sized problems 5–30× faster. Integer matrices get Hermite and Smith normal forms with unimodular transforms, and ℤ-bases of integer kernels.
use ;
use hermite_normal_form_with_transform;
let ctx = new;
// max 5x + 4y s.t. 6x + 4y ≤ 24, x + 2y ≤ 6, x, y ≥ 0
let sol = maximize
.le
.le
.solve.unwrap;
sol.status; // Optimal
sol.x; // [3, 3/2]
sol.objective; // Some(21)
sol.duals; // [3/4, 1/2] shadow prices: yᵀb = 21 = cᵀx*
// x + y ≤ 1 and x + y ≥ 2 cannot both hold — here is the proof
let bad = minimize
.le
.ge
.solve.unwrap;
bad.status; // Infeasible
bad.farkas; // Some([1, -1]) Aᵀy = 0, yᵀb = −1 < 0
// "Is there μ ≥ 0 with Aμ = b?", exactly (Farkas / Carathéodory searches)
feasible_nonneg; // Ok(Some([21/10, 21/10]))
// Integer normal forms: H = U·A (row style), S = U·A·V, ℤ-basis of the kernel
let a = matrix!;
let = hermite_normal_form_with_transform.unwrap;
h; // [[2, 4, 4], [0, 6, 0], [0, 0, 12]]
.eval == h; // true (det U = −1)
a.smith_normal_form.unwrap; // [[2, 0, 0], [0, 6, 0], [0, 0, 12]]
matrix!.integer_nullspace.unwrap; // [(1, 0, −2)ᵀ, (0, 1, −1)ᵀ] — generates every integer solution
Also: linprog (SciPy-shaped), linprog_matrix (from Matrix data), per-variable bounds and free variables, nonneg_combination / feasible_nonneg_certified (cone membership with the separating Farkas vector on failure), column_hermite_normal_form (SymPy's convention), smith_normal_form_with_transforms, is_unimodular, lattice_determinant.
Certified Inequalities and Lean Export
certificates::prove_nonnegative_on_box proves goal ≥ 0 on a box with a Handelman certificate — an exact identity goal = Σ λₖ·Π(xᵢ − lᵢ)^a(uᵢ − xᵢ)^b with λ ≥ 0 found by the exact LP and re-verified with exact polynomial arithmetic — or refutes the claim with an exact counterexample. The certificate exports as a Lean 4 / Mathlib theorem whose proof is nlinarith over exactly those products; Ex::to_lean() renders any elementary expression in Mathlib syntax.
use ;
let ctx = new;
syms!;
let square = ;
let cert = match prove_nonnegative_on_box.unwrap ;
cert.to_string; // -x*y + 1 = -y + y*(-x + 1) + 1, 0 ≤ x ≤ 1, 0 ≤ y ≤ 1
cert.verify; // true — exact re-check, independent of the LP
cert.to_lean.unwrap;
// theorem one_minus_xy (x y : ℝ) (_h_x_lo : (0 : ℝ) ≤ x) (h_x_hi : x ≤ (1 : ℝ)) (h_y_lo : (0 : ℝ) ≤ y)
// (h_y_hi : y ≤ (1 : ℝ)) :
// 0 ≤ -(x * y) + 1 := by
// nlinarith [sub_nonneg.mpr h_y_hi, mul_nonneg (sub_nonneg.mpr h_x_hi) (sub_nonneg.mpr h_y_lo)]
prove_nonnegative_on_box.unwrap;
// Refuted { point: [0, 0], value: -1/2 }
.to_lean.unwrap; // "(x - 1) / (2 * x)"
x.sqrt.gt.to_lean.unwrap; // "0 < Real.sqrt x"
// 0.4: a polyhedron whose facets depend on a parameter j ≥ j₀. On { t ≥ r, t + j·r ≥ j + 1 }
// the goal t − 1 ≥ 0 needs the multiplier λ(j) = 1 + j: (j + 1)(t − 1) = j·h₀ + h₁.
use ;
syms!;
let hyps = ;
let out = prove_nonnegative_on_polyhedron.unwrap;
out.certificate.unwrap.to_string; // (j + 1)*(t - 1) = j*h0 + h1; h0 = -r + t, h1 = j*r - j + t - 1; j ≥ 0
out.certificate.unwrap.to_lean.unwrap;
// … have h0J := mul_nonneg hJ0 h0
// have hg : (0 : ℝ) ≤ (j + 1) * (t - 1) := by linarith only [h0J, h1]
// have hg' := nonneg_of_mul_nonneg_right hg (by linarith only [hJ0])
// linarith only [hg']
Also: prove_polyhedron_empty (the same identity with goal −1: a cell is empty for every j), prove_nonnegative_on_halfline / prove_nonnegative_on_reals (univariate, Pólya multipliers and square factors), lean_steps / lean_hints for dropping a proof into an existing skeleton, LeanOpts::prefer_subtraction, and symplex::polytope::Polytope for the exact geometry of the cells (vertices, volume, cuts, redundancy).
Transforms
let ctx = new;
syms!;
let a = ctx.symbol_with;
.exp.fourier_transform; // Ok(2*a/(a^2 + w^2))
.exp.fourier_transform; // Ok(sqrt(pi)*exp(-1/4*w^2))
.mellin_transform; // Ok((pi/sin(s*pi), re(s) > 0 & 1 > re(s)))
t.sin.laplace; // 1/(s^2 + 1)
.inverse_laplace; // H(t - 2)
x.sign.fourier_series_on.unwrap.truncate;
// 4*sin(x)/pi + 4*sin(3*x)/(3*pi) + 4*sin(5*x)/(5*pi)
Number Theory and Combinatorics
Pollard–Brent rho + ECM factorization, BPSW primality, modular square roots and discrete logarithms, continued fractions, Diophantine equations, and integer sequences.
use *;
use diophantine;
use *;
isprime; // false (Carmichael number)
factorint; // [(1048583, 1), (1048589, 1)] — ~1 ms
sqrt_mod; // Some(3)
discrete_log; // Some(4)
primepi; // Some(78498)
continued_fraction_periodic; // Some(([4], [1, 3, 1, 8]))
pell; // Some((1766319049, 226153980))
sum_of_two_squares; // Some((4, 7))
stirling2; // Some(34105)
partition_count; // Some(190569292)
crt_i64; // Some(23)
igcd; // 6 (gcd_many / lcm_many take BigInt slices)
ilcm; // 60
Numerical Toolbox
Deterministic, budgeted f64 routines — bracketed roots, derivative-free minimisation, global search in a box, least-squares fits — usable on plain closures or directly on expressions (which are compiled first). Bad input is Err(InvalidArgument), a non-finite value is Err(ComputationFailed); nothing panics.
use ;
let ctx = new;
syms!;
// Bracketed roots (Brent–Dekker), on a closure or on a compiled expression
brent_root.unwrap; // 1.41421356237…
.find_root_bracket.unwrap; // 0.739085133215…
// Nelder–Mead: local minimum from a starting point
nelder_mead.unwrap; // x ≈ [1, −2]
let rosen = .powi + 100 * .powi;
let r = rosen.minimize_numeric.unwrap; // r.x ≈ [1, 1], r.fun ≈ 1e-18, r.converged
// Differential evolution: global minimum in a box, deterministic for a given seed
let himmelblau = .powi + .powi;
himmelblau.minimize_global_numeric.unwrap; // fun < 1e-8
// Brent scalar minimisation, and least-squares fits (f64 via Householder QR, or exact rational)
.minimize_scalar_numeric.unwrap; // (0.36787944…, −0.36787944…) = (1/e, −1/e)
poly_fit.unwrap; // ≈ [1, 0, 2] (ascending: 1 + 2x²)
let pts = ;
poly_fit_points.unwrap; // 7/10*x + 7/10 (exact least-squares line)
Also: bisect, newton_root, golden_section, minimize_scalar, poly_fit_exact, linear_fit, trapezoid, eval_poly; RootOpts/MinimizeOpts/DeOpts for tolerances, budgets and seeds.
Code Generation: Rust, C99 and Compiled Closures
Symbolic expressions compile to optimized Rust or C functions with common subexpression elimination, mul_add/fma, integer powers as multiplications, optional domain assertions, and a self-contained special-function runtime.
let ctx = new;
syms!;
let f = &x.sin.powi + &.exp * 3;
f.to_rust_fn.unwrap;
// pub fn f(x: f64, y: f64) -> f64 { 3_f64.mul_add(2_f64.mul_add(x, y).exp(), x.sin().powi(2)) }
f.to_c_fn.unwrap;
// #include <math.h>
// double f(double x, double y) { return fma(3.0, exp(fma(2.0, x, y)), pow(sin(x), 2.0)); }
// Special functions embed only the helpers they need (Rust: `mod symplex_rt`; C: `static inline`)
x.lambertw.to_c_fn.unwrap; // contains symplex_lambert_w0
let opts = CodegenOptions ;
x.ln.to_c_fn_with_options.unwrap; // float g(float x) { return assert(x > 0.0f), logf(x); }
// Compiled closures: Result, arity-checked, Send + Sync; gradients share one CSE pass
let cf = f.compile.unwrap;
cf;
let grad = compile_many.unwrap;
grad.call_vec;
f.to_latex; // \sin^{2}\left(x\right) + 3\exp\left(2x + y\right)
For build.rs pipelines and no_std targets see symplex-build; for the browser see symplex-wasm.
Compile-Time Dimensional Analysis
Physical quantity types are checked at compile time. Adding a Mass to a Length is a compiler error. Differentiation respects dimensions: d(Length)/d(Time) produces Velocity.
use *;
let ctx = new;
let m = symbol;
let a = symbol;
// dim! macro: the `: Force` annotation is a compile-time assertion
let force = dim!; // F = m·a [N]
// Typed calculus: d(Length)/d(Time) → Velocity
syms!; // raw symbols for expr!
let t_var = symbol;
let position = from_ex;
let velocity: Velocity = position.diff_wrt; // g·t [m/s]
// 30 named quantity types, ~100 unit conversions (all exact rationals)
// Mass + Length → compile error
Design Principles
-
Exact by default. Every number is
Ratio<BigInt>. No floating-point contamination.0.1 + 0.2 == 3/10, not0.30000000000000004. Floats only appear on expliciteval_f64(),compile(), orintegrate_numeric().Context::from_f64converts a float to its exact dyadic rational;from_f64_approxto the nearest bounded-denominator rational. -
Explicit contexts. Every expression belongs to a
Context. No hidden global state. Mixing expressions from different contexts is caught immediately (compiler-enforced private field + runtime guard). -
Type-safe expressions.
Ex(numeric),BoolEx(boolean),SetEx(set-valued) are distinct types.sin(bool_expr)is a compile error. -
Thread-safe.
ContextisClone(Arc-based),ExisSend + Sync. Multiple threads can share a context safely. -
No recursion. All tree traversals use explicit stacks. Deep expressions don't blow the call stack.
-
Never silently wrong. Numerical evaluation returns
Result. Operations that can't produce a closed form return unevaluated symbolic nodes —∫x^x dxreturnsIntegral(x^x, x), not garbage.∫₋₁¹ dx/x²isErr(Divergent), not−2.re(z)staysre(z)unlesszis known to be real.
The API Model
Every symbolic operation that might not produce a closed-form result has two entry points:
| Intent | Method | Returns | When to use |
|---|---|---|---|
| Give me math | integrate(&x) |
Ex (always — may contain Integral nodes) |
Interactive exploration, chaining |
| Fail if you can't | try_integrate(&x) |
Result<Ex> |
Pipelines, codegen, safety-critical |
try_ twins exist for diff, integrate, integrate_definite, limit, limit_left/right/dir, series, series_at_infinity, summation, product_over, laplace, inverse_laplace, residue, gosper_sum, solve_ode, solve_gt/ge/lt/le. Check any expression for unevaluated forms:
let anti = hard_expr.integrate;
if anti.has_unevaluated
(RootOf and RootSum are not unevaluated: they are complete algebraic answers.)
Operations that always succeed (simplify, expand, eval, factor, subs, rewrite) return Ex with no try_ variant — "unchanged" is a valid answer.
Result boundaries. Crossing from symbols to numbers (eval_f64, eval_decimal, compile, to_rust_fn, to_c_fn, integrate_numeric) always returns Result. So do operations with structural preconditions (Matrix::inv, cholesky, lu, minor, matmul) and solvers whose failure is a mathematical fact: solve returns Err(InfiniteSolutions) for identities and Err(NoSolution) for contradictions, try_integrate_definite returns Err(Divergent), laplace_final_value returns Err(Divergent) for unstable poles. Transform APIs without an unevaluated node (fourier_transform, mellin_transform, z_transform) are Result-only.
Three-valued queries. is_positive, equals, is_convergent, SetEx::contains, is_subset, Matrix::is_symmetric, is_diagonalizable, is_positive_definite, BoolEx::is_tautology, vector::is_conservative … return Option<bool>: yes, no, or unknown. degree, resultant, discriminant, hypergeometric_ratio return Option<T>.
Comparison with SymPy
| Feature | symplex 0.3 | SymPy |
|---|---|---|
| Arithmetic | Exact Ratio<BigInt> |
Exact (similar) |
| Differentiation | Complete, incl. Bessel/orthogonal/polygamma | Complete |
| Indefinite integration | 15+ strategies incl. Risch + LRT log-to-real | Risch + heurisch + Meijer G (broader) |
| Definite / improper integration | Singularity detection, ~30-entry improper table, divergence reported as Err |
Meijer G-based; much broader table |
| Numeric integration | Adaptive G7/K15 quadrature | via mpmath (more algorithms) |
| Summation | Faulhaber, Gosper, telescoping, binomial, p-series, power-series recognition | + Zeilberger, hypergeometric closed forms (broader) |
| Polynomial solving | Through quartic + RootOf |
Through quartic + CRootOf |
| General solutions | solve_general (periodic families) |
solveset with ImageSet |
| Linear systems | linsolve (unique / parametric / inconsistent, symbolic) |
linsolve (similar) |
| Polynomial systems | Gröbner + FGLM, algebraic solutions | Gröbner, more strategies |
| Polynomial views | Poly over explicit generators with symbolic coefficients, lex terms, coefficient matrices |
Poly with domains, factoring, gcd, resultants (broader) |
| Rational simplification | ratsimp/cancel: heuristic multivariate GCD, opaque subexpressions as indeterminates |
cancel, ratsimp, together, apart (similar; more GCD algorithms) |
| Linear programming | Exact two-phase simplex; duals and Farkas certificates in the public result | sympy.solvers.simplex (lpmin/lpmax/linprog, exact; optimum and argmin only) |
| Integer normal forms | Row and column HNF with transform, SNF with transforms, integer nullspace, lattice index | hermite_normal_form, smith_normal_form (no transforms returned) |
| Numerical optimisation | Brent, bisection, Newton, Nelder–Mead, differential evolution, QR least squares, exact rational fits | Defers to SciPy / mpmath (nsolve, findroot); far broader via SciPy |
| Series expansion | Taylor / Laurent / at ∞ / formal power series with general terms | + O() notation, Puiseux |
| Limits | Gruntz with work budget, one-sided | Gruntz (more mature) |
| Simplification | Multi-strategy fixpoint + public rule engine with AC matching, tracing | More strategies; replace/Wild patterns |
| Factoring | Berlekamp–Zassenhaus (any degree), multivariate via Kronecker | Zassenhaus + Wang (faster multivariate), algebraic extensions |
| Matrices | Eigen/Jordan/exp/sqrt/pow, QR, Cholesky, LDL, LU, RootOf eigenvalues |
More decompositions (SVD, Schur), sparse |
| ODE solving | 16 classes, IVPs, systems | More classes, hints, series solutions |
| Recurrences | Linear constant-coefficient, first-order | rsolve (poly/rational/hyper) |
| Transforms | Laplace, Fourier (3 conventions), Mellin (with strip), Z, Fourier series | Broader tables, Hankel, cosine/sine |
| Sets & logic | Interval algebra, three-valued queries, NNF/CNF/DNF, DPLL | Richer set types (ImageSet, ConditionSet), satisfiable |
| Number theory | rho/ECM, BPSW, sqrt_mod, dlog, Pell, two squares | Broader (quadratic forms, general Diophantine) |
| Combinatorics | Stirling, Bell, partitions, derangements, multinomial | Broader (permutation groups, etc.) |
| Special functions | Γ, ψ⁽ⁿ⁾, erf, B, W, Bessel, Si/Ci/Ei/li, ζ, orthogonal polys | Many more (hypergeometric, elliptic, Meijer G) |
| Algebraic numbers | ℚ(α) field with exact zero/sign testing |
AlgebraicNumber + ANP |
| Code generation | Rust and C99 with CSE, fma, embedded special-function runtime |
Python / C / Fortran / Rust / Julia via codegen |
| Dimensional analysis | Compile-time type checking | Runtime physics.units |
| Thread safety | Send + Sync, no GIL |
GIL-bound |
| Expression type safety | Ex / BoolEx / SetEx at compile time |
Runtime only |
| Language | Rust (compiled, ~174K lines, 92 node types) | Python (interpreted) |
Where SymPy is stronger: geometry, statistics, tensor algebra, quantum mechanics, general Diophantine equations, PDE solving, hypergeometric/Meijer-G machinery, and 30 years of community contributions and testing.
Where symplex is different: compile-time dimensional analysis, thread safety, Rust and C code generation with an embedded runtime, exact arithmetic without Python overhead, algebraic number field arithmetic with exact zero/sign testing, construction-time radical simplification, and exact certificates — LP duals and Farkas vectors, unimodular HNF/SNF transforms, Sturm-verified polynomial signs — as first-class results. Operations that can't complete return honest unevaluated forms or Err rather than guessing.
Migrating from 0.1
The CHANGELOG lists every breaking change with its replacement. The ones most likely to touch your code:
| 0.1 | 0.2 |
|---|---|
expr.compile(&["x"]) → Option<…> |
expr.compile(&["x"])? → Result<CompiledFn> (arity(), try_call()) |
expr.definite_integral(&x, &a, &b) |
expr.integrate_definite(&x, &a, &b) / try_integrate_definite (Err(Divergent)) |
solve returned Ok(vec![]) for identities |
Err(InfiniteSolutions) / Err(NoSolution); roots are eval'd |
m.eigenvals(&lam), m.jordan_form(&lam), m.matrix_exp(&t) |
m.eigenvals(), m.jordan_form(), m.matrix_exp() / m.matrix_exp_t(&t) |
m.cholesky() → Option, m.lu() → tuple |
both Result |
m.minor(i, j) → sub-matrix |
m.minor(i, j) → Result<Ex>; sub-matrix is m.minor_matrix(i, j) |
m.is_symmetric() → bool |
Option<bool> (also is_diagonalizable, vector::is_conservative, …) |
ctx.solve_system(…) → Vec |
Result<LinearSolution> (Unique / Parametric / Inconsistent) |
z.re() assumed z real |
re(z) stays symbolic; declare Assumption::Real |
has_unevaluated() true for RootOf |
RootOf/RootSum are answers, not unevaluated forms |
expr.textplot(…) → String |
Result<String> (all plotting methods) |
iter.sum::<Ex>() on empty iterator → 0 |
panics; use ctx.sum(iter) or Option<Ex> |
The book's migration guide has worked examples.
Migrating from 0.2
No signatures changed between 0.2 and 0.3. A few operations now return a different but equivalent form, which matters only if you compare printed output or match on structure:
| 0.2 | 0.3 |
|---|---|
degree / coeffs / coeff / leading_coeff returned None (and is_polynomial false) for symbolic coefficients (a*x^2 + x) |
They succeed: Some(2), Some([0, 1, a]), Some(a), true |
solve on a linear or quadratic equation with parametric fractional coefficients returned a fraction of fractions |
The root (and the quadratic discriminant) is ratsimp'd: ((3r−1)/(j+1) − (r+1)/(2j)).solve(&r) → (3*j + 1)/(5*j - 1) |
simplify_rational = together + per-symbol cancel; could leave nested fractions uncancelled |
simplify_rational is ratsimp: one fraction, all variables at once, integer-primitive numerator and denominator |
Everything else in the CHANGELOG is additive.
Examples
Every example is self-contained and runs in a few seconds; CI runs all of them.
Getting started:
cargo run --example quickstart # Tour of core operations
cargo run --example repl # Interactive expression evaluation
cargo run --example readme_snippets # Every code block in this README, executed
New in 0.3:
cargo run --example polynomials # Poly views with symbolic coefficients, ratsimp, linear certificates, exact sign on an interval
cargo run --example certificates_to_lean # Handelman certificates on a box (exact LP, exactly re-verified) exported as Mathlib theorems
cargo run --example exact_lp # Exact simplex: optima, shadow prices, Farkas certificates, feasible_nonneg, linprog_matrix
cargo run --example integer_lattices # Row/column HNF with transforms, Smith normal form, integer nullspace, unimodularity, lattice index
cargo run --example numeric_optimization # Brent/Newton roots, Nelder–Mead, differential evolution, polynomial fits (f64 and exact)
New in 0.2:
cargo run --example definite_integration # Improper integrals, divergence detection, quadrature, residues
cargo run --example summation_and_series # Σ/Π closed forms, convergence, formal power series, finite differences
cargo run --example complex_analysis # re/im/conjugate/arg, new constants, Si/Ci/Ei/ζ/polygamma
cargo run --example rule_engine # Custom rewrite rules, tracing, subs_algebraic, targeted simplifiers
cargo run --example linear_systems_and_ivp # solve semantics, linsolve, solve_general, ODE IVPs, rsolve
cargo run --example sets_and_logic # Interval algebra, reduce_inequalities, CNF/DNF, tautology
cargo run --example factoring_and_ntheory # Zassenhaus factoring, factorint, sqrt_mod, Pell, continued fractions
cargo run --example transforms # Fourier, Mellin, Laplace, Fourier series, Z, one-sided limits
cargo run --example matrix_decompositions # QR, Cholesky, LDL, Jordan, matrix_exp_t, RootOf eigenvalues
cargo run --example c_codegen # C99 backend, embedded runtime, compile_many (compiles the C if cc exists)
Engineering workflows:
cargo run --example pid_controller # PID design → stability → Rust codegen
cargo run --example robotics_codegen # DH parameters → Jacobian → optimized Rust
cargo run --example control_system # State-space, transfer functions, pole placement, ZOH
cargo run --example signal_filter # Bilinear transform → digital filter → codegen
cargo run --example dynamics # Lagrangian mechanics, equations of motion
cargo run --example inverse_kinematics # 2-DOF IK via Gröbner bases
Mathematics and science:
cargo run --example calculus # Differentiation, integration, limits, series
cargo run --example equation_solving # Polynomial, transcendental, system solving
cargo run --example matrix_algebra # Eigenvalues, Jordan form, codegen
cargo run --example ode_solving # ODE classification and solving
cargo run --example combinatorics_counting # Stirling numbers, partitions, multinomials
cargo run --example optimization # Gradient, Hessian, critical points
cargo run --example complex_numbers # Euler's formula, complex roots
cargo run --example laplace_transforms # Forward, inverse, z-transforms
Applied problems:
cargo run --example gradient_descent # Symbolic gradient → compiled optimization loop
cargo run --example crypto_rsa # RSA with number theory primitives
cargo run --example number_theory # Primality, factorization, CRT
Dimensional analysis:
cargo run --example units_physics # Compile-time unit checking
cargo run --example units_electrical # Circuit analysis with units
cargo run --example units_engineering # Motor design, imperial conversions
cargo run --example units_kinematics # Kinematics with typed quantities
cargo run --example units_lagrangian # Lagrangian mechanics with units
Output:
cargo run --example latex_output # LaTeX rendering
cargo run --example physics_constants # Physical constants (symbolic + exact)
Dependencies
All MIT or Apache-2.0 licensed. No C bindings. No LGPL.
Core: num-bigint, num-rational, num-traits, num-integer, smallvec, rustc-hash, bitflags, parking_lot, thiserror, astro-float, serde, serde_json, tracing, typenum.
Proc macros: syn, quote, proc-macro2.
Companion crates: symplex-build (build-time codegen for no_std firmware), symplex-wasm (browser bindings).
Requirements
Rust 1.93+ (Edition 2024). No optional features; pure Rust on every platform Rust targets, including wasm32-unknown-unknown.
License
Dual-licensed under MIT and Apache 2.0.