Polynomial system solving

Alkahest solves systems of polynomial — and rational — 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.

import alkahest as ak
from alkahest import ExprPool, Domain, solve, solve_proof_status, solve_side_conditions, nsolve
from alkahest import sqrt, exp, sin, cos, tan

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:

circle = solve([x**2 + y**2 - pool.integer(1), y - x], [x, y])
for sol in circle:
    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 numbers directly: solve(eqs, vars, numeric=True). The domain is ℂ, so a real root is a float and a non-real one a Python complex (solve([s**2 + 1], [s], numeric=True) gives ±1j; domain="real" drops them). Where the symbolic path refuses with E-SOLVE-002 (see below), the numeric path falls back to complex homotopy continuation and returns every finite root or raises E-HOMOTOPY-004. A solution that depends on a free parameter has no numeric value, so it is returned as its exact Expr (solve([x - y], [x], numeric=True) → [{x: y}]) and solve_side_conditions() says so; where there is no symbolic solution to fall back to (the homotopy path) it raises E-SOLVE-006.

#Symbolic parameters: generic solutions, checked

A solution in the parameters is the answer for generic parameter values, and every tuple is checked as a function of them: substituted back into the equations, it must leave a residual that is identically zero (exactly, for rational functions of the parameters; in ball arithmetic at sample points otherwise), or it is dropped. An equation that factors is solved as the union of its factors, so equilibrium systems come out complete:

a, b, c, d = (pool.symbol(n) for n in "abcd")

solve([a*x - b*x*y, -c*y + d*x*y], [x, y])
# → [{x: 0, y: 0}, {x: c/d, y: a/b}]
solve_side_conditions()     # → ['b ≠ 0', 'd ≠ 0']   (empty when b, d are Positive)
solve([x - a, x - b], [x])  # → []   with 'a - b ≠ 0' stated
solve([x*(y - a), x*(y - b)], [x, y])   # → GroebnerBasis: the line x = 0

The side conditions are exactly the non-vanishing hypotheses of the branches that were kept, factored; a coefficient that is identically zero on a branch is never divided by.

#The unknown's declared domain

solve honours a Positive, NonNegative, NonZero or Integer unknown: values provably outside the domain are dropped, and a membership that depends on a parameter is stated rather than assumed.

xp = pool.symbol("xp", Domain.Positive)
solve([xp + 3], [xp])        # → []
solve([xp**2 - 4], [xp])     # → [{xp: 2}]
solve([xp - a], [xp])        # → [{xp: a}];  solve_side_conditions() → ['a > 0']

Domain.Real is the default for every symbol, so it is not treated as a declaration: use domain="real" to drop complex roots.

#Degree three and above: factor, then RootOf

A univariate of degree ≥ 3 — a single equation, or the Lex tail of a system — is factored first. With rational coefficients it is factored over ℚ (FLINT), and each irreducible factor is solved by the first rule that applies:

s = pool.symbol("s")
solve([s**3 + 6*s**2 + 11*s + 6], [s])   # → s = -1, -2, -3
solve([x*(x**2 - 5)], [x])               # → 0, ±√5
solve([x**5 - x - 1], [x])               # → RootOf(x^5 - x - 1, k), k = 0..4

RootOf(f, k) is the k-th root of the irreducible integer polynomial f in a fixed order: real roots ascending, then non-real roots by real part, ties by imaginary part (a conjugate pair lower half-plane first). It is exact and has no free symbols — it stores f by its coefficients — so diff treats it as a constant and subs cannot reach into it. It evaluates through eval_expr (real roots), evaluate(..., mode="complex"), and to any precision through interval_eval / evaluate(..., precision_bits=...) (real roots), using FLINT's certified root isolation. Build one directly with root_of(poly, var, index), which reduces to the irreducible factor the root belongs to and raises E-SOLVE-008 on a non-polynomial input or an out-of-range index. A printed RootOf(x^5 - x - 1, 2) parses back to the same node (parse; the polynomial's variable is local to the call).

With symbolic parameters the univariate is factored over ℚ[params]: s**3 - a*s**2 - s + a gives a, ±1. A factor of degree ≥ 3 in the unknown that still involves a parameter (and is not even over a quadratic) has no closed form here and raises E-SOLVE-002, as does a degree-≥3 step whose coefficients contain algebraic numbers from an earlier step. Repeated roots are listed once.

