Alkahest

Alkahest is a high-performance computer algebra system with a Rust core and Python API.

Install: pip install alkahest from PyPI (Python 3.9–3.13). Optional RL environments: pip install "alkahest[rl]" (Python ≥ 3.10). Default wheels include vendored egglog (egraph), the Gröbner solver (groebner — alkahest.solve and related APIs work out of the box), the Cranelift JIT (cranelift) and the multi-core paths (parallel). Only the LLVM JIT is omitted; Linux linux_x86_64 +jit / +full wheels ship on GitHub Releases. The Getting started chapter has full install notes (venv, LD_LIBRARY_PATH, building from source). See Reinforcement learning for verifiable RL environments.

What it is

A general-purpose symbolic math library designed around three axes:

Performance. The Rust kernel uses hash-consed directed acyclic graphs so structural equality is a pointer comparison and subexpression sharing is automatic. FLINT backs polynomial arithmetic. An LLVM JIT compiles symbolic expressions to native or GPU code at runtime. Common operations run orders of magnitude faster than SymPy.

Correctness. Every simplification and transformation produces a derivation log — the exact sequence of rewrite rules applied, with arguments and side conditions. A subset of operations can export Lean 4 proof terms verifiable by an independent checker.

Ergonomics. The Python API uses operator overloading for natural expression construction. Results are rich objects with .value, .steps, and .certificate attributes. Error messages carry structured codes, location information, and suggested remediations.

Agent loops. Budgets and cancellation, batch APIs that never abort on one bad candidate, versioned compact result envelopes, and session-level claim graphs are first-class — see Autoresearch / agent loops.

Design principles

Explicit representations. The type system distinguishes UniPoly (FLINT-backed univariate polynomial), MultiPoly (sparse multivariate), RationalFunction, and the generic Expr tree. Converting between them is an explicit call. There are no silent representation changes hiding performance cliffs.

Stateless by design. No global assumption contexts. No hidden caches that change behavior. All context (domains, simplification policy, precision, budgets) is passed explicitly or scoped through context(...). This makes results deterministic and parallelism safe.

Composable transformations. trace, grad, jit, and certify operate on a shared traced representation and stack freely: jit(grad(f)) compiles a derivative, jit(grad(grad(f))) compiles a second derivative.

A small primitive set. Each primitive (sin, exp, add, mul, …) registers a full bundle: simplification rule, forward- and reverse-mode differentiation, MLIR lowering, Lean theorem tag, numerical evaluation. New operations are added by registering a primitive, not by adding code paths across the system.

Compared to alternatives

This guide

The guide covers the Rust-level design concepts. For the generated Python API reference (Sphinx), see the API documentation.

Getting started

Install

PyPI (default)

Alkahest is on the Python Package Index. Supported interpreters are Python 3.9 through 3.13 (requires-python on PyPI).

python -m pip install -U pip
pip install alkahest

Use a virtual environment when you also build from source or test multiple Python versions:

python3 -m venv .venv && source .venv/bin/activate   # Windows: .venv\Scripts\activate
python -m pip install -U pip
pip install alkahest

Default PyPI wheels are built with egraph, groebner, cranelift and parallel: the vendored egglog e-graph backend, the Gröbner solver (so alkahest.solve, Diophantine and homotopy work out of the box), the pure-Rust Cranelift JIT (so alkahest.jit_is_available() is True without a system LLVM), and the Rayon-backed multi-core paths. They do not include the LLVM JIT (llvm_jit) or cuda.

parallel is on in the default wheel — but still check it

capabilities()["features"]["parallel"] is True on every wheel published to PyPI, on Linux, macOS and Windows alike. The parallel Cargo feature is what powers the sharded ExprPool, the parallel F4 reduction, and every *_par entry point (numpy_eval_par, simplify_par), so those are genuinely multi-core out of the box.

It is still the one feature worth probing rather than assuming, because its absence is silent. parallel is not a Cargo default: a source build that does not pass --features parallel still gets working numpy_eval_par and simplify_par — they transparently fall back to their single-threaded counterparts, with no error, no warning, and no speedup. Benchmarking numpy_eval_par against numpy_eval on such a build measures the same number twice.

Never infer parallelism from the function existing; ask:

import alkahest as ak
if not ak.capabilities()["features"]["parallel"]:
    ...  # a source build without --features parallel; *_par is a no-op alias here

There is no pip install alkahest[jit] / alkahest[full] that swaps the native extension: pip extras only add Python dependencies, not alternate binaries.

For native LLVM CPU JIT use an opt-in +jit or +full Linux wheel from GitHub Releases (below), or build from source with --features jit. See the repository README.md for the same policy in short form.

Optional Linux wheels (+jit / +full)

Tagged releases attach linux_x86_64 wheels on GitHub Releases (CI builds them on ubuntu-22.04; these are not the manylinux wheels published as the default PyPI binaries). Pick the .whl whose tags match your Python (cp311, etc.) and linux_x86_64.

Example direct installs (replace <version> and the wheel name using the release asset list):

pip install "https://github.com/alkahest-cas/alkahest/releases/download/v<version>/alkahest-<version>+full-cp311-cp311-linux_x86_64.whl"
pip install "https://github.com/alkahest-cas/alkahest/releases/download/v<version>/alkahest-<version>+jit-cp311-cp311-linux_x86_64.whl"

These wheels vendor LLVM and related .so files under site-packages/alkahest.libs/. If import alkahest fails with a missing libLLVM-*.so or libffi-*.so, prepend that directory to LD_LIBRARY_PATH (or install matching system packages).

If your downloader rejects + in the URL, percent-encode it in the filename segment (e.g. 2.0.2%2Bfull).

After +jit or +full, alkahest.jit_is_available() should be True. Gröbner-backed APIs such as alkahest.solve are available in all wheels (including the default PyPI wheel) since groebner became a default Cargo feature in 2.3.1.

macOS and Windows +jit / +full wheels are not produced in CI yet; use building from source there.

Roadmap: a small PEP 503 extras index URL hosting only +jit / +full wheels (PyTorch-style --extra-index-url). Until then, use PyPI for the default wheel or direct URLs / asset downloads from Releases.

Build-profile verification

Every published wheel runs an import-and-capability smoke test in release CI. After installing, inspect the exact native build rather than inferring features from available Python functions:

import alkahest as ak

features = ak.capabilities()["features"]
print(features)

parallel is called out separately because it is the one feature whose absence is silent: numpy_eval_par and simplify_par exist and work in every build, and simply stop being parallel when it is off. Every published wheel now enables it, so the row that can surprise you is the last one — parallel is not a Cargo default, so a source build has to ask for it.

jit and cranelift remain compatibility names in this mapping. Prefer llvm_jit and cranelift_jit when selecting a backend explicitly. cuda indicates that the extension was compiled with NVPTX codegen — it guarantees that ak.compile_cuda and ak.CudaCompiledFn exist, but not that a usable GPU is present at runtime, and it is in no published wheel. There is no groebner_cuda bit: the GPU Gröbner kernel has no Python binding, so the bit was unfalsifiable from Python and was removed in contract v3. Read GPU support before branching on cuda.

Optional: RL environments (alkahest[rl])

Reinforcement-learning environments (symbolic integration, Prime Intellect Hub) are an optional extra. Requires Python ≥ 3.10 (verifiers does not support 3.9).

pip install "alkahest[rl]"

This adds verifiers and datasets. Environment code ships in the main wheel under alkahest.rl. See the RL guide for API details, veRL integration, and Environments Hub publishing.

From source

For optional Cargo features (jit, parallel, cuda, …), GPU/NVPTX, or development, build the PyO3 extension with maturin. The groebner and egraph features are default and included automatically.

Prerequisites (typical): Rust stable (≥ 1.76) and nightly, LLVM 15 (only for --features jit), FLINT (≥ 2.9, 3.x recommended; pulls in GMP/MPFR). See the repository README for distro-specific package names.

FLINT is a hard requirement of every source build. There is no FLINT-free configuration: UniPoly is a FLINT polynomial, and factorization, resultants, normal forms and number_theory call FLINT directly with no pure-Rust fallback. The flint3 Cargo feature selects which FLINT version’s API to use — it does not make FLINT optional. Without it the build stops early with an install hint (sudo apt-get install libflint-dev, dnf install flint-devel, brew install flint, conda install -c conda-forge libflint, …).

Without root, build FLINT into a user-local prefix and point the build at it:

FLINT_LIB_DIR=$PREFIX/lib FLINT_INCLUDE_DIR=$PREFIX/include \
  maturin develop --manifest-path alkahest-py/Cargo.toml --release
export LD_LIBRARY_PATH=$PREFIX/lib      # DYLD_LIBRARY_PATH on macOS

Both variables also feed FLINT version detection, so a locally built FLINT 3 is recognised as FLINT 3. ALKAHEST_SKIP_FLINT_CHECK=1 bypasses the presence probe. If you only need a working install, the PyPI wheels already bundle FLINT.

pip install maturin
git clone https://github.com/alkahest-cas/alkahest.git
cd alkahest
maturin develop --manifest-path alkahest-py/Cargo.toml --release

The default build already includes egraph and groebner. Additional optional features:

# LLVM JIT for native compiled evaluation
maturin develop --manifest-path alkahest-py/Cargo.toml --release --features jit

# Pure-Rust Cranelift JIT (fast compile, no system LLVM required)
maturin develop --manifest-path alkahest-py/Cargo.toml --release --features cranelift

# Parallel simplification and parallel F4 (sharded ExprPool + numpy_eval_par)
maturin develop --manifest-path alkahest-py/Cargo.toml --release --features parallel

# CUDA / NVPTX codegen (requires CUDA toolkit and LLVM with NVPTX target)
maturin develop --manifest-path alkahest-py/Cargo.toml --release --features cuda

# Full native build (JIT + parallel; egraph and groebner are already default)
maturin develop --manifest-path alkahest-py/Cargo.toml --release --features "parallel cranelift jit"

Rust crate

alkahest-cas is also published on crates.io (docs.rs) for use directly from Rust:

[dependencies]
alkahest-cas = "2"

# groebner is included by default; add other optional features as needed:
# alkahest-cas = { version = "2", features = ["parallel", "egraph"] }

System prerequisites (same libraries as the Python build — must be installed before cargo build):

# Debian / Ubuntu
sudo apt-get install -y libflint-dev libgmp-dev libmpfr-dev

# macOS (Homebrew)
brew install flint

The jit feature additionally requires LLVM 15 dev headers (llvm-15-dev / brew install llvm@15).

A self-contained runnable example is in examples/rust_quickstart/.

First steps

Every computation starts with an ExprPool. It owns all expressions; you create symbols and integers from it.

import alkahest
from alkahest import ExprPool, diff, simplify, integrate, sin, exp, cos

pool = ExprPool()
x = pool.symbol("x")
y = pool.symbol("y")

Building expressions

Python operators build expression trees:

expr = x**2 + pool.integer(2) * x + pool.integer(1)
print(expr)  # x^2 + 2*x + 1

Math functions accept expressions:

f = sin(x**2) + exp(x * y)

Parsing expressions from strings

Use parse when the expression comes from user input or a config file:

from alkahest import parse

e = parse("x^2 + 2*x + 1", pool, {"x": x})
print(e)   # x^2 + 2*x + 1

Identifiers not in the symbols dict are auto-created as symbols in pool. Both ^ and ** denote exponentiation. See Parsing from strings for the full syntax reference.

Simplification

r = simplify(x + pool.integer(0))
print(r.value)  # x
print(r.steps)  # [RewriteStep(rule='add_zero', ...)]

Differentiation

dr = diff(sin(x**2), x)
print(dr.value)  # 2*x*cos(x^2)

Integration

r = integrate(exp(x), x)
print(r.value)   # exp(x)

r = integrate(sin(x), x)
print(r.value)   # -cos(x)

Polynomial arithmetic

from alkahest import UniPoly, RationalFunction

