Python API

Numerics and evaluation API

Compiled evaluation, batch evaluation, and ball arithmetic.

Compiled evaluation

function

compile_expr(expr, inputs, expected_evals=None)

#source

Compile a symbolic expression to a fast native function.

Returns a callable Python object (PyCompiledFn).

expected_evals, when given, is the number of evaluations the caller plans (e.g. the length of the arrays it will pass to numpy_eval); it lets a small expression that would otherwise get the interpreter be compiled natively up front. Without it, the first large batch call upgrades the function on its own (see PyCompiledFn::batch_fn).

Given a list (or tuple) of expressions instead of one, returns a CompiledVecFn evaluating all of them in one pass per point, with the subexpressions they share computed once — the shape of a gradient or a Jacobian. numpy_eval on it returns an array with one row per output.

function

eval_expr(expr, bindings)

#source

Numerically evaluate expr; accepts DerivedResult as expr.

Raises DomainError (a ValueError) if expr is mathematically undefined at bindings — e.g. a removable singularity such as (x**2 - 1) / (x - 1) evaluated as written at x = 1 reduces to 0/0 and raises, rather than silently returning the limit (2). Call cancel() first to evaluate the simplified, everywhere-defined form instead.

function

numpy_eval(compiled_fn, *arrays)

#source

Vectorised evaluation of a CompiledFn over arrays.

Phase 25 / PA-8 — NumPy / JAX / PyTorch / DLPack array evaluation.

Evaluation runs entirely natively: input arrays are read via the buffer protocol (a single bulk copy per array, not a Python float per element), the compiled function is called with the GIL released, and the result is written directly into a preallocated output array. No intermediate Python list of floats (.tolist()) is ever created, which is what previously made this path dramatically slower than it needed to be for large arrays.

Calling convention

Two things are easy to get wrong, and both now raise an error that says so:

the first argument is the compiled function, not the expression — f = compile_expr(expr, [x, y]), then numpy_eval(f, xs, ys);
the arrays are separate positional arguments, one per input variable — numpy_eval(f, xs, ys), not numpy_eval(f, [xs, ys]).

(CompiledFn.__call__() goes the other way: it evaluates a single point, and takes that point as one sequence — f([1.0, 2.0]).)

Parameters

compiled_fnCompiledFn: A compiled function returned by compile_expr() or CompileCache. Not an Expr, and not a traced Python function.
*arraysarray-like: One array per input variable, passed as separate positional arguments. All arrays must have the same number of elements. Accepts NumPy arrays, PyTorch CPU tensors, JAX arrays, CuPy arrays (host copy made), or anything with __dlpack__ or __array__.

Returns

numpy.ndarray: Output values with the same shape as the first input array.

Raises

TypeError: compiled_fn is not a compiled function — most often an Expr that has not been through compile_expr().
ValueError: The number of arrays does not match the compiled function's arity, including the common case of one list holding all of them.

Examples

Python

>>> import numpy as np
>>> import alkahest
>>> p = alkahest.ExprPool()
>>> x, y = p.symbol("x"), p.symbol("y")
>>> f = alkahest.compile_expr(x ** 2 + y, [x, y])
>>> xs = np.linspace(0, 1, 1_000_000)
>>> ys = np.linspace(1, 2, 1_000_000)
>>> out = alkahest.numpy_eval(f, xs, ys)     # one argument per variable

function

numpy_eval_par(compiled_fn, *arrays)

#source

Parallel vectorised evaluation of a CompiledFn over arrays.

Identical to numpy_eval() but distributes the N evaluation points across all available CPU cores using Rayon. The Python GIL is released during evaluation so other threads are not blocked.

Requires the parallel build feature (maturin develop --features parallel). Falls back silently to the sequential numpy_eval() when the feature is not compiled in.

Parameters

compiled_fnCompiledFn: A compiled function returned by compile_expr() or CompileCache — not an Expr.
*arraysarray-like: One array per input variable, passed as separate positional arguments (numpy_eval_par(f, xs, ys), not numpy_eval_par(f, [xs, ys])). All arrays must have the same number of elements. Accepts NumPy arrays, PyTorch CPU tensors, JAX arrays, or anything with __dlpack__ or __array__.

Returns