#Rational equations

An equation may be a ratio of polynomials in the declared unknowns. This is the form nodal analysis produces — every admittance term is a reciprocal — so writing the circuit equation directly is the point:

Vo, Vin, R1, C, s = (pool.symbol(n) for n in ("Vo", "Vin", "R1", "C", "s"))

solve([(Vo - Vin) / R1 + Vo * s * C], [Vo])
# → [{Vo: Vin * (1 + s*R1*C)**-1}]

solve_side_conditions()
# → ['(1 + (s * R1 * C)) ≠ 0', 'R1 ≠ 0']

Each equation is put over a common denominator and the numerator system is solved.

#Spurious roots

Clearing a denominator is not an equivalence:

N/D = 0     ⟺     N = 0  and  D ≠ 0

Multiplying up drops the second conjunct, and the root it leaves behind is not a near-miss — it is the one point where the equation has no value at all. x/(x-1) = 1/(x-1) clears to (x-1)**2 = 0, whose only root is x = 1, which is exactly where both sides read 1/0:

# not x = 1
solve([x / (x - 1) - 1 / (x - 1)], [x])   # → []
# a reciprocal is never zero
solve([1 / x], [x])                       # → []
# 0/0 at x = 0
solve([x * x / x], [x])                   # → []
# 1/x has no value at x = 0
solve([1 / (1 / x - 1)], [x])             # → []

So every candidate root is re-tested against the condition that says where the equation has a value at all, and one that makes it vanish is dropped. The test is exact where it can be — rational arithmetic, then the ambient assumptions, then rigorous ball arithmetic — and surviving roots are additionally substituted into the equations as you wrote them, before anything was cleared.

Two things that condition is deliberately not:

• It is not the cancelled denominator. x**2/x reduces to x, whose root x = 0 is precisely the point the original expression is 0/0 at, so the uncancelled form is what gets carried.
• It is not the product of the denominators either. A reciprocal swaps the two halves — (n/d)**-1 = d/n — so an inner denominator becomes an outer numerator and would drop out of such a product. In 1/(1/x - 1) the requirement x != 0 disappears exactly that way, leaving 1 - x (non-zero at x = 0) as the only recorded condition and x = 0 looking like a root. What is carried instead is the domain: everything that has to be non-zero for the expression to denote a number, which for that example is x*(1 - x) — both of the points it is undefined at.

#When the exclusion is undecidable

With symbolic parameters, whether a denominator vanishes at a root may be undecidable — R1 ≠ 0 is not something the solver gets to assume. The root is returned under a stated hypothesis instead, listed by solve_side_conditions() alongside the non-vanishing leading coefficients the back-substitution divided by. Declaring the symbol positive discharges the condition rather than reporting it:

R1 = pool.symbol("R1", "positive")
solve([(Vo - Vin) / R1 + Vo * s * C], [Vo])
solve_side_conditions()          # ['(1 + (s * C * R1)) ≠ 0']   — no 'R1 ≠ 0'

An empty list means every divisor was proved non-zero, not that none was looked at.

When no finite solution list can be produced and a GroebnerBasis comes back instead, that basis is the numerator ideal — no root of it has been tested against the denominators — and the cleared denominators are reported as side conditions so the distinction is visible.

An equation whose denominator is identically zero (1/(x - x)) denotes no function at all and is refused rather than cleared. From Python the refusal comes even earlier: building 1/(x - x) already raises DivisionByZeroError (E-DOMAIN-001). An equation that reaches the solver with such a denominator (from Rust, say) raises SolverError with .code == "E-SOLVE-005".

method="homotopy" and solve_numerical are not covered by this: they take polynomial systems only, and a rational equation still raises E-SOLVE-001 there.

#Beyond polynomials

#Non-polynomial coefficients

Only the unknowns have to appear polynomially. A sub-expression free of them — sin(q1), cos(q1 - q2), sqrt(a), exp(-k*T) — is a coefficient, whatever it is. That is the shape of every multibody equation of motion solved for the accelerations:

q1, q2, qd1, qd2, qdd1, qdd2, g = (
    pool.symbol(n) for n in ("q1", "q2", "qd1", "qd2", "qdd1", "qdd2", "g"))
# Euler–Lagrange equations of a double pendulum (unit masses and lengths)
eom = [
    2*qdd1 + cos(q1 - q2)*qdd2 + sin(q1 - q2)*qd2**2 + 2*g*sin(q1),
    cos(q1 - q2)*qdd1 + qdd2 - sin(q1 - q2)*qd1**2 + g*sin(q2),
]
solve(eom, [qdd1, qdd2])         # [{qdd1: ..., qdd2: ...}], with sin/cos coefficients
solve_side_conditions()          # the mass-matrix determinant: ['(2 + (-1 * cos((q1 + (q2 * -1)))^2)) ≠ 0']

Each such sub-expression is replaced by an opaque parameter, the system is solved, and the sub-expressions are substituted back. Treating sin q and cos q as independent loses identities such as sin² q + cos² q = 1, so nothing is returned unchecked: every hypothesis is zero-tested after substitution, every solution is substituted into the equations as written — in ball arithmetic when it has no free symbol, otherwise at several generic parameter values — and a solution that only exists because an identity was invisible ((sin² a + cos² a − 1)·x = 1) is refused with E-SOLVE-001. An unknown inside a function (sin(x) while solving for x) is not a coefficient; a single such equation goes to the transcendental slice below, a system still refuses.

#Transcendental equations in one unknown

One equation in one unknown with exp, log, powers with the unknown in the exponent, or sin/cos/tan is solved in closed form when it has one of these shapes. The slice reports real solutions (exp(x) = -1 has none), and its hypotheses are listed by solve_side_conditions().

ka, ke, t = pool.symbol("ka"), pool.symbol("ke"), pool.symbol("t")
solve([ka*exp(-ka*t) - ke*exp(-ke*t)], [t])   # → [{t: (log((ke * ka^-1)) * (ke + (-1 * ka))^-1)}]
solve_side_conditions()
# → ['ka ≠ 0', '(ke + (-1 * ka)) ≠ 0', '(ke * ka^-1) > 0']

Declaring ka and ke positive (pool.symbol("ka", "positive")) discharges the sign hypothesis.

Lambert W has both real branches: lambert_w (W₀, values ≥ −1) and lambert_wm1 (W₋₁, real on [−1/e, 0), values ≤ −1; also lambert_w(x, k=-1)). Where the real solution set has a root on each — x·eˣ = c for −1/e < c < 0, x**x = c for e^{−1/e} < c < 1, A'/B < 0 in the log shape — both are returned:

solve([x*exp(x) + pool.rational(1, 4)], [x])
# → [{x: lambert_w(-1/4)}, {x: lambert_wm1(-1/4)}]      (−0.3574…, −2.1533…)

Below −1/e there is no real root ([]). A symbolic argument known to be negative returns both branches and records c + exp(-1) > 0; one of unknown sign records c > 0 and returns the principal root.

Periodic families. A trigonometric equation has infinitely many roots; they are returned as families in an integer symbol (named k, or n, m, j, k_1, … when k is taken), which solve_side_conditions() lists as k ∈ Integer:

half = pool.rational(1, 2)
solve([sin(x) - half], [x])
# → [{x: ((k * pi * 2) + (pi * (1/6)))}, {x: ((k * pi * 2) + (pi * (5/6)))}]
solve_side_conditions()       # → ['k ∈ Integer']
solve([tan(x) - 2], [x])      # → [{x: ((k * pi) + atan(2))}]
solve([sin(x) - half], [x], numeric=True)   # → [{x: 0.5235987755982989}, {x: 2.617993877991494}]

With numeric=True the k = 0 member of each family is returned, and solve_side_conditions() says so; it is one representative per family, not the solution set.

Trigonometric polynomials. An equation polynomial in sin, cos and tan of integer multiples of one affine angle θ is rewritten in S = sin θ, C = cos θ (multiple angles by de Moivre, tan = sin/cos with its denominator cleared), factored over ℚ, and each factor is solved on the circle S² + C² = 1: through sin θ = r or cos θ = r when it reduces to one of them, through tan θ = r when all its terms have degrees of one parity (homogenise with S² + C²; θ = π/2 is checked separately), and otherwise through t = tan(θ/2) — with θ = π, the point that substitution cannot reach, checked directly. Families are merged (π/4 + 2πk and 5π/4 + 2πk become π/4 + πk), offsets that are rational multiples of π are written as such when substitution proves it, a family on which a tan of the equation is undefined is dropped, and every family is substituted back:

solve([2*sin(x)**2 - 1], [x])       # → π/4 + πk, -π/4 + πk
solve([sin(x) - cos(x)], [x])       # → π/4 + πk
solve([sin(2*x) - cos(x)], [x])     # → π/2 + πk, π/6 + 2πk, 5π/6 + 2πk
try:
    solve([sin(x)**2 + cos(x)**2 - 1], [x])
except ak.SolverError as e:
    print(e.code)                   # E-SOLVE-001: an identity

Unknowns inside invertible functions. In a system, an unknown that occurs only inside one exp, log, sin, cos, tan of an affine argument, or a power b**(a*x + c) with a constant base (exp(2x) and exp(x) count as one kernel), is replaced by the function value, the polynomial system is solved, and each value is inverted: exp(f) = u needs u > 0 over the reals (a value proved otherwise drops the tuple, an undecided one is recorded as u > 0), the trigonometric ones give periodic families. solve([exp(x) + y - 3, exp(x) - y - 1], [x, y]) → x = log 2, y = 1.

Complex solutions. Transcendental equations are solved over the reals. With domain="complex" — or an unknown declared complex — exp(f) = c returns the family log(c) + 2πik for any c ≠ 0, sin(x) = 2 returns its two complex families, and a Lambert-W equation is refused (its complex roots lie on infinitely many branches). From Rust, wrap the call in solver::with_complex_domain.

Anything else — two kernels that do not reduce (sin(x) + exp(x)), a kernel with a non-affine argument, an unknown both inside and outside a function — keeps the coded refusal E-SOLVE-001; with numeric=True the message points at the numeric route below.

#Numeric roots with no closed form: nsolve and interval=

For a single equation with no closed form — an implied volatility BS(σ) = price, cos x = x, sin x = x/10 — nsolve finds real roots numerically and proves what it can about them:

r = nsolve(cos(x) - x, x, bracket=(0, 1))     # Brent on a sign change
# → NumericRoot(0.7390851332151607, enclosure=[0.7390851332151603, 0.739085133215161],
#               status='verified')
nsolve(cos(x) - x, x, x0=1.0)                 # damped Newton from a guess
iso = nsolve(sin(x) - x/10, x, all_in=(-12, 12))
# → RootIsolation(roots=[-8.4232…, …, 0.0, …, 8.4232…], unresolved=[], complete=True)

solve([sin(x) - x/10], [x], numeric=True, interval=(-12, 12))   # the same 7 roots

The search runs in f64; each enclosure is then checked in 128-bit outward-rounded ball arithmetic. "verified" means f changes sign across [lo, hi] and an enclosure of f′ there excludes zero, so exactly one root lies inside; "exists" drops the uniqueness half; "heuristic" is the bare iterate (a double root, a non-smooth equation). all_in subdivides the interval, discarding a piece only when an enclosure of f excludes zero and accepting one only by the monotonicity test, so its root list is complete exactly when unresolved is empty; a double root such as (x − 1)² is listed there instead of being guessed at. solve(..., numeric=True, interval=(a, b)) uses the same isolation, lists each enclosure in solve_side_conditions(), and refuses with E-SOLVE-012 when the list cannot be proved complete. Every symbol but the unknown must have a value (E-SOLVE-009 otherwise), and the loops honour an active Budget.

#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.

#Engines

GroebnerBasis.compute picks the engine by order; every engine returns the same reduced (monic, interreduced) basis, which is unique.

