Triangularise a polynomial system instead of eliminating it in radicals

BREAKING-CHANGES.md

@@ -61,6 +61,7 @@ read first.
 | loud | implicit `Entity[]` to `Entity` | made a `FiniteSet`, discarding order and repeats | removed |
 | **silent** | a `MathS.Settings` scope across an `await`, or inside a task | lost, or somebody else's | follows the call |
 | **silent** | a `RewriteRecording` across an `await`, or work started under it | lost, or somebody else's | follows the call |
+| loud | a polynomial system with more equations than unknowns | `WrongNumberOfArgumentsException` | solved |
 
 ---
 
@@ -261,6 +262,36 @@ single-threaded case, which is what `Simplify` is, is unaffected.
 
 Being off is still free: no recording open still costs one ambient read per rule set and
 allocates nothing, which `RewriteAllocationTest` continues to hold to.
+### An over-determined polynomial system is solved rather than refused
+
+`EquationSystem.Solve` insisted on as many equations as unknowns and threw otherwise:
+
+```csharp
+MathS.Equations("x^2 + y^2 - 25", "x + y - 7", "x*y - 12").Solve("x", "y");
+// was: WrongNumberOfArgumentsException
+// is:  the two solutions, (3, 4) and (4, 3)
+```
+
+The count was a consequence of how the old solver worked — it eliminated one variable per
+equation — and not of the problem. A Gröbner basis has no use for the equality, and an
+extra equation that happens to be a consequence of the others is not an error to report.
+
+An inconsistent system now reports itself as one, which it also could not do before:
+
+```csharp
+MathS.Equations("x^2 + y^2 - 25", "x + y - 7", "x*y - 99").Solve("x", "y");
+// was: WrongNumberOfArgumentsException
+// is:  null — no solutions
+```
+
+**What breaks.** Code that catches `WrongNumberOfArgumentsException` to detect a
+malformed system will no longer see it for the polynomial case. **Fewer** equations than
+unknowns still throws: a free variable means infinitely many solutions, which this does
+not enumerate.
+
+The relaxation only applies where the system is a polynomial one over `Q` in at most eight
+variables and its solutions are rational. Everything else reaches the previous solver
+exactly as before, including the equation-count check.
 
 ### `Minusf`'s two operands exchanged names
 

Sources/.editorconfig