numpy.ndarray: Output values with the same shape as the first input array.

Notes

For small N (< ~~1 000 points) thread-scheduling overhead may exceed the computation time; [numpy_eval()](~~/api/numerics/#numpy_eval) is faster in that regime.

Examples

Python

>>> import numpy as np
>>> import alkahest
>>> p = alkahest.ExprPool()
>>> x = p.symbol("x")
>>> f = alkahest.compile_expr(x ** 2, [x])
>>> xs = np.linspace(0, 1, 10_000_000)
>>> ys = alkahest.numpy_eval_par(f, xs)   # multi-core evaluation

class

class CompiledVecFn

#source

Several outputs compiled together: compile_expr([e1, e2, …], inputs).

Calling it on one point returns a list with one value per output; numpy_eval(f, xs, ys) returns an array of shape (n_outputs, *xs.shape). Subexpressions the outputs share are evaluated once per point, and every output equals what compile_expr of that expression alone returns.

Properties

n_inputs

#source

Number of input variables.

n_outputs

#source

Number of outputs.

tier

#source

Backend evaluating the outputs: "interpreter", "cranelift" or "llvm".

Methods

__call__(point)

#source

Evaluate every output at one point, given as a single sequence of floats (one per input variable). Returns a list of n_outputs floats.

class

class CompiledFn

#source

A compiled function returned by compile_expr().

Properties

n_inputs

#source

Number of input variables.

tier

#source

Backend that evaluates single-point calls: "interpreter", "cranelift" or "llvm".

Methods

__call__(*args)

#source

Evaluate at one point, given as a single sequence of floats.

f([1.0, 2.0]) — one argument holding every input value, not one argument per input. To evaluate over arrays of points instead, use alkahest.numpy_eval(f, xs, ys), which takes the arrays as separate positional arguments.

call_batch_raw(inputs_flat, n_vars, n_points)

#source

Batch-evaluate over N points (Phase 25 — NumPy/JAX array evaluation).

inputs_flat is a flat list of length n_vars * n_points laid out var-major: [x0[0], x0[1], …, x1[0], x1[1], …].

Returns a flat list of N outputs. In Python this is wrapped by the alkahest.numpy_eval helper that handles the buffer-protocol conversion.

call_batch_raw_par(inputs_flat, n_vars, n_points)

#source

Parallel batch evaluation using Rayon (requires --features parallel).

Identical to call_batch_raw() but distributes the N points across all available CPU cores. The GIL is released during evaluation so other Python threads are not blocked.

Use numpy_eval_par() for a NumPy-friendly wrapper.

call_batch_buffer(inputs, output)

#source

Fast-path batch evaluation directly against buffer-protocol objects (NumPy float64 arrays, array.array('d', …), memoryviews, …).

Unlike call_batch_raw(), this never materialises a Python list or float object per element: each input array's contents are copied into a native Vec<f64> in one bulk memcpy-style operation (via PyO3's buffer protocol), the compiled function runs with the GIL released, and the result is written back into the caller-supplied output buffer the same way.

inputs must contain exactly n_inputs buffers, each of length n_points, C-contiguous, and holding float64 elements. output must be a writable buffer of the same shape (e.g. numpy.empty).

This is the method used internally by numpy_eval(); most callers should use that helper rather than calling this directly.

call_batch_buffer_par(inputs, output)

#source

Parallel counterpart to call_batch_buffer() (requires --features parallel). Distributes the n_points evaluations across all available CPU cores via Rayon; the GIL is released for the duration of the native call.

Used internally by numpy_eval_par().

function

emit_c(expr, var, var_name="x", fn_name="eval_poly")

#source

Emit a C function that evaluates the Horner form of a polynomial.

Returns a str containing a complete C function definition.

Parameters

exprExpr: A univariate polynomial in var.
varExpr or list[Expr]: The polynomial variable. emit_c only supports univariate polynomials, so a list/tuple is accepted as a convenience but must contain exactly one Expr (e.g. [x] is equivalent to x).
var_namestr: The C variable name (default "x").
fn_namestr: The C function name (default "eval_poly").

Ball arithmetic

class

class ArbBall(mid, rad=0.0, prec=128)

#source

Create a real ball [mid ± rad].