• "grevlex" — multimodular. The basis is computed modulo random word-size primes (below 2³¹) by a matrix F4 (all critical pairs of minimal sugar degree at once, symbolic preprocessing, row reduction of a sparse matrix in machine words), the images are combined by Chinese remaindering and lifted by rational reconstruction, and a lift is accepted only after one more prime agrees with it and it passes verification over ℚ: every input reduces to zero and every critical pair surviving the Gebauer–Möller criteria reduces to zero, i.e. the result is a Gröbner basis of an ideal containing the input. For homogeneous input that proves it is the basis; for inhomogeneous input it does not, and proof_status says so. A prime whose image has a different shape (an unlucky prime) is outvoted. If the path is not conclusive, Buchberger over ℚ is the fallback. A single generator, a constant among the generators, or a system of degree at most one goes straight to Buchberger, which is just as fast there and is a proof.
• "lex" — for a zero-dimensional ideal, the grevlex basis above is converted by FGLM, itself run modulo primes, lifted, and then verified exactly: each lifted polynomial is shown to lie in the ideal of the grevlex basis through exact normal forms, and its staircase has exactly as many monomials as that ideal's quotient has dimensions, so the lex basis is as sure as the grevlex basis it came from. Positive-dimensional ideals use Buchberger over ℚ directly.
• "grlex" — Buchberger over ℚ (sugar selection, Gebauer–Möller criteria, fraction-free reduction).
• GroebnerBasis.compute_f5 — a signature-based algorithm of the F5 family (incremental, position-over-term signatures, syzygy and rewrite criteria) in exact rational arithmetic. Under "lex" a multi-variable zero-dimensional ideal goes through grevlex and FGLM like compute; anything else runs F5 directly under lex.

With the parallel feature (in the default wheels) the multimodular engines compute several primes at once and spread the verification over the cores. Set RAYON_NUM_THREADS=1 for a single-threaded run; an active Budget is honoured on every worker.

#How sure is this answer?

GroebnerBasis.proof_status is "proved" or "probabilistic", set from the path the computation actually took.

gb = GroebnerBasis.compute([x**2 - y**2, x*y], [x, y], order="grevlex")  # homogeneous
gb.proof_status                                  # 'proved'

gb = GroebnerBasis.compute([x**2 + y**2 - pool.integer(1), x*y - pool.integer(1)], [x, y], order="grevlex")
gb.proof_status                                  # 'probabilistic'

The failure model. Since I ⊆ ⟨B⟩, equal leading terms would give equal ideals, so a wrong answer needs a majority of the primes used, at least two and each drawn independently from the roughly 5·10⁷ primes in [2³⁰, 2³¹), to give leading terms that I does not have. Such primes divide a fixed non-zero integer that depends only on the input, so if k primes in that range are unlucky for your input, the chance is of order (k / 5·10⁷)². The primes are drawn after the input is fixed, so no input can arrange it in advance, though enormous coefficients do raise k.

What a probabilistic basis still proves. gb.contains(p) == False is a proof that p is not in I. gb.contains(p) == True is only as sure as the basis. Every zero of B is a zero of I, but I could have more.

solve. A returned tuple has been substituted back into your equations whichever status applies, so it is a solution. The status is about whether the list is complete:

sols = solve([x**2 + y**2 - pool.integer(1), x - y], [x, y])
solve_proof_status()                             # 'probabilistic'

solve([x**2 + y**2 - pool.integer(1), x**2 + y**2 - pool.integer(2)], [x, y])  # []
solve_proof_status()                             # 'proved'

An empty answer from an inconsistent system is proved. {1} passes the verification over ℚ without testing anything, so it gets its own proof: the system is homogenised with a new variable and run through the same engine, where homogeneous input is proved, and 1 is in the ideal exactly when the result contains a pure power of that variable. On cyclic-5 and 6 and katsura-5 to 7 plus a generic quadric this took 0.02 s to 2.4 s, where Buchberger over ℚ took over 180 s. If that run is inconclusive, solve falls back to Buchberger, and an active Budget can stop it with BudgetExceededError. GroebnerBasis.compute(..., order="grevlex") does not fall back: it reports "probabilistic" for {1} in that rare case.

Needing a proof. compute_f5 and order="grlex" are proved, at their cost. From Rust, GroebnerBasis::compute_exact runs Buchberger under any order. Primary decomposition, triangularize and rosenfeld_groebner compute their bases exactly and are not affected by any of this.

#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 on the CPU engines above 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

one = pool.integer(1)

# 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, a] — 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
try:
    gb.specialize([-1])