@@ -24,6 +24,12 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed
 [Tests/UnitTests/Core/Transformations/*.cs]
 file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
 
+[AngouriMath/Functions/Algebra/Groebner/*.cs]
+file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
+
+[Tests/UnitTests/Algebra/GroebnerSystemTest.cs]
+file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
+
 [Tests/UnitTests/Common/ListArgumentOverloadTest.cs]
 file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
 

Sources/AngouriMath/Functions/Algebra/Groebner/Buchberger.cs

@@ -0,0 +1,321 @@
+//
+// Copyright (c) 2019-2026 Angouri.
+// AngouriMath is licensed under MIT.
+// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
+// Website: https://am.angouri.org.
+//
+
+using System;
+using System.Diagnostics;
+using PeterO.Numbers;
+
+namespace AngouriMath.Functions.Algebra.Groebner
+{
+    /// <summary>
+    /// What a Gröbner computation is allowed to spend before it gives up.
+    /// </summary>
+    /// <remarks>
+    /// Buchberger is doubly exponential in the worst case, so declining has to be reachable:
+    /// a caller waiting forever is a worse answer than "not this one". Four separate
+    /// ceilings because the ways it runs away are not the same — a system can blow up in the
+    /// number of pairs, in the size of one polynomial, in the width of the rationals while
+    /// everything else stays small, or in none of those while simply taking too long.
+    /// Coefficient width is here because it is the one that actually fires: a system has
+    /// been seen to give up with 53 pairs and 64 terms and coefficients 188 digits wide,
+    /// which no count of pairs or terms would have caught.
+    /// </remarks>
+    internal sealed class GroebnerBudget
+    {
+        internal int MaxPairs { get; init; } = 20000;
+        internal int MaxBasisSize { get; init; } = 500;
+        internal int MaxTerms { get; init; } = 20000;
+        internal int MaxCoefficientDigits { get; init; } = 400;
+        internal int MaxQuotientDimension { get; init; } = 512;
+        internal TimeSpan Limit { get; init; } = TimeSpan.FromSeconds(5);
+
+        internal string? Exceeded { get; private set; }
+
+        private readonly Stopwatch stopwatch = Stopwatch.StartNew();
+
+        internal bool Spend(string what)
+        {
+            if (Exceeded is not null)
+                return false;
+            if (stopwatch.Elapsed > Limit)
+            {
+                Exceeded = "time";
+                return false;
+            }
+            _ = what;
+            return true;
+        }
+
+        internal bool Allow(bool within, string what)
+        {
+            if (Exceeded is not null)
+                return false;
+            if (!within)
+            {
+                Exceeded = what;
+                return false;
+            }
+            return true;
+        }
+
+        internal bool CheckPolynomial(MultivariatePolynomial polynomial)
+            => Allow(polynomial.TermCount <= MaxTerms, "terms")
+            && Allow(polynomial.MaxCoefficientDigits() <= MaxCoefficientDigits, "coefficients");
+    }
+
+    /// <summary>Buchberger's algorithm, with the Gebauer–Möller pair criteria.</summary>
+    internal static class Buchberger
+    {
+        /// <summary>
+        /// A Gröbner basis of the ideal generated by <paramref name="generators"/>, or
+        /// <see langword="null"/> where <paramref name="budget"/> ran out.
+        /// </summary>
+        internal static List<MultivariatePolynomial>? Compute(
+            IReadOnlyList<MultivariatePolynomial> generators, MonomialOrder order, GroebnerBudget budget)
+        {
+            var basis = new List<MultivariatePolynomial>();
+            foreach (var generator in generators)
+                if (!generator.IsZero)
+                    basis.Add(generator.MakeMonic(order));
+
+            var pairs = new List<(int Left, int Right)>();
+            for (var i = 0; i < basis.Count; i++)
+                for (var j = i + 1; j < basis.Count; j++)
+                    pairs.Add((i, j));
+
+            var considered = 0;
+            while (pairs.Count > 0)
+            {
+                if (!budget.Spend("time")) return null;
+                if (!budget.Allow(++considered <= budget.MaxPairs, "pairs")) return null;
+                if (!budget.Allow(basis.Count <= budget.MaxBasisSize, "basis")) return null;
+
+                var chosen = ChooseNormalStrategy(basis, pairs, order);
+                var (left, right) = pairs[chosen];
+                pairs.RemoveAt(chosen);
+
+                var first = basis[left];
+                var second = basis[right];
+                var firstLeading = first.LeadingMonomial(order);
+                var secondLeading = second.LeadingMonomial(order);
+                var lcm = MultivariatePolynomial.MonomialLcm(firstLeading, secondLeading);
+
+                if (Coprime(firstLeading, secondLeading))
+                    continue;
+                if (IsRedundantByChain(basis, pairs, order, left, right, lcm))
+                    continue;
+
+                var s = SPolynomial(first, second, order, lcm);
+                if (s is null) { _ = budget.Allow(false, "degree"); return null; }
+                if (!budget.CheckPolynomial(s)) return null;
+
+                var remainder = TopReduce(s, basis, order, budget);
+                if (remainder is null) return null;
+                if (remainder.IsZero)
+                    continue;
+
+                basis.Add(remainder.MakeMonic(order));
+                for (var k = 0; k < basis.Count - 1; k++)
+                    pairs.Add((k, basis.Count - 1));
+            }
+
+            return Reduced(basis, order, budget);
+        }
+
+        static int ChooseNormalStrategy(
+            List<MultivariatePolynomial> basis, List<(int Left, int Right)> pairs, MonomialOrder order)
+        {
+            var best = 0;
+            var bestLcm = MultivariatePolynomial.MonomialLcm(
+                basis[pairs[0].Left].LeadingMonomial(order), basis[pairs[0].Right].LeadingMonomial(order));
+            for (var k = 1; k < pairs.Count; k++)
+            {
+                var lcm = MultivariatePolynomial.MonomialLcm(
+                    basis[pairs[k].Left].LeadingMonomial(order), basis[pairs[k].Right].LeadingMonomial(order));
+                var degree = MultivariatePolynomial.TotalDegree(lcm);
+                var bestDegree = MultivariatePolynomial.TotalDegree(bestLcm);
+                if (degree < bestDegree || (degree == bestDegree && MultivariatePolynomial.Greater(order, bestLcm, lcm)))
+                {
+                    best = k;
+                    bestLcm = lcm;
+                }
+            }
+            return best;
+        }
+
+        /// <summary>
+        /// Buchberger's first criterion: leading monomials sharing no variable give an
+        /// S-polynomial that always reduces to zero, so it need not be built.
+        /// </summary>
+        static bool Coprime(ulong left, ulong right)
+        {
+            for (var variable = 0; variable < MultivariatePolynomial.MaxVariables; variable++)
+                if (MultivariatePolynomial.PowerOfMonomial(left, variable) > 0
+                    && MultivariatePolynomial.PowerOfMonomial(right, variable) > 0)
+                    return false;
+            return true;
+        }
+
+        /// <summary>
+        /// The chain criterion: a third basis element whose leading monomial divides this
+        /// pair's lcm, and whose own two pairs are already dealt with, makes this one
+        /// redundant.
+        /// </summary>
+        static bool IsRedundantByChain(
+            List<MultivariatePolynomial> basis, List<(int Left, int Right)> pairs,
+            MonomialOrder order, int left, int right, ulong lcm)
+        {
+            for (var k = 0; k < basis.Count; k++)
+            {
+                if (k == left || k == right)
+                    continue;
+                if (!MultivariatePolynomial.MonomialDivides(basis[k].LeadingMonomial(order), lcm))
+                    continue;
+                var withLeft = (Math.Min(left, k), Math.Max(left, k));
+                var withRight = (Math.Min(right, k), Math.Max(right, k));
+                if (!pairs.Contains(withLeft) && !pairs.Contains(withRight))
+                    return true;
+            }
+            return false;
+        }
+
+        static MultivariatePolynomial? SPolynomial(
+            MultivariatePolynomial first, MultivariatePolynomial second, MonomialOrder order, ulong lcm)
+        {
+            var firstLeading = first.LeadingMonomial(order);
+            var secondLeading = second.LeadingMonomial(order);
+            var fromFirst = first.TimesTerm(
+                MultivariatePolynomial.MonomialQuotient(lcm, firstLeading),
+                ERational.One.Divide(first.LeadingCoefficient(order)));
+            var fromSecond = second.TimesTerm(
+                MultivariatePolynomial.MonomialQuotient(lcm, secondLeading),
+                ERational.One.Divide(second.LeadingCoefficient(order)));
+            if (fromFirst is null || fromSecond is null)
+                return null;
+            return fromFirst.Subtract(fromSecond);
+        }
+
+        /// <summary>
+        /// Reduces only while the leading monomial is divisible. Enough for Buchberger — a
+        /// remainder that is nonzero and no longer top-reducible is a correct thing to add —
+        /// and cheaper than reducing the tail nobody looks at. <see cref="FullyReduce"/> is
+        /// the one for callers who need every term standard.
+        /// </summary>
+        static MultivariatePolynomial? TopReduce(
+            MultivariatePolynomial polynomial, List<MultivariatePolynomial> basis,
+            MonomialOrder order, GroebnerBudget budget)
+        {
+            var remainder = polynomial;
+            var reducing = true;
+            while (reducing && !remainder.IsZero)
+            {
+                if (!budget.Spend("time")) return null;
+                if (!budget.CheckPolynomial(remainder)) return null;
+                reducing = false;
+                var leading = remainder.LeadingMonomial(order);
+                foreach (var divisor in basis)
+                {
+                    if (divisor.IsZero)
+                        continue;
+                    var divisorLeading = divisor.LeadingMonomial(order);
+                    if (!MultivariatePolynomial.MonomialDivides(divisorLeading, leading))
+                        continue;
+                    var scaled = divisor.TimesTerm(
+                        MultivariatePolynomial.MonomialQuotient(leading, divisorLeading),
+                        remainder.CoefficientOf(leading).Divide(divisor.LeadingCoefficient(order)));
+                    if (scaled is null) { _ = budget.Allow(false, "degree"); return null; }
+                    remainder = remainder.Subtract(scaled);
+                    reducing = true;
+                    break;
+                }
+            }
+            return remainder;
+        }
+
+        /// <summary>
+        /// The true normal form: every term divided out, not only the leading one. FGLM
+        /// reads a reduction as a vector over the standard monomials, so a tail left
+        /// unreduced is a term with no coordinate to sit in.
+        /// </summary>
+        internal static MultivariatePolynomial? FullyReduce(
+            MultivariatePolynomial polynomial, IReadOnlyList<MultivariatePolynomial> basis,
+            MonomialOrder order, GroebnerBudget budget)
+        {
+            var remainder = MultivariatePolynomial.Zero(polynomial.VariableCount);
+            var work = polynomial;
+            while (!work.IsZero)
+            {
+                if (!budget.Spend("time")) return null;
+                if (!budget.CheckPolynomial(work)) return null;
+
+                var leading = work.LeadingMonomial(order);
+                var coefficient = work.CoefficientOf(leading);
+                MultivariatePolynomial? divisor = null;
+                foreach (var candidate in basis)
+                    if (!candidate.IsZero
+                        && MultivariatePolynomial.MonomialDivides(candidate.LeadingMonomial(order), leading))
+                    {
+                        divisor = candidate;
+                        break;
+                    }
+
+                if (divisor is null)
+                {
+                    var standing = MultivariatePolynomial.Term(work.VariableCount, leading, coefficient);
+                    remainder = remainder.Add(standing);
+                    work = work.Subtract(standing);
+                    continue;
+                }
+
+                var scaled = divisor.TimesTerm(
+                    MultivariatePolynomial.MonomialQuotient(leading, divisor.LeadingMonomial(order)),
+                    coefficient.Divide(divisor.LeadingCoefficient(order)));
+                if (scaled is null) { _ = budget.Allow(false, "degree"); return null; }
+                work = work.Subtract(scaled);
+            }
+            return remainder;
+        }
+
+        /// <summary>Drops what other elements already cover, and reduces what is left.</summary>
+        static List<MultivariatePolynomial>? Reduced(
+            List<MultivariatePolynomial> basis, MonomialOrder order, GroebnerBudget budget)
+        {
+            var kept = new List<MultivariatePolynomial>();
+            for (var i = 0; i < basis.Count; i++)
+            {
+                var leading = basis[i].LeadingMonomial(order);
+                var covered = false;
+                for (var j = 0; j < basis.Count && !covered; j++)
+                {
+                    if (i == j)
+                        continue;
+                    var otherLeading = basis[j].LeadingMonomial(order);
+                    if (MultivariatePolynomial.MonomialDivides(otherLeading, leading)
+                        && (otherLeading != leading || j < i))
+                        covered = true;
+                }
+                if (!covered)
+                    kept.Add(basis[i]);
+            }
+
+            var result = new List<MultivariatePolynomial>(kept.Count);
+            for (var i = 0; i < kept.Count; i++)
+            {
+                var others = new List<MultivariatePolynomial>(kept.Count - 1);
+                for (var j = 0; j < kept.Count; j++)
+                    if (i != j)
+                        others.Add(kept[j]);
+                var reduced = FullyReduce(kept[i], others, order, budget);
+                if (reduced is null)
+                    return null;
+                if (!reduced.IsZero)
+                    result.Add(reduced.MakeMonic(order));
+            }
+            return result;
+        }
+    }
+}