Properties

mid

#source

Midpoint rounded to the nearest f64.

This is a lossy readout: past ~16 digits mid ± rad read through the f64 accessors is not an enclosure. Use mid_str, mid_exact or to_mpf for the full-precision midpoint, and lo / hi for outward-rounded f64 bounds.

prec

#source

Working precision of the ball, in bits.

rad

#source

Radius, rounded up to f64.

Nearest-rounding here would report a radius smaller than the true one and turn a valid enclosure into an invalid one at the Python boundary.

lo

#source

Lower endpoint, rounded down to f64.

The ball is computed at prec bits — far finer than f64 — so a nearest-rounded endpoint can land inside the true interval. A caller writing lo <= v <= hi would then get False for a value the ball genuinely encloses, which defeats the entire purpose of rigorous arithmetic. Rounding outward keeps the f64 view a valid enclosure of the exact one, at the cost of being at most one ulp wider.

hi

#source

Upper endpoint, rounded up to f64. See [lo].

is_exact

#source

True if the ball has zero radius (a single exact value).

Methods

mid_str() / rad_str()

#source

Midpoint as a decimal string carrying every digit of the stored binary value (round-trips exactly), e.g. "1.2020569...e0".

rad_str()

#source

Radius as a decimal string (full precision, no f64 underflow).

mid_exact() / rad_exact()

#source

Midpoint as an exact fractions.Fraction.

rad_exact()

#source

Radius as an exact fractions.Fraction.

to_mpf()

#source

(mid, rad) as mpmath.mpf values, exact conversions of the stored binary floats (requires mpmath).

contains(v)

#source

sin()

#source

cos()

#source

exp()

#source

log()

#source

sqrt()

#source

function

interval_eval(expr, bindings, prec=None)

#source

Evaluate a symbolic expression using rigorous interval (ball) arithmetic.

bindings is a dict mapping Expr → ArbBall.

GPU evaluation

function

compile_cuda(expr, inputs)

#source

alkahest.compile_cuda(expr, inputs) -> CudaCompiledFn

Compile a symbolic expression to a CUDA GPU kernel targeting the NVPTX backend. Returns a CudaCompiledFn object whose .ptx attribute contains the generated PTX assembly.

Evaluation

function

evaluate(expr, bindings, *, mode="auto", precision_bits=None)

#source

No docstring.

class

class EvaluationResult

#source

Structured result returned by the unified evaluator (stable top-level).

Properties

value

#source

status

#source

backend

#source

requested_mode

#source

requested_precision_bits

#source

achieved_precision_bits

#source

enclosure

#source

reason

#source

is_enclosure

#source

class

class ModularEvaluation

#source

The result of evaluating a ModularRecurrence, with the evidence that makes its residues trustworthy.

singular_indices is the field to read when a residue is surprising: it lists the steps where the recurrence's leading coefficient was not a unit mod p and the working precision had to absorb the loss.

Properties

prime

#source

The prime the evaluation ran at.

precision

#source

k — the precision that was asked for, and delivered.

working_precision

#source

K >= k — the precision the forward pass actually ran at.

working_precision - precision is the total p-adic precision lost to singular steps, and is 0 for a recurrence whose leading coefficient is a unit throughout.

modulus

#source

p**k.

n_singular

#source

How many singular steps there were, in total.

steps

#source

How many forward steps the evaluation took.

Methods

residues()

#source

The residues, one per requested index, each in [0, p**k).

singular_indices()

#source

Indices n where p divides the leading coefficient a_J(n).

Truncated to the first 64; n_singular is the full count.

function

guess_relation(
    constants,
    precision_bits=664,
    max_abs_coeff=None,
    check_precision=True,
    *,
    digits=None
)

#source

Search for integers aᵢ with Σ aᵢ·constantsᵢ ≈ 0.

Returns the coefficients as a list of Python int (of any size), or None when no relation is found.

Precision travels with the input. The search runs at the smaller of precision_bits and the precision the constants carry, and its detection threshold is derived from that effective width:

================================== ================================== Input Search precision it allows ================================== ================================== float 53 bits (the 64-bit engine minimum) mpmath.mpf its mantissa, passed exactly — never through float decimal str / Decimal its significant digits ArbBall what its radius leaves of 53 bits int, Fraction unbounded (exact) digits= a further cap ================================== ==================================