except ak.ParamGroebnerError as e:
    print(e.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.

#specialize(values, verify=True)

Most of what conditions() lists are leading coefficients that were inverted somewhere inside the Buchberger loop and then cancelled, so a large share of the points it excludes are ordinary. Measured over eight parametric systems on the small-integer box {-2..2}: 646 refusals, of which 239 are genuine poles, 73 are necessary, and 334 are unnecessary. The fraction is box-dependent and dominated by parameters equal to exactly 0 — {-2..2} gives 52%, {-2,-1,1,2} gives 25%, {-3..3} gives 58% — so the honest headline is a quarter to a half of refusals on small-integer grids, not a single number.

gb = GroebnerBasis.compute([a*x + b*y - one, c*x + d*y - one], [x, y],
                           params=[a, b, c, d])
gb.conditions()                          # [c, a*d - b*c, a]

try:
    gb.specialize([0, 1, 1, 1])
except ak.ParamGroebnerError as e:
    print(e.code)                        # E-PARAMGB-004: a = 0 is on the locus
gb.specialize([0, 1, 1, 1], verify=True).to_exprs()   # [x, y - 1] — the refusal was unnecessary
try:
    gb.specialize([1, 1, 2, 2], verify=True)
except ak.ParamGroebnerError as e:
    print(e.code)                        # still E-PARAMGB-004: a*d - b*c = 0

verify=True re-solves the specialised system over ℚ and compares it with the specialised generic basis, instead of deciding on the recorded conditions alone. It is strictly more complete and never less sound: the refusals were separately checked to be sound (1,930 regular points across three sweeps, zero disagreements between specialize and the basis computed directly at the point), so this closes a completeness gap, not a correctness one. It is off by default because on the locus it pays for a second Gröbner basis over ℚ.

#Feeding a parametric basis back in

reduce and contains accept input that is rational in the parameters. A den**-1 factor is an ordinary element of Q(params), not a non-polynomial, so a basis accepts its own to_exprs() output — only a denominator in a ring variable is refused:

gb = GroebnerBasis.compute([a*x - one], [x], params=[a])
gb.to_exprs()                            # ['(x + (-1 * a^-1))']
gb.contains(gb.to_exprs()[0])            # True

equals_ideal answers the question that needs: do two parametric bases generate the same ideal of Q(params)[vars]?

g1 = GroebnerBasis.compute([a*x - one], [x], params=[a])
g2 = GroebnerBasis.compute([a*a*x - a], [x], params=[a])
g1.equals_ideal(g2)                      # True
g1.contains_ideal(g2)                    # True  — one direction only

This is exact, not generic: reduction over Q(params) only ever divides by non-zero field elements, so neither basis's conditions() enters the answer. Those conditions still bound what each basis says about a specialised parameter point — equal ideals over the fraction field can specialise differently on the locus. Bases over different numbers of variables or parameters compare False; there is no shared ring.

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, contains, contains_ideal and equals_ideal. 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.

#Differential elimination with the parameters in Q(params)

rosenfeld_groebner(dae, params=[...]) runs the prolongation loop of Rosenfeld–Gröbner over Q(params) instead of putting the model parameters in the ring, and returns a ParametricRosenfeldGroebnerResult whose final_basis() is a ParametricGroebnerBasis. That is the input–output relation of an ODE model computed end to end, rather than prolonged by hand:

from alkahest import DAE, rosenfeld_groebner

S, I, R, y0 = (pool.symbol(n) for n in ("S", "I", "R", "y0"))
b, g = pool.symbol("b"), pool.symbol("g")
dS, dI, dR, dy0 = (pool.symbol(n) for n in ("dS", "dI", "dR", "dy0"))

# SIR:  S' = -b*S*I,  I' = b*S*I - g*I,  R' = g*I,  y0 = I
dae = DAE.new([dS + b*S*I, dI - b*S*I + g*I, dR - g*I, y0 - I],
              [S, I, R, y0], [dS, dI, dR, dy0], t)

r = rosenfeld_groebner(dae, params=[b, g], eliminate=[S, I, R],
                       minimal=True, order="lex")
r.minimal_prolongation_rounds            # 2
states = ("S", "I", "R")
io = r.final_basis().eliminate([v for v in r.variables()
                                if str(v).lstrip("d").split("/")[0] in states])
io.to_exprs()                            # y0*y0'' - y0'^2 + b*y0^2*y0' + b*g*y0^3, over Q(b,g) (printed monic)

eliminate=[...] names the variables the caller intends to eliminate — typically the unobserved states. Their whole jet chain goes with them (x also names dx/dt, d2x/dt2, …), and the ranking is built with those jets first so the Lex elimination contract holds.

One prolongation too many is expensive out of all proportion. On the SIR model above, stopping at the first informative round costs 0.03 s and gives one 4-term relation; prolonging one round further does not finish in ten minutes. The over-supplied answer is correct, just enormously more expensive, so nothing else signals it. minimal=True stops at the first informative round; without it, a run that prolonged past one emits a UserWarning naming the round it could have stopped at, and minimal_prolongation_rounds records it either way.

Scope of "minimal". minimal_prolongation_rounds is "the first round at which eliminating those variables leaves a generator", not a theorem about the differential ideal. For a single-output model it coincides with the jet order the input–output relation needs; for multi-output models the criterion is known to be wrong, because one output can become informative several rounds before the others and the truncated basis then misses their relations. Treat it as a cost signal, not a certificate — and note that a minimal=True result is always flagged truncated.

#Which structural identifiability this decides

Where this pipeline is used for structural identifiability, it answers multi-experiment identifiability, not single-experiment: by Ovchinnikov–Pillay–Pogudin–Scanlon, Computing all identifiable functions of parameters for ODE models, Theorem 19, a function of the parameters is multi-experiment identifiable iff it is input–output identifiable, for general rational systems. That is why an IO-elimination route can call a model globally identifiable where a single-experiment tool such as SIAN reports "locally, not globally" — both verdicts are right, about different questions.

Hypothesis, not a guarantee. Theorem 19 requires the input–output equations to be the characteristic presentation of I_Σ ∩ C(θ){y,u}. An algebraic lex elimination at a hand-picked finite jet order — what the surface above does — is not guaranteed to produce one. So "this is a multi-experiment tool" holds under that standard; where it fails, the computed field can be a proper subfield of the identifiable one.

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

#Performance

Wall time in seconds for the standard Gröbner benchmarks, one process per run, on an AMD Ryzen Threadripper PRO 3955WX (shared machine, 8-core CPU quota). "default" is GroebnerBasis.compute with the parallel feature; "1 thread" is the same with RAYON_NUM_THREADS=1; both measured 2026-11-04 with a 120 s cap, as is compute_f5. SymPy is 1.14 groebner(..., method="buchberger") / method="f5b" with its cache cleared before each run, measured 2026-10-25 on the same machine. Basis sizes agree with SymPy's wherever SymPy finishes.

The "default" and "1 thread" columns on these systems are "probabilistic" and compute_f5 is "proved"; see How sure is this answer?. Lex bases of these systems are large: the katsura-6 lex basis has coefficients of about 6 500 bits over 6 500-bit denominators, so lifting it takes several hundred primes, and that, not the modular arithmetic, is what the lex rows measure.

#Generalized Pell equations: every solution class

diophantine(x**2 - D*y**2 - N, [x, y]) returns one particular solution of x² − D·y² = N plus the unit. But N can have several inequivalent fundamental solutions (x² − 10y² = 9 has (3, 0), (7, 2) and (13, 4)), and the one-class answer omits the rest. pell_classes(D, N) returns all of them:

from alkahest import pell_classes

r = pell_classes(10, 9)
r["unit"]      # (19, 6): fundamental solution of x² − 10y² = 1
r["classes"]   # [(3, 0), (7, 2), (13, 4)]

Every integer solution is x + y√D = ±(x_c + y_c√D)·(ux + uy√D)^k for exactly one class c and one k ∈ ℤ. Each representative satisfies √|N| ≤ x + y√D < √|N|·u, so a class and its conjugate are listed separately unless they coincide. negative_unit is the fundamental solution of x² − D·y² = −1 when one exists, and classes is empty when the equation has no solution. The search is complete (Nagell's bound, or the convergents of √D when |N| < √D); an instance too large to decide raises DiophantineError (E-DIOPH-003) rather than returning a partial list, and the active Budget is honoured. D must be a positive non-square and N ≠ 0.