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.

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 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 and raises E-SOLVE-006.

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:

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.

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 + R1*C*s)**-1}]

solve_side_conditions()
# → ['(1 + (R1 * C * s)) ≠ 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:

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

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 + (C * s * 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 raises SolverError with .code == "E-SOLVE-005" rather than being cleared.

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:

eom = [...]                      # Euler–Lagrange equations of a double pendulum
solve(eom, [qdd1, qdd2])         # → [{qdd1: ..., qdd2: ...}], with sin/cos coefficients
solve_side_conditions()          # → the mass-matrix determinant ≠ 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().

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

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

Lambert W is available on its principal branch only. Where the real solution set also has a root on W₋₁ — x·eˣ = c for −1/e < c < 0, x**x = c for e^{−1/e} < c < 1, A'/B < 0 in the log shape — the call is refused rather than answered with half the roots; for a symbolic c the condition that makes W₀ complete (c > 0) is recorded.

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:

solve([sin(x) - half], [x])
# → [{x: ((k * pi * 2) + asin(1/2))}, {x: (pi + (k * pi * 2) + (-1 * asin(1/2)))}]
solve_side_conditions()       # → ['k ∈ Integer']
solve([tan(x) - 2], [x])      # → [{x: ((k * pi) + atan(2))}]
solve([sin(x) - half], [x], numeric=True)   # → [{x: 0.5236}, {x: 2.618}]  (k = 0)

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.

Anything else — two kernels that do not reduce (sin(x) + exp(x)), a kernel with a non-affine argument, a system with an unknown inside a function — keeps the coded refusal E-SOLVE-001. There is no numeric root-finder behind numeric=True for those.

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.

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]

gb.specialize([0, 1, 1, 1])              # E-PARAMGB-004: a = 0 is on the locus
gb.specialize([0, 1, 1, 1], verify=True) # [x, y - 1] — the refusal was unnecessary
gb.specialize([1, 1, 2, 2], verify=True) # 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:

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

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

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.

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.