Searching wider than the data used to zero-pad it into exact rationals, which both hid true relations — 90·ζ(4) = π⁴ from 115-digit strings came back None under the 664-bit default — and manufactured false ones among the padded values. An mpf went through float(), so 120-digit inputs were searched as doubles.

Raises PslqError (E-PSLQ-004) when the relation found is larger than the inputs' precision can justify, or does not hold to the precision they carry, rather than returning a number that is an artifact of how the inputs were written; and E-PSLQ-005 when the inputs are exact rationals and the relation is simply false for them:

Python

guess_relation([float(pi), float(e), float(log(2))])
# was: [-60771139, 67263243, 11653676]   ← confident and meaningless
# now: None, or E-PSLQ-004 — never that relation

The test is on the result, not the input, so small relations among floats still come back: [1.0, 2.0, 3.0] needs almost no precision to pin down and is returned normally. See relation_confidence() for the numbers.

A relation among decimal strings or mpf values that passes those tests is returned, but unendorsed unless digits= declares how many digits are trustworthy — relation_confidence() reports credible: None for it until then. digits= also caps the search width:

Python

guess_relation([mpf(x) for x in floats], digits=16)   # raises E-PSLQ-004
guess_relation(two_hundred_digit_strings, digits=200) # judged, and passes

Pass check_precision=False to accept a relation among the supplied values themselves (the search width is still capped at their precision).

function

relation_confidence(
    constants,
    coeffs,
    *,
    digits=None,
    precision_bits=None,
    margin_digits=None
) → dict

#source

Judge a found relation against the precision its inputs actually carry.

Pinning down n coefficients bounded by H selects one of the (2H+1)ⁿ integer vectors in that box, so it takes about n·log10(2H+1) digits of agreement. When that exceeds the digits the inputs have, the relation was purchasable from the available precision — the search would have found something whatever the constants were — and it is evidence of nothing. Because PSLQ returns the smallest relation the precision can buy, a purchased one lands just under the bound rather than over it, so a relation must clear it by margin_digits (default 10) to count as credible.

The relation is also evaluated. Σ aᵢ·cᵢ is computed exactly on the values supplied (a float or mpmath.mpf is an exact binary rational, a decimal string an exact decimal one). It must cancel to the precision the inputs carry — residual_digits — and the digits it holds to beyond what the coefficient size explains — excess_digits — must reach margin_digits. A relation that fails either test is credible: False whatever the declared precision: a cheap relation that does not hold is not rescued by being cheap.