# Convert to FLINT-backed univariate polynomial
p = UniPoly.from_symbolic(x**3 + pool.integer(-1), x)
q = UniPoly.from_symbolic(x + pool.integer(-1), x)
print(p.gcd(q))          # x - 1
print(p // q)            # x^2 + x + 1

Compiled evaluation

from alkahest import compile_expr, eval_expr

# Scalar evaluation via a dict binding
result = eval_expr(x**2 + y, {x: 3.0, y: 1.0})
print(result)  # 10.0

# JIT-compiled callable
f = compile_expr(x**2 + pool.integer(1), [x])
print(f([3.0]))  # 10.0

Vectorized evaluation over NumPy arrays

import numpy as np
from alkahest import compile_expr, numpy_eval

f = compile_expr(sin(x) * exp(pool.integer(-1) * x), [x])
xs = np.linspace(0, 10, 1_000_000)
ys = numpy_eval(f, xs)  # vectorised; much faster than a Python loop

Context manager

with alkahest.context(pool=pool, simplify=True):
    z = alkahest.symbol("z")  # uses the active pool
    expr = z**2 + alkahest.sin(z)

Agent / autoresearch loops

For a fan-out of candidates under a wall-clock or step budget, with results that survive context compaction:

import alkahest as ak

pool = ak.ExprPool()
x = pool.symbol("x")

with ak.context(pool=pool, budget=ak.Budget(wall_ms=100, seed=1)):
    outs = ak.integrate_many([x**2, ak.sin(x)], x)
    for item in outs:
        if item.ok:
            print(item.value.to_dict(mode="compact")["verification"]["status"])
        else:
            print(item.error["code"])  # e.g. E-INT-001 or E-BUDGET-001

Full picture: Autoresearch / agent loops, Budgets, Batch, Claim graphs.

Running the examples

The examples/ directory in the Git repository has runnable end-to-end scripts. With alkahest installed (pip install alkahest or maturin develop as above), from the repository root run:

python examples/calculus.py
python examples/polynomials.py
python examples/jit_eval.py
python examples/ball_arithmetic.py
python examples/ode_modeling.py

If you are developing without installing the extension into the active environment, set PYTHONPATH=python so the pure-Python package is importable alongside your build.

Python API reference

The Sphinx-generated Python API is published alongside this guide:

Open the Python API documentation

It includes ExprPool, simplification, calculus, polynomials, numerics, transforms, matrices, ODE/DAE, solvers, codegen, error types, and the search / workload surface (Budget, batch_map, DerivedResult.to_dict, …).

Conceptual chapters for agent-facing plumbing:

Submodules

Not everything lives on the top-level namespace as a function. These are reached as alkahest.<name> and documented in their own chapters:

alkahest.ansatz, alkahest.crosscheck and alkahest.smt are new in 3.8. They are in alkahest.__all__ and resolve on attribute access without a separate import, as do their error classes AnsatzError, CrossCheckError and SmtError.

Kernel design

The expression kernel is the foundation everything else builds on. It lives in alkahest-core/src/kernel/.

Hash-consed DAG

Every expression is represented as a directed acyclic graph stored in an ExprPool. Nodes are interned: before inserting a new node, the pool checks whether a structurally identical node already exists. If it does, the existing ExprId is returned instead of allocating a new node.

This gives three properties:

1. Structural equality is a pointer comparison. id_a == id_b iff the expressions are structurally identical. No tree traversal required.
2. Automatic subexpression sharing. If sin(x²) appears in ten different expressions, there is only one sin(x²) node in memory.
3. Hash-based memoization is cheap. Caching the result of a transformation keyed by ExprId is O(1) and correct. Hot recursive paths (simplify, differentiation, integration guards, JIT interpreter) use per-call HashMap<ExprId, T> memo tables so shared DAG nodes are processed once, not once per tree occurrence.

ExprPool

ExprPool is the intern table. It owns all expressions in a session.

pool = ExprPool()
x = pool.symbol("x")       # intern a Symbol node
n = pool.integer(42)       # intern an Integer node

Multiple pools are independent. An ExprId from one pool must not be mixed into another — the pool validates this in debug builds.

Persistent pool (V1-14). A pool can be serialized to disk and reopened, preserving all ExprIds across sessions:

pool.save_to("session.alkp")
pool2 = ExprPool.load_from("session.alkp")

Sharded pool. With --features parallel, the intern table uses a sharded concurrent hashmap (DashMap), allowing multiple threads to insert expressions without contention.

ExprData variants

Each interned node is one of:

Add and Mul are n-ary: a + b + c is one Add node with three children, not two nested Add nodes. Children are sorted at construction time so that commutativity is structural — a + b and b + a produce the same interned node.

Domains

Every symbol carries a domain as part of its structural identity:

x_real = pool.symbol("x", "real")
x_complex = pool.symbol("x", "complex")
# x_real and x_complex are distinct expressions — different ExprIds

The domain is not a global assumption; it is part of what the symbol is. Simplification rules can query a symbol’s domain to decide whether a rewrite is valid (e.g. sqrt(x²) → x requires x to be non-negative).

Available domains: real, positive, nonnegative, integer, complex. The default when no domain is specified is real.

ExprId and memory

ExprId is a 32-bit index into the pool’s internal arena. It is Copy, Send, and Sync. Cloning an ExprId is free. No reference counting is needed because the pool owns all nodes; expressions are not freed until the pool is dropped.

That last clause is a hard limit, not an implementation detail. The arena is append-only: there is no clear, no truncate, no refcount and no GC, so nothing is ever reclaimed while the pool is alive, and the storage cannot shrink. A distinct expression costs roughly 200 bytes of resident memory per node, permanently. A loop that builds a module-scope pool once and then calls into it forever grows linearly and without bound, at flat per-call latency — so it OOMs with no slowdown to warn you first. Every Expr, Matrix, Series and DerivedResult holds a strong reference to its pool, so retaining one result retains the whole history.

The supported pattern for unattended work is one pool per problem, dropped when the problem is done: see Budgets → ExprPool never reclaims.

The kernel is designed with parallelism as a first-class property. All kernel types are Send + Sync. The simplification and differentiation passes can run concurrently on disjoint ExprIds from the same pool.

Interning cost model

Interning a new node requires:

1. Hash the ExprData.
2. Look up in the concurrent hash map.
3. On miss: allocate the node in the arena and insert into the map.
4. On hit: return the existing ExprId.

Step 4 (the common case in a running computation) is a single hash lookup plus a pointer load. The arena uses bump allocation, so step 3 is also fast.

The memory benchmark group in alkahest-core/benches/alkahest_bench.rs verifies that rebuilding an identical expression tree does not grow the pool. Note the scope of that guarantee: it covers repeated work, not new work. A stream of distinct inputs grows the pool by every node it interns, and none of it comes back — see the warning under ExprId and memory.

One documented exception to “identical input does not grow the pool”: Matrix.eigenvals() interns a fresh gensym per call, so it grows by about 1.9 KB per call even on the same matrix. Cache its result rather than recomputing.

Expression representations

Alkahest exposes multiple representation types rather than hiding everything behind a single Expr. This is a deliberate design decision: the representation is visible, conversions are explicit, and performance characteristics are predictable.

Choosing a representation

Conversion to a specialized type is always an explicit opt-in:

expr = x**3 + pool.integer(-2) * x + pool.integer(1)
p = UniPoly.from_symbolic(expr, x)   # explicit conversion

If the expression cannot be represented in the target type (e.g. sin(x) as a polynomial), a ConversionError is raised with a remediation hint.

Expr

The generic symbolic expression. All other types convert to and from Expr. Built by operator overloading on the Python side:

expr = x**2 + pool.integer(3) * x * y - pool.integer(1)

Operations like diff, simplify, and integrate work on Expr and return DerivedResult objects wrapping an Expr.

UniPoly

Dense univariate polynomial backed by FLINT. Coefficients are exact integers or rationals stored in a FLINT polynomial object.

from alkahest import UniPoly

# x^3 - 2x + 1
p = UniPoly.from_symbolic(x**3 + pool.integer(-2) * x + pool.integer(1), x)

print(p.degree)          # 3
print(p.coefficients())  # [1, -2, 0, 1]  (constant first)
print(p.leading_coeff)   # 1  (a property, and an exact Python int)

# Arithmetic — all FLINT-backed, exact
q = UniPoly.from_symbolic(x + pool.integer(-1), x)
print(p * q)             # x^4-x^3-2*x^2+3*x-1
print(p.gcd(q))          # x-1
print(p // q)            # x^2+x-1
print(p % q)             # 0

# Powers
r = UniPoly.from_symbolic(x + pool.integer(1), x)
print(r ** 3)            # x^3+3*x^2+3*x+1

degree, is_zero and leading_coeff are properties (zero-argument O(1) accessors); coefficients() is a method because it allocates a list. Note that coefficients() is i64 and truncates coefficients that do not fit, whereas leading_coeff is exact at any size.

UniPoly is the right choice when you are doing heavy univariate polynomial arithmetic (GCD chains, resultants, factorization) because FLINT applies highly optimized algorithms with exact arithmetic.

MultiPoly

Sparse multivariate polynomial over ℤ (integers). Terms are stored as a map from exponent vectors to coefficients.

from alkahest import MultiPoly

# x^2*y + x*y^2 - 1
expr = x**2 * y + x * y**2 + pool.integer(-1)
mp = MultiPoly.from_symbolic(expr, [x, y])

print(mp.total_degree)       # 3
print(mp.integer_content())  # 1

# Arithmetic
mp2 = MultiPoly.from_symbolic(x * y, [x, y])
print(mp + mp2)              # -1 + xy + xy^2 + x^2y
print(mp * mp2)              # -xy + x^2y^3 + x^3y^2

Terms print in the polynomial’s internal (ascending exponent-vector) order, not in descending degree.

Variable order matters for the exponent-vector key. Pass variables in a consistent order when constructing MultiPoly objects that will be combined.

MultiPolyFp

Sparse multivariate polynomial over 𝔽ₚ = ℤ/pℤ. Used as the working type for black-box sparse interpolation and sparse modular GCD.

from alkahest import sparse_interp_univariate, sparse_interp, gcd_sparse, MultiPoly

p = 32749  # prime

# Recover a sparse univariate from 2·term_bound black-box evaluations
# (Ben-Or/Tiwari).  Returns a list of (coefficient, exponent) pairs.
f = sparse_interp_univariate(lambda v: (v**5 + 3*v**3 + 7) % p, term_bound=3, prime=p)
print(f)   # [(7, 0), (3, 3), (1, 5)]   i.e. x^5 + 3*x^3 + 7

# Recover a sparse multivariate via Zippel's algorithm — this one is a MultiPolyFp,
# printed over positional variables x0, x1, ... in the order given by `vars`
f2 = sparse_interp(
    lambda vals: (vals[0]**3 * vals[1]**2 + vals[0] * vals[1]**4) % p,
    vars=[x, y], term_bound=2, degree_bound=5, prime=p,
)
print(f2)  # 1*x0*x1^4 + 1*x0^3*x1^2 (mod 32749)

# Sparse modular GCD over ℤ[x₁,...,xₙ] — substrate for exact GCD algorithms
a = MultiPoly.from_symbolic((x + y) * (x - y), [x, y])
b = MultiPoly.from_symbolic((x + y) * (x + pool.integer(1)), [x, y])
h = gcd_sparse(a, b, term_bound=4, degree_bound=4)
print(h)   # y + x

sparse_interp_univariate uses Berlekamp–Massey + BSGS discrete-log + Vandermonde solve and requires exactly 2 * term_bound oracle calls. sparse_interp uses Zippel’s variable-by-variable algorithm with batched Vandermonde lifting.

RationalFunction

Quotient of two MultiPoly objects, automatically reduced by their GCD.

from alkahest import RationalFunction

# (x^2 - 1) / (x - 1) → normalized to x + 1
numer = x**2 + pool.integer(-1)
denom = x + pool.integer(-1)
rf = RationalFunction.from_symbolic(numer, denom, [x])
print(rf)   # 1 + x

# Arithmetic preserves the rational form
rf_x = RationalFunction.from_symbolic(x, pool.integer(1), [x])
rf_inv = RationalFunction.from_symbolic(pool.integer(1), x, [x])
print(rf_x + rf_inv)   # (1 + x^2) / (x)

GCD normalization runs at construction, so every RationalFunction is in lowest terms.

ArbBall

A real interval [midpoint ± radius] backed by FLINT’s Arb library. Arithmetic on ArbBall values produces guaranteed enclosures of the true result.

from alkahest import ArbBall, ExprPool, interval_eval, sin

# ArbBall(midpoint, radius, precision_bits=53)
a = ArbBall(2.0, 0.5)    # [1.5, 2.5]
b = ArbBall(3.0, 0.0)    # exactly 3

# A ball prints as midpoint ± radius, not as an interval
print(a + b)   # ArbBall(5.000000 ± 5.00e-1)   i.e. [4.5, 5.5]
print(a * b)   # ArbBall(6.000000 ± 1.50e0)    i.e. [4.5, 7.5]

# Evaluate a symbolic expression rigorously
pool = ExprPool()
x = pool.symbol("x")
result = interval_eval(sin(x), {x: ArbBall(1.0, 1e-10)})
print(result.lo, result.hi)   # tight enclosure of sin(1)

The output ball is guaranteed to contain the true value for any input in the input balls. This is useful for:

• Certified numerical evaluation
• Proving bounds on symbolic expressions
• Verification workflows alongside Lean certificate export

See Ball arithmetic for more detail.

Converting back to Expr

Conversion back to Expr is per-type, not universal. GbPoly — the Gröbner representation — round-trips:

from alkahest import expr_to_gbpoly, diff

g = expr_to_gbpoly(x**2 + pool.integer(1), [x])
expr_again = g.to_expr()   # (x^2 + 1)
dr = diff(expr_again, x)   # DerivedResult(value=(x * 2))

UniPoly, MultiPoly and RationalFunction do not currently expose a symbolic conversion. Read their coefficients (UniPoly.coefficients(), MultiPolyFp.terms) or keep the original Expr alongside the polynomial — from_symbolic does not consume it.

Parsing expressions from strings

alkahest.parse converts a human-readable math string into an Expr node using a Pratt (top-down operator precedence) recursive-descent parser.

import alkahest
from alkahest import ExprPool, parse, diff, simplify

pool = ExprPool()
x = pool.symbol("x")

e = parse("x^2 + 2*x + 1", pool, {"x": x})
print(e)                    # x^2 + 2*x + 1

dr = diff(e, x)
print(dr.value)             # 2*x + 2

Syntax

Whitespace (spaces, tabs, newlines) is ignored everywhere.

Operator precedence

From lowest to highest:

So -x^2 parses as -(x^2), not (-x)^2, and x^2^3 parses as x^(2^3) = x^8.

Supported functions

abs, acos, asin, atan, atan2, ceil, cos, cosh, erf, erfc, exp, floor, gamma, log, round, sign, sin, sinh, sqrt, tan, tanh

The symbols map

By default, every new identifier is interned as a fresh pool.symbol(name). Pass a pre-built symbols dict to bind identifiers to existing Expr objects, or to collect the symbols that were created:

# Pre-bind x to an existing symbol
x = pool.symbol("x")
e = parse("sin(x)^2 + cos(x)^2", pool, {"x": x})

# Collect auto-created symbols after parsing
sym_map: dict = {}
e = parse("a*x^2 + b*x + c", pool, sym_map)
print(sym_map.keys())   # dict_keys(['a', 'x', 'b', 'c'])

Identifiers not in the map are created and then added to the map, so the same string name always resolves to the same Expr within a single parse call.

Error handling

parse raises ParseError (code E-PARSE-001) on any lexical or syntax error. The exception’s .span attribute gives the (start, end) byte range of the offending token, and .remediation provides a hint:

from alkahest import ParseError

try:
    parse("sin(x) @ 2", pool, {"x": x})
except ParseError as e:
    print(e)           # unexpected character '@' at offset 7
    print(e.span)      # (7, 8)

try:
    parse("zeta(x)", pool, {"x": x})
except ParseError as e:
    print(e.remediation)  # known functions: abs, acos, asin, ...

Round-trip with pretty-printing

parse is the inverse of str() for expressions built from the operators and functions listed above:

from alkahest import latex, unicode_str

e = parse("sin(x)^2 + cos(x)^2", pool, {"x": x})
print(latex(e))        # \sin\!\left(x\right)^{2} + \cos\!\left(x\right)^{2}
print(unicode_str(e))  # sin(x)² + cos(x)²

Simplification

Alkahest provides two complementary simplification engines that operate on the same expression pool.

Rule-based simplification

simplify applies a fixed set of algebraic rewrite rules until no more apply (fixpoint). It is fast, predictable, and always terminates.

from alkahest import simplify

r = simplify(x + pool.integer(0))   # → x
r = simplify(x * pool.integer(1))   # → x
r = simplify(pool.integer(2) * pool.integer(3))  # → 6  (constant folding)

The default rule set covers:

• Identity and absorbing elements (x + 0 → x, x * 1 → x, x * 0 → 0)
• Constant folding (integer and rational arithmetic)
• Basic polynomial simplification (x + x → 2*x, x² * x → x³)
• Commutativity and associativity (normalized at construction)

Domain-specific rule sets

from alkahest import simplify_trig, simplify_log_exp, simplify_expanded

# Pythagorean identity and double-angle formulas
r = simplify_trig(sin(x)**2 + cos(x)**2)  # → 1

# Conservatively leaves branch-sensitive identities unchanged
r = simplify_log_exp(exp(log(x)))

# Expand products and collect like terms
r = simplify_expanded((x + pool.integer(1))**3)

Customizing the rule set

from alkahest import simplify_with, make_rule

# Add a custom rule: sin²(x) → 1 - cos²(x)
my_rule = make_rule("sin_sq_to_cos", lhs=sin(x)**2, rhs=pool.integer(1) - cos(x)**2)
r = simplify_with(expr, rules=[my_rule])

Conditional simplification (colored e-graphs)

Branch-sensitive rewrites are opt-in. In Python, create an explicit experimental context tied to one expression pool; x > 0 then enables sqrt(x²) → x:

from alkahest.experimental import Assumptions

assumptions = Assumptions(pool)
assumptions.refine(pool.gt(x, pool.integer(0)))
r = assumptions.simplify(sqrt(x**2))  # → x

The current fact language recognizes conjunctions of positive and non-zero predicates. Unsupported predicates are retained for contradiction detection but do not authorize a rewrite; a definitive contradiction raises AssumptionError with code E-SIMPLIFY-001. Contexts never modify the thread-local context() helper or global simplifier state.

Without a proven fact, branch-cut identities such as sqrt(x²) → x, exp(log(x)) → x, and log(a*b) → log(a) + log(b) remain unchanged. Algebraic cancelations such as x/x → 1 and x^0 → 1 still fire in the default simplifier and record NonZero side conditions in the derivation log. The colored pass runs after ordinary rule simplification and preserves repeated terms and factors.

The literal-zero carve-out

b · b⁻¹ → 1 is a documented convention for a symbolic base, but it is not a convention anybody holds when the base is the literal integer 0: 0⁻¹ is division by zero, so 0 · 0⁻¹ is the indeterminate form 0 · ∞ and has no value. Through 3.7 simplify returned 1 for it, simplify_egraph returned 0, and simplify(5 · 0⁻¹ · 0) returned 0 — three answers that were their own proof that at least two were wrong. As of 3.8 all of them decline:

import alkahest as ak

pool = ak.ExprPool()
Z = pool.integer(0)
undefined = Z * Z ** pool.integer(-1)

ak.simplify(undefined).value          # (0 * 0^-1) — left alone
ak.simplify_egraph(undefined).value   # (0 * 0^-1) — same
ak.simplify(pool.integer(5) * Z ** pool.integer(-1) * Z).value   # (0 * 0^-1 * 5)

The guard tests for a literal zero base. Because the rule engine normalises strictly bottom-up, that also covers every base the simplifier can itself reduce to zero — x - x included, so diff(2/(x - x), x) no longer reports 1 for a function whose domain is empty. A base that is zero but not provably so keeps the b · b⁻¹ → 1 convention: deciding it would need a three-valued zero test on the Mul rewrite path, which costs several 128-bit ball evaluations per node.

Unaffected controls, so you can see the boundary: x · x⁻¹ → 1, 0 · x → 0, and 2x − 2x → 0 all still fire.

Parallel simplification

from alkahest import simplify_auto, simplify_par, simplify_redex, simplify_strategy

simplify_par(expr)      # fork-join: best on wide expressions
simplify_redex(expr)    # level-scheduled: best on deep ones, deterministic log
simplify_auto(expr)     # picks one of the two from the expression's shape
simplify_strategy(expr) # "fork_join" | "level_scheduled" | "sequential"

Each takes a single expression and returns the same result as simplify; only the schedule differs. Published wheels are built with --features parallel, so all three schedulers are real there. parallel is not a Cargo default, though: a source build that omits it falls back to the sequential path, so all four remain callable and simplify_strategy then reports "sequential" — which is the way to tell the two situations apart.

Neither parallel scheduler dominates:

Measured on 32 cores, best time over 1–32 threads:

Fork-join keeps each chain on one worker and wins on wide expressions through cache locality. Level scheduling wins on deep ones, where fork-join finds no wide node to fork on and runs essentially sequentially. At one thread the level scheduler is faster on every shape measured. Reproduce with cargo run --release --features parallel --example simplify_three_way.

Both are experimental: alkahest_cas::experimental::{simplify_par, simplify_redex}.

E-graph simplification

simplify_egraph uses equality saturation via egglog to explore many equivalent forms simultaneously before committing to the best one via a cost function.

from alkahest import simplify_egraph

# The e-graph can discover non-obvious equivalences
r = simplify_egraph(x * x - pool.integer(1))  # may factor or simplify

E-graph saturation is more powerful than rule-based simplification for some inputs but slower and has non-deterministic performance for complex expressions. See E-graph saturation for configuration options (EgraphConfig, including disjoint_schedule for match-disjoint rule groups).

User-defined PatternRule sets can use discrimination-net indexing on the Rust side (PatternRuleSet, simplify_with_pattern_rules); the built-in RewriteRule engine still scans rules linearly.

Choosing between the two

For most workflows: use simplify (or a domain-specific variant) first. Reach for simplify_egraph when you need the system to discover a non-obvious equivalence.

Collect and normalize

Two utility passes that sit between the two engines:

from alkahest import collect_like_terms, poly_normal

# 2*x + 3*x → 5*x
r = collect_like_terms(pool.integer(2) * x + pool.integer(3) * x)

# Normalize to canonical polynomial form over given variables
r = poly_normal(x**2 + pool.integer(2) * x * y + y**2, [x, y])

Complex constructors (experimental)

alkahest.experimental exposes symbolic conjugate, re, im, and principal arg. These are symbolic-only: they are not registered for f64, ball, or JIT evaluation.

Safe simplifications include involution of conjugate, real/integer literals for re/im, arg of a strictly positive literal or Domain.Positive symbol, and exact arg(I) / arg(-I). Branch-sensitive cases such as conjugate(log(z)), arg(0), negative reals, and generic complex inputs stay unevaluated. Principal Arg uses the conventional range (−π, π] with a cut on the negative real axis; do not rewrite through atan2, log, or sqrt yet.

Rule engine

The rule engine underlies both simplify and the e-graph backend. Rules are the atomic units of algebraic knowledge.

Anatomy of a rule

A RewriteRule has:

• A name — stable string identifier (used in derivation logs and Lean certificate output)
• A LHS pattern — an expression template with pattern variables
• A RHS template — the replacement
• Optional side conditions — predicates that must hold for the rule to fire

Pattern syntax

Patterns are regular expressions with a subset of ExprData nodes used as wildcards. From the Python side, pattern variables are Expr objects whose names start with ?:

from alkahest import make_rule, match_pattern

pool = ExprPool()
x = pool.symbol("x")

# Pattern variable — matches any subexpression
pv = pool.symbol("?a")

# Rule: ?a + 0 → ?a
add_zero = make_rule("add_zero", lhs=pv + pool.integer(0), rhs=pv)

Pattern variables capture any subexpression and must bind consistently: if ?a appears twice in the LHS it must match the same expression in both positions.

Matching

match_pattern applies a pattern to an expression and returns all match substitutions:

matches = match_pattern(sin(x)**2 + cos(x)**2, pattern)
for subst in matches:
    print(subst)  # dict mapping pattern variable → matched expr

The matcher is associative-commutative (AC): a + b matches b + a, and a + b + c matches any ordering.

Built-in rule sets

The rule sets loaded by simplify and the domain-specific simplifiers are:

Defining custom rules

from alkahest import make_rule, simplify_with

pool = ExprPool()
x = pool.symbol("x")
a = pool.symbol("?a")
b = pool.symbol("?b")

# Commutativity of subtraction rewrite: a - a → 0
self_cancel = make_rule(
    "self_cancel",
    lhs=a + pool.integer(-1) * a,
    rhs=pool.integer(0),
)

# Apply the custom rule alongside the default set
r = simplify_with(expr, rules=[self_cancel])

Custom rules are recorded in derivation logs with the name you provide. If you tag the rule with a Lean theorem name (via the Rust PrimitiveRegistry API), the corresponding step can be exported as a Lean proof term.

Rule execution model

simplify applies rules in a fixpoint loop:

1. For each node in the expression (post-order traversal):
   • Try each rule in the rule set.
   • If a rule matches and its side conditions are satisfied, apply it, emit a RewriteStep, and restart the loop for the modified subtree.
2. Repeat until no rule fires in a full pass.

This is an inner-outer loop strategy rather than exhaustive bottom-up application. It is fast but not complete — some sequences of rewrites require rules to be applied in a specific order. The e-graph engine removes this ordering dependency.

Side conditions

A rule can carry a side condition checked against the matched substitution:

# sqrt(x^2) → x  only when x is non-negative
sqrt_sq = make_rule(
    "sqrt_sq_nonneg",
    lhs=sqrt(a**2),
    rhs=a,
    condition="nonnegative",   # checked against the domain of ?a
)

Side conditions that reference symbol domains are sound: sqrt_sq_nonneg will only fire when ?a is bound to a symbol with domain positive or nonneg. They propagate into the derivation log as SideCondition entries and into Lean output as assumptions.

E-graph saturation

The e-graph backend exposes a fundamentally different approach to simplification: rather than applying rules one at a time in a fixed order, it builds a structure that represents many equivalent expressions simultaneously, then extracts the best one.

What is an e-graph?

An e-graph partitions expressions into equivalence classes (e-classes). When a rewrite rule fires, it does not replace the LHS — it adds the RHS to the same e-class as the LHS. At the end of saturation, an extraction step picks the cheapest representative from each e-class according to a cost function.

This eliminates the phase-ordering problem: rules can fire in any order without risk of committing to a suboptimal form. The e-graph remembers all explored forms and chooses among them at the end.

Using the e-graph

from alkahest import simplify_egraph, simplify_egraph_with

# Default configuration
r = simplify_egraph(expr)

# With explicit config
from alkahest import EgraphConfig, simplify_egraph_with

cfg = EgraphConfig(node_limit=10_000, iter_limit=20)
r = simplify_egraph_with(expr, cfg)

Cost functions

The extraction step minimizes a cost function over e-class representatives. Three built-in cost functions:

StabilityCost is aware of numerical stability issues: it penalizes subtractive cancellation patterns and prefers numerically stable rearrangements.

Configuration

The e-graph runs until saturation (no new e-class merges) or until a limit is hit:

• node_limit — maximum number of e-nodes. Once reached, saturation stops and extraction runs on the current state.
• iter_limit — maximum number of saturation rounds.

For large or complex expressions, saturation can be expensive. The rule-based simplify is often sufficient and should be preferred on hot paths.

Rule sets in the e-graph

The e-graph uses the same RewriteRule objects as the rule-based engine. By default it loads the arithmetic rules. Domain-specific rules (trig, log/exp) are kept separate to avoid e-class explosions on expressions that do not involve those operations.

The default e-graph rule set includes trig identities (sin²+cos²→1) and safe log/exp cancellation. Disable per domain via EgraphConfig(include_trig_rules=False) or include_log_exp_rules=False.

Match-disjoint scheduling

By default (disjoint_schedule=True), shrink and explore rules are split into match-disjoint egglog rulesets (shrink-add, shrink-mul, shrink-pow, explore-trig, explore-log, explore-mul) and run as separate (run …) steps within each phase. This reduces cross-rule interference during saturation. Set disjoint_schedule=False to use the legacy single-ruleset schedule.

When e-graphs help

The e-graph is especially powerful when:

• Multiple non-obvious rewrites must be combined in a specific order that is hard to predict.
• The “right” form is not syntactically similar to the input (e.g. factoring followed by cancellation).
• You want the globally cheapest form under a custom cost function, not just any simplified form.

It is less useful when:

• The expression is already in near-canonical form and only identity cleanup is needed.
• You need predictable performance on a hot path.
• The expression is large and associative-commutative, where the e-graph can grow combinatorially.

Colored e-graphs (conditional rewrites)

Separate from the egglog saturation backend, Alkahest implements a native colored e-graph (simplify/colored_egraph.rs) for conditional simplification under explicit assumptions (e.g. x > 0 ⊢ sqrt(x²) → x). When SimplifyConfig::assumptions is non-empty, simplify_with runs this pass before the rule engine. See Simplification — conditional simplification.

AC matching in the e-graph

The egglog backend handles associativity and commutativity structurally: Add and Mul children are sorted at pool-insertion time, so there is a single canonical ordering. The e-graph does not need to enumerate permutations.

This is more efficient than classical AC-completion but requires that the canonical ordering is established at construction, which the kernel enforces.

Calculus

Alkahest supports symbolic differentiation and integration with full derivation logging.

Differentiation

diff(expr, var) computes the symbolic derivative of expr with respect to var.

from alkahest import diff, sin, cos, exp, log

pool = ExprPool()
x = pool.symbol("x")

# Polynomial
dr = diff(x**3 + pool.integer(2) * x, x)
print(dr.value)   # 3*x^2 + 2

# Chain rule
dr = diff(sin(x**2), x)
print(dr.value)   # 2*x*cos(x^2)

# Product rule
dr = diff(x * exp(x), x)
print(dr.value)   # exp(x) + x*exp(x)

# Logarithm
dr = diff(log(x**2 + pool.integer(1)), x)
print(dr.value)   # 2*x / (x^2 + 1)

Registered primitives

Every primitive in the registry has a differentiation rule. The 23 currently registered primitives include:

sin, cos, tan, asin, acos, atan, atan2, sinh, cosh, tanh, exp, log, sqrt, abs, sign, erf, erfc, gamma, floor, ceil, round, min, max

Derivation log

The DerivedResult returned by diff records every rule application:

dr = diff(sin(x**2), x)
for step in dr.steps:
    print(f"  {step['rule']:25s}  {step['before']}  →  {step['after']}")

Forward-mode automatic differentiation

diff_forward computes the derivative using forward-mode AD (dual numbers). It produces the same result as diff but through a different computational path:

from alkahest import diff, diff_forward

sym = diff(x**3, x)
fwd = diff_forward(x**3, x)
# fwd.value == sym.value

Forward mode is useful for checking that the symbolic rules agree with dual-number evaluation.

Symbolic gradient (symbolic_grad)

symbolic_grad(expr, vars) returns a list of Expr — one partial derivative per variable. It does not use @trace and is not composable with jit directly.

from alkahest import symbolic_grad

pool = ExprPool()
x = pool.symbol("x")
y = pool.symbol("y")

expr = x**2 * y + sin(x * y)
grads = symbolic_grad(expr, [x, y])
# grads[0] = ∂/∂x = 2*x*y + y*cos(x*y)
# grads[1] = ∂/∂y = x^2 + x*cos(x*y)

For the JAX-style gradient of a traced Python function (compose with jit), use alkahest.grad — not symbolic_grad. See Transformations.

Integration

integrate(expr, var) computes the symbolic antiderivative of expr with respect to var.

from alkahest import integrate, sin, cos, exp

# Polynomials
r = integrate(x**3, x)
print(r.value)    # x^4/4

# Known functions
r = integrate(sin(x), x)
print(r.value)    # -cos(x)

r = integrate(exp(x), x)
print(r.value)    # exp(x)

r = integrate(x**pool.integer(-1), x)
print(r.value)    # log(x)

Integration rules

The integration engine applies a rule table for common forms, then escalates to the Risch decision procedure for harder cases:

Rule table (fast path)

• Power rule: ∫ xⁿ dx = xⁿ⁺¹/(n+1) for integer n ≠ -1
• Logarithm: ∫ 1/x dx = log(x)
• Exponential tower: ∫ exp(a*x + b) dx, ∫ xⁿ·exp(x) dx (poly × exp)
• Linear substitution: ∫ f(a*x + b) dx
• Trigonometric: ∫ sin(x) dx, ∫ cos(x) dx, etc.
• Standard table entries for erf, inverse trig, etc.

Risch algorithm (escalation)

• Rational functions A(x)/D(x): Hermite reduction (repeated factors → rational part), then Rothstein–Trager (rational residues → log), irreducible quadratics (negative discriminant → arctan; positive discriminant → log with √Δ coefficients), and irreducible factors of degree ≥ 3 via a RootSum node (Lazard–Rioboo–Trager over the number field ℚ[t]/Q(t)).
• Exp tower with rational coefficient: ∫ f(x)·exp(η) dx where f ∈ ℚ(x) — solved via the rational Risch DE (Bronstein §6.1).
• Polynomial × exp / log towers: poly-RDE and known-table rules.

# Rational functions
r = integrate(pool.integer(1) / (x**2 - pool.integer(1)), x)  # → ½·log((x-1)/(x+1))

r = integrate(pool.integer(1) / (x**2 + pool.integer(1)), x)  # → arctan(x)

# Rational coefficient × exp
r = integrate((x - pool.integer(1)) / x**2 * exp(x), x)       # → exp(x)/x

# Degree-≥3 denominator → RootSum
r = integrate(pool.integer(1) / (x**3 - pool.integer(3)*x + pool.integer(1)), x)
# r.value contains a RootSum node (sum over algebraic residues)

Non-elementary certification: when the integrand is provably non-elementary (Liouville’s theorem — e.g. sin(x)/x, exp(x)/x, exp(x²)), integrate raises IntegrationError with code E-INT-004 (NonElementary) rather than a generic “not implemented”:

from alkahest import IntegrationError

try:
    integrate(exp(x) / x, x)
except IntegrationError as e:
    print(e.code)         # E-INT-004
    print(e.remediation)  # "no elementary antiderivative (NonElementary)"

For integrands outside the supported classes (e.g. sqrt(P(x)), mixed algebraic+transcendental), integrate raises IntegrationError with code E-INT-001 (NotImplemented).

Verification

A common pattern is to verify an antiderivative by differentiating it back:

antideriv = integrate(expr, x).value
check = simplify(diff(antideriv, x).value)
# check.value should equal expr

Higher derivatives

Chain calls to diff:

d2 = diff(diff(sin(x), x).value, x)
print(d2.value)   # -sin(x)

The derivation log of the outer diff does not include the inner steps. If you need the full trace, concatenate dr1.steps + dr2.steps.

Asymptotics at scale

series and Gruntz limit expand a function; asymptotic_expand handles the power/log/exp scales of a function at infinity. What a loop working on combinatorics and analysis needs on top of that is the asymptotics of a sum — how Σ_{k=a}^{n} f(k) behaves as n → ∞. That is how a conjecture about a growth rate actually gets settled.

from alkahest import ExprPool
from alkahest.experimental import euler_maclaurin

pool = ExprPool()
k, n = pool.symbol("k"), pool.symbol("n")

r = euler_maclaurin(pool.integer(1) / k, k, 1, n, corrections=2)
r.leading           # log(n)
r.terms            # [log(n), γ, 1/(2n), −1/(12n²)] — most significant first
r.rigor            # "numerically_consistent"

The formula

For a smooth summand,

Σ_{k=a}^{n} f(k) = ∫_a^n f(t) dt + (f(a) + f(n))/2
                   + Σ_{j=1}^{m} B_{2j}/(2j)! · (f^{(2j-1)}(n) − f^{(2j-1)}(a))
                   + R_m

Read as n → ∞, the a-endpoint pieces collapse into one additive constant and what remains is an asymptotic expansion in n. For f(k) = 1/k that is the classical H_n ~ log n + γ + 1/(2n) − 1/(12n²) + ….

The constant is fitted, and the report says so

Euler–Maclaurin does not determine the additive constant from the n-side terms. For the harmonic numbers that constant is Euler’s γ, and no amount of boundary algebra at k = a produces it.

This module therefore obtains it numerically, from the exactly computed sum at large n, and is explicit about that: rigor comes back as "numerically_consistent" rather than "proved", and the fitted constant appears in hypotheses with status "assumed". The shape of the expansion is derived symbolically — only that one scalar is empirical.

r.all_hypotheses_checked      # False
[s for s, _ in r.hypotheses]  # ["checked", "assumed", "assumed"]

This distinction is the whole reason the result is an AsymptoticReport rather than a bare list of terms. A loop that records “γ ≈ 0.5772 was fitted, not proved” can revisit it; one handed an undifferentiated expansion cannot.

Ordering and the verification gate

Terms are returned most-significant first, and the ordering is by magnitude at the check points, not by where they came from in the formula. This matters: the constant sits below every growing term and above every decaying one, so for Σ k it lands after n²/2 and n/2, while for H_n it lands right after log n.

Every term then goes through the same o()-gate as asymptotic_expand: the truncated expansion is compared against the exactly computed sum at increasing n, and a term is kept only if it genuinely refines its predecessor against the true value. Terms that fail are dropped; if nothing survives, the call refuses rather than emitting an unverified expansion. r.verification carries the evidence and r.max_relative_error summarises it.

Refusals

Singularity analysis of generating functions

For a rational f(z), the growth of [zⁿ] f(z) is governed entirely by the singularity of smallest modulus:

from alkahest.experimental import coefficient_asymptotics

gf = pool.integer(1) / (pool.integer(1) - z - z*z)   # Fibonacci
r = coefficient_asymptotics(gf, z, n)
r.terms[0]      # ~ C·φⁿ

A pole ρ of multiplicity m contributes C · n^{m-1} · ρ^{-n}. The shape — which power of n, which exponential base — is exact; the single leading constant is obtained from the exact power series by Richardson extrapolation, and the report says so.

That extrapolation is not a nicety. Reading the constant off a single finite index absorbs the subleading term: for 1/(1-z)², where [zⁿ] = n+1 exactly, taking C = a₄₈/48 = 49/48 leaves a permanent 2% bias, so the relative error stops shrinking as n grows — which is precisely what an asymptotic statement must not do. One Richardson step cancels the 1/k term and recovers C = 1.

What it refuses

The transfer theorem needs a unique dominant singularity. 1/(1-z²) has poles at both +1 and −1; its coefficients oscillate (1,0,1,0,…) and no single power-law term describes them. Rather than reporting one pole as if it won, the routine declines — as it does for a complex dominant pole (necessarily one of a conjugate pair) and for non-rational input.

From a recurrence: Poincaré–Perron

A certified recurrence already determines how fast its sequence grows, so after zeilberger or guess_holonomic the growth law is one call away:

from alkahest import ExprPool
from alkahest.experimental import asymptotics_from_recurrence

pool = ExprPool()
n = pool.symbol("n")

# (n+1)·u(n+1) − (4n+2)·u(n) = 0 — the central binomial coefficients.
r = asymptotics_from_recurrence([(-2, -4), (1, 1)], n, terms=[1])

r.growth_rate_exact          # 4          — derived
r.polynomial_exponent_exact  # -1/2       — derived
r.connection_constant        # 0.5641895… — fitted; this is 1/√π
r.verdict                    # "single_dominant_root"

rec may be a ZeilbergerCertificate, a GuessedRecurrence, or a plain list of coefficient polynomials [p_0, …, p_J] for Σ_i p_i(n)·u(n+i) = 0; each p_i is an Expr in n or a tuple of ascending integer coefficients.

What is derived and what is fitted

Write D = max_i deg p_i, take the coefficient of n^D in each p_i to build the characteristic polynomial χ(t) = Σ_i a_i tⁱ, and the coefficient of n^{D-1} to build χ₁. Poincaré’s theorem says the sequence grows like a root of χ; Perron’s refinement pins the polynomial factor:

u(n) ~ C · ρⁿ · n^α,     α = −χ₁(ρ) / (ρ · χ'(ρ))

ρ and α are functions of the recurrence and of nothing else, and when ρ is rational both are available exactly. C is not. It is determined by the initial conditions, so it is extrapolated numerically from the exact terms — run forward from the recurrence in exact rational arithmetic — and reported on its own, with connection_constant_converged and connection_constant_drift from a second extrapolation over a smaller range of indices. evidence() splits the two:

r.evidence()["derived"]["growth_rate"]         # 4.0
r.evidence()["fitted"]["connection_constant"]  # 0.5641895…
r.evidence()["fitted"]["relative_drift"]       # 2.5e-10

This is the same discipline euler_maclaurin applies to its additive constant, for the same reason: no amount of algebra on the recurrence produces 1/√π.

What it refuses to answer

The theorem needs the roots of χ to have distinct moduli and the leading coefficient to be eventually non-zero. When they do not hold, verdict says so and growth_rate is None — a growth rate is never invented:

Multiplicity is exact — it comes from the squarefree decomposition of χ over ℚ, not from clustering the numeric roots. That is not fussiness: A359643’s characteristic polynomial is (t−1)³·(27t−283), whose triple root is real and sits well below the dominant one, and a tolerance that merged them would refuse a case the theory handles perfectly.

One more hypothesis is easy to miss. Poincaré’s conclusion is that u(n+1)/u(n) tends to some root, not necessarily the largest. u(n+2) = 3u(n+1) − 2u(n) with u(0) = u(1) = 1 is the constant sequence, and its component along the dominant root 2 is zero. With terms supplied that is detected and reported as follows_dominant_root == False; without them it is an explicitly assumed hypothesis in report().hypotheses.

Scope of this release

Shipped: the Euler–Maclaurin route for Σ_{k=a}^{n} f(k), with Bernoulli corrections, magnitude ordering, the numeric gate, and the checked-versus- assumed hypothesis ledger in AsymptoticReport.

Also shipped: singularity analysis for rational generating functions, and Poincaré–Perron growth from a P-recursive recurrence (asymptotics_from_recurrence).

Not shipped, and tracked as follow-up (the shared scaffolding — AsymptoticReport, the gate, exact Bernoulli numbers, rational-function extraction and a complex root finder — is already in place for them):

• Algebraic and log-type generating functions — the transfer theorem beyond poles (√(1-4z) for the Catalan numbers, log singularities). Only the rational case ships here.
• Sequence asymptotics from a closed form — Stirling-based expansions of a ratio of factorials. The recurrence route ships (above); the closed-form one does not.
• Full Birkhoff–Trjitzinsky asymptotics — the cases asymptotics_from_recurrence reports and declines to answer: equal-modulus roots, a repeated dominant root, and the degenerate leading coefficient that produces ρⁿ·n^{cn} growth.
• Laplace / saddle-point / stationary-phase asymptotics of parameter integrals.

Note that Σ log k gives Stirling’s formula for log n! through this route already, so the factorial case is reachable today.

Transformations

Transformations are higher-order operations that take a function and return a new function. They compose freely and operate on traced symbolic representations.

Tracing

trace symbolically executes a Python function by replacing its inputs with symbolic variables and recording the computation as an expression DAG.

import alkahest
from alkahest import ExprPool

pool = ExprPool()

@alkahest.trace(pool)
def f(x, y):
    return x**2 + alkahest.sin(y)

print(f.expr)    # x^2 + sin(y)
print(f.symbols) # [x, y]

The decorator takes the pool as an argument. Variable names are inferred from the function signature.

Numeric evaluation

TracedFn objects are callable with numeric values:

print(f(3.0, 0.0))   # 9.0

import numpy as np
xs = np.linspace(0, 1, 1000)
ys = np.zeros(1000)
result = f(xs, ys)   # vectorised automatically

Gradient (grad, not symbolic_grad)

alkahest.grad applies to a TracedFn from @trace. For partial derivatives of a bare Expr, use symbolic_grad instead.

grad differentiates a traced function symbolically with respect to all (or a subset of) its inputs:

df = alkahest.grad(f)
# df(x, y) returns [∂f/∂x, ∂f/∂y] = [2*x, cos(y)]

grads = df(1.0, 0.0)   # [2.0, 1.0]

Differentiate with respect to a subset:

df_x = alkahest.grad(f, wrt=[f.symbols[0]])  # ∂f/∂x only

JIT compilation

jit wraps a traced function in the LLVM JIT backend. The first call triggers compilation; subsequent calls run the compiled code directly.

fast_f = alkahest.jit(f)
print(fast_f(3.0, 0.0))  # 9.0, via LLVM-compiled code

Vectorised evaluation is automatic when array inputs are detected:

xs = np.linspace(0, 10, 1_000_000)
ys = np.zeros_like(xs)
result = fast_f(xs, ys)   # zero-copy batch path

Composing transformations

Transformations stack:

# Compiled gradient
fast_df = alkahest.jit(alkahest.grad(f))
grads = fast_df(xs, ys)   # compiled, vectorised gradient

# Second derivative: grad of grad
d2f = alkahest.grad(alkahest.grad(f))

Note that grad returns a GradTracedFn, not a TracedFn. jit can be applied to GradTracedFn when it wraps a single scalar output. For multi-output cases, compile each gradient expression individually with compile_expr.

trace_fn

Functional (non-decorator) version of trace:

from alkahest import trace_fn

fn = trace_fn(lambda x, y: x * alkahest.exp(y), pool)

PyTrees

Transformations work over nested data structures (lists, dicts, tuples, dataclasses). The Python layer flattens and unflattens them automatically:

from alkahest import flatten_exprs, unflatten_exprs, map_exprs

# A system of equations as a list
eqs = [x**2 + y - pool.integer(1), x - y**2 + pool.integer(1)]

# Map simplification over the list
simplified = map_exprs(simplify, eqs)

# Flatten to a list of ExprIds and the structure descriptor
flat, treedef = flatten_exprs(eqs)
restored = unflatten_exprs(flat, treedef)

This follows the JAX pytree pattern. The Rust kernel sees only flat sequences; structure reconstruction is a Python-layer concern.

Context manager

alkahest.context sets a default pool and configuration for a block:

with alkahest.context(pool=pool, simplify=True):
    z = alkahest.symbol("z")          # uses the active pool
    expr = z**2 + alkahest.sin(z)     # auto-simplified

Inside the context, alkahest.symbol(name) creates a symbol in the active pool without passing the pool explicitly. This is a convenience wrapper — the pool is still explicit at the structural level.

Code generation

Alkahest can compile symbolic expressions to fast native or GPU code. Compiled code bypasses Python entirely during evaluation.

The compilation pipeline

Expressions lower through multiple IR levels:

ExprPool (hash-consed DAG)
    ↓  e-graph extraction + canonicalization
Canonical expression form
    ↓  alkahest MLIR dialect
High-level MLIR (math-aware ops: horner, poly_eval, interval_eval)
    ↓  lowering passes
Standard MLIR (arith, math, linalg, gpu)
    ↓
LLVM IR / PTX / StableHLO (depending on target)
    ↓
Native machine code / GPU kernel / XLA

The custom alkahest MLIR dialect is where math-aware optimizations happen: Horner’s method for polynomials, fused multiply-add emission, numerically stable rearrangements via StabilityCost.

compile_expr

compile_expr produces a callable from a symbolic expression and a list of input variables:

from alkahest import ExprPool, compile_expr, sin, cos

pool = ExprPool()
x = pool.symbol("x")
y = pool.symbol("y")

f = compile_expr(x**2 + sin(y), [x, y])
print(f([3.0, 0.0]))   # 9.0

The callable takes a list of floats (one per variable) and returns a float. For batch evaluation see numpy_eval below.

CPU compilation uses a three-tier dispatch chosen by CompileConfig (expression size and expected_evals, not feature-flag order):

1. Interpreter — small DAGs and few planned evaluations (zero compile latency)
2. Cranelift (--features cranelift) — pure Rust, ~10× faster compile than LLVM; no system LLVM required
3. LLVM (--features jit) — inkwell / LLVM 15 MCJIT; best throughput for large batch sweeps

Default PyPI wheels use the interpreter only. +jit / +full release wheels enable LLVM. Add cranelift to a from-source build for the fast-compile tier without LLVM.

Native JIT backends also emit a bulk entry point (alkahest_eval_bulk) for column-major batch evaluation; CompiledFn.call_batch / call_bulk use the same layout as numpy_eval.

CompileCache

Repeated compilation of the same expression is expensive. CompileCache memoizes by (ExprId, input variables) — hash-consing makes ExprId a stable content key:

from alkahest import ExprPool, compile_expr, CompileCache

pool = ExprPool()
x = pool.symbol("x")
cache = CompileCache()
f = cache.compile(x**2, [x], pool)   # JIT compiles on first call
g = cache.compile(x**2, [x], pool)   # cache hit — O(1)
print(cache.stats())                 # hits, compiles, hit_rate

eval_expr

For one-off evaluation without compiling:

from alkahest import eval_expr

result = eval_expr(x**2 + sin(y), {x: 3.0, y: 0.0})
print(result)   # 9.0

eval_expr is slower than a compiled function for repeated evaluation but has no compilation overhead.

numpy_eval

numpy_eval vectorises a compiled function over NumPy arrays via the batch path:

import numpy as np
from alkahest import numpy_eval

f = compile_expr(sin(x) * cos(x), [x])
xs = np.linspace(0, 2 * 3.14159, 1_000_000)
ys = numpy_eval(f, xs)   # vectorised, zero-copy

Also accepts PyTorch CPU tensors and JAX arrays via DLPack.

Parallel batch evaluation

With --features parallel, numpy_eval_par distributes evaluation across CPU cores (Rayon), releasing the GIL during computation:

from alkahest import numpy_eval_par

ys = numpy_eval_par(f, xs)   # same API as numpy_eval; multi-core

numpy_eval_par is only parallel when the extension was built with parallel. Every published wheel is — but parallel is not a Cargo default, so a source build that omits it gets a silent no-op alias: numpy_eval_par falls back to numpy_eval with the same correct result, the same single-threaded throughput, and no warning. Timing the two against each other on such a build measures one code path twice.

import alkahest as ak
ak.capabilities()["features"]["parallel"]   # True on PyPI wheels

If that prints False, rebuild with --features parallel — see Getting started.

A CompiledFn belongs to the thread that compiled it

numpy_eval_par fans out internally, so you do not need threads of your own to use it. If you add them anyway, compile inside each thread: the object returned by compile_expr owns JIT code pages and a backend handle, and is pinned to its creating thread. Touching it from another one is refused:

f = compile_expr(expr, [x])          # thread A
threading.Thread(target=lambda: numpy_eval(f, xs)).start()
# pyo3_runtime.PanicException: alkahest::PyCompiledFn is unsendable,
#                              but sent to another thread

Two things to know about that refusal. It has nothing to do with parallel — plain numpy_eval is refused identically. And PanicException derives from BaseException, not Exception, so a worker wrapped in a bare except Exception: will not catch it and the thread will die silently.

ExprPool and Expr have no such restriction: share the pool, build expressions anywhere, and call compile_expr per thread.

Horner-form emission

horner rewrites a polynomial expression into Horner’s form, which is numerically better conditioned and faster to evaluate:

from alkahest import horner

# x^3 + 2x^2 + 3x + 4 → x*(x*(x + 2) + 3) + 4
h = horner(x**3 + pool.integer(2)*x**2 + pool.integer(3)*x + pool.integer(4), x)

emit_c emits a C function string for embedding in other projects:

from alkahest import emit_c

c_code = emit_c(expr, [x, y], fn_name="f")
# → "double f(double x, double y) { return ...; }"

MLIR dialect

The alkahest-mlir crate exposes the custom MLIR dialect. The dialect ops are:

Three lowering targets are available:

• ArithMath — lowers to arith + math MLIR dialects; uses math.fma for Horner chains
• StableHlo — lowers to StableHLO ops for XLA/JAX integration
• Llvm — lowers to llvm dialect for LLVM IR / PTX emission

from alkahest import to_stablehlo

# Emit textual MLIR in the StableHLO dialect
mlir_text = to_stablehlo(expr, [x, y], fn_name="my_fn")
print(mlir_text)  # valid input to mlir-opt / XLA

GPU codegen (NVPTX)

With --features cuda and an LLVM 15 installation with NVPTX support — not in any published wheel, so pip install alkahest never has this. Full detail, including build prerequisites, the supported node set, error codes and the state of testing, is in GPU support (CUDA).

from alkahest import compile_cuda

f_gpu = compile_cuda(expr, [x, y])
result = f_gpu.call_batch(inputs)   # runs on CUDA device 0

Guard on the capability bit before reaching for it, since the name does not exist at all without the feature:

if alkahest.capabilities()["features"]["cuda"]:
    ...

The GPU compiler:

1. Lowers the expression through inkwell to NVPTX LLVM IR for sm_86 (Ampere)
2. Links libdevice.10.bc for transcendental functions (__nv_sin, etc.)
3. Emits PTX via LLVM’s target machine
4. Loads the PTX via the CUDA driver (cudarc)

The benchmark nvptx/nvptx_polynomial_1M shows 16.2× speedup over the CPU JIT on a 1M-point polynomial evaluation on an RTX 3090.

Upcoming (v1.1): AMD ROCm / amdgcn target (hardware-blocked pending RDNA3 availability).

Caching

Use CompileCache for explicit per-session memoization of compiled functions (see above). The persistent ExprPool (V1-14) can serialize expression DAGs across sessions; combine with CompileCache to avoid recompilation after reload.

Tier dispatch always falls back to the interpreter when native JIT features are unavailable or compilation fails.

GPU support (CUDA)

Alkahest has two independent CUDA features. Neither is in the wheel published to PyPI, so pip install alkahest gives you no GPU support at all — a source build is required, and each feature has different build prerequisites.

The two do not imply each other. cuda = ["jit", "dep:cudarc"] and groebner-cuda = ["groebner", "dep:cudarc"] (alkahest-core/Cargo.toml).

Building

# NVPTX expression codegen. Needs LLVM 15 built with NVPTX:
#   llc --version | grep nvptx     # must list nvptx64
maturin develop --manifest-path alkahest-py/Cargo.toml --release --features cuda

# GPU Gröbner kernel (Rust-only; nothing changes at the Python surface)
maturin develop --manifest-path alkahest-py/Cargo.toml --release --features groebner-cuda

# Both
maturin develop --manifest-path alkahest-py/Cargo.toml --release \
    --features "cuda groebner-cuda"

cudarc uses dynamic loading, so the extension builds on a machine with no CUDA installed; the driver is only needed when a kernel actually launches. LLVM 15 with NVPTX, by contrast, is needed at build time for cuda — a build without it fails or produces E-CUDA-001 at compile time.

libdevice.10.bc (from the CUDA toolkit) is linked into every generated module so that sin, cos, … resolve to __nv_*. If it is not found automatically, point at it explicitly:

export ALKAHEST_LIBDEVICE_PATH=/usr/local/cuda/nvvm/libdevice/libdevice.10.bc

What capabilities() reports

import alkahest as ak

features = ak.capabilities()["features"]
features["cuda"]           # `--features cuda` was compiled in
features["llvm_jit"]       # True on any `cuda` build: `cuda` implies the LLVM backend

Read these bits precisely — each says what was linked, and nothing more:

• cuda == True guarantees ak.compile_cuda and ak.CudaCompiledFn exist and that PTX can be emitted on the host. It does not promise a GPU: the driver is loaded lazily, so a machine with no device compiles happily and fails at call_batch with E-CUDA-003. Ask ak.cuda_device_count() — it reports 0 when there is no usable device, and never raises.
• llvm_jit == True on a cuda build even when alkahest-py’s own jit feature was never named, because alkahest-core’s cuda feature turns on jit. Cranelift and LLVM are not mutually exclusive; a CUDA build can link both.

There is no groebner_cuda bit (contract v3 and later — capabilities()["features"] raises KeyError for it). It was removed rather than wired up because it was unfalsifiable: no Python observation distinguished True from False. See below.

ak.CudaError is importable on every build, CUDA or not — it is an exception class, not an entry point, and code that writes except ak.CudaError around a compile step must keep working when moved between wheels. compile_cuda and CudaCompiledFn genuinely do not exist without the feature, and are appended to __all__ only when they do.

compile_cuda

import alkahest as ak

pool = ak.ExprPool()
x, y = pool.symbol("x"), pool.symbol("y")
expr = ak.sin(x) * ak.cos(y) + (x * x + y * y) * pool.rational(1, 100)

fn = ak.compile_cuda(expr, [x, y])   # -> CudaCompiledFn
fn.n_inputs                          # 2
fn.ptx                               # generated PTX assembly (str), `.target sm_86`

out = fn.call_batch([xs, ys])        # list[float], one output per point

call_batch takes one column per symbolic input (structure-of-arrays: xs is every x value, not the first point), all of equal length, and returns a Python list with one float per point. It copies host → device, launches on device 0, and copies back; a mismatched column count or ragged columns raise ValueError before anything touches the GPU. fn.call_batch_on(ordinal, inputs) is the same call on a chosen device.

Pipeline: expression → LLVM IR via inkwell → link libdevice.10.bc → internalize and DCE → PTX for sm_86 (Ampere) → loaded through the CUDA driver by cudarc.

Discovering the valid device ordinals

n = ak.cuda_device_count()          # 0 when there is no GPU here
for dev in range(n):
    fn.call_batch_on(dev, [[0.0]] * fn.n_inputs)

The valid arguments to call_batch_on are exactly range(ak.cuda_device_count()). Both directions are pinned by tests that run on hardware (tests/test_cuda.py::test_device_count_agrees_with_the_ordinals_that_launch and nvptx_gpu::cuda_device_count_matches_the_ordinals_that_launch): every ordinal below the count launches, and the ordinal at the count is refused.

cuda_device_count() never raises. Every “no GPU here” shape — no driver, no device, driver too old — reports 0, because that is the single answer a caller acts on. This matters more than it looks: cudarc panics rather than returning Err when libcuda.so cannot be dlopen’d at all, so a naive binding would abort the process on precisely the machines a capability probe exists to report on. The same bug bit groebner_cuda.rs::gpu_available and nvptx_gpu.rs::device_available.

Earlier releases documented a workaround here — launch on an ordinal and catch E-CUDA-003 — because a cuda_device_count binding could not be verified by anything on an ordinary dev box: cuda implies LLVM 15 with NVPTX, so it could not even be compiled, and no CI job built the Python extension with the feature. It shipped once both could be done on real hardware, which is the standard the capability overclaims this page documents were failing.

Limits worth knowing before you reach for it

• sm_86 is hard-coded. Newer or older architectures rely on the driver’s PTX JIT.
• f64 only, one output value per point. There is no vector or complex return.
• Supported nodes: integer/rational/float constants, +, *, **, and the unary functions sin, cos, tan, exp, log, sqrt, abs. Integer exponents in 0..=16 are unrolled to multiplies; anything else goes through __nv_pow. Any other function — atan, sinh, erf, … — is refused with E-CUDA-002 rather than approximated, as is any symbol you forgot to pass in inputs.
• Host lists in, host list out. The zero-copy device-pointer entry point (call_device_ptrs) exists in the Rust crate only; it has no PyO3 binding, so a CuPy or Torch CUDA tensor is round-tripped through host memory today.

Errors

All are CudaError, a subclass of AlkahestError, each carrying .code and .remediation (see Error handling).

groebner-cuda is not reachable from Python

The feature compiles a real, tested CUDA kernel — MacaulayMatrix::reduce_gpu plus a multi-prime CRT lift — and exports compute_groebner_basis_gpu from the Rust crate. But no shipped code path calls it. GroebnerBasis.compute, solve, and triangularize all go through compute_buchberger_basis on the CPU, and alkahest-py never references the GPU entry point at all. So on a --features groebner-cuda build:

• no Python name appears or disappears,
• no Python call gets faster,
• and, since 3.8, no capability bit claims otherwise.

This is deliberate rather than accidental — the crossover policy in docs/symbolic-gpu-benchmarks.md says production dispatch must not prefer the GPU until the benchmark harness says it wins, and that wiring does not exist yet. Rust users can call alkahest_cas::poly::groebner::compute_groebner_basis_gpu directly.

Why the bit was removed rather than wired up

capabilities()["features"]["groebner_cuda"] used to report True here. It was the only occurrence of the string groebner_cuda anywhere in alkahest-py — there was no binding to go with it, no *gpu* name in the public or the private module, and GroebnerBasis exposing only CPU methods. That made it strictly worse than the cuda overclaim fixed in d139a46, which at least had a private route in.

An unreachable True is the same class of defect as a silent wrong answer: it makes a caller trust something it should not. The two ways out are to add a binding or to drop the claim, and dropping it was the right one days before a release — a binding would have been new public API that no CI job can build (no job builds the Python extension with either CUDA feature; see below) and that nobody without a GPU can run. Adding unverifiable surface is how the original defect got in. The bit is gone; the kernel is unchanged and still Rust-reachable. If dispatch ever prefers the GPU, the binding lands first and a bit follows it.

compute_groebner_basis_gpu now reports where it ran

The Rust entry point falls back to CPU row reduction when no device is present, and it used to say nothing about having done so — a device_id: None run, a run whose driver calls all failed, and a real GPU run returned identical, indistinguishable values. A function named ..._gpu that quietly runs on the CPU is a footgun of exactly the kind this release has spent its time eliminating, so both it and reduce_batch now return a GpuBackendReport alongside the polynomials:

#![allow(unused)]
fn main() {
use alkahest_cas::poly::groebner::{compute_groebner_basis_gpu, MonomialOrder};

let (basis, backend) = compute_groebner_basis_gpu(gens, MonomialOrder::Lex, Some(0))?;
assert!(backend.ran_on_gpu(), "fell back to the CPU: {backend:?}");
}

ran_on_gpu() is true only when at least one mod-p row reduction executed on a device and none fell back; fell_back_to_cpu() is its counterpart; reductions_on_gpu, reductions_on_cpu and first_gpu_error carry the detail. The stderr warning on fallback is still emitted, but it is no longer the only channel. This is a breaking change to the Rust signature — a compile error on upgrade, which is the correct failure mode for a caller who was reading a result as a GPU result.

State of testing — read this before trusting the feature

Rust, on hardware. alkahest-core/tests/nvptx_gpu.rs and alkahest-core/tests/groebner_cuda.rs run under .github/workflows/cuda_nightly.yml on a self-hosted dual-RTX-3090 runner: 17 CUDA-gated tests, plus compute-sanitizer memcheck and racecheck. The last full run was green and both sanitizers clean. --target-processes all is load-bearing in that workflow — without it the sanitizer instruments cargo, a process that makes no CUDA calls, and reports success having inspected nothing.

Python. tests/test_cuda.py covers the binding: the capability/namespace contract, PTX emission, the CudaError refusals, call_batch argument validation, and — the point of the exercise — GPU-versus-CPU numerical agreement on polynomial, transcendental, compile_expr and numpy_eval comparisons. Only the contract tier runs without the feature; everything else skips, which is what happens in CI and on the wheel. Setting ALKAHEST_GPU_TESTS=1 (as the nightly does for Rust) turns those skips into a hard error, so a job that promises hardware cannot quietly report success without reaching it.

The honest gap: no CI job has ever built the Python extension with cuda or groebner-cuda. The nightly runs cargo, never maturin, so the Python tier is only exercised when someone builds with the feature and runs pytest by hand on a GPU box. Until a maturin develop --features cuda + pytest tests/test_cuda.py step is added to the nightly, treat the Python GPU surface as verified by hand and not by CI — which is precisely how the compile_cuda export gap survived three releases while capabilities() advertised the feature.

A second gate that was inspecting nothing. cargo test --features groebner-cuda could not pass on a machine with no NVIDIA driver at all, contradicting the header comment of alkahest-core/tests/groebner_cuda.rs. cudarc panics rather than returning Err when libcuda.so cannot be dlopened, so gpu_available() — whose whole job is to answer “should the GPU tier run?” — aborted the three GPU tests instead of skipping them. It now treats a missing library and a missing device alike (both mean not available), while still failing hard when ALKAHEST_GPU_TESTS=1 asserted a device that is not usable. The ALKAHEST_GPU_TESTS=1 tier additionally asserts GpuBackendReport::ran_on_gpu(), so a “GPU test” that silently reduced every matrix on the CPU now fails rather than passing on identical results.

See also

• Code generation — the CPU JIT tiers, emit_c, StableHLO
• Error handling — the full code registry
• examples/gpu_batch_eval.py — CPU/GPU batch comparison, degrades cleanly to CPU-only on a wheel without the feature

Ball arithmetic

Unified evaluation (experimental)

alkahest.experimental.evaluate provides one result contract for exact rational, f64, and rigorous interval evaluation. It returns an EvaluationResult; mathematically unsupported inputs return status == "unsupported" and a stable E-EVAL-* reason code rather than raising. Invalid API input, such as an invalid mode or zero precision, still raises ValueError.

from fractions import Fraction
import alkahest as ak
from alkahest.experimental import evaluate

p = ak.ExprPool()
x = p.symbol("x")
result = evaluate(x + p.rational(1, 3), {x: Fraction(1, 6)})
assert result.value == Fraction(1, 2)
assert result.backend == "exact_rational"

Use mode="f64" for ordinary floating-point evaluation and mode="interval" with ArbBall bindings for an enclosure. In interval mode, result.value and result.enclosure are the same ArbBall; its lo and hi bound the true result. mode="auto" selects intervals for ArbBall bindings (or an explicit precision), exact rationals when possible, and otherwise falls back to f64.

Ball arithmetic provides rigorous enclosures: every operation produces an interval guaranteed to contain the true result. Alkahest uses FLINT’s Arb library as the backend.

ArbBall

An ArbBall represents the real interval [midpoint ± radius]:

from alkahest import ArbBall

a = ArbBall(2.0, 0.5)    # [1.5, 2.5]
b = ArbBall(3.0, 0.0)    # exactly 3.0

print(a.mid)   # 2.0
print(a.rad)   # 0.5
print(a.lo)    # 1.5
print(a.hi)    # 2.5

An ArbBall can also carry a precision (in bits) for the midpoint:

a = ArbBall(2.0, 1e-30, prec=128)  # 128-bit midpoint

Ball arithmetic operations

All arithmetic on ArbBall values produces a guaranteed enclosure. The radius grows to account for rounding and operation error:

a = ArbBall(2.0, 0.1)
b = ArbBall(3.0, 0.1)

print(a + b)    # [4.8, 5.2]  — radius grows by sum of input radii
print(a * b)    # guaranteed enclosure of [1.9, 2.1] * [2.9, 3.1]
print(a ** 2)   # [3.24, 4.41]  (squares the interval)

interval_eval

interval_eval evaluates a symbolic expression with ArbBall inputs:

from alkahest import ExprPool, ArbBall, interval_eval, sin, exp

pool = ExprPool()
x = pool.symbol("x")

# sin(1 ± 1e-10) — guaranteed enclosure
result = interval_eval(sin(x), {x: ArbBall(1.0, 1e-10)})
print(result.lo, result.hi)

# Multivariate
y = pool.symbol("y")
expr = sin(x) * exp(y)
result = interval_eval(expr, {
    x: ArbBall(1.0, 0.01),
    y: ArbBall(0.0, 0.01),
})

interval_eval guarantees that the output ball contains the true value for any input in the given input balls, accounting for all rounding in the intermediate computation.

AcbBall

Complex ball arithmetic for expressions over ℂ:

from alkahest import AcbBall

z = AcbBall(1.0, 0.0, 1.0, 0.0)  # 1 + i, exact

Use cases

Certified numerical evaluation. Compute a value and prove it lies within a tight bound without symbolic proof:

# Prove sin(1) ∈ [0.841, 0.842]
r = interval_eval(sin(x), {x: ArbBall(1.0, 0.0)})
assert r.lo > 0.841 and r.hi < 0.842

Numerical verification of symbolic results. After deriving a symbolic simplification, verify it numerically with rigorous bounds:

# Verify sin²(x) + cos²(x) = 1 at x = 1
lhs = sin(x)**pool.integer(2) + cos(x)**pool.integer(2)
r = interval_eval(lhs, {x: ArbBall(1.0, 0.0)})
assert 1.0 in r  # ball contains 1

Sensitivity analysis. Pass an input ball representing parameter uncertainty and observe how the output uncertainty grows:

# x = 1 ± 0.1 (10% uncertainty)
r = interval_eval(x**pool.integer(3), {x: ArbBall(1.0, 0.1)})
print(r)  # output uncertainty

Relationship to Lean certificates

Ball arithmetic and Lean certificate export are complementary:

• Ball arithmetic gives numerical certainty within floating-point computation.
• Lean certificates give symbolic/logical certainty for the rewrite steps applied.

Combining them: certify(interval(differentiate(f))) gives a derivative, evaluated with rigorous interval bounds, with a machine-checkable proof of the symbolic differentiation step.

Rigorous global bounds (Taylor models)

Ball arithmetic gives rigorous pointwise enclosures: evaluate f at a ball and the true value is inside the ball you get back. A research loop routinely needs rigorous global statements instead:

• the maximum of f on [a,b]×[c,d] is at most M;
• ∫_a^b f dx lies in [I₁, I₂];
• f has no root anywhere in this box.

These turn a numeric observation into a theorem, which is exactly the step that takes a candidate from “survived the sweep” to “established”.

import alkahest as ak

pool = ak.ExprPool()
x = pool.symbol("x")

r = ak.bound_on_box(x * (pool.integer(1) - x), [(x, 0.0, 1.0)])
r.lower, r.upper        # encloses the true range [0, 1/4]
r.budget_exhausted      # False — converged within the work budget

ak.verified_integral(ak.sin(x), x, 0.0, 3.14159)   # Enclosure containing ~2
ak.verified_no_roots(x*x + pool.integer(1), [(x, -5.0, 5.0)])   # "true"
ak.verified_sign(ak.exp(x), [(x, -5.0, 5.0)], "positive")       # "true"

The soundness contract

Every returned enclosure is a rigorous outer bound. Rounding is outward at every step, and any operation whose remainder cannot be bounded rigorously refuses rather than returning something plausible.

The consequence worth internalising: an enclosure may be wide, but it is never wrong. A wide-but-true bound is a fine answer for a loop — it just means “not settled yet”. A tight-but-false one is a poisoned lemma that everything downstream inherits.

Running out of work budget is therefore not an error. The enclosure is returned anyway, with budget_exhausted = True, still sound.

Why not just use interval arithmetic

Naive interval evaluation is rigorous but suffers the dependency problem: it forgets that the two xs in x - x are the same number.

A Taylor model carries a polynomial in the box’s normalised coordinates plus a rigorously enclosing remainder interval, so cancellation happens symbolically in the polynomial part and only the genuinely uncertain residue stays in the interval.

Branch-and-bound

bound_on_box runs Moore–Skelboe branch-and-bound: one pass for the minimum, one for the maximum. Each keeps a rigorous bound on the extremum and prunes sub-boxes whose enclosure proves they cannot contain it, so work concentrates where the extremum actually is.

This pruning is not an optimisation detail, it is what makes the module usable. A stopping rule that instead required every sub-box’s enclosure to be narrower than tol would be demanding a pointwise-tight model everywhere, which no finite budget achieves for a function with a wide range — exp on [-5,5] would exhaust the budget and still return an enclosure loose enough to straddle zero, so even exp > 0 would come back undecided.

Three-valued predicates

verified_no_roots and verified_sign return one of "true", "false", "undecided". The third is never collapsed into the other two.

A sign predicate is a universally quantified claim, so a single point where it provably fails disproves it. verified_sign uses that: it evaluates the expression rigorously at the box centre, the per-axis endpoints and (in low dimension) the corners, and reports "false" when one of those point enclosures lies strictly on the wrong side. Without it, x > 0 on [-1,1] — plainly false — could only ever be "undecided", since the range enclosure straddles zero by construction.

Proving a root exists

verified_no_roots returns "false" only with a proof in hand, and the proof is the intermediate value theorem. The full-box enclosure succeeding is already a continuity certificate — a Taylor model is only ever built where every elementary step stayed strictly inside its domain — and a box is convex. So if two points of the box can be found where f is rigorously proven to have opposite signs, the segment between them stays in the box and f vanishes somewhere on it.

The two points do not have to be the box’s own endpoints, and that is what makes the test usable: the search subdivides the box, records the sign of any sub-box whose enclosure has a determined one, and samples the centres of the rest.

The last row is the honest limit. A double root never changes sign, so no witness pair exists; "undecided" is the answer, and it is not upgraded to "false" on the strength of an enclosure that merely touches zero.

Inequalities that are tight at an endpoint

The interesting inequalities are usually the sharp ones, and sharp means the margin goes to zero somewhere. Subdivision alone cannot certify those: where the margin vanishes, every enclosure of the range straddles zero however fine the boxes get.

Two separate things are done about it. tol is an absolute width, so it is the wrong stopping rule for a sign question — an expression whose minimum is 10⁻¹³ meets a 1e-9 tolerance while its enclosure still straddles zero. verified_sign therefore re-runs the search with the sign itself as the goal, refining while the bound straddles zero rather than to a fixed width. And where the margin genuinely reaches zero, the box is split: a collar [a, a+δ] at the endpoint is handled by a truncated Taylor expansion with a proven Lagrange remainder, the rest by ordinary branch-and-bound. The pieces are closed and share the join point, so their union is the original box.

x = pool.symbol("x")
# Cusa–Huygens, denominator cleared: x(2 + cos x) − 3 sin x ≥ 0, tight at x = 0
f = x * (pool.integer(2) + ak.cos(x)) - pool.integer(3) * ak.sin(x)
ak.verified_sign(f, [(x, 0.0, 1.5)], "nonnegative")   # "true"

Mitrinović–Adamović, Wilker, Huygens and Jordan’s inequality behave the same way. The remainder is proven, not assumed: a Taylor coefficient counts as zero only when substitution and simplify land on a literal integer 0 — no numeric enclosure can prove a value is zero — and the tail is bounded by sup|g⁽ᵐ⁾|/m! enclosed over the whole collar, with analyticity certified by requiring every derivative up to g⁽ᵐ⁾ to enclose successfully there.

The limits are worth knowing:

(x − 7/10)²(x + 1) on [0, 3/2] is non-negative and touches zero in the middle; it stays "undecided" rather than being upgraded on the strength of an enclosure that merely touches zero.

Which functions are covered — ask before you build the workload

Taylor models reach the elementary fragment: exp, log, sqrt, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, abs, plus arithmetic and integer/rational powers — and, since 3.9.0, erf and erfc. Outside it are bessel_j0, bessel_j1, digamma, lambert_w, gamma, the elliptic integrals, floor and ceil, and so is any two-argument function such as atan2.

The three inverse hyperbolics carry the domain restriction their branch has: acosh needs the whole box strictly above 1 and atanh needs it strictly inside (-1, 1). A box that only touches the boundary is refused with E-VALIDATED-003, because the derivative is unbounded there and no Taylor remainder exists — that is a statement about the box, not about coverage, so bounds_supported still answers True.

That boundary is queryable, so a search loop can choose a certifiable route instead of discovering it by hitting E-VALIDATED-001:

ak.bounds_supported(ak.sin(x) * ak.exp(x))     # truthy
answer = ak.bounds_supported(ak.bessel_j0(x))
bool(answer), answer.functions                  # (False, ['bessel_j0'])
answer.blocker                                  # "function `bessel_j0`"

# Per primitive, in the agent contract:
{row["name"] for row in ak.capabilities()["primitives"] if row["taylor_model"]}

numeric_ball is not this flag. It reports pointwise ball arithmetic, which bessel_j0, digamma, lambert_w and floor all have; a Taylor model additionally needs a rule with a rigorous Lagrange remainder, which they do not. Both bits are honest — they answer different questions. taylor_model and bounds_supported are derived by running the Taylor evaluator, not from a maintained list, so neither can drift from what bound_on_box accepts.

A True answer means “will not be refused with E-VALIDATED-001”. It is not a promise of success: a covered function can still hit a domain violation or an infinite enclosure on a particular box, which is a property of the box and not of the expression.

Refusals

For E-VALIDATED-003 the search first tries bisecting away from the trouble, since a domain violation on a wide box is often a boundary effect rather than a genuine interior pole. Only after the box has been bisected far enough for that explanation to be exhausted does it refuse — which is the right answer for a real interior singularity, where the range is not a bounded interval at all.

Removable singularities in verified_integral

∫₀¹ ln(1+x)/x dx = π²/12 has nothing singular about it — only the expression is singular at x = 0, and the integrand extends continuously to 1 there. A Taylor model still refuses, because the reciprocal’s enclosure contains zero.

verified_integral recognises this shape. If the integrand splits as N(x)/D(x) and there is a point p of the offending sub-interval at which N and D both vanish, it enclosures that piece with Cauchy’s mean value theorem instead:

N(p) = D(p) = 0,  D' ≠ 0 on J   ⟹   ∀ x ∈ J\{p} :  N(x)/D(x) = N'(ξ)/D'(ξ)  for some ξ ∈ J
                                ⟹   ∫_J N/D dx ∈ |J| · range(N'/D' on J)

so the piece is bounded by an enclosure of N'/D', which is perfectly regular. The number returned is the integral of the continuous extension.

Three guards keep this from swallowing a genuine pole:

• N(p) = 0 and D(p) = 0 are established symbolically (substitute the exact rational p, simplify, require a literal zero). No numeric enclosure can prove a value is exactly zero, so none is asked to.
• D' must be certified non-vanishing on the sub-interval. That is what fails for sin(x)/x², where the denominator has a double zero and the integral does not converge.
• N and D must each have a successful enclosure over the whole sub-interval, which certifies they are analytic — and hence that the symbolic derivatives really are their derivatives.

ak.verified_integral(ak.log(pool.integer(1) + x) / x, x, 0.0, 1.0)  # ≈ π²/12
ak.verified_integral(ak.sin(x) / x, x, -1.0, 1.0)                   # ≈ 1.8921661
ak.verified_integral(pool.integer(1) / x, x, -1.0, 1.0)             # refuses: N(0) ≠ 0

What is still refused

An integrable but non-removable singularity is refused, and the message says so rather than implying the integral does not exist:

These need an integrable-tail bound or a singularity-removing substitution, neither of which can be derived rigorously from the expression alone today. The refusal is the honest answer; widening an enclosure to make them pass would break the contract that makes the module worth using.

Relation to the rest of the stack

This is the slow, certifying half of falsify fast, certify slow:

verified_integral complements the symbolic integrate_definite: use the symbolic path when you want a closed form, and this one when you want a guaranteed numeric interval — including for integrands with no elementary antiderivative.

Scope of this release

Shipped: Taylor model arithmetic over a box (arithmetic, powers, division, and the elementary functions with rigorous remainders), range enclosure by branch-and-bound, verified 1-D definite integrals including removable singularities, root absence, root existence and sign predicates.

Not shipped: multivariate verified quadrature (verified_integral is 1-D), improper integrals, integrable-but-not-removable singularities, and Taylor-model-based ODE enclosures.

ODE and DAE modeling

Alkahest provides symbolic infrastructure for ordinary differential equations (ODEs) and differential-algebraic equations (DAEs), including structural analysis and automatic index reduction.

Every snippet on this page is executed by tests/test_docs_ode_dae.py, so it runs against the version of Alkahest it ships with.

ODE

ODE represents an explicit first-order system d(state_vars)/dt = rhs. The constructor takes three positional arguments — state variables, right-hand sides, and the independent variable — and ODE(...) and ODE.new(...) are the same thing:

from alkahest import ExprPool, ODE

pool = ExprPool()
t = pool.symbol("t")
x = pool.symbol("x")
v = pool.symbol("v")

# Simple harmonic oscillator x'' + x = 0, written as a first-order system:
#   x' = v,  v' = -x
ode = ODE.new([x, v], [v, pool.integer(-1) * x], t)

ode.order            # 2 — number of state variables
ode.state_vars()     # [x, v]
ode.rhs()            # [v, (x * -1)]
ode.is_autonomous()  # True — t does not appear on the right

Initial conditions are attached one at a time, and with_ic returns a new ODE rather than mutating:

ode_with_ic = ode.with_ic(x, pool.integer(1)).with_ic(v, pool.integer(0))

Lowering to first order

lower_to_first_order takes the pieces of a scalar higher-order equation — the unknown, the right-hand side, the order, and the independent variable — not an existing ODE:

from alkahest import lower_to_first_order

# x'' = -4x  →  [x' = x_1, x_1' = -4x]
ode = lower_to_first_order(x, pool.integer(-4) * x, 2, t)

ode.state_vars()   # [x, x_1]
ode.rhs()          # [x_1, (x * -4)]

The auxiliary states are named x_1, x_2, … for the successive derivatives.

DAE

DAE represents a system of implicit equations mixing differential and algebraic constraints, F(t, variables, derivatives) = 0. Build one with DAE.new(equations, variables, derivatives, time_var) — four positional arguments:

• equations — each expression means g = 0. Write x' = f as dx - f.
• variables — the dependent variables.
• derivatives — a separate symbol standing for the time derivative of variables[i], e.g. pool.symbol("dx/dt"). Alkahest does not parse the name; the positional pairing is what makes it a derivative.
• time_var — the independent variable.

A variable with no entry in derivatives is purely algebraic. That is how you declare a Lagrange multiplier, and it is what makes the system high-index:

from alkahest import DAE, ExprPool

pool = ExprPool()
t = pool.symbol("t")
x, y, u, w = (pool.symbol(n) for n in ("x", "y", "u", "w"))
lam = pool.symbol("lam")               # Lagrange multiplier — algebraic
dx, dy, du, dw = (pool.symbol(n) for n in ("dx/dt", "dy/dt", "du/dt", "dw/dt"))
one, two = pool.integer(1), pool.integer(2)

# Cartesian pendulum (index 3): x' = u, y' = w, u' = -lam*x, w' = -lam*y - 1,
# subject to x**2 + y**2 = 1.
dae = DAE.new(
    [dx - u, dy - w, du + lam * x, dw + lam * y + one, x**two + y**two - one],
    [x, y, u, w, lam],   # five variables …
    [dx, dy, du, dw],    # … but only four have derivatives
    t,
)

dae.n_equations   # 5
dae.n_variables   # 5
dae.equations()   # the five expressions, each meaning "= 0"
dae.derivatives() # [dx/dt, dy/dt, du/dt, dw/dt]
dae.time_var      # t

Pantelides algorithm

The Pantelides algorithm performs structural index reduction: it repeatedly differentiates the equations that a maximum bipartite matching leaves unmatched, until every equation is matched to a variable. It returns the reduced DAE — not a separate report object:

from alkahest import pantelides

reduced = pantelides(dae)

reduced.index        # 1 — differentiation rounds used
reduced.n_equations  # 6 — one more than the input
reduced.equations()[-1]
# ((x * dx/dt * 2) + (y * dy/dt * 2)) — the differentiated constraint

reduced.index is the number of rounds, so 0 means the input already had a perfect structural matching and nothing was differentiated. The equations appended by differentiation are visible in reduced.equations(), and the higher jets they introduce (ddx/dt/dt, …) in reduced.derivatives().

pantelides raises ValueError (E-DAE-002) above index 10. dae_index_reduce(dae) runs Pantelides first and falls back to rosenfeld_groebner when it hits that cap.

Because the algorithm looks only at which variables occur in which equations, it can reduce a system whose coefficients make it unsolvable; structural regularity is not numerical regularity.

Differential elimination

rosenfeld_groebner prolongs the system — differentiating each equation and treating the new jets as fresh indeterminates — and computes a Gröbner basis after each round. The basis holds the algebraic consequences of the differential system, which is what elimination-based work (input–output equations, structural identifiability) needs:

from alkahest import DAE, ExprPool, rosenfeld_groebner

pool = ExprPool()
t, x, dx = pool.symbol("t"), pool.symbol("x"), pool.symbol("dx/dt")

dae = DAE.new([dx - x], [x], [dx], t)          # x' = x
result = rosenfeld_groebner(dae, max_prolong_rounds=1)

result.consistent    # True — the unit ideal was never reached
result.truncated     # True — stopped at the round budget, not at saturation
result.variables()   # [t, x, dx/dt, ddx/dt/dt] — jets, in exponent-slot order

[str(e) for e in result.final_basis().to_exprs()]
# ['(x + (-1 * ddx/dt/dt))', '(dx/dt + (-1 * ddx/dt/dt))']

final_basis() returns a GroebnerBasis that knows its variable ordering, so to_exprs() reads the relations back as Expr (each meaning = 0). See Solving for the rest of the GroebnerBasis surface.

truncated=True means prolongation stopped because max_prolong_rounds ran out, not because the chain saturated. A truncated basis is a sound set of consequences but need not be complete — “not in the basis” then does not mean “not a consequence”. Nonlinear jets often do not saturate in finitely many algebraic steps, so this is the common case.

final_basis() returns None when consistent is False.

Sensitivity analysis

Sensitivity analysis computes how solutions depend on parameters. sensitivity_system(ode, params) augments the state with ∂x/∂p; adjoint_system(ode, objective_grad) takes the gradient of the objective with respect to the state, as a list parallel to ode.state_vars():

from alkahest import ExprPool, ODE, adjoint_system, sensitivity_system

pool = ExprPool()
t, y, k = pool.symbol("t"), pool.symbol("y"), pool.symbol("k")

ode = ODE.new([y], [pool.integer(-1) * k * y], t)   # y' = -k*y

sens = sensitivity_system(ode, [k])
sens.original_dim                # 1
sens.n_params                    # 1
sens.extended_ode.state_vars()   # [y, dS_y_k]

# Objective J = y(T)**2  →  dJ/dy = 2y
adj = adjoint_system(ode, [pool.integer(2) * y])
adj.state_vars()                 # [lambda_y]
adj.rhs()                        # [(k * lambda_y)]

Acausal modeling

Acausal component modeling lets you describe physical systems by their component equations without manually choosing which direction information flows:

from alkahest import AcausalSystem, ExprPool, capacitor, resistor, voltage_source

pool = ExprPool()
t = pool.symbol("t")

# Component constructors return {"name", "n_equations", "n_ports", "component"}.
src = voltage_source("V1", pool.symbol("Vs"))["component"]
res = resistor("R1", pool.symbol("R"))["component"]
cap = capacitor("C1", pool.symbol("C"))["component"]

circuit = AcausalSystem(pool)
circuit.add_component(src)
circuit.add_component(res)
circuit.add_component(cap)

# Wire the loop: Vs.p -> R.p, R.n -> C.p, C.n -> Vs.n
circuit.connect(src.port("V1.p"), res.port("R1.p"))
circuit.connect(res.port("R1.n"), cap.port("C1.p"))
circuit.connect(cap.port("C1.n"), src.port("V1.n"))

# Flatten into a DAE
dae = circuit.flatten(t)
dae.n_equations   # 10

Built-in components (resistor, capacitor, voltage_source, and others registered via the component API) generate their constitutive equations automatically. AcausalSystem.flatten assembles them, plus the Kirchhoff-style connection equations, into a DAE that Pantelides can reduce. See examples/acausal_and_laplace.py for a runnable end-to-end example.

Laplace transform

The Laplace transform lives in alkahest.experimental (the calculus/ODE/transform surface is not yet semver-frozen):

from alkahest import ExprPool
from alkahest.experimental import inverse_laplace_transform, laplace_transform

pool = ExprPool()
s, t = pool.symbol("s"), pool.symbol("t")
F = laplace_transform(pool.integer(1), t, s)     # s^-1
f = inverse_laplace_transform(F, s, t)           # back to 1

Hybrid systems

HybridODE adds event handling to an ODE: at a crossing event, the state is reset and integration resumes. Event.new(name, condition, reset_map) takes the reset map as a list of (variable, new_value) pairs, and both HybridODE.new and add_event return new objects:

from alkahest import Event, ExprPool, HybridODE, ODE

pool = ExprPool()
t, x, v = pool.symbol("t"), pool.symbol("x"), pool.symbol("v")

# Bouncing ball: x' = v, v' = -1; velocity reverses at floor contact
base_ode = ODE.new([x, v], [v, pool.integer(-1)], t)
bounce = Event.new("bounce", x, [(v, pool.integer(-1) * v)])

hybrid = HybridODE.new(base_ode).add_event(bounce)
hybrid.n_events   # 1
hybrid.guards()   # [x] — fires when x crosses 0

Event.rising() and Event.falling() restrict an event to one crossing direction; by default it fires in both.

Polynomial system solving

Alkahest solves systems of polynomial equations symbolically using Gröbner bases.

solve

solve finds the solutions of a system of polynomial equations in a list of variables. It uses the groebner Cargo feature, which is included in all PyPI wheels (default feature since 2.3.1) and in all source builds.

from alkahest import ExprPool, solve, sqrt

pool = ExprPool()
x = pool.symbol("x")
y = pool.symbol("y")

# Linear system
solutions = solve([x + y - pool.integer(1), x - y], [x, y])
# → [{x: 1/2, y: 1/2}]

# Circle intersected with a line: irrational solutions
solutions = solve(
    [x**2 + y**2 - pool.integer(1), y - x],
    [x, y]
)
# → [{x: sqrt(2)/2, y: sqrt(2)/2}, {x: -sqrt(2)/2, y: -sqrt(2)/2}]

# Parametric solve: free symbols omitted from `vars` stay as parameters
solutions = solve([x**2 - y], [x])
# → [{x: sqrt(y)}, {x: -sqrt(y)}]

Solutions are symbolic: irrational roots are returned as Expr trees (e.g. sqrt(2)/2) rather than floats. Quadratic elimination produces exact symbolic answers. Free symbols that appear in the equations but are not listed in vars are treated as parameters, so solutions may depend on those symbols.

Solution types

The return value is a list of dicts mapping Expr variable → Expr solution:

for sol in solutions:
    for var, val in sol.items():
        print(f"{var} = {val}")
        # Evaluate numerically if needed
        from alkahest import eval_expr
        numeric = eval_expr(val, {})

solve returns an empty list for inconsistent systems and a GroebnerBasis handle for parametric families (infinite solution sets).

Pass numeric=True to return float values directly: solve(eqs, vars, numeric=True).

GroebnerBasis

A GroebnerBasis can be constructed directly for ideal-theoretic operations:

from alkahest import GroebnerBasis

# Compute a Gröbner basis (lex by default)
polys = [x**2 + y**2 - pool.integer(1), x - y]
gb = GroebnerBasis.compute(polys, [x, y])

# Check ideal membership
print(gb.contains(x - pool.rational(1, 2)))  # False

# Reduce a polynomial modulo the ideal — the remainder is a GbPoly
reduced = gb.reduce(x**3 + y**3)
print(reduced.to_expr())                     # y

Reading a basis

A GroebnerBasis is a sequence of GbPoly, and each GbPoly converts back to an Expr. This is how you read an elimination result — the generators of a Lex basis that are free of the eliminated variables are the eliminated relations:

len(gb)                       # 2 — number of generators
gb.order                      # "lex"
gb.variables()                # [x, y] — what exponent slots 0, 1 refer to

for g in gb:
    print(g.to_expr(), "= 0")
# (y^2 + -1/2) = 0
# (x + (y * -1)) = 0

gb.to_exprs()                 # the same list in one call
gb[0].terms()                 # [((0, 0), Fraction(-1, 2)), ((0, 2), 1)]

terms() gives (exponent tuple, coefficient) pairs with the coefficient as an exact Python int or fractions.Fraction; the exponent tuple is parallel to variables().

The conversion runs the other way with expr_to_gbpoly, which is what reduce and contains accept alongside plain Expr:

from alkahest import expr_to_gbpoly

p = expr_to_gbpoly(x**2 + y**2 - pool.integer(1), [x, y])
p.n_terms                     # 3
gb.contains(p)                # True
GroebnerBasis.compute_raw([p])

A GbPoly stores exponent vectors, not names, so converting one back needs the variable list its slots refer to. Every basis Alkahest hands out carries that list — including the ones from solve (solve variables followed by the free parameters), triangularize and rosenfeld_groebner — so to_expr() normally takes no arguments. Naming too few variables raises ValueError rather than quietly misreading the exponent slots.

Monomial orders

Supported orders: Lex (lexicographic), GrLex (graded lexicographic), GRevLex (graded reverse lexicographic). GRevLex is generally fastest for basis computation; Lex is required for elimination.

Parallel F4

With --features "groebner parallel", Gröbner basis computation uses Rayon for parallel S-polynomial reduction via the F4 algorithm.

GPU-accelerated Macaulay matrix (groebner-cuda) — Rust only, not wired into the solver

--features "groebner-cuda" compiles a CUDA kernel for the mod-p row reduction of the Macaulay matrix, with multi-prime CRT lifts reconstructing rational coefficients, and falls back to pure-Rust row reduction when no CUDA device is present.

It does not accelerate anything on this page. GroebnerBasis.compute, solve and triangularize run Buchberger/F4 on the CPU regardless; the GPU routine is reachable only as the Rust function alkahest_cas::poly::groebner::compute_groebner_basis_gpu, and production dispatch deliberately does not prefer it until the benchmark harness says it wins. There is correspondingly no capabilities()["features"]["groebner_cuda"] bit — it used to exist and report that the kernel had been compiled in, which no Python observation could confirm or refute, so 3.8 removed it. See GPU support.

Because the Rust entry point falls back to CPU row reduction when no device is present, it returns a GpuBackendReport alongside the basis: let (basis, backend) = compute_groebner_basis_gpu(gens, order, Some(0))?; and backend.ran_on_gpu() is the only way to tell a real GPU run from a fallback, since the basis is identical either way.

Elimination ideals

GroebnerBasis.eliminate computes the elimination ideal I ∩ k[remaining vars] by dropping every generator whose support mentions one of the given variables. Under a lex basis with the eliminated variables ordered first, what is left is a Gröbner basis for that ideal:

# Implicitize the parametric curve (t, t**2): eliminate the parameter t.
gb = GroebnerBasis.compute([x - t, y - t**2], [t, x, y])
gb.to_exprs()                      # [(t + (x * -1)), ((y * -1) + x^2)]

implicit = gb.eliminate([t])
implicit.to_exprs()                # [((y * -1) + x^2)]  —  y = x**2

Note the variable order passed to compute: t comes first, so lex eliminates it. eliminate requires the basis to know its variables (gb.variables()), and rejects a variable it is not written over.

Coefficient fields: Q(params) instead of Q[vars, params]

GroebnerBasis.compute(polys, vars, params=[...]) moves the listed symbols into the coefficient field Q(params) instead of the polynomial ring. They never enter the monomial order and never generate S-pairs, which is the difference between eliminating states from Q[states, Y, params] and from Q(params)[states, Y] — the parameter count no longer inflates the staircase.

from alkahest import GroebnerBasis

# a lives in the coefficient field Q(a), not the ring Q[x, y, a]
gb = GroebnerBasis.compute([a*x - y, x + y - one], [x, y], params=[a])
type(gb)                      # ParametricGroebnerBasis
[g.to_expr() for g in gb]     # coefficients are rational functions of a

Measured on a catenary compartmental ODE model (a linear chain of n states, output the first compartment, eliminating the states from the jet equations down to the input–output relation): at n = 4 states / 7 rate constants the parametric route runs in 0.27s against 4.2s putting the rate constants in the ring (Lex, --release, ~15×) and leaves 5 total basis generators against 25; at n = 5 states / 9 rate constants the parametric route finishes in 6.9s while the direct computation had not finished after 240s. These are wall-clock numbers on one machine, illustrating the shape of the difference (S-pairs among the parameters are exactly what the ring route pays for and the coefficient-field route never generates) rather than a promised ratio.

The result is generic. A leading coefficient in Q(params) can be a non-zero rational function of the parameters and still vanish at a specific parameter point, and there the basis this computation built is not the one the same algorithm would build over ℚ at that point:

gb.conditions()               # [a + 1] — the basis says nothing at a = -1
gb.is_regular_at([3])         # True
gb.is_regular_at([-1])        # False

gb.specialize([3])            # an ordinary GroebnerBasis over Q
gb.specialize([-1])           # raises ParamGroebnerError, code "E-PARAMGB-004"

conditions() lists the hypersurfaces the computation assumed non-vanishing — every leading-coefficient inversion contributes its numerator and denominator, every input coefficient contributes its denominator — factored into irreducible, primitive pieces so the report is a list of conditions rather than one opaque polynomial in many parameters. The list is sufficient, not necessary: it can flag a point that turns out fine (a removable coincidence the bookkeeping cannot see), but it never misses a point where the generic basis is genuinely wrong. specialize refuses on the flagged locus with ParamGroebnerError (E-PARAMGB-004) rather than returning something that is not a basis; check is_regular_at first if a degenerate point is a normal outcome for your caller.

The read path matches GroebnerBasis: the object is a sequence of ParametricGbPoly, each with to_expr() / terms(), and the basis itself has to_exprs(), eliminate(vars) (same Lex-with-eliminated-variables-first contract, refuses to eliminate a coefficient-field parameter since there is nothing to eliminate), reduce, and contains. GroebnerBasis.compute(..., params=None) or params=[] is the unmodified Q[vars] path; ParametricGroebnerBasis.compute(polys, vars, params, order=None) is the equivalent direct constructor in alkahest.experimental.

This surface is experimental (alkahest.experimental.ParametricGroebnerBasis / ParametricGbPoly) and requires --features groebner.

Performance

On the solve_circle_line benchmark (2-variable quadratic system), Alkahest is approximately 40× faster than SymPy due to the FLINT-backed polynomial arithmetic and the compiled F4 core.

Upcoming (v2.0): F5 / signature-based Gröbner basis, real root isolation, primary decomposition, and other advanced algorithms.

Positivity certificates (SOS / Positivstellensatz)

decide answers real-algebraic questions by CAD, and pays doubly-exponential cost for it. Most positivity questions that actually arise — is this bound valid, is this Lyapunov candidate non-negative, is this inequality true on a box — do not need a decision procedure at all. They need a certificate: a short algebraic identity that makes the answer checkable by anyone, including a proof assistant.

decide is not complete in this implementation: on some sentences it refuses with E-CAD-001 rather than answering. See decide refuses rather than guessing below — this changed in 3.8 and it changed because the alternative was answering wrongly.

import alkahest as ak

pool = ak.ExprPool()
x, y = pool.symbol("x"), pool.symbol("y")

cert = ak.sos_decompose(x*x - pool.integer(2)*x*y + y*y, [x, y])
cert.kind        # "sos"
cert.identity    # p = 1*(x - y)^2
cert.verify()    # True — re-expands to the target exactly
cert.to_lean()   # Lean 4 rendering, or None

Constrained, on a basic semialgebraic set:

g1, g2 = x, pool.integer(1) - x          # the box 0 ≤ x ≤ 1
cert = ak.prove_nonneg(x - x*x, [x], constraints=[g1, g2])
cert.kind        # "handelman"  —  x − x² = x·(1 − x)

Three outcomes, deliberately kept apart

This is the part that matters for a search loop, and the part most CAS get wrong by collapsing:

E-SOS-002 must be recorded as unknown, never as “not SOS”. It reports that this search — the LP-representable subcone described below, at this basis_degree (or level) — found nothing. Three different worlds produce it and the error cannot tell them apart: p is SOS but its Gram matrix lies outside the subcone; p is SOS only at a higher basis degree; p is non-negative but genuinely not SOS. A loop that maps E-SOS-002 to “the conjecture is false” or “p is not a sum of squares” closes a branch on evidence that does not support it, and a wrongly closed branch is invisible — nothing downstream will ever contradict it. Only E-SOS-003 is a refutation, and it carries a witness point.

E-SOS-002 is not a claim that the polynomial is not a sum of squares, and certainly not that it is negative. The canonical illustration is the Motzkin polynomial x⁴y² + x²y⁴ − 3x²y² + 1, which is non-negative everywhere but provably not a sum of squares itself — asked to decompose it directly (no multiplier), this module refuses with E-SOS-002, correctly: it does not report it as negative, and it does not invent a decomposition. sos_decompose’s full pipeline does not stop there, though (see “What the search actually covers” below) — it also tries multiplying by a power of x²+y² before giving up, and that succeeds for Motzkin, so the end-to-end call returns a certificate, not a refusal. The homogeneous 3-variable form of Motzkin still refuses even through the full pipeline (multiplier search included), and — unlike the affine case above — this is not just “not yet reached”: the classical fact needs multiplier power N = 2, not N = 1, for this specific homogeneous ternary form, and N = 2 has now been attempted to closure by a genuinely harder search (symmetry reduction plus an exact algebraic zero-vector restriction, cutting 165 free parameters down to 16) without succeeding — see “What’s still open” below for the numbers. That is the actual reachable illustration of a genuine E-SOS-002 from this module today: the refusal is a property of the search, not of the polynomial, and which polynomials it applies to shifts as the search grows more complete.

The three-way branch a loop should write:

try:
    cert = ak.sos_decompose(p, [x, y])
    verdict, evidence = "nonneg", cert    # proved, with a checkable identity
except ak.SosError as e:
    if e.code == "E-SOS-003":
        verdict, evidence = "negative", str(e)   # refuted; witness point in the message
    else:
        verdict, evidence = "unknown", str(e)    # E-SOS-002 lands here — leave the branch open

What the search actually covers

A sum-of-squares decomposition p = zᵀQz over the monomial basis z exists iff there is a positive semidefinite Gram matrix Q matching p’s coefficients. Deciding general PSD feasibility is a semidefinite programme, and no floating-point SDP solver is ever trusted with a certificate here — a rounded Q is not a proof.

The search tries three things, in order, before refusing:

1. The diagonally-dominant (DSOS) subcone — non-negative combinations of squares of a fixed generator set, (e_i)² and (a·e_i ± b·e_j)² for small coprime (a, b), solved with the exact rational simplex in real::sos::lp (Bland’s rule, so termination is unconditional and there are no epsilon tolerances). Every generator is literally a square, so a feasible point is sound by construction — but the cone is a strict subset of the SOS cone: diagonal dominance is not invariant under scaling the basis, so a perfect square as ordinary as (x/2 + 1/3)² has a Gram matrix — its only one — that is PSD but not DD.
2. The full PSD Gram cone, when DSOS fails (real::sos::psd::psd_search). This subsumes DSOS but is not free: it leans on a floating-point search (Jacobi eigendecomposition, PSD-cone projection, an annealed schedule of shrinking eigenvalue floors with several random restarts — real::sos::sdp) to propose a Gram matrix, which is then rounded to nearby rationals and re-expanded to check it equals p exactly before anything is returned. A Some here is always sound regardless of what the numeric search converged to; a None means only “the search did not turn up a certificate”, never “not SOS”.
3. A Reznick multiplier search, when even step 2 fails on p itself: tries (x_1²+…+x_n²)^N·p for N = 1..4 and reruns step 2 on the product. Some positive-definite forms are not SOS at all (Hilbert 1888 — this is what Motzkin’s polynomial witnesses), but Reznick’s theorem guarantees (Σxᵢ²)^N·p is SOS for some N; the search does not know N in advance and reports budget exhaustion honestly rather than a disproof.

E-SOS-002 at the end of all three is phrased as a statement about the search, not the polynomial — the search’s incompleteness, at any step.

Step 3’s search has to work harder than plain alternating projection, because the multiplier certificates it exists for are frequently tight — Motzkin’s polynomial and Robinson’s form (the textbook PSD-not-SOS examples) both have witnessing Gram matrices that are singular, sitting exactly on the boundary of the PSD cone rather than its interior. A first version of this search (annealed alternating projection with several random restarts) converged toward that boundary monotonically (confirmed by a diagnostic trajectory) but never reliably closed the last, asymptotically slow stretch — the textbook behaviour of alternating projection at a tangential (non-transversal) set intersection. The search now also tries Douglas–Rachford splitting with over-relaxation and a facial-reduction step — both standard escapes for exactly this stall — and with them, both (x²+y²)·Motzkin(x,y) and (x²+y²+z²)·Robinson(x,y,z) are found and exactly re-verified. What’s still open: the homogeneous 3-variable form of Motzkin at multiplier power N = 2, (x²+y²+z²)²·(x⁴y²+x²y⁴−3x²y²z²+z⁶) (N = 1 is not classically expected to work for this specific homogeneous ternary form at all), is not — and this has now been pushed well past a budget skip. A new fallback, symmetry reduction (real::sos::psd::symmetry_reduced_search), restricts the search to the subspace fixed by the target’s own signed-permutation symmetry (here, order 16: swap x, y, and independently flip the sign of each variable, since every exponent is even) whenever that subspace is genuinely smaller — cutting this case’s 165 free parameters to 26. An exact, non-numeric restriction on top of that — Motzkin’s known zero at (1,1,1) forces Q·z(1,1,1) = 0 on any witnessing Gram matrix Q, and z(1,1,1) is literally the all-ones vector, no rounding involved — cuts it again to 16. Even on that 16-parameter family, though, Douglas–Rachford converges only very slowly: 6,000,000 iterations bring the minimum eigenvalue to roughly −1.4·10⁻⁸, and rational rounding still fails at every stage even with denominators up to roughly 10⁹. This is now a genuine, quantified numerical-hardness finding, not an unexplored avenue or an under-tuned budget — checked to be a real search limitation and not a bug in the machinery the same way as before: an independent sanity check confirms the affine Gram-matrix family (and each reduced family) is constructed correctly, and a synthetic planted example with a singular Gram matrix of the same size is found and exactly re-verified. The remaining test for the 3-variable Motzkin form records undecided rather than a false certificate.

Constrained certificates

With constraints g_i ≥ 0, prove_nonneg searches for a Handelman certificate

p = Σ_α c_α · Π_i g_i^{α_i},    c_α ≥ 0 rational,   Σ_i α_i ≤ level

which is again an exact LP in the weights c_α. level is a user-visible parameter; exceeding it is a refusal (E-SOS-002), not a failure. Handelman is complete for polytopes given a high enough level, but the level needed is not known in advance — so raising it on refusal is a meaningful retry.

Verification

Every certificate is re-expanded in exact rational arithmetic and compared against the target identically before it is returned. A candidate that fails is refused (E-SOS-005), never returned with a caveat. verify() runs the same check on demand so a downstream consumer never has to trust the search that produced the certificate.

When to use this versus decide

The intended pattern is: try the certificate route first because it is cheap and its output is citable; fall back to decide on E-SOS-002 when you need a verdict rather than a certificate and can afford the cost. Note that neither route is complete, so “both refused” is a real and expected outcome — it means undecided by these methods, not false.

decide refuses rather than guessing

decide implements CAD over a bounded fragment: purely polynomial bodies over ℚ in one or two real variables, with a quantifier prefix of at most two. Outside that fragment it raises CadError (E-CAD-001). Inside it, there is one further refusal, and it is the important one.

The CAD sample set is built from rational points — bracket endpoints, refined brackets, midpoints. For a strict atom (<, >) that is complete: strict solution sets are open, so if a solution exists, a whole interval of rational points solves it too. For a non-strict atom (=, ≠, ≤, ≥) the solution set can be a single boundary point, and if that point is irrational it is never in the sample set. Concluding “no sample satisfied it, therefore unsatisfiable” would then be a claim about a point that was never tested — and via ∀x. φ ≡ ¬∃x. ¬φ, that fabricated false becomes a machine-checked-looking proof of a false universal theorem.

So when a boundary root has not been shown rational and the body has a non-strict atom, decide refuses:

import alkahest as ak

pool = ak.ExprPool()
x = pool.symbol("x")

# Rational double root: found exactly, so the verdict is real.
body = pool.gt((pool.integer(3) * x + pool.integer(2)) ** pool.integer(2), pool.integer(0))
ak.decide(ak.Forall(x, body))        # (False, None) — false at x = -2/3

# Irrational double root at ±sqrt(2): refuses instead of answering.
irr = pool.gt((x ** pool.integer(2) - pool.integer(2)) ** pool.integer(2), pool.integer(0))
try:
    ak.decide(ak.Forall(x, irr))
except ak.CadError as e:
    print(e.code)                     # E-CAD-001

Three consequences worth planning for:

• E-CAD-001 is “I did not establish this”, not “false”. A search loop must not record it as a closed branch. It is the same class of answer as E-SOS-002.
• Witnesses are verified. When decide reports (True, {...}) for an existential, the point is substituted back and checked; if it does not satisfy the sentence the witness is reported as None rather than as a certificate that fails. ∃x. 3x − 2 = 0 gives (True, {'x': '2/3'}); ∃x. x² = 2 gives (True, None), because no rational witness exists and a midpoint of the isolating interval is not one.
• Mixed-alternation sentences refuse more often than same-flavour ones. ∀x∃y. p > 0 is decided through ¬∃x∀y. p ≤ 0, and De Morgan turns a strict body into a non-strict one, so it can land in the refusal case even though the original body was strict.

If you need an answer where decide refuses, the routes are: a positivity certificate (above), alkahest.smt with a nonlinear-real solver, or rigorous numerics (validated bounds) if a quantified-over-a-box statement is good enough.

Scope of this release

Shipped: exact rational SOS over the DSOS generator cone, a general PSD Gram search (floating-point proposal, exact rational verification) for cases DSOS alone refuses — with Douglas–Rachford splitting and a facial-reduction step alongside the original annealed alternating projection, specifically so boundary-only (singular Gram matrix) certificates are reachable — a Reznick multiplier search ((Σxᵢ²)^N·p for N ≤ 4) on top of that (finds both Motzkin’s polynomial and Robinson’s form), Handelman certificates on basic semialgebraic sets, exact verification, and Lean export.

Not yet shipped: reliable certification of every boundary-case example — the homogeneous 3-variable form of Motzkin specifically is still out of reach (see above), so this is a real but narrower gap than “Motzkin doesn’t certify” was in the prior release — a proper interior-point solver that would close it more systematically, and Putinar-style certificates with genuine SOS — rather than non-negative constant — multipliers on the constraints. CertificateKind::Putinar exists in the certificate type so those can be added without a shape change.

Interoperability

Alkahest integrates with the Python numerical ecosystem at well-defined boundaries.

NumPy

Batch evaluation

numpy_eval vectorises a compiled function over NumPy arrays with zero unnecessary copies:

import numpy as np
from alkahest import ExprPool, compile_expr, numpy_eval, sin

pool = ExprPool()
x = pool.symbol("x")

f = compile_expr(sin(x) ** 2 + x, [x])
xs = np.linspace(0, 2 * np.pi, 1_000_000)
ys = numpy_eval(f, xs)   # returns a NumPy array, shape (1_000_000,)

Two argument conventions meet here, and mixing them up is the most common way this call fails:

f = compile_expr(x ** 2 + y, [x, y])   # a CompiledFn — numpy_eval takes this, not the Expr

numpy_eval(f, xs, ys)                  # ✓ one positional argument per input variable
numpy_eval(f, [xs, ys])                # ✗ ValueError — that is one argument, not two
f([1.0, 2.0])                          # ✓ CompiledFn evaluates ONE point, given as one sequence
f(1.0, 2.0)                            # ✗ TypeError

Both wrong spellings raise an error naming the right one. numpy_eval(expr, xs) — the expression instead of the compiled function — is a TypeError telling you to call compile_expr first.

Inputs are converted to f64 arrays via DLPack or __array__. The call is vectorised through CompiledFn.call_batch_buffer in Rust: each array is read via the buffer protocol in one bulk copy (no per-element Python float boxing, no .tolist()), the native call runs with the GIL released, and the result is written directly into a preallocated output array — no Python loop, and no intermediate Python list on either side. numpy_eval_par additionally distributes points across CPU cores via Rayon — but only when the extension was built with --features parallel. Every published wheel is; parallel is not a Cargo default, so a source build that omits it gets a silent alias for numpy_eval: same result, no speedup, no warning. Check capabilities()["features"]["parallel"] before attributing a timing to it; see Getting started.

Array protocol

CompiledFn objects implement __array__ for direct NumPy coercion:

result = np.asarray(f([1.0]))  # scalar result as a 0-d array

PyTorch

PyTorch CPU tensors are accepted wherever NumPy arrays are (via __dlpack__):

import torch
xs = torch.linspace(0, 1, 10_000)
ys = numpy_eval(f, xs)   # returns a NumPy array

For GPU tensors, use the compile_cuda path (requires --features cuda; see GPU support). Note that its Python call_batch takes and returns host sequences: a CUDA tensor is copied to the host and back. The zero-copy device-pointer entry point (call_device_ptrs) exists in the Rust crate only and has no PyO3 binding.

JAX

numpy_eval with JAX arrays

JAX arrays implement __dlpack__ and are accepted by numpy_eval:

import jax.numpy as jnp
xs = jnp.linspace(0, 1, 10_000)
ys = numpy_eval(f, xs)

JAX primitive source (to_jax)

to_jax registers a symbolic expression as a JAX primitive, making it callable inside JAX computations including jax.jit, jax.grad, and jax.vmap:

from alkahest import to_jax, ExprPool, sin

pool = ExprPool()
x = pool.symbol("x")

jax_fn = to_jax(sin(x) ** 2, [x])

import jax
import jax.numpy as jnp

# Use inside jax.jit / jax.grad
jit_fn = jax.jit(jax_fn)
grad_fn = jax.grad(lambda x: jax_fn(x).sum())

The primitive registers:

• A concrete def_impl that calls the Rust evaluator
• An abstract evaluation rule for shape/dtype propagation
• A JVP (forward-mode) rule derived from the symbolic gradient
• A vmap batching rule

StableHLO / XLA

to_stablehlo emits textual MLIR in the StableHLO dialect, which XLA and JAX’s XLA backend can compile:

from alkahest import to_stablehlo

mlir_text = to_stablehlo(expr, [x, y], fn_name="my_kernel")
# Pass to xla_client.compile() or save to .mlir file

SymPy interop

Alkahest’s kernel does not import SymPy. Two supported bridges exist on top of it: alkahest.crosscheck.to_sympy translates an Expr into a SymPy expression, and alkahest.crosscheck drives SymPy as a differential-testing oracle. The test oracle in tests/test_oracle.py uses SymPy as a ground-truth reference. For ad-hoc mixed workflows, converting through the string representation is fine.

The interop trap: casus-irreducibilis cube roots

Read this before round-tripping a symbolic result into another CAS. An expression can be correct in Alkahest and evaluate to a wrong number somewhere else, because the two systems do not agree on which branch a cube root denotes.

Matrix.eigenvals() on a 3×3 with an irreducible cubic characteristic polynomial and three real roots returns the Cardano form, and in the casus irreducibilis one of the cube-root radicands is negative:

import alkahest as ak

pool = ak.ExprPool()
I = pool.integer
M = ak.Matrix.from_rows([[I(2), I(0), I(-2)], [I(2), I(0), I(-1)], [I(1), I(1), I(2)]])

for value in M.eigenvals():
    print(value)
# two conjugate-looking siblings, then:
# (4/3 + (sqrt(298/27) + -89/27)^(1/3) + (-89/27 + (-1 * sqrt(298/27)))^(1/3))

That expression denotes the eigenvalue under the real cube-root convention. Alkahest is consistent about this and honest at the boundary: eval_expr on it refuses with E-EVAL-009, and interval_eval returns ArbBall(1.629231 ± inf) — an enclosure that is true and useless, rather than a number that is neither.

Hand the same expression to a principal-branch evaluator — SymPy, NumPy, most calculators — and (negative)^(1/3) takes the principal complex root instead. You get a confident number back, and it is not an eigenvalue. In one sweep of 720 random integer matrices, 14 produced eigenvalues of this shape.

So, when a loop exports symbolic results to another tool:

• Prefer transporting a verified numeric enclosure (refine_root, interval_eval, bound_on_box) rather than a radical expression, whenever the consumer only needs a number.
• If you must transport the expression, evaluate it in Alkahest first. A refusal (E-EVAL-009, or an infinite ball) is the signal that the expression is branch-sensitive and must not be handed over as-is.
• Never treat “the other tool produced a float” as confirmation. Substitute the value back into the characteristic polynomial (or whatever defined it) and check the residual.

This is the general shape of the hazard, not a quirk of eigenvals: an honest refusal inside Alkahest becomes somebody else’s silent error the moment the expression crosses the boundary. alkahest.crosscheck reports exactly this situation as incomparable rather than diverge, for the same reason.

DLPack

All DLPack-compatible arrays (NumPy, PyTorch, JAX, CuPy) are accepted at the numpy_eval boundary. The DLPack conversion is zero-copy for CPU arrays with matching dtypes; a device array is copied to the host first. There is no device-pointer boundary exposed to Python — call_device_ptrs is a Rust-crate API.

Exporting C code

emit_c generates a standalone C function for embedding in other projects:

from alkahest import emit_c

c_code = emit_c(
    sin(x) * exp(pool.integer(-1) * x),
    [x],
    var_name="x",
    fn_name="damped_sin",
)
print(c_code)
# double damped_sin(double x) { return sin(x) * exp(-x); }

The emitted code uses only standard <math.h> functions and has no Alkahest dependency.

Reinforcement learning environments

alkahest.rl turns the Alkahest CAS into verifiable RL environments: generators produce tasks, verifiers re-derive correctness (never storing reference answers in dataset rows), and optional curriculum schedulers advance Risch difficulty tiers.

The package has two layers:

Install

The RL stack is an optional extra on the main PyPI package. It requires Python ≥ 3.10 because verifiers does not support 3.9.

pip install "alkahest[rl]"

This pulls verifiers and datasets. The integration environment itself ships inside the alkahest wheel — no separate install is required for local use.

From source (development):

maturin develop --manifest-path alkahest-py/Cargo.toml --release --features egraph
pip install "verifiers>=0.1.5" datasets

Build with the egraph feature when possible: the integration verifier uses e-graph simplification as its second verification layer and falls back gracefully when it is unavailable.

Quick start — symbolic integration

from alkahest.rl.envs.integration import IntegrationVerifier, load_environment

# Standalone verifier (any trainer)
verifier = IntegrationVerifier()
reward = verifier.verify(
    "x^2",
    {"f_expr": pool_expr, "is_elementary": True, "pool": pool},
)

# Prime Intellect verifiers environment (after pip install "alkahest[rl]")
env = load_environment(
    difficulty_tier=0,
    n_train=1000,
    n_eval=100,
    hard_negative_fraction=0.25,
    adaptive=True,
)
# env.curriculum.record(reward)  # inside your training loop

Dataset rows store a parseable f_str integrand (not live Expr objects) so HuggingFace / Arrow serialization works. The async reward function reconstructs expressions at scoring time.

Risch tiers

Hard negatives (hard_negative_fraction) inject integrands certified NonElementary (e.g. exp(x²), sin(x)/x) so models learn to refuse honestly.

Core API

BaseGenerator

Produces (prompt, metadata) dicts. Never include a reference_answer key.

BaseVerifier

def verify(self, completion: str, metadata: dict) -> float:
    """Return reward in [-1, 1]."""

CurriculumScheduler

Tracks rolling pass rate at the current tier; advances when advance_threshold (default 0.70) is met over a sliding window (default 256 rewards).

Rubric

Framework-agnostic weighted reward functions (Rubric.score(**kwargs)). Prime Intellect’s vf.Rubric is used inside load_environment() for Hub compatibility.

veRL adapter

from recipes.verl_integration_reward import compute_score

# ground_truth = row dict from the integration dataset builder
score = compute_score(solution_str, ground_truth)

Pass compute_score to veRL’s reward_fn config. The ground_truth dict should include at least f_str, is_elementary, and fields the verifier expects after reconstruction (see IntegrationVerifier.verify).

Adding a new domain

1. Create python/alkahest/rl/envs/<domain>/ mirroring integration/.
2. Subclass BaseGenerator — emit unsolved tasks + metadata, no reference answer.
3. Subclass BaseVerifier — check the model output with Alkahest primitives.
4. Add load_environment() returning vf.SingleTurnEnv (optional verifiers import).
5. Add a Hub pyproject.toml + README.md if publishing independently.

Publishing to the Prime Intellect Environments Hub

The integration environment includes Hub metadata at python/alkahest/rl/envs/integration/. See Hub checklist below.

After install, users run:

prime env install alkahest/alkahest-symbolic-integration
prime eval run alkahest/alkahest-symbolic-integration -m <model>

Hub checklist

1. Alkahest on PyPI — the Hub package depends on alkahest>=3.7.0 (includes alkahest.rl). For local development before PyPI catches up, use maturin develop or a release wheel from GitHub.

2. Install the Prime CLI and log in:

   uv tool install prime
   prime login
   

3. Smoke-test locally from the Hub package directory:

   cd python/alkahest/rl/envs/integration
   pip install -e .                    # pulls verifiers + alkahest + datasets
   python -c "from alkahest.rl.envs.integration.env import load_environment; load_environment(n_train=4, n_eval=2)"
   

4. Push to the Hub (from the same directory):

   prime env push
   # or: prime env push --team alkahest-cas --auto-bump
   

   Hub CI validates the package, load_environment() entrypoint, and dependencies.

5. Run a Hub eval against a hosted model:

   prime eval run alkahest/alkahest-symbolic-integration -m <model>
   

6. Iterate — bump version in pyproject.toml (or use --auto-bump) for each publish; tag releases in git when verifier tiers or reward logic change.

Monorepo note

Implementation code lives in the main alkahest package; the Hub directory is a thin installable manifest that pins dependencies and declares the [tool.verifiers] entrypoint. This avoids duplicating environment code while still letting prime env install resolve everything from PyPI.

Further reading: Prime Intellect — Create & Upload Environment, Verifiers Environments Hub.

Derivation logs

Most transformations in Alkahest return a DerivedResult that records the exact sequence of rewrite steps applied. This log is the foundation for both human inspection and Lean proof export.

DerivedResult

DerivedResult is the return type of diff, simplify, integrate, sum_*, and most other transforming operations. Notable exceptions: limit returns a bare Expr, and series returns a Series.

from alkahest import diff, sin

pool = ExprPool()
x = pool.symbol("x")

dr = diff(sin(x**2), x)

Attributes

Methods

Rewrite steps

Each step in .steps is a dict with:

for step in dr.steps:
    print(f"  {step['rule']:25s}  {step['before']}  →  {step['after']}")

Side conditions

A side condition is a predicate that must hold for a rewrite to be sound:

• Positive(x) — x must be positive (e.g. for sqrt(x²) → x)
• NonZero(x) — x must be non-zero (e.g. for x/x → 1). For a symbolic x the rewrite fires and the condition is recorded; for a literal zero base it does not fire at all, since 0 · 0⁻¹ has no value (literal-zero carve-out)
• Integer(n) — n must be an integer (e.g. for some power rules)
• BranchCut(f, x) — records that f may have a branch cut at x

Side conditions propagate into the derivation log as SideCondition entries and are aggregated in dr.verification["side_conditions"]. A generated Lean source artifact is evidence that can be checked; it is not a claim that the project has checked the artifact with Lean.

For antiderivatives, exactly_verified means that the in-kernel symbolic residual d/dx(F) - f simplified to zero. numerically_checked means only that the integration soundness gate found agreement at several floating-point samples; it is useful evidence, but it is not an exact proof. lean_checked remains reserved for an actual completed external Lean check.

evidence = dr.verification
if evidence["status"] == "certificate_available":
    assert not evidence["externally_verified"]
    lean_source = dr.certificate
    # Invoke a pinned Lean/Mathlib checker before treating this as lean_checked.

Inspecting a derivation

dr = diff(sin(x**2), x)

print(f"Result: {dr.value}")
print(f"Steps ({len(dr.steps)}):")
for step in dr.steps[:5]:
    rule = step['rule']
    before = step['before']
    after = step['after']
    print(f"  [{rule}]: {before} → {after}")
    for condition in step["side_conditions"]:
        print(f"    side condition: {condition}")

DerivationLog overhead

Logging is always on and is cheap — a Vec<RewriteStep> appended to during traversal. The benchmark group log_overhead in alkahest-core/benches/alkahest_bench.rs measures logging cost separately from computation.

For production workloads where you only need .value, the steps list is still populated but you can ignore it. There is no way to disable logging in the current API (disabling it would compromise the Lean certificate pipeline).

Combining logs

When you chain operations, the logs are separate:

simplified = simplify(expr)
derived = diff(simplified.value, x)

# Full derivation: simplify steps first, then diff steps
all_steps = simplified.steps + derived.steps

For operations like integrate that internally call simplify, the log includes the simplification sub-steps interleaved with the integration steps.

Machine-parseable output: to_dict / to_json

Agents pay for every character a call returns. .steps, .verification, and .certificate_status are convenient to poke at interactively, but stitching them into one payload for logging, RPC, or a context window means writing that glue yourself, on every call site, forever. DerivedResult.to_dict() and DerivedResult.to_json() give you the stitched, versioned envelope directly:

dr = diff(sin(x**2), x)

full = dr.to_dict()                     # mode="full" is the default
compact = dr.to_dict(mode="compact")    # short keys, token-efficient
json_str = dr.to_json(mode="compact")   # json.dumps(dr.to_dict(mode="compact"))

Envelope shape (mode="full")

{
  "kind": "alkahest.derived_result",
  "schema_version": 1,
  "steps_schema_version": 1,
  "value": "<display string>",
  "verification": { "status": "...", "evidence": "...", "externally_verified": false, "artifact_format": "...", "side_conditions": [...], "method": "..." },
  "certificate_status": { "certifiable": true, "reason": "...", "blocking_steps": [] },
  "steps": [ {"rule": "...", "before": "...", "after": "...", "side_conditions": [...]}, ... ],
  "has_certificate": true
}

verification and certificate_status are exactly the dicts returned by the .verification and .certificate_status getters; steps is exactly .steps. kind is a stable discriminator string — useful when logs or RPC payloads mix DerivedResult envelopes with other structured outputs (e.g. error envelopes carrying E-SUBSYSTEM-NNN codes).

Schema versions

Two independent version constants, both starting at 1:

Also available as DerivedResult.SCHEMA_VERSION / DerivedResult.STEPS_SCHEMA_VERSION class attributes, and documented alongside the field-name contract in alkahest._result_schema (STEP_FIELDS, STEP_FIELDS_COMPACT). Either constant is bumped independently if its shape ever changes, so pinning schema_version/steps_schema_version in your own parsing code is safe across upgrades that don’t touch the piece you depend on.

Compact mode

mode="compact" keeps the same top-level envelope shape but shrinks the biggest token costs:

• Steps use short keys — r for rule, s for side_conditions — and omit before/after entirely. Those two expression strings are usually the largest part of a multi-step derivation and the single biggest win for token budget. s is itself omitted from a step’s dict when that step has no side conditions (the common case).
• verification is pruned to status and externally_verified only. These are the two fields that carry the honesty signal — whether the result is verified, and whether that verification happened out-of-process — so they are never renamed, abbreviated, or dropped in compact mode.
• certificate_status is pruned to certifiable and reason; the blocking_steps diagnostic list (which repeats before/after text) is dropped.
• No mode ever includes Lean certificate source text. has_certificate (bool) plus certificate_status["reason"] is enough to know whether a certificate exists and, if not, why — without paying for the source. Use the .certificate getter when you actually need the Lean source.

dr.to_dict(mode="compact")
# {
#   "kind": "alkahest.derived_result",
#   "schema_version": 1,
#   "steps_schema_version": 1,
#   "value": "...",
#   "verification": {"status": "certificate_available", "externally_verified": false},
#   "certificate_status": {"certifiable": true, "reason": "emitted"},
#   "steps": [{"r": "diff_sin"}, {"r": "sqrt_of_square_positive", "s": ["x > 0"]}],
#   "has_certificate": true
# }

Prefer to_dict(mode="compact") / to_json(mode="compact") over reading .steps directly in hot loops — batch derivations, autoresearch search plumbing, or anywhere you’re serialising many DerivedResults and only need the rule names, side conditions, and verification status rather than full before/after expression text.

An invalid mode (anything other than "full"/"compact") raises ValueError.

Beyond one call

DerivedResult is per-call. To accumulate many results into a citable, serialisable, re-verifiable artifact — a DAG of claims with stable IDs, hypotheses, and certificate status — see claim graphs.

See also

• Autoresearch / agent loops
• Batch — produce many DerivedResults without aborting on one failure
• Budgets — bound the call that produced the derivation
• Certificate coverage — certificate_status in the envelope

Claim graphs

DerivedResult (derivation logs) is a per-call object. A research loop that runs for days produces thousands of them and has nowhere to put them, so its output degrades into a transcript — and transcripts do not survive context compaction.

alkahest.research supplies the session-level artifact that is missing: a directed acyclic graph of claims, serialisable to disk, diffable, re-verifiable, and renderable into a document a human referee can read.

import alkahest as ak

pool = ak.ExprPool()
x = pool.symbol("x")

with ak.research.session(title="Worked example", pool=pool, capture=True) as s:
    integrand = x / (x**pool.integer(2) + pool.integer(1))
    definite = ak.integrate(integrand, x, pool.integer(0), pool.integer(1))
    ak.simplify(pool.integer(2) * definite.value - ak.log(pool.integer(2)))

graph = s.graph
graph.save("claims.json")
print(graph.to_markdown())

What a claim carries

Stable identity

The ID is sha256 over a canonical JSON encoding of the normalised statement, the sorted hypotheses, and the method — nothing else. It deliberately excludes timestamps, dependency edges, and the library version, so the same claim derived twice in different sessions receives the same ID and citing across runs is meaningful:

from alkahest.research import claim_id

claim_id("(1/2 * log(2))", ["x > 0"], "integrate")   # 'clm_…', deterministic

Expression statements are put into the kernel’s normal form (via simplify) before hashing, so 1 + x**2 and x**2 + 1 produce the same ID. Pass normalize=False to session() in very hot loops; IDs then depend on how the expression was built.

Because IDs are content-addressed, re-deriving the same claim is not an error — the stored claim keeps its status and derivation and merely gains the new dependency edges.

Hypotheses travel with the claim

This is the whole point of Assumptions being first-class. The predicates in the active assumption context — whether set by session(assumptions=…) or by an enclosing alkahest.context(assumptions=…) — are recorded on every claim and are part of its identity. A claim proved under x > 0 and the same claim proved unconditionally are different claims with different IDs.

Recording

Automatic capture

session(capture=True) wraps the DerivedResult-producing functions in the alkahest module namespace, so every result computed inside the block is recorded, and dependency edges are inferred by looking for previously-recorded values among the subexpressions of each call’s arguments.

Automatic capture has a boundary, and it is stated rather than hidden:

• it sees calls made through the module namespace — alkahest.integrate(...);
• it does not see a name bound before the hooks were installed (from alkahest import integrate) or methods on objects (Assumptions.simplify).

alkahest.research.captured_operations() lists exactly what is hooked, and session.capture_report() reports the mode, the hooked operations, the claim count, and any error raised inside the hook. Hooks are installed once and never removed — removing them when one session exits would silently disable capture for a session still running on another thread.

Explicit recording

For anything outside that boundary, record() is one line and always complete:

with ak.research.session(pool=pool) as s:
    result = ak.integrate(integrand, x)
    s.record(result, method="integrate", label="Antiderivative", sources=[integrand])

    # or, equivalently, letting the session make the call:
    s.run(ak.integrate, integrand, x)

Conjectures

conjecture() records a claim that is not proved. Its status is hard-wired to "unverified"; there is deliberately no parameter to say otherwise.

s.conjecture(
    "-2 * integral(x/(x^2+1), dx, 0, 1) + 1 * log(2) = 0",
    evidence="integer relation found by guess_relation at 60 digits",
)

The honesty invariant

The recording layer never upgrades a claim’s status.

• A claim’s status is whatever DerivedResult.verification["status"] said.
• conjecture() always produces "unverified".
• verify() may only lower confidence: a failed re-check sets "refuted"; a successful one appends an audit entry and promotes nothing.

The renderers follow the same rule. An emitted-but-unchecked Lean certificate is marked [CERT ONLY, UNCHECKED], never as a proof, and every document opens with the exact machine-checkable fraction:

Machine-checkable subset: 1 of 4 claims (25%). Only claims marked verified were checked by a checker. Everything else is recorded evidence and must not be read as proved.

Claim.machine_checked is true only for exactly_verified and lean_checked.

Querying

graph.by_status("unverified")
graph.machine_checkable()
graph.dependencies(cid)     # direct
graph.ancestors(cid)        # transitive, what it rests on
graph.dependents(cid)       # direct
graph.impact(cid)           # transitive — "what dies if this claim is false?"
graph.roots(), graph.leaves(), graph.summary(), graph.topological_order()

impact() is the query a loop needs when a lemma turns out to be wrong: it returns every claim that transitively cited it.

Edges always point at claims already in the graph, so cycles are impossible by construction; graphs loaded from JSON are topologically checked and raise CycleError.

Serialisation

to_json() / from_json() round-trip losslessly, with an explicit schema_version. Output uses sorted keys, so a byte diff is a content diff. Reading a document whose schema_version is newer than this build raises ClaimGraphError rather than guessing.

graph.save("claims.json")
same = ak.research.ClaimGraph.load("claims.json")

graph.to_json(stable=True)   # drops timestamps: byte-identical across runs
graph.digest()               # sha256 of the stable form

Volatile data (recorded_at, started_at, finished_at, audit timestamps) is confined to explicitly-marked fields and never contributes to identity or to the digest, so two runs of the same computation produce identical stable=True documents and identical digests.

Rendering

to_markdown() emits a document with a verification summary table, then one section per claim carrying its statement (LaTeX where available), hypotheses, method, dependency links, derivation table, certificate status, and last re-check.

to_latex() emits an article (or just the body, with standalone=False) in which every claim carries \label{clm:<id>} and dependencies are \hyperref links — so a loop that ran for a week emits a writeup with every claim linked to its derivation and its certificate status. Both build on latex() and the DerivedResult machinery rather than reimplementing them.

Re-verification

graph.verify() walks the graph and re-checks whatever carries a recipe, re-parsing each claim’s recorded text into a fresh pool — so a graph loaded from disk can be revalidated against a newer library build rather than trusted blindly.

Outcomes are ok (exact), numeric_ok (the exact check was inconclusive but a numeric residual passed — evidence, not a proof), failed, inconclusive, or skipped.

report = graph.verify()
report.ok            # False if anything was refuted
report.summary()     # {'ok': 2, 'numeric_ok': 1, 'skipped': 1}
print(report.to_markdown())

A long session accumulates, by design

Two things grow monotonically for as long as a session is open, and neither is a leak — both are the feature working.

• The claim graph holds one claim per captured operation, plus one string per capture failure. That is the point of a provenance record, but it means capture=True around a million-call sweep builds a million-entry graph in memory. Snapshot with to_dict() / to_markdown() and start a new session periodically rather than running one session for the life of the process.
• The pool never reclaims, and a session pins one for its whole lifetime. See ExprPool never reclaims — a session scoped to one problem, with its own pool, is the pattern that survives a multi-day run.

A complete loop

examples/pslq_research_loop.py runs the full experimental-mathematics loop — high-precision quadrature, guess_relation for the integer relation, a recorded conjecture, symbolic proof, Lean certificate, JSON round trip, and re-verification — and prints the rendered document.

See also

• Autoresearch / agent loops — how claim graphs fit next to budgets, batch APIs, and compact result envelopes
• Budgets — bound one candidate so a hard instance cannot stall the session
• Batch — fan out many DerivedResult producers without aborting on one failure
• Derivation logs — to_dict(mode="compact") for token-efficient claim payloads

Lean certificates

Alkahest can emit Lean 4 source for derivations. Generated source becomes a machine-checked proof only after it typechecks with the pinned Lean/Mathlib toolchain and without admitted placeholders.

Three levels of evidence

Derivation logs — always on, always cheap. Records every rewrite rule applied, with rule name and arguments. Human-readable; machine-parseable; forms the basis for Lean export.

Lean certificate export — for computations expressible as sequences of rewrites tagged with Lean theorem names. The library emits a .lean file containing a proof term. Emission makes a certificate available; it does not by itself mean Lean has checked the source.

In the agent contract, certificate_available has this same meaning. Only a corpus artifact compiled by pinned Lean/Mathlib without admissions can be described as lean_checked — and because that compilation happens out of process, in CI, no DerivedResult ever reports lean_checked as its status. It is therefore not listed in capabilities()["verification"]["statuses"], which enumerates only the statuses the library can actually emit.

Algorithmic certificates — planned evidence for operations where rewrite sequences do not work. Do not treat a derivation log as an independently verified witness.

Theorem mapping

Every primitive in the registry is tagged with a Lean 4 / Mathlib theorem name:

The full mapping lives in alkahest-core/src/lean/.

Exporting a certificate

from alkahest import diff, sin

pool = ExprPool()
x = pool.symbol("x", "real")

dr = diff(sin(x**2), x)

# The certificate is in dr.certificate when Lean export is enabled
if dr.certificate:
    with open("proof.lean", "w") as f:
        f.write(dr.certificate)

The emitted .lean file imports Mathlib and contains a proof term that Lean can verify:

import Mathlib.Analysis.SpecialFunctions.Trigonometric.Deriv

-- Alkahest certificate: d/dx sin(x²) = 2*x*cos(x²)
theorem alkahest_diff_sin_sq (x : ℝ) :
    HasDerivAt (fun x => Real.sin (x ^ 2)) (2 * x * Real.cos (x ^ 2)) x := by
  have h1 : HasDerivAt (fun x => x ^ 2) (2 * x) x := ...
  exact (Real.hasDerivAt_sin _).comp x h1

Strict Lean CI

The CI pipeline (.github/workflows/lean.yml) generates a deliberately small, strict corpus of basic arithmetic rewrites and d/dx x³. Every corpus entry must have the expected non-empty derivation log, must contain no sorry, admit, or axiom, and must typecheck with warnings treated as errors. The pinned Lean 4.9 compiler does not provide a --no-sorries command-line flag, so the source admission check is explicit in both the generator and CI.

1. Generates proof files via tests/lean_corpus.py
2. Compiles them with the pinned Lean/Mathlib toolchain (with Mathlib cached)
3. Fails the build if a proof contains an admission or does not typecheck

Knowing where the line is, before you compute

The boundary above is not something to discover by running an operation and finding .certificate is None. It is tabulated in Certificate coverage — a table generated by running a corpus and recording what actually emitted, checked into the repo, and pinned by a CI drift check. Three ways to use it:

import alkahest as ak

p = ak.ExprPool()
x = p.symbol("x")

# 1. Ask before you commit to a route. Truthy iff a certificate is produced;
#    carries the reason when it is not.
answer = ak.certifiable("integrate", ak.log(x), x)
bool(answer)        # False
answer.reason       # 'class_withheld'

# `mode="ledger"` answers from the table alone, running nothing at all —
# for scoring many candidate routes. The default `mode="verify"` never says
# True without the certificate in hand, and hands the computed result back on
# `answer.result` so the check costs nothing extra.

# 2. Read the whole boundary.
[r["shape"] for r in ak.certificate_coverage("diff") if r["verdict"] == "certified"]

# 3. Refuse to degrade silently.
with ak.context(require_certificate=True):
    ak.integrate(ak.log(x), x)   # raises CertificateUnavailableError (E-CERT-001)

capabilities()["verification"]["coverage"] summarises the same ledger, and each primitive’s lean_theorem bit is read from it — so the agent contract and the coverage table cannot disagree.

Coverage

The strict CI corpus currently covers:

• Basic arithmetic rewrites (add_zero, mul_one, mul_zero, constant folding, and pow_one)
• The polynomial differentiation fast path for d/dx x³
• Indefinite integrals of sin, cos, exp, and xⁿ, certified via the FTC derivative relation deriv (fun x => F) x = f
• Definite integrals ∫ x in a..b, f x = F b - F a of the same base family (sin, cos, exp, xⁿ), plus finite sums and numeric-literal constant multiples of those terms (∫ (sin x + cos x), ∫ 3·cos x, ∫ -exp x, ∫ (x² + sin x + 3·cos x), …), certified via Mathlib’s interval-FTC lemma (intervalIntegral.integral_eq_sub_of_hasDerivAt) composed with HasDerivAt.add/.const_mul/.mul_const and the matching IntervalIntegrable combinators. A symbolic (non-literal) coefficient, or any addend outside the base family, withholds the whole certificate.

Other exports are generated source, not CI-qualified Lean proofs. In particular, non-polynomial differentiation, conditional logarithm/power rewrites, integration, limits, and unsupported expression forms can require side conditions or currently use a placeholder tactic. They must remain unverified until their proof encoding and strict corpus coverage are added.

Planned algorithmic certificates include:

• Polynomial factoring
• Polynomial GCD

Side conditions in proofs

Side conditions (domain constraints and branch-cut restrictions) are recorded in the derivation log. They are not yet translated into Lean hypotheses by the exporter, so conditional rewrites are excluded from the strict corpus.

Certificate coverage

DerivedResult.certificate returns Lean 4 source only where the emitter can prove the statement without sorry. This page is the map of that boundary. It is generated by running a corpus and recording what actually emitted — never hand-maintained — from tests/lean_corpus.py (the strict, CI-typechecked corpus) and tests/textbook_gate/ (first-course calculus and algebra).

Corpus: 240 observations over 156 shape classes — 66 certified, 2 conditional, 0 partial, 88 withheld.

Reading a row

A shape class is an operation plus a structural fingerprint of its arguments: which primitive functions occur (funcs), whether they are applied to the bare variable, to an integer power of it, or to something else (fn_arg), what kinds of exponents appear (pow, pow_base), how factors combine (mul), and the top-level form of the expression (form).

certifiable() answers True only for certified classes, and in its default verify mode confirms by running the operation before saying so. It under-claims by construction: a shape the corpus has never reached answers False with reason unknown_shape.

diff

63 shape classes — 26 certified, 37 withheld.

integrate

41 shape classes — 14 certified, 27 withheld.

simplify

16 shape classes — 5 certified, 2 conditional, 9 withheld.

simplify_log_exp

6 shape classes — 6 certified.

simplify_trig

21 shape classes — 15 certified, 6 withheld.

sum_definite

8 shape classes — 8 withheld.

sum_indefinite

1 shape classes — 1 withheld.

Querying it

import alkahest as ak

p = ak.ExprPool()
x = p.symbol("x")

# Will this route give me something a referee can check?
answer = ak.certifiable("integrate", ak.log(x), x)
bool(answer)      # False
answer.reason     # 'class_withheld'

# Plan across candidate routes without computing any of them.
ak.certifiable("diff", ak.sin(x), x, mode="ledger")

# Fail loudly instead of degrading silently.
with ak.context(require_certificate=True):
    ak.integrate(ak.log(x), x)   # raises E-CERT-001

Error handling

Alkahest uses a structured exception hierarchy. Every error carries a stable diagnostic code, a human-readable message, an optional source span, and an optional remediation hint.

Exception hierarchy

AlkahestError (base)
├── ConversionError   (E-POLY-*)   — expression → polynomial/rational conversion
├── DomainError       (E-DOMAIN-*) — mathematical side conditions violated
├── DiffError         (E-DIFF-*)   — differentiation failed
├── IntegrationError  (E-INT-*)    — integration failed
├── MatrixError       (E-MAT-*)    — matrix shape / singularity / undecidable determinant
│   ├── LinearAlgebraError (E-LINALG-*) — elimination, decompositions, canonical forms
│   └── EigenError         (E-EIGEN-*)  — eigenvalues, eigenvectors, Jordan form
├── CadError          (E-CAD-*)    — real quantifier elimination, see [Positivity](./positivity.md#decide-refuses-rather-than-guessing)
├── OdeError          (E-ODE-*)    — ODE construction or lowering
├── DaeError          (E-DAE-*)    — DAE structural analysis
├── SolverError       (E-SOLVE-*)  — polynomial system solving
├── JitError          (E-JIT-*)    — LLVM/JIT codegen
├── CudaError         (E-CUDA-*)   — NVPTX compile, kernel launch, or driver, see [GPU support](./gpu.md)
├── PoolError         (E-POOL-*)   — ExprPool misuse
├── AnsatzError       (E-ANSATZ-*) — ansatz family construction or fitting, see [Ansatz families](./ansatz.md)
├── CrossCheckError   (E-XCHECK-*) — cross-CAS check could not be posed, see [Cross-CAS testing](./crosscheck.md)
├── SmtError          (E-SMT-*)    — SMT-LIB export, solver run, or model lift, see [SMT bridge](./smt.md)
└── BudgetExceededError (E-BUDGET-*) — budget/cancellation trip, see [Budgets](./budgets.md)

Error attributes

Every exception instance exposes:

import alkahest
from alkahest import ExprPool, UniPoly, ConversionError

pool = ExprPool()
x = pool.symbol("x")

try:
    # sin(x) cannot be represented as a polynomial
    p = UniPoly.from_symbolic(alkahest.sin(x), x)
except ConversionError as e:
    print(e.code)          # E-POLY-001
    print(e.message)       # "expression contains non-polynomial term: sin(x)"
    print(e.remediation)   # "Use Expr directly, or expand sin(x) as a series first"

Common errors and remediations

ConversionError (E-POLY-*)

Raised when an expression cannot be converted to a polynomial or rational function.

DomainError (E-DOMAIN-*)

Raised when a mathematical side condition is violated.

IntegrationError (E-INT-*)

SolverError (E-SOLVE-*)

PrimaryDecompositionError (E-IDEAL-*)

Refusals: when Alkahest declines to answer

A refusal is not a malfunction. These codes all mean “I could not establish this, and the alternative to saying so is a confident wrong answer” — the outcome an unattended loop must record as undecided, never as a negative result.

E-SERIES-003 travels out of band for the same reason (SeriesError is exhaustive) but is wired into the bindings: series returns SeriesError::InvalidOrder with calculus::series::take_series_refusal() pending, and the Python layer raises SeriesError with .code == "E-SERIES-003" — or BudgetExceededError when a budget was what stopped it.

E-IDEAL-005, E-IDEAL-006 and E-SOLVE-004 are new in 3.8 and travel out of band: PrimaryDecompositionError and SolverError are public exhaustive enums that cannot gain a variant in a patch release, so the refusal is returned inside an existing variant and the real code is available from ideal::take_ideal_refusal() / solver::regular_chains::take_triangularize_refusal(). The Python bindings consult both, so radical and primary_decomposition raise AlkahestError with .code == "E-IDEAL-005" / "E-IDEAL-006", and triangularize raises SolverError with .code == "E-SOLVE-004". AlkahestError subclasses ValueError, so code that catches ValueError is unaffected.

The takers are consuming, which is what keeps the carrier variant honest: a genuinely non-polynomial equation still reports E-SOLVE-001, because no refusal is pending for it. Both readings of the shared variant stay distinguishable.

The three-valued zero test behind E-LINALG-010 / E-MAT-004 is new in 3.8. Before it, “could not prove det ≠ 0” was silently read as “det = 0”, and Matrix.nullspace() returned a confident wrong basis for any 2×2 with a symbolic determinant.

import alkahest as ak

pool = ak.ExprPool()
a = pool.symbol("a")
zero, one = pool.integer(0), pool.integer(1)

# `mystery` has no evaluation rule, so its vanishing is genuinely undecidable.
opaque = pool.func("mystery", [a])
m = ak.Matrix([[opaque, zero], [zero, one]])

try:
    m.inverse()
except ak.MatrixError as e:
    print(e.code)          # E-MAT-004
    print(e.remediation)   # substitute concrete values for the parameters

try:
    ak.Matrix([[opaque, zero], [zero, zero]]).nullspace()
except ak.LinearAlgebraError as e:
    print(e.code)          # E-LINALG-010

LinearAlgebraError and EigenError are both subclasses of MatrixError, so except ak.MatrixError catches all three families; catch the subclass when you want to distinguish them. Note that eigenvects() raises EigenError — with code E-LINALG-010, because the code identifies what could not be decided, not which wrapper it surfaced through.

Catching errors by code

For programmatic error handling:

try:
    result = alkahest.integrate(expr, x)
except alkahest.AlkahestError as e:
    if e.code.startswith("E-INT-"):
        print(f"Integration failed: {e.remediation}")
    else:
        raise

Error taxonomy

Every error is classified on two independent axes: subsystem (determines the code prefix and exception class) and cause (informs the remediation hint).

Subsystem axis

E-RESIDUE-005 is raised only at the Python boundary — the Rust residue takes an already-parsed point and cannot reach that state — so it is deliberately absent from alkahest-core’s REGISTRY, on the same footing as E-SMT-001/003/004 in alkahest/smt.py and E-BATCH-001 in alkahest/_batch.py. It exists because residue(f, z, a) with a symbolic a reads perfectly well and used to escape as a bare AttributeError naming an attribute of the implementation, which is not an AlkahestError and so was invisible to except ak.AlkahestError.

Three of these describe outcomes that are results rather than malfunctions, and the wording of each is deliberate. E-ANSATZ-003 means no member of this family satisfies the constraints — for a search loop that is a closed branch worth recording, not a failure. E-XCHECK-002 means no oracle is installed, and exists so that a missing oracle can never be mistaken for agreement. E-SMT-003 refuses a model containing an algebraic number that cannot be lifted exactly, rather than truncating it to a float — a float witness recorded as an exact one is precisely the silent-error shape these subsystems exist to prevent.

E-CERT-* — certificate policy

This one is unusual: the computation succeeded. What is missing is the machine-checkable evidence, so it is a policy failure rather than a mathematical one. It is raised only when you ask for it, by alkahest.require_certificate(result) or ambiently inside with alkahest.context(require_certificate=True):. The remediation names the blocking rewrite rules where they can be identified.

import alkahest as ak

p = ak.ExprPool()
x = p.symbol("x")

with ak.context(require_certificate=True):
    ak.diff(ak.sin(x), x)        # fine — certifies
    ak.integrate(ak.log(x), x)   # raises E-CERT-001

Cause axis

1. User-input — the expression or argument is outside the supported fragment. Always has a remediation; carries a span once parsing lands.
2. Domain — input is syntactically fine but violates a mathematical side condition. Remediation is “substitute a different value,” not “reformulate.”
3. Unsupported — the operation is not implemented for this case. Must name the missing capability so users can file a feature request.
4. Resource/environment — CUDA device absent, out-of-memory, JIT target mismatch, pool closed. Typically no span; remediation references the environment, not the expression.
5. Internal invariant — a bug. Should never reach users in release; in debug it carries a backtrace. Use E-INTERNAL-001.

Adding a new error code

1. Does it fit an existing subsystem? Add a variant and a code one higher than the current max for that prefix.
2. Does it name a new subsystem? Add a prefix, a class, and an entry in REGISTRY in the same PR. Do not reuse prefixes across unrelated subsystems.
3. Write the remediation before the message — if you cannot say what the user should do, the taxonomy is telling you this is an internal bug, not a user error.

Users match on subsystem (the exception class); triagers filter on cause (the code suffix and remediation text).

Autoresearch / agent loops

Alkahest is useful inside unsupervised or lightly-supervised math search loops because it is designed to be called many times under a budget, with results that stay auditable. The pieces below are the search plumbing that sits next to the mathematics:

A minimal loop shape:

import alkahest as ak

pool = ak.ExprPool()
x = pool.symbol("x")
candidates = [x**2, ak.sin(x), ak.log(ak.log(x))]

with ak.research.session(title="Sweep", pool=pool, capture=True) as s:
    with ak.context(pool=pool, budget=ak.Budget(wall_ms=200, max_steps=50_000, seed=7)):
        for item in ak.integrate_many(candidates, x, parallel=True):
            if not item.ok:
                # E-BUDGET-* → deprioritize; E-INT-* → record and move on
                continue
            # Token-cheap record for the next iteration / a human referee
            _ = item.value.to_dict(mode="compact")

print(s.graph.to_markdown())

That snippet uses one pool for the whole sweep, which is right for a sweep that ends. It is wrong for a loop that runs for days: see running for days without dying below.

Honesty rules that matter in a loop:

• A budget trip is a fine answer, not a crash — catch BudgetExceededError (E-BUDGET-*) and deprioritize that candidate.
• A refusal is not a negative result. E-CAD-001, E-SOS-002, E-LINALG-010, E-MAT-004, E-SMT-003 and E-ANSATZ-003 all mean undecided by this route. Recording any of them as “proved false” or “no such object exists” is the most expensive mistake a search loop can make, because it closes a branch permanently. The only codes that are genuine mathematical verdicts are the ones documented as such — e.g. E-INT-004 (proven non-elementary). See Refusals.
• A batch never drops a slot — failures become BatchItem(ok=False, error=…).
• Compact mode never hides verification status — verification["status"] stays readable; Lean source is omitted on purpose.
• Certificates are withheld rather than lied about — see certificate coverage.

Running for days without dying

Four limits bound an unattended run. None of them is a bug you can wait out; all four are properties of the design, and a loop has to be written around them.

1. Memory is not budgeted, and ExprPool never reclaims. A pool created once at startup grows linearly and forever — roughly 200 bytes per interned node, on the order of 2–3.5 KB per integrate — at flat per-call latency, so the run dies by OOM with no slowdown to warn you. Use one pool per problem, drop it when the problem is done, and carry to_dict() envelopes rather than live Expr handles between iterations (holding any Expr, Matrix or DerivedResult pins its entire pool). Full treatment.

2. wall_ms is cooperative, and its granularity is one primitive operation. The call stops at the first checkpoint after the deadline. On a high-degree integrand that operation is a FLINT call, which nothing short of an OS-level kill interrupts — a 300 ms budget can return after ~2 s there. Details.

3. run_with_wall_fallback does not bound wall time for an uncooperative callee. It joins its worker before raising, so it returns when the callee returns. run_with_wall_fallback(time.sleep, 3.0, budget=Budget(wall_ms=50)) raises after 3000 ms. The only hard bound is an outer process timeout. Details.

4. Some questions get refused, not answered. decide is not complete; it raises E-CAD-001 rather than fabricate a verdict it cannot justify, and the linear-algebra zero test refuses with E-LINALG-010 / E-MAT-004 rather than pick a branch. Budget for refusals in the loop’s control flow, not just for failures.

The skeleton that respects all four:

import alkahest as ak

def run_one(problem):
    pool = ak.ExprPool()                  # (1) fresh pool per problem
    x = pool.symbol("x")
    with ak.context(pool=pool, budget=ak.Budget(wall_ms=500, seed=7)):   # (2)
        try:
            result = ak.integrate(build(pool, problem), x)
        except ak.BudgetExceededError:
            return {"status": "undecided", "why": "budget"}
        except ak.AlkahestError as e:
            kind = "verdict" if e.code == "E-INT-004" else "undecided"   # (4)
            return {"status": kind, "code": e.code}
    return {"status": "ok", "result": result.to_dict(mode="compact")}    # no live Expr escapes

Run the driver itself under an OS-level timeout (3), not run_with_wall_fallback.

See also the runnable experimental-mathematics demo examples/pslq_research_loop.py.

Budgets, cancellation, and determinism

A fan-out loop trying thousands of candidate rewrites/integrals/Gröbner bases cannot afford one pathological candidate to hang the whole batch — and an orchestrator that decides a candidate isn’t worth more time needs a way to stop it now, not wait for an OS-level kill (SIGKILL, a process timeout). alkahest.Budget and alkahest.context(budget=...) give heavy engines a cheap, structured way to bail out honestly — raising BudgetExceededError — instead of running unbounded.

import alkahest as ak

p = ak.ExprPool()
x = p.symbol("x", "real")

with ak.context(pool=p, budget=ak.Budget(wall_ms=50, max_steps=10_000, seed=7)):
    try:
        ak.integrate(hard_expr, x)
    except ak.BudgetExceededError as e:
        assert e.code.startswith("E-BUDGET-")
        # ... deprioritize this candidate and move on to the next one ...

Model

A Budget is an immutable (wall_ms, max_steps, seed) triple. Every field is optional; Budget() never trips a check on its own — only alkahest.request_cancel() can stop a call entered with a bare Budget().

context(budget=...) pushes the budget into a thread-local stack on the Rust side (alkahest_cas::budget) for the scope of the with block, and pops it on exit — including on an exception, matching every other resource the context manager owns. Budgets nest like every other context(...) key: only the innermost frame is consulted, so a nested context(budget=...) shadows the outer one rather than combining limits with it. A nested context(...) that omits budget= leaves the outer one active (nothing is pushed, so nothing shadows it):

with ak.context(pool=p, budget=ak.Budget(seed=1, max_steps=1000)):
    ak.budget_seed()   # 1
    with ak.context(pool=p):
        ak.budget_seed()  # 1 — no budget= here, outer frame still active
    with ak.context(pool=p, budget=ak.Budget(max_steps=2)):
        ak.budget_seed()  # None — this frame set no seed; it does not inherit
    ak.budget_seed()   # 1 — back to the outer frame

What checks the budget today

The Rust engines call a cheap cooperative checkpoint (alkahest_cas::budget::check) at a handful of strategic points — not blanket-inserted into every loop:

• alkahest.integrate — at the top-level entry (covers every route: algebraic, Risch/transcendental, rational-function, log-derivative); at the integrate_inner recursion boundary; at every integrate_raw entry, which the sum rule and the constant-multiple rule recurse through, so a long sum is bounded between summands; once per candidate of the derivative-divides u-substitution search (each surviving candidate runs a full recursive integrate, and there are up to twelve); at the stage boundaries of the rational-function route (normalisation, Hermite reduction, Rothstein–Trager, the partial-fraction pass and each of its irreducible factors); and inside the two Euclidean loops that dominate a hard rational integrand — the ℚ[x] GCD used to reduce A/D to lowest terms and the number-field GCD of Lazard–Rioboo–Trager. A trip raises BudgetExceededError — integration has a Result return type with a natural place to signal it.

  Those last few are not decoration. See how tightly wall_ms binds — before they existed a 300 ms budget on ∫ cos x·sin¹²x/(sin⁹x + sin x + 1) dx returned after 3.4 s, and the same integrand at degree 40 never returned at all.

• alkahest.limit — at every limit_inner recursion boundary, in the Gruntz comparability sweep, in the pole-clearing loop of the x ↦ 1/t substitution, and between Taylor coefficients of the local expansion (the loop that can grow without bound on nested radicals). LimitError is an exhaustive public enum and cannot grow a Budget variant without a major semver break, so a trip is reported internally as LimitError::DepthExceeded and the E-BUDGET-* cause is recovered out-of-band (alkahest_cas::calculus::limits::last_budget_trip); the Python binding raises BudgetExceededError exactly as integrate does. With no budget active the same paths are bounded by an internal work ceiling, so an unsolvable limit refuses with LimitError / E-LIMIT-004 instead of running unboundedly.

• alkahest.simplify (and simplify_with, simplify_batch) — once per full bottom-up rewrite pass. simplify has no error channel (DerivedExpr isn’t a Result), so a trip here stops further passes early and returns the best value simplified so far — exactly like running out of the existing max_iterations cap already does, silently. If you need a hard raise on a simplify call specifically, wrap it in alkahest.run_with_wall_fallback (below).

Other heavy primitives (Gröbner bases, homotopy continuation, …) do not yet check the budget; wiring them is a follow-up, not part of this cut. Calling check() is cheap when no budget is active and cancellation has not been requested (an atomic load, and — only if a budget is active — an Instant::now()), so it is safe to sprinkle at more call sites over time without a performance concern gating it.

How tightly wall_ms binds

wall_ms is cooperative: the call stops at the first checkpoint after the deadline, so it always overshoots by however long the engine had left in the stretch it was in. That makes the useful question “how long is the longest stretch”, not “is it exact” — and a budget whose overshoot grows without bound is not a budget at all.

Measured on ∫ cos x·sinⁿx/(sin^d x + sin x + 1) dx, the family that goes through the Weierstrass half-angle substitution and then Rothstein–Trager (elapsed until the trip, wall_ms=300):

and across budget sizes on the worst of them (n=40, d=17): 53 ms for wall_ms=50, 106 ms for 100, 305 ms for 300, 1071 ms for 1000, 3158 ms for 3000 — the overshoot is a small additive term, not a multiple of the budget and not a function of the problem size.

What is left, honestly. The residual granularity is one primitive polynomial operation, and past a certain degree that operation is a FLINT call — factorisation over ℤ, or a bivariate resultant. Those are single foreign-function calls: nothing short of an OS-level kill stops one part-way, and adding checkpoints around them cannot help. On a degree-62 integrand (1/(sin³¹x + sin x + 1)) one such call measured about 2 s, so a 300 ms budget there returns after roughly that long. That is the honest floor, and only an outer process timeout goes below it — not run_with_wall_fallback, which joins the same call rather than preempting it (below).

(The pure-Rust loops that used to dominate — the ℚ[x] and number-field Euclidean GCDs — are now checkpointed, which is what removed the growth. A per-step check was also tried on the ℚ long division underneath them and measured no further improvement, so it was dropped rather than kept for the look of it: it would only have made max_steps count faster for nothing.)

So the guarantee worth relying on is: the budget is checked between operations, and one operation on a high-degree integrand can take seconds. It is not a hard real-time deadline, and no cooperative mechanism can make it one.

Cancellation

alkahest.request_cancel() sets a single process-wide flag — deliberately not scoped to a thread or a Budget frame. It models “the orchestrator wants the current heavy operation to stop right now”, e.g. because a fan-out loop decided a candidate has used enough wall time, and the operation might be running on a different thread than the one that decided to give up on it. alkahest.is_cancelled() reads it; alkahest.clear_cancel() resets it — call this before starting the next candidate, or every subsequent call trips E-BUDGET-003 immediately.

import threading

def watchdog():
    time.sleep(0.05)
    ak.request_cancel()

threading.Thread(target=watchdog, daemon=True).start()
try:
    ak.integrate(hard_expr, x)
except ak.BudgetExceededError as e:
    assert e.code == "E-BUDGET-003"
finally:
    ak.clear_cancel()

The watchdog runs while the call runs

For that example to mean anything, the watchdog thread has to be able to execute during the call it is trying to cancel. alkahest.integrate and alkahest.limit — the two budget-honouring engines — therefore release the GIL around their core call (py.allow_threads, the same idiom simplify_par uses for its Rayon workers). Without that the flag was only ever observed if it had been set before the call: the watchdog could not run a single bytecode until the operation it wanted to stop had already finished, which is the opposite of what a fan-out search loop needs.

Two things that follow, and one that does not:

• Cancellation is cooperative, not preemptive. The flag is observed at the checkpoints listed above, so the call stops at the next one — not instantly. An engine stretch with no checkpoint runs to its end.
• Other calls still hold the GIL. Only integrate and limit release it (plus the parallel simplifiers and the compiled-function batch paths, for unrelated reasons). request_cancel() cannot reach a running Gröbner basis or homotopy continuation, because those do not check the budget at all yet.
• Nothing about pool safety changes. ExprPool is Send + Sync and interns through a lock-free index; releasing the GIL around a call that holds only a shared &ExprPool is strictly weaker than the concurrent Rayon access simplify_par already performs on the same structure.

Determinism seed

Budget(seed=...) doesn’t do anything by itself — it makes the seed available via alkahest.budget_seed() (Rust: alkahest_cas::budget::seed()) to any RNG-consuming sampler that chooses to consult it, instead of threading an explicit seed parameter through every call in a pipeline. Two runs entering Budget(seed=7) observe the same budget_seed() at every call site that reads it, so a search loop that seeds its own sampling from the ambient budget is reproducible run-to-run.

Budgets and threads

The budget frame is thread-local; the cancellation flag is process-wide. Every surprise in this area follows from that pair, so it is worth stating plainly:

• A worker thread does not inherit the budget its parent entered. Handing work to a concurrent.futures.ThreadPoolExecutor yourself runs it unbudgeted unless you re-enter the budget inside the worker.
• alkahest.batch_map / batch_map_iter / the *_many helpers do that for you under parallel=True: the active budget is snapshotted on the calling thread and re-entered in each worker task, so a trip is reported as E-BUDGET-00x on the item that tripped, exactly as it would be sequentially. wall_ms stays a single sweep-wide deadline; max_steps becomes per-item (the Rust step counter is not readable from Python). See Batch and streaming evaluation.
• alkahest.run_with_wall_fallback likewise enters its budget argument on the worker thread it spawns, so cooperative call sites there actually observe it.
• request_cancel() needs no propagation, and that cuts both ways: it stops every in-flight cooperative call in the process, not just the one you had in mind. A single candidate’s budget trip therefore never sets it — nothing in batch_map touches the flag.

The Python-layer wall-clock fallback

Because simplify cannot raise through its own return type, context(budget=...) alone only bounds it the same way max_iterations already does — silently, by returning early. alkahest.run_with_wall_fallback turns that silent truncation into a raised, coded error: it runs the call on a worker thread (with budget entered on that thread) and raises BudgetExceededError (E-BUDGET-001) when the call overruns wall_ms.

result = ak.run_with_wall_fallback(ak.simplify, big_expr, budget=ak.Budget(wall_ms=200))

It does not bound wall time for an uncooperative callee

Read this before putting it in a loop. run_with_wall_fallback joins its worker before the exception propagates, so it returns control when the callee returns — not at wall_ms. Measured: run_with_wall_fallback(time.sleep, 3.0, budget=Budget(wall_ms=50)) raises E-BUDGET-001 after 3000 ms. The error message reports the real elapsed time (“returned control after 3000 ms”) precisely so this shows up in a log instead of being inferred later.

For a callee that does reach a cooperative checkpoint the wait is short, because the worker now runs inside the budget and stops on it (integrate on a hard integrand: wall_ms=300 returns in about 320 ms) — but that is the case where context(budget=...) alone would already have bounded it. The uncooperative case, the one this function looks like it exists for, is the one it cannot bound.

Why not return at the deadline and let the worker run on? Python cannot kill a thread, so “return early” means leaking a live thread that still takes the GIL in bursts, still allocates into the pool, and can only be asked to stop through the process-wide cancel flag — which aborts every unrelated in-flight call, and which nobody can then clear safely (clearing it before the orphan observes it is a no-op; leaving it set makes every later cooperative call fail with E-BUDGET-003). In a multi-day loop that trades a bounded stall for unbounded orphan accumulation plus collateral cancellation. Joining is the lesser evil, so it is what the function does.

What actually bounds wall time, in order of preference:

1. context(budget=...) around an engine that checks the cooperative budget — integrate and limit today. This is the real mechanism; run_with_wall_fallback is a reporting shim over it.
2. An OS-level bound for anything else: run the work in a subprocess with a timeout, or put a process-level watchdog around the loop. Nothing inside one Python process can preempt a thread mid-FLINT-call — see what is left, honestly.

So reach for run_with_wall_fallback to get a raise out of a cooperatively-budgeted call that would otherwise hand back a silently truncated answer. Do not reach for it to contain an unknown callee.

ExprPool never reclaims

Budget bounds time and steps. Nothing bounds memory, and the shape of the memory growth is the single most likely way a multi-day loop dies. This section is as important as everything above it.

The mechanism

ExprPool is an append-only hash-consed arena. It has no clear, no truncate, no refcount and no garbage collector; the underlying storage cannot shrink. The only way to reclaim interned nodes is to drop the entire pool. And every Expr, Matrix, Series and DerivedResult holds a strong reference to the pool it came from, so keeping one interesting result alive keeps every node ever interned alive with it — which is exactly the usage pattern a research loop has.

Measured on this machine, 20 000 integrate calls with a distinct integrand each time:

one shared pool for the whole loop ......  1 992 bytes/call, forever, linear
a fresh pool per iteration ..............      0 bytes/call

Two properties make this nastier than an ordinary leak:

• Time stays flat. Per-call latency does not degrade as the pool grows, so there is no early warning — the loop runs at full speed until the OOM killer arrives. Growth is O(n) in memory with O(1) time.
• You cannot measure it from Python. ExprPool exposes no __len__ and no stats(), so a loop cannot watch its own footprint and decide to recycle.

Per-call cost depends on the operation. As a rough guide, roughly 200 bytes of resident memory per interned node, and on the order of 0.8 KB/call for diff or simplify, 2–3.5 KB for integrate, ~8 KB for a crosscheck.check, ~12.5 KB for a series of order 6. At one integrate per second on one pool that is gigabytes within a day.

The supported pattern: one pool per problem

import alkahest as ak

for problem in problems:
    pool = ak.ExprPool()             # fresh pool per iteration
    x = pool.symbol("x")
    with ak.context(pool=pool, budget=ak.Budget(wall_ms=500)):
        result = ak.integrate(build(pool, problem), x)
        record(str(result.value))    # keep a *string* / dict, not the Expr
    del pool, x, result              # dropping the pool reclaims everything

The critical line is record(str(result.value)). Holding the Expr (or the DerivedResult, or a Matrix derived from it) pins the pool and defeats the whole scheme. DerivedResult.to_dict() / .to_json() exist partly for this: they give you a plain-Python envelope that outlives the pool. Do not carry live Expr handles between iterations of an unattended loop; re-parse or rebuild them in the new pool if you need them again.

One operation grows even on identical input

Matrix.eigenvals() mints a fresh gensym (__eigen_lambda_N) into the pool on every call, so re-asking the same eigenvalue question keeps allocating for no new information — measured at about 1.9 KB/call on the same 2×2 integer matrix over 20 000 calls, where simplify on identical input is exactly 0. Cache eigenvalue results yourself rather than recomputing them in a loop. (Every other Python-facing entry point measured is flat on repeated input.)

If you enable the LLVM JIT

The jit (LLVM) feature leaks a whole LLVM Context per compile — a true leak with no pool to drop, on the error paths as well as the success path. Cranelift (the default wheel’s JIT) is unaffected. Do not compile in a loop under a +jit / +full build.

Error codes

All three are Cause::Resource in the Rust registry (alkahest_cas::errors::codes) — a budget/cancellation trip is an environment/policy limit, not a statement about the mathematics, so it is never conflated with e.g. IntegrationError::NonElementary (a proof that no elementary antiderivative exists). alkahest.integrate and alkahest.limit raise BudgetExceededError directly rather than wrapping it in IntegrationError / LimitError, so callers can catch it uniformly regardless of which engine tripped it:

try:
    ak.integrate(hard_expr, x)
except ak.BudgetExceededError as e:
    ...  # deprioritize and move on
except ak.IntegrationError as e:
    ...  # a genuine "no elementary antiderivative" or "not implemented" verdict

API reference

On the Rust side (alkahest_cas::budget): Budget, enter, BudgetGuard, check, seed, is_active, request_cancel, clear_cancel, is_cancelled, BudgetError.

See also

• Autoresearch / agent loops
• Batch and streaming evaluation — budgets compose with *_many / batch_map, including under parallel=True; a trip becomes one failed BatchItem carrying E-BUDGET-00x, not a killed process and not a mathematical verdict
• Error handling — E-BUDGET-* in the exception hierarchy
• Claim graphs — session-level provenance around budgeted calls

Batch and streaming evaluation

Search loops are embarrassingly parallel at the candidate level: try to integrate a hundred generated integrands, simplify a thousand rewrite targets, differentiate every entry in a lookup table. Every Alkahest entry point is one-call-one-answer, so today that fan-out is written by hand at every call site — and the first candidate that raises aborts the whole batch unless the caller remembers try/except around every single call.

alkahest.batch_map (and the *_many convenience wrappers over integrate, simplify, and diff) do that fan-out once. They never raise for a single bad element — the exception is caught and turned into a structured BatchItem carrying the failing exception’s stable E-* diagnostic code, so a loop can tell “this candidate has no elementary antiderivative” (a fine, expected answer) from “the whole batch process crashed”.

import alkahest as ak

pool = ak.ExprPool()
x = pool.symbol("x")

outs = ak.integrate_many([x**2, ak.log(ak.log(x)), ak.sin(x)], x)
for item in outs:
    if item.ok:
        print(item.index, "=>", item.value.value)
    else:
        print(item.index, "FAILED", item.error["code"], item.error["message"])

0 => (x^3 * 1/3)
1 FAILED E-INT-001 [E-INT-001] integrate: not implemented: ...
2 => (-1 * cos(x))

Honesty invariant

batch_map always returns exactly one BatchItem per input, in input order — a batch of 100 items yields a list of 100 items, full stop. Nothing in this module silently drops a failing candidate; a failure is recorded as ok=False with its error, never as a missing slot.

BatchItem

Exactly one of value / error is populated: ok=True implies error is None.

error["code"] is the raised exception’s .code when it is an AlkahestError-like exception — including the native error types, which expose the same attribute — otherwise alkahest._batch.UNEXPECTED_ERROR_CODE ("E-BATCH-001"), the fallback for a failure whose exception carries no diagnostic code of its own (e.g. a plain ValueError raised by caller code passed to batch_map).

batch_map and batch_map_iter

def batch_map(fn, items, *, parallel=False, max_workers=None, **kwargs) -> list[BatchItem]: ...
def batch_map_iter(fn, items, *, parallel=False, max_workers=None, **kwargs) -> Iterator[BatchItem]: ...

Both call fn(item, **kwargs) once per item. parallel=True fans the calls out over a concurrent.futures.ThreadPoolExecutor; some Alkahest hot paths (integrate, limit, the parallel simplifiers, NumPy evaluation) release the GIL for their native work, so a thread pool can genuinely overlap them. For calls that hold the GIL throughout, parallel=True mainly helps when fn itself does I/O or otherwise yields the GIL — it never makes anything incorrect, only sometimes not faster.

Order guarantees

• batch_map always returns results in input order, whether or not parallel=True. This is the guarantee to reach for when you need zip(items, batch_map(...)) to line up.
• batch_map_iter documents two different behaviours by design:
  • parallel=False streams in input order — item i is fully computed and yielded before item i + 1 starts.
  • parallel=True streams in completion order, not input order. This is the whole point of streaming under fan-out: a fast failure surfaces immediately instead of waiting behind a slow item that happened to be submitted first. Every yielded BatchItem still carries its original index, so a caller that needs input order can sort by it, or just use batch_map.

# Streaming: react to failures as they arrive, without waiting for the slowest item.
for item in ak.batch_map_iter(ak.simplify, candidates, parallel=True):
    if not item.ok:
        log.warning("candidate %d failed: %s", item.index, item.error["code"])

integrate_many / simplify_many / diff_many

Thin batch_map wrappers over the three most common derivation entry points:

def integrate_many(exprs, var, *bounds, parallel=False, max_workers=None, **kwargs) -> list[BatchItem]: ...
def simplify_many(exprs, *, parallel=False, max_workers=None, **kwargs) -> list[BatchItem]: ...
def diff_many(exprs, var, *, parallel=False, max_workers=None, **kwargs) -> list[BatchItem]: ...

integrate_many accepts optional trailing bounds (a, b) for a batch of definite integrals, exactly like alkahest.integrate. **kwargs on every helper is forwarded to the underlying call (e.g. assumptions= for simplify_many).

outs = ak.simplify_many(candidates, parallel=True)
ok = [o.value for o in outs if o.ok]
failed = [(o.index, o.error) for o in outs if not o.ok]

Never raises — except for real interpreter signals

batch_map and batch_map_iter catch Exception, not BaseException: a KeyboardInterrupt or SystemExit still propagates and stops the batch, since swallowing those would make the process unkillable. Everything else — including every Alkahest E-* error and any exception your own fn raises — is captured.

Combining with budgets

Wrap the batch in context(budget=…) so each candidate inherits the same cooperative wall/step limit (and optional seed). A trip surfaces as BatchItem(ok=False, error={"code": "E-BUDGET-00x", …}) rather than aborting the rest of the batch — see Budgets.

with ak.context(pool=pool, budget=ak.Budget(wall_ms=50, max_steps=10_000, seed=7)):
    outs = ak.integrate_many(candidates, x, parallel=True)

This works under parallel=True as well as parallel=False, but the two are not identical field-for-field, because a Rust budget frame lives on a thread-local stack and a worker thread does not inherit its parent’s. batch_map therefore snapshots the active budget on the calling thread and re-enters it inside every worker task:

The wall_ms deadline is captured when batch_map is called, not when context(budget=…) was entered — Python cannot read the frame’s start instant — so a batch launched partway through a budgeted block gets the full wall_ms again. That is one budget’s worth of slack for the whole fan-out, not per item.

A budget trip is not a mathematical verdict

This is the reason the propagation matters more than the speed-up. Before it, a fanned-out sweep ran completely unbudgeted, and the candidates a sequential sweep reported as E-BUDGET-001 came back as E-INT-001 instead — the integrator’s verdict that no elementary antiderivative exists. A research loop records that as a permanently closed branch, when in truth nothing was decided and the machine merely ran out of the time it was given. E-BUDGET-00x is Cause::Resource; keep the two apart when you interpret a BatchItem:

for item in outs:
    if item.ok:
        accept(item.value)
    elif item.error["code"].startswith("E-BUDGET-"):
        requeue(item.index)          # ran out of budget — undecided, try again with more
    else:
        close(item.index, item.error)  # a real verdict about the mathematics

Cancellation across a batch

request_cancel() needs no propagation — the flag is process-wide, so every worker already sees it and a caller can abort a whole in-flight sweep with it (each item then reports E-BUDGET-003). The converse is deliberate: one item tripping its budget never cancels its siblings. batch_map never sets the flag itself; the trip is recorded on the item that tripped, and the rest of the sweep runs out the shared deadline.

One pool for the batch, not for the process

A batch shares one ExprPool across all its items, which is right — the whole point is that the items are related. What is not right is reusing that pool for the next batch, and the next: ExprPool never reclaims, so a driver that keeps one module-scope pool and runs batch_map in a loop grows linearly and forever at flat latency. Construct the pool per batch and drop it, and carry item.value.to_dict(mode="compact") forward rather than the DerivedResult itself (holding one pins the whole pool). See ExprPool never reclaims.

See also

• Autoresearch / agent loops
• Budgets, cancellation, and determinism
• Derivation logs — compact envelopes
• Error handling

Creative telescoping (Zeilberger’s algorithm)

zeilberger decides combinatorial identities over a well-defined class and hands back a certificate you can re-check independently. That combination — a decision procedure whose output is a short, verifiable algebraic object — is why it matters for an autoresearch loop: it turns “this identity survived a numeric sweep” into “this identity is proved”, without a human in the loop.

import alkahest as ak

pool = ak.ExprPool()
n, k = pool.symbol("n"), pool.symbol("k")
one = pool.integer(1)

# F(n, k) = C(n, k), written as the ratio of gammas the engine recognises.
F = ak.gamma(n + one) / (ak.gamma(k + one) * ak.gamma(n - k + one))

cert = ak.zeilberger(F, n, k)
cert.order          # 1
cert.coeffs         # [a_0(n), a_1(n)] — here proportional to [-2, 1]
cert.certificate    # R(n, k)
cert.boundary       # "vanishes" — so the recurrence holds for the *sum*

The result says: with S(n) = Σ_{k=k_lo}^{k_hi} F(n,k),

Σ_i a_i(n)·S(n+i) = b(n)

where b(n) is a boundary term that the summand-level identity does not by itself make zero. Here it is zero — cert.boundary says so, and says it after computing it — so the recurrence reads S(n+1) − 2·S(n) = 0, which together with S(0) = 1 gives Σ_k C(n,k) = 2^n. When it is not zero the inhomogeneous recurrence is returned instead, and when neither can be proved nothing is claimed about the sum at all. That three-way verdict is the subject of its own section below, and it is the difference between a certificate and a theorem.

What the certificate asserts

The returned certificate is a rational function R(n, k) such that, with G(n,k) = R(n,k)·F(n,k),

Σ_i a_i(n)·F(n+i, k) = G(n, k+1) − G(n, k)

holds identically. That identity in k is the whole of what is verified. Because it is a rational-function identity, a reader — or a referee, or another CAS — can check it by clearing denominators and expanding, with no reference to how it was found.

The boundary hypothesis, and the verdict on it

Summing that identity over k = k_lo .. k_hi telescopes the right-hand side to G(n, k_hi+1) − G(n, k_lo) — a boundary difference, not zero. The familiar homogeneous recurrence for S(n) therefore holds only when that difference vanishes: the natural boundary hypothesis, which Zeilberger’s algorithm does not establish.

It holds in the usual case, where F vanishes outside 0 ≤ k ≤ n, and that covers every classical identity in this chapter. It fails, for instance, for F(n,k) = C(n,k)/(k+1): there G(n,0) = −1, and the true relation is (n+2)·S(n+1) − (2n+2)·S(n) = 1, not 0. Reading the homogeneous recurrence off the certificate there gives a false lemma.

zeilberger decides it rather than leaving it stated. cert.boundary is one of three values:

cert = ak.zeilberger(F, n, k)          # limits default to k = 0..n
cert.boundary            # "vanishes"
cert.boundary_rhs        # None — the right-hand side is 0
cert.limits              # (0, n), echoed back so the assumption is on the record
cert.boundary_reason     # why the verdict came out this way
cert.implies_sum_recurrence   # True for "vanishes" and "nonzero"

"nonzero" is a result, not a refusal. OEIS A279013, a(n) = Σ_{k=0}^{n} C(2k,k)/(k+1)·C(2n−1,n−k), gets a verified order-2 certificate in a tenth of a second; the homogeneous recurrence read off it fails against the sequence at the very first term. The engine now returns the inhomogeneous recurrence, which does hold:

cert.boundary       # "nonzero"
cert.boundary_rhs   # b(n), an explicit hypergeometric term in n
# Σ_i a_i(n)·S(n+i) = b(n) — checked exactly against 2, 8, 35, 161, 768, 3773

The summation range is part of the claim

The verdict is about the range in cert.limits, and it changes with it — truncating a sum by one term generally turns "vanishes" into "nonzero". Pass limits=(k_lo, k_hi) (each an Expr or an int) to say what you are summing:

cert = ak.zeilberger(F, n, k, limits=(0, n - pool.integer(1)))
cert.boundary                         # "nonzero" — the k = n term was dropped

# Or ask the same certificate about another range, without re-running the search:
cert.boundary_at(0, n)["boundary"]    # "vanishes"

Two things this design does deliberately. The default (0, n) is stated and echoed back, not inferred from the summand, so a caller summing over something else can see the mismatch. And a range the analysis cannot place — endpoints that are not integer-affine in n — is "unknown", never "vanishes".

What "vanishes" is worth

It is a proof, not a numeric check. Each endpoint of G is evaluated by exact order counting in Q(n): the multiplicity of the endpoint as a root of the certificate’s numerator and denominator, plus −e for every Γ(a·n+b·k+c)^e factor whose argument lands on a non-positive integer there (a pole of Γ, or a zero of 1/Γ). A strictly positive total order is an exact zero. A negative one means G is unbounded at the endpoint, and that is reported as "unknown". Nothing that merely looks like zero can produce "vanishes".

The verdict also accounts for a subtlety that is easy to miss: when the limits move with n, Σ_{k=0}^{n} F(n+i,k) is not S(n+i). For Σ_{k=0}^{n} C(n,k) the telescoped difference alone is −1, and it is the missing term C(n+1,n+1) = 1 that cancels it. The full statement is

b(n) = G(n, k_hi+1) − G(n, k_lo) + Σ_i a_i(n)·D_i(n)

with D_i the finitely many values of F between the range at n and the range at n+i, and it is b(n) that the verdict is about.

Symmetrically, "nonzero" needs a witness: an integer n₀ at which b(n₀) is nonzero in exact rational arithmetic. Sampling that finds only zeros proves nothing and yields "unknown".

cert.boundary_term still returns G(n,k) = R(n,k)·F(n,k) if you want to discharge the hypothesis yourself.

side_conditions is a list[str], the same shape as DerivedResult.verification["side_conditions"]: things the result depends on that were assumed rather than proved. It tracks the verdict — a discharged hypothesis, a refuted one and an open one read differently, so a loop that only reads this list still cannot mistake the three. It is never empty: even a proved boundary is a statement about the n at which everything involved is defined, and a permanent record of which range was assumed.

Verification is not optional

Every certificate is re-checked as an exact identity in Q(n)(k) before it is returned. A candidate that fails verification is discarded and the search continues; it is never returned with a caveat. This is the same withhold-rather-than-lie discipline as the Lean certificate exporter: in a loop, one unverified certificate becomes a false lemma that every downstream derivation inherits.

Verification runs in exact rational arithmetic — no floating point is involved at any stage of the algorithm.

The class it decides, and where it refuses

The supported class is proper hypergeometric terms:

F(n, k) = R(n, k) · z^k · w^n · ∏_j Γ(a_j·n + b_j·k + c_j)^(e_j)

with R a rational function, z, w nonzero rationals, a_j, b_j integers and c_j rational. Factorials, binomials and Pochhammer symbols are recognised and normalised into this form.

Everything else is refused, with a structured error rather than a guess:

E-HOLO-002 is worth dwelling on. It does not mean “no recurrence exists” — it means none was found within the bounds you set. The distinction matters: an agent can raise the bounds and retry, whereas E-HOLO-001 is a permanent answer about the input and the branch can be closed for good.

Cost and bounds

max_order (default 4) and max_degree (default 16) bound the search. Solving the linear system over Q(n) at one (J, d) pair gets rapidly more expensive as either grows — measured on Σ (−1)^k C(n,k)³ at order 1, a single probe goes from 0.7 ms at d = 0 to 0.6 s at d = 7 to 84 s at d = 12, and one extra order costs about what three extra degrees cost.

Both are upper bounds, not starting points. The (J, d) pairs are visited by iterative deepening, cheapest estimated candidate first, and the first relation that passes exact verification is returned. Raising a bound therefore widens the reach without moving where the search starts — Σ (−1)^k C(n,k)³ is decided at order 2 in 0.8 s at the defaults and 0.6 s at max_order=6, max_degree=64 (a hand-tuned max_order=2, max_degree=4 costs 0.2 s, because a tight bound also truncates the cheap probes the deepening would interleave first):

# Same answer, same order of magnitude — the bounds are not a starting point.
cert = ak.zeilberger(F, n, k)                            # defaults
cert = ak.zeilberger(F, n, k, max_order=6, max_degree=64)

What the bounds do control is the price of a refusal: a term with no recurrence inside them pays the full grid before raising E-HOLO-002. Set them to the largest search you are willing to wait through when the answer is “no”.

cert = ak.zeilberger(F, n, k, max_order=2, max_degree=6)

Is the order minimal?

Usually the interesting part of a certified recurrence is not that it exists but that it is short: an order-4 relation where the literature records a guessed order-5 is a result, and an order-4 relation that might have been order 3 is a coincidence. So the question has to be asked explicitly, and the answer is on the certificate:

cert.order_is_minimal   # True only when the search established it

Cheapest-first is not order-ascending. The deepening above orders probes by estimated cost, 3·(J−1) + d, so it can reach a cheap order-2 probe long before an expensive order-1 one — which is exactly what makes Dixon, Franel and Apéry decidable at the default bounds. A returned order 2 therefore does not establish that no order-1 relation exists, and order_is_minimal is False to say so. False means not established, never “a lower order exists”: a lower-order relation that had been found would have been the one returned.

It is True for free at order 1, and True whenever the cost-ordered plan happened to spend every lower order before the probe that succeeded — which it does at narrow max_degree, since the plan interleaves less there. To get the claim in general, ask for it:

cert = ak.zeilberger(F, n, k, max_degree=6, minimal=True)
cert.order_is_minimal   # True — every degree ≤ 6 at every lower order was refused

minimal=True walks the grid order-major: every degree 0..max_degree at order J is probed and refused before order J+1 is tried at all. Same bounds, same exact verification, same certificate — only what was ruled out along the way differs. The flag is computed from the probes that actually happened rather than from the mode, so it cannot drift away from what the search did.

The price is the whole hopeless low-order sweep the default plan exists to avoid, and it is charged against max_degree because that is the bound minimality is claimed relative to. Measured on this machine at max_order=4:

The default column is flat in max_degree — that is the deepening doing its job — and the minimal=True column is not, because it is the column that has to sweep.

The default is unchanged, deliberately — the cost-ordered plan is what makes these terms decidable at all, and minimal=True is an opt-in for when minimality is the result you intend to publish. Note the shape of the table: the sweep grows like 3^d, so the honest move is usually to claim minimality against the smallest max_degree you are willing to state rather than against the default 16.

Guessing the recurrence first

The other half of the loop is guess_holonomic, which fits a P-recursive recurrence to the first terms of a sequence in exact rational arithmetic — the guess in guess-then-prove, with zeilberger supplying the proof:

motzkin = [1, 1, 2, 4, 9, 21, 51, 127, 323, 835, 2188, 5798, 15511,
           41835, 113634, 310572, 853467, 2356779, 6536382, 18199284, 50852019]

guess = ak.guess_holonomic(motzkin)
guess.order, guess.degree   # (2, 1)
guess.surplus_terms         # 14 equations confirmed it without being needed
guess.confirmed             # True

A fitted recurrence is a conjecture, and the number that says how much of one is surplus_terms. See Guessing recurrences for the guard and what it refuses.

The q-analogue (alkahest.experimental.q_zeilberger)

q-hypergeometric sums — Gaussian binomials [n;k]_q, q-Pochhammer symbols (a;q)_n — are not proper hypergeometric terms in (n,k), so zeilberger refuses them (correctly) with E-HOLO-001. q_zeilberger is the q-shifted twin of the same algorithm, and the same discipline: the certificate is re-checked as an exact identity in Q(q)(qⁿ)(q^k) before it is returned.

import alkahest as ak
from alkahest.experimental import q_zeilberger, qbinomial

pool = ak.ExprPool()
q, n, k = pool.symbol("q"), pool.symbol("n"), pool.symbol("k")

# Σ_k [n;k]_q² · q^{k²} = [2n;n]_q  — the q-analogue of Σ_k C(n,k)² = C(2n,n).
b = qbinomial(pool, n, k)
cert = q_zeilberger(b * b * q ** (k * k), q, n, k)

cert.order            # 1
cert.boundary         # "vanishes"
cert.support          # ("0", "n") — where the summand is proved to live
cert.sum_term(3)      # the exact q-series value S(3), a polynomial in q

sum_term(n0) is the part worth reaching for. It evaluates the sum from the definition of the q-Pochhammer symbol, not through the shift quotients the search used, so checking Σ_i a_i(qⁿ)·S(n+i) = 0 against it is an independent check of the returned recurrence rather than a restatement of the certificate.

What it accepts

F(n,k) = R(qⁿ, q^k) · z^k · w^n · q^{A·k² + B·n·k + C·n² + D·k + E·n}
         · Π_j (q^{u_j}; q^{d_j})_{v_j}^{e_j}

with u_j, v_j integer-affine in n, k. Written as an expression: the heads qbinomial(N, K) and qpochhammer(u, d, v) (meaning (q^u; q^d)_v), powers of q whose exponent is a degree-≤2 polynomial in n and k, powers with a base free of n and k, and any rational function of q, qⁿ, q^k. Half-integer quadratic coefficients are fine — q^{k(k−1)/2} is not rational in q^k but all of its shift quotients are, which is the property the algorithm needs.

E-HOLO-024 is the interesting one. (q^k; q²)_n shifted in k moves its first argument by 1, which the base q² does not divide, so the shift quotient is an infinite product and no algorithm in this family applies. That is a permanent answer about the input, like E-HOLO-020, not a budget problem.

The boundary verdict is two-valued here

cert.boundary is "vanishes" or "unknown" — there is no "nonzero" arm. The sum it is about is S(n) = Σ_{k ∈ Z} F(n,k), which the analysis also proves is a finite sum, over the window in cert.support. Fixing the range at all of Z is what makes the proof short: the range does not move with n, so there are no D_i correction terms, and "vanishes" follows from two structural facts about the summand alone — that it vanishes outside an affine window in k, and that it is finite at every integer k.

The certificate is not evaluated at an endpoint, and that is deliberate rather than lucky: R genuinely has poles at integer k — on the summand above it has a double pole exactly where the summand has a double zero, and G(n, n+1) is a finite non-zero limit of 0·∞. What the proof does instead is find one k far to the right that is past both the window and the (finitely many) poles, where G = R·0 = 0 with no indeterminacy, and then induct downwards on G(n,k) = G(n,k+1) − Σ_i a_i(qⁿ)·F(n+i,k), whose right-hand side the support analysis has already shown is finite everywhere. That gives every G a finite value, poles included, without evaluating the product at one; G is then constant and zero beyond the window at both ends, and the sum over Z telescopes to zero. (Read analytically at generic q with 0 < |q| < 1; the conclusion is an identity between rational functions of q that holds on an open set, so it holds in Q(q).)

What is not implemented is the inhomogeneous arm: computing b(n) for a q-sum needs endpoint values of G that are not rational in qⁿ, so a summand whose support the analysis cannot bound gets "unknown" and no claim about its sum, not a guessed inhomogeneity.

One more caveat, and it is on every verdict’s side_conditions: q is treated as transcendental. Everything here is an identity in Q(q). Specialising q to a root of unity — which is what the q-supercongruence literature does — is a separate step with its own hypotheses.

Specialising at a root of unity (specialize_at_root_of_unity)

A proved Q(q) recurrence does not, by itself, license setting q = ζ_d for a primitive d-th root of unity: a coefficient or a sum value can have a pole there, and specialising anyway is exactly the q-analogue of the A279013 failure mode — a certificate that re-checks perfectly while the specialised claim is false. QZeilbergerCertificate.specialize_at_root_of_unity(d, n) takes that step as a decision, not an assumption:

from alkahest.experimental import cyclotomic_polynomial, q_zeilberger, qbinomial

pool = ak.ExprPool()
q, n, k = pool.symbol("q"), pool.symbol("n"), pool.symbol("k")
b = qbinomial(pool, n, k)
cert = q_zeilberger(b * b * q ** (k * k), q, n, k)

spec = cert.specialize_at_root_of_unity(3, 2)  # q = zeta_3, at n = 2
spec.status                       # "specializes" / "obstructed" / "unknown"
spec.sum_value(0)                 # S_zeta(2), the canonical rep in Q[q]/(Phi_3)
spec.sum_valuation(0)             # the exact Phi_3-adic valuation of S(2)
spec.modulus()                    # Phi_3(q) = q^2 + q + 1, exposed for a by-hand check

The hypotheses — no pole in any coefficient a_i(qⁿ) or sum value S(n+i) at ζ_d — are decided exactly, by polynomial divisibility by Φ_d(q) over Q in the cyclotomic field Q(ζ_d) = Q[q]/(Φ_d(q)); nothing is evaluated numerically at any stage. Φ_d is irreducible over Q, so “does p vanish at ζ_d” is exactly “does Φ_d divide p”, which is what makes the valuation — and therefore the decision — exact rather than approximate. cyclotomic_polynomial(pool, d) returns Φ_d(q) directly, so a caller can redo the whole check by hand.

status is three-valued, and the three are not interchangeable:

• "specializes" — proved: every coefficient and every sum value has non-negative Φ_d-adic valuation, so the specialisation map is defined on all of them, and the specialised identity was re-checked as an exact statement in Q(ζ_d) before being returned. Three further things are reported on this verdict rather than folded into it, because each of them makes a true verdict mean less than it looks:
  • is_vacuous — every coefficient died at ζ_d (the q → 1 limit at d = 1 is always like this), so the recurrence is 0 = 0. Still a theorem; it constrains nothing.
  • leading_coefficient_survives — False means the specialised recurrence no longer determines the last value from the earlier ones, even though it is not vacuous.
  • support_shrinks / effective_support — the q-Lucas phenomenon: [2;1]_q = 1 + q is non-zero in Q(q) and zero at ζ_2, so the surviving window at a root of unity can be a strict subset of the generic one. It can never grow.
• "obstructed" — a pole at ζ_d was exhibited: some coefficient or sum value has negative Φ_d-adic valuation (available via sum_valuation even on this verdict, since a negative valuation is the obstruction). Nothing is offered — sum_value and coefficient raise — and this is not a claim that the specialised identity is false, only that this route to it is blocked.
• "unknown" — the generic boundary verdict was already "unknown", so there is no proved Q(q) statement to specialise in the first place.