Coupled polynomial systems are unsolvable: eliminate by Groebner basis, not by radicals

[Rafael-SOWNet]

EquationSystem.Solve cannot solve coupled polynomial systems. It dies at four coupled variables:

system	current
2 coupled equations, 2 vars	2 ms
3 coupled equations, 3 vars	2 500 ms
4 coupled equations, 4 vars	no answer in 300 s

This is not a size problem. Four uncoupled quartics producing 256 solutions take 17 ms, and a degree-8 univariate takes 12 ms.

Root cause

EquationSolver.InSolveSystem eliminates one variable at a time by calling SolveEquation, which applies closed-form radical formulas (quadratic, Cardano, Ferrari). With numeric coefficients that is cheap; with symbolic ones it is not:

polynomial	result size
x⁴ - 10x³ + 35x² - 50x + 24 (numeric coefficients)	5 nodes
x⁴ + ax³ + bx² + cx + d (symbolic coefficients)	5661 nodes

Each elimination turns the next step's coefficients into nested radicals, so size compounds multiplicatively, and it repeats per branch. Traced on the 4-variable case:

eq[0..3] (deg 1,2,3,4)                    19 nodes each
eliminate x_3 via the linear equation  →  25 nodes each
eliminate x_2 via the quadratic        →  109 nodes   (coefficients now radicals)
eliminate x_1 via the resulting cubic  →  3698 nodes, 14.1 s   for ONE branch

Gröbner bases avoid this: they triangularise so the last variable satisfies a univariate polynomial with rational coefficients, and back-substitution then works with numbers rather than symbols.

Measured with a prototype

A throwaway lex Buchberger over the same packed representation MultivariatePolynomial already uses (one byte per variable in a ulong, ERational coefficients — so lex monomial comparison is just ulong comparison, and LeadingMonomial() already exists):

system	current	prototype
4 coupled variables	>300 s	0 ms
5 coupled variables	>60 s	3 ms
6 coupled variables	>60 s	14 ms

The 4-variable case emits exactly x⁴ - 10x³ + 35x² - 50x + 24, which the existing PolynomialSolver already solves in 3 ms — so the end-to-end path is milliseconds.

Over-determined systems work natively, including inconsistency, which the current API cannot even express because it throws on the equation-count mismatch:

system	prototype
circle + line (2 eq, 2 var)	basis 2, correct
circle + line + a redundant consequence (3 eq, 2 var)	basis 2, same answer
circle + line + an inconsistent equation	basis {1} — the textbook "no solutions"

Sketch

• GroebnerSystemSolver in front of EquationSolver.SolveSystem. If every equation passes MultivariatePolynomial.TryParse (polynomial, rational coefficients, ≤ 8 variables) take the Gröbner path; otherwise fall through to today's InSolveSystem unchanged, so nothing that works now changes.
• Buchberger over lex, with the Gebauer–Möller pair criteria. Lex because the packed ulong is the lex comparison and lex directly yields the triangular form. The prototype's criteria skipped 90–96 % of pairs.
• Back-substitution hands the univariate to the existing PolynomialSolver, inheriting its rational → radical → numeric convention rather than introducing a RootOf node.
• A hard budget on pairs, basis size, intermediate terms and wall time, falling back rather than hanging. Buchberger is doubly exponential in the worst case.
• Drop #eqs == #vars on the Gröbner path.

Two things the prototype disproved

Encoding dominates, and the implementation must pass the user's equations through unmodified rather than trying to be clever:

encoding	prototype
degree-4 membership generators + 2 power sums, 5 vars	gave up after 100 s
the same system as low-degree power sums	3 ms
sum-of-squares compression into fewer equations (degree 8)	gave up after 167 s
the same constraints uncompressed	0 ms

Known gaps

• Back-substitution degrades to floating point once a root is not expressible in radicals, and deciding whether a residual is zero then becomes a tolerance question — the PrecisionErrorZeroRange trap from Precision is being lost somewhere #602. The real fix is algebraic-number arithmetic, which the library does not have; that is arguably the more foundational piece and worth its own issue.
• The 8-variable cap of MultivariatePolynomial stays, so this covers systems users actually type, not large combinatorial ones.
• Coefficient growth was not measured by the prototype.

Additive, so 2.1 rather than 2.0.

[Rafael-SOWNet]

Correcting the sketch above: lex is the wrong order, and the prototype now says so with numbers.

The lex recommendation rested on the packed ulong exponent vector being the lex comparison — free, and directly triangular. That held on everything measured at the time, but everything measured at the time was sparse and highly structured: power sums and a Sudoku encoding. That is the easy case for an elimination order, and it does not generalise.

On cyclic-n and Katsura-n, the two families a Gröbner engine is normally benchmarked with, lex dies at five variables:

system	vars	lex	DegRevLex
cyclic-4	4	0 ms	1 ms
cyclic-5	5	gave up, 87-digit coefficients	400 ms, 6-digit
cyclic-6	6	gave up, 196-digit	gave up, 123-digit
katsura-3	4	8 ms, 15-digit	0 ms, 4-digit
katsura-4	5	gave up after 143 s, 188-digit	9 ms, 5-digit
katsura-5	6	gave up, 160-digit	258 ms, 7-digit
katsura-6	7	gave up, 81-digit	16.9 s, 15-digit

The failure mode is not the one the sketch guarded against

Katsura-4 in lex gave up after considering 53 pairs, never exceeding 64 terms, with coefficients at 188 digits. DegRevLex finishes the same system in 9 ms with 5-digit coefficients.

The term count was never the problem — the rationals were. "Coefficient growth was not measured" was filed above as a minor gap; it is in fact the entire failure mode. A budget that bounds pairs, basis size and terms — as the sketch proposed — would not have noticed this system dying, because none of those three ever got large.

Consequences for the plan

• DegRevLex, then FGLM to convert to lex for back-substitution. That is a second algorithm, which the sketch explicitly wanted to avoid; the measurement says it is the price of the feature working on anything dense.
• The order abstraction is needed after all. DegRevLex cannot be a single integer compare: total degree, then a scan for the last differing variable. It costs about 4× on the sparse systems where lex was already fine (7 coupled variables: 76 ms → 288 ms). Losing three orders of magnitude on dense input to save that factor of four is not a trade worth making.
• The budget must bound coefficient size, not just terms and pairs.
• The budget must be checked inside the arithmetic. Katsura-4 ran 143 s against a 45 s limit — 3.2× overshoot — because the check sits at loop heads and one operation on 188-digit rationals runs straight through it. That was a hypothesis in the sketch; it is now a measurement.

Still unmeasured

• FGLM's cost. It is now on the critical path for back-substitution and has not been prototyped at all. Largest remaining unknown.
• Multi-modular arithmetic. DegRevLex still cannot do cyclic-6. Running Buchberger modulo several primes and reconstructing attacks swell directly rather than avoiding it, and would likely lift the remaining ceiling. A follow-up, not a precondition.

None of this changes the case for the issue. Today's solver dies at four coupled variables; the prototype does seven sparse in 76 ms and an over-determined 4×4 Sudoku in ~100 ms. What changed is the claimed reach, not the direction.

The prototype is a harness now rather than scratch code, with both orders behind a flag:

dotnet run --project work/groebner -c release dense           # lex
dotnet run --project work/groebner -c release dense degrevlex