The input's precision has to be knowable, and often it is not. A float is 53 bits however it is printed, so it is judged. int and Fraction are exactly the rationals they spell, so precision cannot be the reason a relation among them is spurious — instead the relation itself is evaluated (see exact_residual below). A decimal string is unknown: the digits "3.14159265358979" are equally the exact rational 314159265358979/10¹⁴ and π truncated to 15 places, and nothing in the string says which. An mpmath.mpf is unknown for the same reason: its accuracy is not recoverable from the object. Unknown is not unbounded, though: a literal's length is a ceiling on what it can carry (a 20-digit string holds at most 20 digits, an mpf at most its mantissa's bits), which is enough to refute a relation. Pass digits= (decimal) or precision_bits= (binary) to say how many digits are trustworthy; a declaration is a cap, so declaring 200 digits for a float still yields ~16, and for a 20-digit string still yields 20.

Returns a dict:

credible True, False, or None for unknown — the precision of the inputs could not be established, so no judgement was made. None is never a pass: treat it as "not yet checked" and re-run with digits=. A relation can still be refuted (False) with the precision unknown. available_digits Digits the inputs are known to carry (inf for exact inputs), or None when unknown. consumed_digits Digits of agreement the relation costs, n·log10(2H+1). spare_digits available - consumed, or None when unknown. margin_digits Spare digits demanded of a credible relation. precision_source Where available_digits came from: "float", "exact", "declared", "content" (a declaration capped by the length of the literals), "ball", or "unknown". residual_digits Digits to which Σ aᵢ·cᵢ cancels relative to Σ |aᵢ·cᵢ| (inf for an exact zero), or None when it cannot be evaluated. excess_digits min(residual_digits, ceiling) - consumed_digits: agreement the coefficient size does not explain. A credible relation has at least margin_digits of it. exact_residual Σ aᵢ·cᵢ as an exact Fraction when every constant is an exact rational and no precision was declared, otherwise None. A nonzero value refutes the relation, and is the one way credible can be False while available_digits is inf: nothing is short of precision, the relation is just false.

Note what credible: True means on the exact branch: the relation holds for the numbers supplied. If those numbers are rounded stand-ins for something else — Fraction(*float(pi).as_integer_ratio()) is exactly the double, not π — then so is the verdict, and no gate reading only the constants can tell. Declare digits= when they are stand-ins; that is what turns them back into approximations with a precision to be judged.

Python

>>> import alkahest as ak
>>> ak.relation_confidence([1.0, 2.0, 3.0], [1, 1, -1])["credible"]
True
>>> ak.relation_confidence([0.1, 0.2, 0.7], [60771139, 67263243, 11653676])["credible"]
False

Exact inputs are judged by evaluating the relation, not by counting digits:

Python

>>> ak.relation_confidence([1, 2, 3], [1, 1, -1])["credible"]
True
>>> ak.relation_confidence([1, 2, 3], [1, 1, 1])["credible"]
False

A decimal string is not endorsed until its precision is declared, but a relation that does not hold is refuted either way:

Python

>>> sqrt2 = "1.41421356237309504880168872420969807856967187537694806841301"
>>> two_sqrt2 = "2.82842712474619009760337744841939615713934375075389613682602"
>>> print(ak.relation_confidence([sqrt2, two_sqrt2], [-2, 1])["credible"])
None
>>> ak.relation_confidence([sqrt2, two_sqrt2], [-2, 1], digits=60)["credible"]
True
>>> ak.relation_confidence([sqrt2, two_sqrt2], [-3, 1], digits=60)["credible"]
False

A ten-digit-per-coefficient relation among 20-digit strings is exactly what 20 digits can buy, and declaring more digits than the strings hold does not change that:

Python

>>> pi_20 = "3.1415926535897932385"
>>> e_20 = "2.7182818284590452354"
>>> big = [5144503108, -5945642943]
>>> ak.relation_confidence([pi_20, e_20], big, digits=60)["credible"]
False

Rigorous bounds

function

bound_on_box(
    expr,
    box,
    *,
    order=6,
    prec=128,
    tol=1e-9,
    max_subdivisions=2048
)

#source

alkahest.bound_on_box(expr, box, *, order=6, prec=128, tol=1e-9, max_subdivisions=2048)

Rigorous enclosure of the range of expr over an axis-aligned box, via Taylor models plus Moore–Skelboe branch-and-bound.

box is a list of (variable, lo, hi). The returned Enclosure is a guaranteed outer bound; when the work budget runs out it is returned anyway with budget_exhausted=True — wide but true.

function

bounds_supported(expr)

#source

alkahest.bounds_supported(expr) -> BoundsSupport

Can the validated-bounds subsystem bound this expression at all?

bound_on_box(), verified_integral(), verified_no_roots() and verified_sign() all evaluate through the same Taylor-model evaluator, and it refuses any construct it has no rigorous rule for with E-VALIDATED-001. This answers that question without running the bound, so a planning loop can pick a certifiable route instead of discovering the boundary by hitting it.

The answer is produced by running the real evaluator on a probe box, so it cannot drift from what bound_on_box() does. The same information per primitive is the taylor_model flag in capabilities()["primitives"].

What it does not promise. True means no E-VALIDATED-001. A supported function can still be refused on a particular box for a domain violation (E-VALIDATED-003, e.g. log on [-2, -1]) or a non-finite enclosure (E-VALIDATED-004) — those depend on the box, not on the expression, so no box-free predicate can rule them out.

Note this is a different question from certifiable(), which asks whether an operation emits a Lean certificate. A rigorous enclosure is not a Lean proof term, and the validated subsystem is not in the certificate ledger; conflating the two under one predicate would make a True mean two different kinds of evidence.

Examples

Python

>>> import alkahest as ak
>>> p = ak.ExprPool(); x = p.symbol("x")
>>> bool(ak.bounds_supported(ak.sin(x) * ak.exp(x)))
True
>>> bool(ak.bounds_supported(ak.bessel_j0(x)))   # covered since 3.9.0
True
>>> answer = ak.bounds_supported(ak.floor(x))
>>> bool(answer), answer.functions
(False, ['floor'])

class

class BoundsSupport

#source

Whether the validated-bounds subsystem can bound an expression at all.

Returned by bounds_supported(). Truthy exactly when every construct in the expression has a rigorous Taylor-model rule, so it drops into if ak.bounds_supported(f):.

Attributes

supportedbool: The verdict. bool(self) is the same value.
blockerstr or None: The evaluator's own description of the first construct it has no rule for ("function `bessel_j0`"), or None when supported.
functionslist of str: Every function in the expression with no Taylor-model rule, sorted. Empty when supported is True — and possibly empty when it is False, if the blocker is a node kind rather than a function.
detailstr: One-sentence human explanation.

Properties

supported

#source

blocker

#source

functions

#source

detail

#source

Methods

as_dict()

#source

The verdict as a plain dict, for logging and JSON.

class

class Enclosure

#source

A rigorous enclosure: the true value is guaranteed to lie in [lower, upper].

Returned by bound_on_box() and verified_integral(). The bound may be wide, but it is never wrong — that is the whole contract. Check budget_exhausted to see whether it is as tight as you asked for.

Properties

lower

#source

upper

#source

budget_exhausted

#source

subdivisions

#source

width

#source

Width of the enclosure.

Methods

contains(v)

#source

True if v is inside the rigorous enclosure.

function

verified_integral(
    expr,
    var,
    a,
    b,
    *,
    order=8,
    prec=128,
    tol=1e-9,
    max_subdivisions=4096
)

#source

alkahest.verified_integral(expr, var, a, b, *, order=8, prec=128, tol=1e-9, max_subdivisions=4096)

Rigorous enclosure of ∫_a^b expr d(var) by Taylor-model quadrature.

Unlike integrate_definite(), this never needs a closed form — and unlike a floating-point quadrature, the answer is a theorem: the true integral is guaranteed to lie in the returned interval. Refuses on singular or improper integrands rather than guessing.

A removable singularity is not a refusal: an integrand written as N(x)/D(x) whose two halves vanish at the same point p, to the same order, with some D**(d)(p) != 0 — log(1+x)/x on [0, 1], sin(x)/x on [-1, 1], (1-cos(x))/(x*x) on [0, 1] — is enclosed via Cauchy's mean value theorem, iterated as many times as the order of the zero requires, and the value returned is the integral of the continuous extension. The zeros are checked symbolically, so a genuine pole is never mistaken for a removable one.

A bounded integrand whose Taylor model runs out of domain rather than out of function is not a refusal either. asin and acos on [0, 1], sqrt(1 - x*x) on [0, 1], sqrt(x) and x**x on [0, 1] are all finite and continuous on the closed interval, but the Taylor rule needs a derivative bound it does not have at the end of the domain, so no model exists on the last panel however far it is bisected. That panel is closed with a width * range bound computed by interval arithmetic, which needs no derivative. The panel is of order 2**-60 of the interval by then, so the crude bound costs essentially nothing.

An unbounded integrand is still refused, whether or not the integral converges (-log(x) on [0, 1], 1/sqrt(1-x*x) on [0, 1]): the integral exists, but no rigorous enclosure of the integrand does, and no bounded range exists either. The ValidatedError message names the Taylor rule that stopped and where, and does not claim the integral fails to exist.

function

verified_no_roots(
    expr,
    box,
    *,
    order=6,
    prec=128,
    tol=1e-9,
    max_subdivisions=2048
)

#source

alkahest.verified_no_roots(expr, box, ...) -> "true" | "false" | "undecided"

Three-valued: "true" proves expr has no root anywhere in the box, "false" proves it does, and "undecided" means neither could be established within the budget. The third is never collapsed into the other two.

"false" is certified by the intermediate value theorem: the box is first proven free of poles and branch cuts (so expr is continuous on it), then the box is subdivided until two points are found at which expr is rigorously proven to have opposite signs. A box is convex, so the segment joining them stays inside it and expr must vanish somewhere along it. Because the two points need not be the box's own endpoints, an even number of roots no longer defeats the test — x*x - 2 on [-2, 2] is "false".

A root that never produces a sign change — a double root such as (x-1)**2 — stays "undecided": no witness exists, and none is invented.

function

verified_sign(
    expr,
    box,
    predicate,
    *,
    order=6,
    prec=128,
    tol=1e-9,
    max_subdivisions=2048
)

#source

alkahest.verified_sign(expr, box, predicate, ...) -> "true" | "false" | "undecided"

predicate is one of "positive", "negative", "nonnegative", "nonpositive". A "false" verdict is itself certified — either the enclosure proves the predicate fails everywhere, or a rigorously evaluated point witnesses the failure.

An inequality that is tight at an endpoint of the box is still decided. Subdivision alone cannot do it — where the margin vanishes, every enclosure of the range straddles zero — so the box is split: a collar at the endpoint is handled by a truncated Taylor expansion with a proven Lagrange remainder, the rest by branch-and-bound. That covers the classical sharp trigonometric inequalities (Cusa–Huygens, Mitrinović–Adamović, Wilker, Huygens, Jordan) on boxes reaching x = 0. Tightness in the interior is not covered and stays "undecided".

"undecided" is not monotone in the box: it is a search verdict, not a logical one, and a box that stops just short of the point an inequality is tight at gets a much weaker collar than one that reaches it. A box that stops short of x = 0 is retried with the collar planted at 0 — a "true" on the larger box implies "true" on the box asked about — so the classical inequalities decide on [1e-12, pi/2] as well as on [0, pi/2]. For a statement tight somewhere else the effect remains: on "undecided", try the larger box whose endpoint is the tight point before concluding the statement is out of reach.

tol sets the tolerance of the enclosure, not of the verdict: it is an absolute width, so it does not bound how close to zero the answer may be. Once the enclosure has been computed, the search is re-run with the sign itself as the stopping rule, which is what decides an inequality whose margin is narrower than tol.

function

verified_ode(
    rhs,
    vars,
    t,
    y0,
    t_span,
    *,
    order=20,
    prec=128,
    tol=1e-16,
    max_steps=20000,
    h_min=None,
    max_rel_width=1.0
)

#source

alkahest.verified_ode(rhs, vars, t, y0, t_span, *, order=20, prec=128, tol=1e-16, max_steps=20000, h_min=None, max_rel_width=1.0)

Validated integration of y' = f(t, y), y(t0) ∈ y0: returns boxes that are guaranteed to contain the true solution (every solution, when y0 is a box), or refuses with a coded ValidatedError.

Method: interval Taylor series of order order (coefficients by automatic differentiation of rhs in ball arithmetic), an a-priori enclosure per step validated by the Picard–Lindelöf operator, a Lagrange remainder bound, and Lohner's QR coordinate frame against the wrapping effect.

Parameters

rhs: list of Expr — f_i(t, y).
vars: list of state symbols y_i (same length as rhs).
t: the time symbol, or None for an autonomous system.
y0: initial values — each a float (a point) or a (lo, hi) pair.
t_span: (t0, t1) with t0 < t1.
order: Taylor order per step (1..60).
prec: working precision in bits.
tol: local error target (relative to max(1, |y|)); only steers the step size, never soundness.
max_steps: refuse with E-VALIDATED-008 beyond this many steps.
h_min: smallest step tried (default 1e-12 * (t1 - t0)).
max_rel_width: refuse with E-VALIDATED-007 when a component's width exceeds this times max(1, |y|).

Returns

: VerifiedOdeSolution.

Raises

Example

y = pool.symbol("y") sol = ak.verified_ode([-y], [y], None, [1.0], (0.0, 1.0)) lo, hi = sol.enclosure(1.0)[0] # contains exp(-1), width ~1e-25

class

class VerifiedOdeSolution

#source

Rigorous enclosures of an ODE solution, returned by verified_ode().

times[i] is an exact grid point and enclosures[i] a list of (lo, hi) intervals, one per state component, that is guaranteed to contain y(times[i]). Bounds are rounded outward to float. enclosure() encloses the solution at any t in the span.

Properties