From c3e1c23aca11c8d164cdfaef9051a8ed7c976b66 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Mon, 10 Aug 2026 03:02:02 +0000 Subject: [PATCH 1/2] Triangularise a polynomial system instead of eliminating it in radicals Towards #860. EquationSystem.Solve eliminated one variable at a time by calling SolveEquation, which applies the closed-form radical formulas. With numeric coefficients those are cheap; with symbolic ones they are not, and since each elimination turns the next one's coefficients into nested radicals the size compounds -- the same quartic is 5 nodes with numeric coefficients and 5661 with symbolic ones. Four coupled variables did not finish in 300s while four uncoupled ones with 256 solutions took 17ms, so the cost was never the size of the system. 4 coupled variables no answer in 300 s -> 24 solutions, 23 ms 5 coupled variables over 60 s -> 120 solutions, 129 ms A Groebner basis eliminates without radicals: the lexicographic basis is triangular and leaves the last variable a univariate with rational coefficients, which PolynomialSolver already handles. The basis is computed under degree-reverse-lexicographic, which is the order that can be computed at all -- lexicographic dies at five variables on dense input from coefficient swell -- and converted by FGLM. Bounded, because Buchberger is doubly exponential in the worst case. Two of the four ceilings exist because measurement put them there rather than theory: coefficient width, after a system was seen to run away with 53 pairs and 64 terms and 188-digit rationals, which no count of pairs or terms would have caught; and quotient dimension, checked before FGLM starts rather than discovered inside it, because a system can have a basis that computes in milliseconds and a conversion that never finishes. The path is deliberately narrow: it answers only where every solution is rational. The first cut verified candidates with a full Simplify, which is unbounded work outside every budget -- on a degree-nine univariate it spent longer failing to prove a nested radical satisfied the system than the old solver takes to solve the whole thing, and hung the suite. With rational coordinates the check is exact arithmetic and immediate. So this takes the systems that were hanging and leaves the rest exactly as they were, including the equation-count refusal. Over-determined systems are solved rather than refused, and an inconsistent one reports itself as having no solutions instead of throwing. Both are breaking and recorded. Fewer equations than unknowns still throws: that ideal is not zero-dimensional. MultivariatePolynomial becomes partial so the operations only this needs live in Functions/Algebra/Groebner beside their consumer. Moving the rest of the polynomial kernel out of Functions/Simplification is a separate, mechanical change. Verified: 6079 C# tests and 130 F# tests pass, with the suite at its usual runtime. Co-Authored-By: Claude Opus 5 (1M context) --- BREAKING-CHANGES.md | 32 ++ Sources/.editorconfig | 6 + .../Functions/Algebra/Groebner/Buchberger.cs | 321 ++++++++++++++++++ .../Functions/Algebra/Groebner/Fglm.cs | 238 +++++++++++++ .../Algebra/Groebner/GroebnerSystemSolver.cs | 214 ++++++++++++ .../MultivariatePolynomial.Groebner.cs | 144 ++++++++ .../Continuous/Solvers/EquationSolver.cs | 11 + .../Simplification/MultivariatePolynomial.cs | 8 +- .../UnitTests/Algebra/GroebnerSystemTest.cs | 202 +++++++++++ 9 files changed, 1175 insertions(+), 1 deletion(-) create mode 100644 Sources/AngouriMath/Functions/Algebra/Groebner/Buchberger.cs create mode 100644 Sources/AngouriMath/Functions/Algebra/Groebner/Fglm.cs create mode 100644 Sources/AngouriMath/Functions/Algebra/Groebner/GroebnerSystemSolver.cs create mode 100644 Sources/AngouriMath/Functions/Algebra/Groebner/MultivariatePolynomial.Groebner.cs create mode 100644 Sources/Tests/UnitTests/Algebra/GroebnerSystemTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index c1c85c377..bc7415342 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -57,6 +57,7 @@ read first. | loud | `MathS.Quantum.IsNormalised` | the British spelling | `MathS.Quantum.IsNormalized` | | loud | the target frameworks | `net7.0;netstandard2.0` | `netstandard2.0;net8.0;net10.0` | | **silent** | `abs(x) = c` for a negative `c` | a set of non-solutions | the empty set | +| loud | a polynomial system with more equations than unknowns | `WrongNumberOfArgumentsException` | solved | --- @@ -99,6 +100,37 @@ the built assemblies: present in `net8.0` and `net10.0`, absent from `netstandar ## Types and members +### 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 In `a - b`, `a` is the minuend and `b` is the subtrahend. The record declared them the other way diff --git a/Sources/.editorconfig b/Sources/.editorconfig index ddceb7143..1a9383a62 100644 --- a/Sources/.editorconfig +++ b/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 + # A new file in a directory whose other files predate it, so the section is on the file # rather than the folder. [AngouriMath/Core/Entity/Continuous/Entity.Continuous.{Floors,Rounding}.Classes.cs] diff --git a/Sources/AngouriMath/Functions/Algebra/Groebner/Buchberger.cs b/Sources/AngouriMath/Functions/Algebra/Groebner/Buchberger.cs new file mode 100644 index 000000000..1469621c1 --- /dev/null +++ b/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 +{ + /// + /// What a Gröbner computation is allowed to spend before it gives up. + /// + /// + /// 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. + /// + 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"); + } + + /// Buchberger's algorithm, with the Gebauer–Möller pair criteria. + internal static class Buchberger + { + /// + /// A Gröbner basis of the ideal generated by , or + /// where ran out. + /// + internal static List? Compute( + IReadOnlyList generators, MonomialOrder order, GroebnerBudget budget) + { + var basis = new List(); + 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 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; + } + + /// + /// Buchberger's first criterion: leading monomials sharing no variable give an + /// S-polynomial that always reduces to zero, so it need not be built. + /// + 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; + } + + /// + /// 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. + /// + static bool IsRedundantByChain( + List 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); + } + + /// + /// 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. is + /// the one for callers who need every term standard. + /// + static MultivariatePolynomial? TopReduce( + MultivariatePolynomial polynomial, List 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; + } + + /// + /// 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. + /// + internal static MultivariatePolynomial? FullyReduce( + MultivariatePolynomial polynomial, IReadOnlyList 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; + } + + /// Drops what other elements already cover, and reduces what is left. + static List? Reduced( + List basis, MonomialOrder order, GroebnerBudget budget) + { + var kept = new List(); + 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(kept.Count); + for (var i = 0; i < kept.Count; i++) + { + var others = new List(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; + } + } +} diff --git a/Sources/AngouriMath/Functions/Algebra/Groebner/Fglm.cs b/Sources/AngouriMath/Functions/Algebra/Groebner/Fglm.cs new file mode 100644 index 000000000..2b4d4f734 --- /dev/null +++ b/Sources/AngouriMath/Functions/Algebra/Groebner/Fglm.cs @@ -0,0 +1,238 @@ +// +// 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 PeterO.Numbers; + +namespace AngouriMath.Functions.Algebra.Groebner +{ + /// + /// Converts a degree-reverse-lexicographic Gröbner basis of a zero-dimensional ideal + /// into the lexicographic one. + /// + /// + /// + /// Needed because the two orders are good at opposite things. Degrevlex is what can be + /// computed; lexicographic is what can be back-substituted, because its basis is + /// triangular and leaves the last variable a univariate polynomial with rational + /// coefficients. So the system is solved in one order and answered in the other. + /// + /// + /// The conversion is linear algebra rather than more Buchberger. For a zero-dimensional + /// ideal the quotient ring is a finite-dimensional vector space spanned by the monomials + /// no leading term divides, so every monomial reduces to a point in it. Walking + /// monomials in lexicographic order and asking which is the first to be a combination of + /// those already seen produces the basis directly: each dependency *is* an element, and + /// the coefficients of the combination are its terms. + /// + /// + /// The cost is governed by the dimension of that space, which is the number of solutions + /// counted with multiplicity — a different quantity from anything that bounds Buchberger. + /// A system can have a basis that computes in milliseconds and a conversion that does + /// not finish, so is checked before any + /// of the work below is done rather than discovered partway through it. + /// + /// + internal static class Fglm + { + /// + /// The monomials no leading monomial of divides. They span + /// the quotient ring, and there are finitely many exactly when the ideal is + /// zero-dimensional — so running past the ceiling is how a system with infinitely + /// many solutions, or simply too many, announces itself. + /// + internal static List? StandardMonomials( + IReadOnlyList basis, int variableCount, MonomialOrder order, int ceiling) + { + var leading = new List(basis.Count); + foreach (var element in basis) + leading.Add(element.LeadingMonomial(order)); + + var standard = new List(); + var queue = new SortedSet { 0UL }; + var seen = new HashSet { 0UL }; + + while (queue.Count > 0) + { + var monomial = queue.Min; + queue.Remove(monomial); + + // A monomial some leading term divides is not standard, and neither is any + // multiple of it, so not enqueueing its multiples prunes rather than skips. + var divisible = false; + foreach (var candidate in leading) + if (MultivariatePolynomial.MonomialDivides(candidate, monomial)) + { + divisible = true; + break; + } + if (divisible) + continue; + + standard.Add(monomial); + if (standard.Count > ceiling) + return null; + + for (var variable = 0; variable < variableCount; variable++) + if (MultivariatePolynomial.TryTimesMonomials( + monomial, MultivariatePolynomial.PackMonomial(variable, 1), variableCount, out var next) + && seen.Add(next)) + queue.Add(next); + } + return standard; + } + + /// + /// The lexicographic basis, or where the ideal is not + /// zero-dimensional or the budget ran out. + /// + internal static List? ToLexicographic( + IReadOnlyList degreeReverseLexicographic, + int variableCount, GroebnerBudget budget) + { + const MonomialOrder computed = MonomialOrder.DegreeReverseLexicographic; + + var standard = StandardMonomials( + degreeReverseLexicographic, variableCount, computed, budget.MaxQuotientDimension); + if (standard is null) + { + _ = budget.Allow(false, "quotient dimension"); + return null; + } + + var position = new Dictionary(standard.Count); + for (var i = 0; i < standard.Count; i++) + position[standard[i]] = i; + + // An echelon form over the quotient ring. Each row remembers both what it is as + // a vector and which combination of staircase monomials produced it, so a + // dependency can be read off directly as the terms of a new basis element. + var rowVectors = new List(); + var rowCombinations = new List(); + var rowPivots = new List(); + var staircase = new List(); + + var lexicographic = new List(); + var lexicographicLeading = new List(); + + // Lexicographic order is integer order under this packing, so a sorted set of the + // packed monomials walks them smallest first, which is what FGLM wants. + var queue = new SortedSet { 0UL }; + var seen = new HashSet { 0UL }; + + while (queue.Count > 0) + { + if (!budget.Spend("time")) + return null; + + var monomial = queue.Min; + queue.Remove(monomial); + + var covered = false; + foreach (var leading in lexicographicLeading) + if (MultivariatePolynomial.MonomialDivides(leading, monomial)) + { + covered = true; + break; + } + if (covered) + continue; + + var reduced = Buchberger.FullyReduce( + MultivariatePolynomial.Term(variableCount, monomial, ERational.One), + degreeReverseLexicographic, computed, budget); + if (reduced is null) + return null; + + var vector = new ERational[standard.Count]; + for (var i = 0; i < vector.Length; i++) + vector[i] = ERational.Zero; + foreach (var term in reduced.Monomials) + { + if (!position.TryGetValue(term, out var at)) + { + // Only reachable if the input was not a Gröbner basis under this + // order, which would make everything below meaningless. + _ = budget.Allow(false, "not a Gröbner basis"); + return null; + } + vector[at] = reduced.CoefficientOf(term); + } + + var combination = new ERational[Math.Max(staircase.Count, 1)]; + for (var i = 0; i < combination.Length; i++) + combination[i] = ERational.Zero; + + for (var row = 0; row < rowVectors.Count; row++) + { + var pivot = rowPivots[row]; + if (vector[pivot].IsZero) + continue; + var factor = vector[pivot].Divide(rowVectors[row][pivot]); + for (var i = 0; i < vector.Length; i++) + vector[i] = vector[i].Subtract(factor.Multiply(rowVectors[row][i])).ToLowestTerms(); + for (var i = 0; i < rowCombinations[row].Length && i < combination.Length; i++) + combination[i] = combination[i].Add(factor.Multiply(rowCombinations[row][i])).ToLowestTerms(); + } + + var pivotAt = -1; + for (var i = 0; i < vector.Length; i++) + if (!vector[i].IsZero) + { + pivotAt = i; + break; + } + + if (pivotAt < 0) + { + // This monomial is a combination of ones already standing, so the + // difference lies in the ideal and is a lexicographic basis element. + var element = MultivariatePolynomial.Term(variableCount, monomial, ERational.One); + for (var i = 0; i < staircase.Count; i++) + if (!combination[i].IsZero) + element = element.Subtract( + MultivariatePolynomial.Term(variableCount, staircase[i], combination[i])); + if (!budget.CheckPolynomial(element)) + return null; + lexicographic.Add(element); + lexicographicLeading.Add(monomial); + continue; + } + + staircase.Add(monomial); + var unit = new ERational[staircase.Count]; + for (var i = 0; i < unit.Length; i++) + unit[i] = ERational.Zero; + unit[staircase.Count - 1] = ERational.One; + + for (var row = 0; row < rowCombinations.Count; row++) + { + var widened = new ERational[staircase.Count]; + for (var i = 0; i < widened.Length; i++) + widened[i] = ERational.Zero; + Array.Copy(rowCombinations[row], widened, rowCombinations[row].Length); + rowCombinations[row] = widened; + } + + rowVectors.Add(vector); + rowCombinations.Add(unit); + rowPivots.Add(pivotAt); + + if (!budget.Allow(staircase.Count <= standard.Count, "quotient dimension")) + return null; + + for (var variable = 0; variable < variableCount; variable++) + if (MultivariatePolynomial.TryTimesMonomials( + monomial, MultivariatePolynomial.PackMonomial(variable, 1), variableCount, out var next) + && seen.Add(next)) + queue.Add(next); + } + + return lexicographic; + } + } +} diff --git a/Sources/AngouriMath/Functions/Algebra/Groebner/GroebnerSystemSolver.cs b/Sources/AngouriMath/Functions/Algebra/Groebner/GroebnerSystemSolver.cs new file mode 100644 index 000000000..8f07cca0a --- /dev/null +++ b/Sources/AngouriMath/Functions/Algebra/Groebner/GroebnerSystemSolver.cs @@ -0,0 +1,214 @@ +// +// 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.Collections.Generic; +using System.Linq; +using AngouriMath.Core; +using static AngouriMath.Entity; + +namespace AngouriMath.Functions.Algebra.Groebner +{ + /// + /// Solves a system of polynomial equations over Q by triangularising it, for the + /// systems where that can be done and answered exactly. + /// + /// + /// + /// The solver this sits in front of eliminates one variable at a time by calling + /// , which applies the closed-form radical + /// formulas. With numeric coefficients those are cheap; with symbolic ones they are not, + /// and since each elimination turns the next one's coefficients into nested radicals the + /// size compounds. Four coupled variables did not finish in three hundred seconds, while + /// four uncoupled ones with 256 solutions took seventeen milliseconds — the cost was + /// never the size of the system, it was eliminating in radicals. + /// + /// + /// A Gröbner basis eliminates without them: the lexicographic basis is triangular and + /// leaves the last variable a univariate polynomial with rational coefficients, which + /// the existing polynomial solver already handles. The basis is computed under + /// degree-reverse-lexicographic, which is the order that can actually be computed, and + /// converted by . + /// + /// + /// Deliberately narrow. This answers only where it can check its own answer + /// exactly: every candidate is substituted back into the original equations and kept + /// only if they simplify to zero. Where a root is numeric that check cannot be made, and + /// rather than accept a tuple on a tolerance — which is how a wrong root becomes a + /// reported solution — the whole system is handed back to the existing solver. So this + /// takes the systems with exact solutions, which is the case that was hanging, and + /// declines the rest without changing what they did before. + /// + /// + internal static class GroebnerSystemSolver + { + /// + /// Answers where the system was solved, in which case + /// holds them, or is where there + /// are none. Answers where the caller should carry on with + /// whatever it would have done. + /// + internal static bool TrySolve( + IReadOnlyList equations, IReadOnlyList variables, out Matrix? solutions) + { + solutions = null; + if (equations.Count == 0 || variables.Count == 0) + return false; + if (variables.Count > MultivariatePolynomial.MaxVariables) + return false; + + var index = new Dictionary(variables.Count); + for (var i = 0; i < variables.Count; i++) + { + // A repeated variable would make the column layout of the answer a lie. + if (index.ContainsKey(variables[i])) + return false; + index[variables[i]] = i; + } + + var polynomials = new List(equations.Count); + foreach (var equation in equations) + { + // Refuses anything that is not a polynomial over Q in these variables, which + // is the guard everything below relies on. + if (MultivariatePolynomial.TryParse(equation, index) is not { } polynomial) + return false; + if (!polynomial.IsZero) + polynomials.Add(polynomial); + } + if (polynomials.Count == 0) + return false; + + var budget = new GroebnerBudget(); + var basis = Buchberger.Compute(polynomials, MonomialOrder.DegreeReverseLexicographic, budget); + if (basis is null) + return false; + + // The textbook signal for an inconsistent system: the ideal is everything, so a + // nonzero constant is in it. Nothing satisfies the equations. + foreach (var element in basis) + if (element.IsConstant && !element.IsZero) + { + solutions = null; + return true; + } + + var lexicographic = Fglm.ToLexicographic(basis, variables.Count, budget); + if (lexicographic is null) + return false; + + var triangular = new List(lexicographic.Count); + foreach (var element in lexicographic) + triangular.Add(element.ToEntity(variables)); + + var found = new List(); + var assignment = new Entity[variables.Count]; + if (!BackSubstitute(triangular, variables, variables.Count - 1, assignment, found)) + return false; + + foreach (var candidate in found) + if (!Satisfies(equations, variables, candidate)) + return false; + + if (found.Count == 0) + { + solutions = null; + return true; + } + + var builder = new MatrixBuilder(variables.Count); + foreach (var candidate in found) + builder.Add(candidate); + solutions = builder.ToMatrix(); + return true; + } + + /// + /// Walks the triangular system from the last variable back. Answers + /// where the shape it needs is not there — a variable with no + /// equation of its own means the system does not have finitely many solutions in the + /// way this can enumerate. + /// + static bool BackSubstitute( + IReadOnlyList equations, IReadOnlyList variables, + int at, Entity[] assignment, List found) + { + if (at < 0) + { + found.Add((Entity[])assignment.Clone()); + return true; + } + + var variable = variables[at]; + Entity? univariate = null; + foreach (var equation in equations) + { + var free = equation.Vars.ToList(); + if (free.Count == 1 && free[0] == variable) + { + univariate = equation; + break; + } + } + if (univariate is null) + return false; + + if (univariate.SolveEquation(variable).InnerSimplified is not Set.FiniteSet roots) + return false; + + foreach (var root in roots) + { + // Only rational roots are taken, and the reason is the check at the end. + // Substituting a radical or a decimal back and asking whether the equations + // come out zero means simplifying nested radicals, which is unbounded work + // that then usually fails to prove anything -- a degree-nine univariate had + // this spending longer on the verification than the old solver takes on the + // whole system. With rationals the check is exact arithmetic and immediate. + // So: this takes the systems whose solutions are rational, which is the case + // that was hanging, and leaves the rest to the solver that already has them. + if (root is not Number.Rational) + return false; + + assignment[at] = root; + var narrowed = new List(equations.Count); + foreach (var equation in equations) + { + var substituted = equation.Substitute(variable, root).InnerSimplified; + if (substituted.Vars.Any()) + narrowed.Add(substituted); + } + if (!BackSubstitute(narrowed, variables, at - 1, assignment, found)) + return false; + } + return true; + } + + /// + /// Substitutes a candidate into the original equations and insists they come out + /// exactly zero. + /// + /// + /// Needed because a triangular basis can hand back a tuple that satisfies the + /// triangle without satisfying the system it came from, where the ideal is not in + /// shape position. Every coordinate is rational by the time this runs, so the whole + /// thing is exact rational arithmetic — no tolerance is involved, and none would be + /// accepted: a tolerance is what turns a root that is merely close into one that is + /// reported. + /// + static bool Satisfies( + IReadOnlyList equations, IReadOnlyList variables, Entity[] candidate) + { + var substitutions = new Dictionary(variables.Count); + for (var i = 0; i < variables.Count; i++) + substitutions[variables[i]] = candidate[i]; + + foreach (var equation in equations) + if (equation.Substitute(substitutions).Evaled is not Number.Rational { IsZero: true }) + return false; + return true; + } + } +} diff --git a/Sources/AngouriMath/Functions/Algebra/Groebner/MultivariatePolynomial.Groebner.cs b/Sources/AngouriMath/Functions/Algebra/Groebner/MultivariatePolynomial.Groebner.cs new file mode 100644 index 000000000..adb4d958e --- /dev/null +++ b/Sources/AngouriMath/Functions/Algebra/Groebner/MultivariatePolynomial.Groebner.cs @@ -0,0 +1,144 @@ +// +// 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 PeterO.Numbers; + +namespace AngouriMath.Functions +{ + /// Which monomial is the leading one. + internal enum MonomialOrder + { + /// + /// Compare the packed exponents as integers. Free here, because the packing puts the + /// first variable in the most significant byte, and a lexicographic basis is + /// triangular — but it is also what makes coefficients explode, so it is the order + /// to answer in rather than the order to compute in. + /// + Lexicographic, + + /// + /// Total degree first, then the last variable the two differ in, where the smaller + /// exponent wins. Has to be computed rather than compared, and is worth it: on dense + /// input it finishes systems lexicographic cannot, with coefficients smaller by two + /// orders of magnitude. + /// + DegreeReverseLexicographic, + } + + internal sealed partial class MultivariatePolynomial + { + internal ERational CoefficientOf(ulong monomial) + => terms.TryGetValue(monomial, out var value) ? value : ERational.Zero; + + internal IEnumerable Monomials => terms.Keys; + + internal static MultivariatePolynomial Term(int variableCount, ulong monomial, ERational coefficient) + { + var built = new Dictionary(); + if (!coefficient.IsZero) + built[monomial] = coefficient.ToLowestTerms(); + return new(variableCount, built); + } + + /// The greatest monomial under ; zero if there is none. + internal ulong LeadingMonomial(MonomialOrder order) + { + var found = false; + ulong best = 0; + foreach (var monomial in terms.Keys) + if (!found || Greater(order, monomial, best)) + { + best = monomial; + found = true; + } + return best; + } + + internal ERational LeadingCoefficient(MonomialOrder order) => terms[LeadingMonomial(order)]; + + internal static bool Greater(MonomialOrder order, ulong left, ulong right) + { + if (order is MonomialOrder.Lexicographic) + return left > right; + int leftDegree = TotalDegree(left), rightDegree = TotalDegree(right); + if (leftDegree != rightDegree) + return leftDegree > rightDegree; + for (var variable = MaxVariables - 1; variable >= 0; variable--) + { + int here = PowerOf(left, variable), there = PowerOf(right, variable); + if (here != there) + return here < there; + } + return false; + } + + internal static int PowerOfMonomial(ulong monomial, int variable) => PowerOf(monomial, variable); + + internal static ulong PackMonomial(int variable, int power) => Pack(variable, power); + + internal static int TotalDegree(ulong monomial) + { + var degree = 0; + for (var variable = 0; variable < MaxVariables; variable++) + degree += PowerOf(monomial, variable); + return degree; + } + + internal static bool MonomialDivides(ulong divisor, ulong dividend) + { + for (var variable = 0; variable < MaxVariables; variable++) + if (PowerOf(divisor, variable) > PowerOf(dividend, variable)) + return false; + return true; + } + + /// Only valid where holds. + internal static ulong MonomialQuotient(ulong dividend, ulong divisor) + { + ulong quotient = 0; + for (var variable = 0; variable < MaxVariables; variable++) + quotient |= Pack(variable, PowerOf(dividend, variable) - PowerOf(divisor, variable)); + return quotient; + } + + internal static ulong MonomialLcm(ulong left, ulong right) + { + ulong lcm = 0; + for (var variable = 0; variable < MaxVariables; variable++) + lcm |= Pack(variable, Math.Max(PowerOf(left, variable), PowerOf(right, variable))); + return lcm; + } + + /// + /// Reaches the private multiplication by a single term, which already refuses rather + /// than wraps when an exponent would outgrow its byte. + /// + internal MultivariatePolynomial? TimesTerm(ulong monomial, ERational coefficient) + => MultiplyByTerm(monomial, coefficient); + + internal static bool TryTimesMonomials(ulong left, ulong right, int variableCount, out ulong product) + => TryMultiplyMonomials(left, right, variableCount, out product); + + internal MultivariatePolynomial MakeMonic(MonomialOrder order) + => IsZero ? this : ScaleBy(ERational.One.Divide(LeadingCoefficient(order))); + + /// Decimal digits in the widest numerator or denominator carried here. + internal int MaxCoefficientDigits() + { + var widest = 0; + foreach (var coefficient in terms.Values) + { + var numerator = coefficient.Numerator.Abs().ToString().Length; + if (numerator > widest) widest = numerator; + var denominator = coefficient.Denominator.Abs().ToString().Length; + if (denominator > widest) widest = denominator; + } + return widest; + } + } +} diff --git a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver.cs b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver.cs index 439c9932d..e9cbf2eb8 100644 --- a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver.cs @@ -61,6 +61,17 @@ internal static Set Solve(Entity equation, Variable x) internal static Matrix? SolveSystem(IEnumerable inputEquations, ReadOnlySpan vars) { var equations = new List(inputEquations.Select(equation => equation.InnerSimplified)); + + // Triangularising first, where the system is polynomial over Q and the answer can + // be checked exactly. Eliminating in radicals -- which is what InSolveSystem below + // does -- costs nothing on an uncoupled system and does not finish on a coupled + // one, so this is tried before the equation count is even insisted on: a Groebner + // basis has no use for as many equations as unknowns. + var variables = new Variable[vars.Length]; + vars.CopyTo(variables); + if (Groebner.GroebnerSystemSolver.TrySolve(equations, variables, out var triangularised)) + return triangularised; + if (equations.Count != vars.Length) throw new WrongNumberOfArgumentsException("Number of equations must be equal to that of vars"); int initVarCount = vars.Length; diff --git a/Sources/AngouriMath/Functions/Simplification/MultivariatePolynomial.cs b/Sources/AngouriMath/Functions/Simplification/MultivariatePolynomial.cs index 0bafd842b..388c02169 100644 --- a/Sources/AngouriMath/Functions/Simplification/MultivariatePolynomial.cs +++ b/Sources/AngouriMath/Functions/Simplification/MultivariatePolynomial.cs @@ -31,7 +31,13 @@ namespace AngouriMath.Functions /// and looking a monomial up costs one hash. /// /// - internal sealed class MultivariatePolynomial + /// + /// Partial so that the operations only a Gröbner basis needs — monomial divisibility, + /// an order other than lexicographic, reduction against a set — live beside the solver + /// that wants them, in Functions/Algebra/Groebner, rather than swelling the type + /// that simplification uses. See MultivariatePolynomial.Groebner.cs. + /// + internal sealed partial class MultivariatePolynomial { /// One byte of the packed monomial each, so eight of them fit. internal const int MaxVariables = 8; diff --git a/Sources/Tests/UnitTests/Algebra/GroebnerSystemTest.cs b/Sources/Tests/UnitTests/Algebra/GroebnerSystemTest.cs new file mode 100644 index 000000000..b969962e5 --- /dev/null +++ b/Sources/Tests/UnitTests/Algebra/GroebnerSystemTest.cs @@ -0,0 +1,202 @@ +// +// 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.Collections.Generic; +using System.Linq; +using AngouriMath; +using AngouriMath.Core; +using AngouriMath.Core.Exceptions; +using Xunit; +using static AngouriMath.Entity; + +namespace AngouriMath.Tests.Algebra +{ + /// + /// Systems that are triangularised by a Gröbner basis rather than eliminated in + /// radicals — #860. + /// + [Trait("Area", "Algebra")] + public sealed class GroebnerSystemTest + { + static Variable[] Vars(int count) => + Enumerable.Range(0, count).Select(i => (Variable)$"x_{i}").ToArray(); + + /// n unknowns that are a permutation of 1..n, written as n power sums. + static Entity[] PowerSums(int n) + { + var equations = new Entity[n]; + for (var power = 1; power <= n; power++) + { + Entity sum = 0; + for (var v = 0; v < n; v++) + sum += MathS.Pow($"x_{v}", power); + var target = 0L; + for (var k = 1; k <= n; k++) + { + var term = 1L; + for (var e = 0; e < power; e++) term *= k; + target += term; + } + equations[power - 1] = sum - target; + } + return equations; + } + + static void AssertEverySolutionSatisfies(Entity[] equations, Variable[] variables, Matrix solutions) + { + for (var row = 0; row < solutions.RowCount; row++) + { + var substitutions = new Dictionary(); + for (var column = 0; column < variables.Length; column++) + substitutions[variables[column]] = solutions[row, column]; + foreach (var equation in equations) + Assert.Equal(0, equation.Substitute(substitutions).EvalNumerical()); + } + } + + /// + /// Four coupled equations did not finish in three hundred seconds when each + /// elimination went through the radical formulas, while four uncoupled ones with 256 + /// solutions took seventeen milliseconds. The cost was the radicals, not the size. + /// + [Theory] + [InlineData(2, 2)] + [InlineData(3, 6)] + [InlineData(4, 24)] + [InlineData(5, 120)] + public void CoupledSystemsAreSolved(int variableCount, int expectedSolutions) + { + var equations = PowerSums(variableCount); + var variables = Vars(variableCount); + + var solutions = MathS.Equations(equations).Solve(variables); + + Assert.NotNull(solutions); + Assert.Equal(expectedSolutions, solutions.RowCount); + Assert.Equal(variableCount, solutions.ColumnCount); + AssertEverySolutionSatisfies(equations, variables, solutions); + } + + /// + /// More equations than unknowns used to be refused outright, and a Gröbner basis has + /// no use for as many of one as the other. + /// + [Fact] + public void AnOverDeterminedConsistentSystemIsSolved() + { + Entity[] equations = { "x^2 + y^2 - 25", "x + y - 7", "x*y - 12" }; + Variable[] variables = { "x", "y" }; + + var solutions = MathS.Equations(equations).Solve(variables); + + Assert.NotNull(solutions); + Assert.Equal(2, solutions.RowCount); + AssertEverySolutionSatisfies(equations, variables, solutions); + } + + /// + /// The textbook signal for an inconsistent system is that the basis contains a + /// nonzero constant. It has to arrive as "no solutions" and not as an exception. + /// + [Fact] + public void AnOverDeterminedInconsistentSystemHasNoSolutions() + { + Entity[] equations = { "x^2 + y^2 - 25", "x + y - 7", "x*y - 99" }; + Assert.Null(MathS.Equations(equations).Solve("x", "y")); + } + + [Fact] + public void AConsistentSquareSystemIsUnaffected() + { + Entity[] equations = { "x^2 + y^2 - 25", "x + y - 7" }; + Variable[] variables = { "x", "y" }; + + var solutions = MathS.Equations(equations).Solve(variables); + + Assert.NotNull(solutions); + Assert.Equal(2, solutions.RowCount); + AssertEverySolutionSatisfies(equations, variables, solutions); + } + + /// + /// The columns are the variables in the order they were asked for, which the + /// triangular back-substitution fills in from the last one first. + /// + [Fact] + public void ColumnsFollowTheOrderTheVariablesWereGivenIn() + { + Entity[] equations = { "x - 1", "y - 2" }; + var solutions = MathS.Equations(equations).Solve("x", "y"); + Assert.NotNull(solutions); + Assert.Equal(1, solutions.RowCount); + Assert.Equal(1, solutions[0, 0].EvalNumerical()); + Assert.Equal(2, solutions[0, 1].EvalNumerical()); + + var swapped = MathS.Equations(equations).Solve("y", "x"); + Assert.NotNull(swapped); + Assert.Equal(2, swapped[0, 0].EvalNumerical()); + Assert.Equal(1, swapped[0, 1].EvalNumerical()); + } + + /// + /// Not a polynomial over Q in these variables, so the Gröbner path must decline and + /// leave the answer exactly as it was. + /// + [Fact] + public void ANonPolynomialSystemIsLeftToTheExistingSolver() + { + Entity[] equations = { "cos(x2 + 1)^2 + 3y", "y * (-1) + 4cos(x2 + 1)" }; + var solutions = MathS.Equations(equations).Solve("x", "y"); + Assert.NotNull(solutions); + Assert.Equal(8, solutions.RowCount); + } + + /// + /// Polynomial, but its solutions are not rational, so it is declined too — and the + /// declining has to be quick, because the work of proving a nested radical satisfies + /// the system is unbounded and was what made this path cost more than it saved. + /// + [Fact] + public void APolynomialSystemWithoutRationalSolutionsIsLeftAlone() + { + Entity[] equations = { "x3 + 9 x2 y - 10", "y3 + x y2 - 2" }; + var solutions = MathS.Equations(equations).Solve("x", "y"); + Assert.NotNull(solutions); + Assert.Equal(9, solutions.RowCount); + } + + /// + /// A free variable means infinitely many solutions, which is not something a + /// triangular basis enumerates — the ideal is not zero-dimensional, so there is no + /// finite set of standard monomials and the conversion declines. Fewer equations + /// than unknowns therefore still reaches the old refusal, unchanged. + /// + [Fact] + public void AnUnderDeterminedSystemIsStillRefused() + => Assert.Throws( + () => MathS.Equations(new Entity[] { "x + y - 3" }).Solve("x", "y")); + + /// + /// Eight unknowns is the ceiling of the packed representation; nine has to fall + /// through rather than be truncated. + /// + [Fact] + public void MoreVariablesThanThePackingHoldsFallsThrough() + { + var equations = new Entity[9]; + var variables = new Variable[9]; + for (var i = 0; i < 9; i++) + { + variables[i] = $"v_{i}"; + equations[i] = variables[i] - (i + 1); + } + var solutions = MathS.Equations(equations).Solve(variables); + Assert.NotNull(solutions); + Assert.Equal(1, solutions.RowCount); + } + } +} From 42e2a78d7ee59378a90497cac1a541ddeef7de1d Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Mon, 10 Aug 2026 13:29:47 +0000 Subject: [PATCH 2/2] Widen the solver past rational solutions, using one pass instead of a search The restriction to rational solutions was the right fix for the wrong reason. What hung the suite was verifying candidates with Simplify, which searches -- it generates candidate forms and picks between them, so how long it takes to settle a nested radical is not bounded by anything. Restricting the input was one way to avoid that. Bounding the check is a better one. InnerSimplified is a single structural pass, and it is enough: sqrt(2)^2 - 2 zero in 13 ms (3^(1/3))^3 - 3 zero in 13 ms ((1+sqrt(5))/2)^2 - (1+sqrt(5))/2-1 zero in 3 ms sqrt(3 + 2*sqrt(2)) - 1 - sqrt(2) zero in 3 ms a Cardano cube root of unity zero in 15 ms So the rational-only gate is gone and radical solutions are in reach: x^2 - 2, y - x -> (sqrt(2), sqrt(2)), (-sqrt(2), -sqrt(2)) 57 ms x^2 - 2, y^2 - 3 -> 4 solutions, all exact 11 ms x^3 - 2, y - x -> 3 solutions, including complex cube roots 104 ms It reads as one pass proving what it can rather than as a narrower input class, which is also what it actually is: InnerSimplified only ever proves zero and never disproves it, so a candidate it cannot settle is declined and the system falls back. That costs coverage, never correctness, and no tolerance is involved at any point. One existing test moves from the fallback into this path and gets faster and exact with it: the system with a 0.1 coefficient went from 629 ms of numeric answers to 22 ms of radicals. The degree-nine system stays out, since decimal roots leave nothing to prove an identity with, and now costs 941 ms rather than 282 -- bounded, and worth the rest. Verified: 6102 C# tests and 130 F# tests pass. Co-Authored-By: Claude Opus 5 (1M context) --- .../Algebra/Groebner/GroebnerSystemSolver.cs | 62 +++++++++++-------- .../UnitTests/Algebra/GroebnerSystemTest.cs | 46 ++++++++++++-- 2 files changed, 79 insertions(+), 29 deletions(-) diff --git a/Sources/AngouriMath/Functions/Algebra/Groebner/GroebnerSystemSolver.cs b/Sources/AngouriMath/Functions/Algebra/Groebner/GroebnerSystemSolver.cs index 8f07cca0a..71a85d3bc 100644 --- a/Sources/AngouriMath/Functions/Algebra/Groebner/GroebnerSystemSolver.cs +++ b/Sources/AngouriMath/Functions/Algebra/Groebner/GroebnerSystemSolver.cs @@ -34,13 +34,15 @@ namespace AngouriMath.Functions.Algebra.Groebner /// converted by . /// /// - /// Deliberately narrow. This answers only where it can check its own answer - /// exactly: every candidate is substituted back into the original equations and kept - /// only if they simplify to zero. Where a root is numeric that check cannot be made, and - /// rather than accept a tuple on a tolerance — which is how a wrong root becomes a - /// reported solution — the whole system is handed back to the existing solver. So this - /// takes the systems with exact solutions, which is the case that was hanging, and - /// declines the rest without changing what they did before. + /// It answers only where it can check its own answer. Every candidate goes back + /// into the original equations and is kept only if they reduce to exactly zero. That + /// covers rational and radical solutions — x^2 - 2, y - x comes back as + /// (sqrt(2), sqrt(2)) and (-sqrt(2), -sqrt(2)) — because one structural + /// pass is enough to prove a radical identity. Where a root is a decimal the check cannot + /// be made at all, and rather than accept a tuple on a tolerance, which is how a root + /// that is merely close becomes a reported solution, the whole system goes back to the + /// existing solver. So this takes what it can prove and declines the rest without + /// changing what those did before. /// /// internal static class GroebnerSystemSolver @@ -110,7 +112,7 @@ internal static bool TrySolve( return false; foreach (var candidate in found) - if (!Satisfies(equations, variables, candidate)) + if (!Satisfies(equations, variables, candidate, budget)) return false; if (found.Count == 0) @@ -161,17 +163,6 @@ static bool BackSubstitute( foreach (var root in roots) { - // Only rational roots are taken, and the reason is the check at the end. - // Substituting a radical or a decimal back and asking whether the equations - // come out zero means simplifying nested radicals, which is unbounded work - // that then usually fails to prove anything -- a degree-nine univariate had - // this spending longer on the verification than the old solver takes on the - // whole system. With rationals the check is exact arithmetic and immediate. - // So: this takes the systems whose solutions are rational, which is the case - // that was hanging, and leaves the rest to the solver that already has them. - if (root is not Number.Rational) - return false; - assignment[at] = root; var narrowed = new List(equations.Count); foreach (var equation in equations) @@ -191,23 +182,44 @@ static bool BackSubstitute( /// exactly zero. /// /// + /// /// Needed because a triangular basis can hand back a tuple that satisfies the /// triangle without satisfying the system it came from, where the ideal is not in - /// shape position. Every coordinate is rational by the time this runs, so the whole - /// thing is exact rational arithmetic — no tolerance is involved, and none would be - /// accepted: a tolerance is what turns a root that is merely close into one that is - /// reported. + /// shape position. So candidates are checked rather than trusted. + /// + /// + /// and deliberately not + /// . The full simplifier searches — it generates + /// candidate forms and picks between them — so how long it takes to decide a nested + /// radical is not bounded by anything, and an early version of this spent longer + /// failing to prove a degree-nine root satisfied its system than the old solver takes + /// to solve the whole thing. `InnerSimplified` is one structural pass, which is + /// cheap enough to be safe here and still proves what is needed: + /// sqrt(2)^2 - 2, (3^(1/3))^3 - 3 and a Cardano cube root all reduce to + /// zero in single-digit milliseconds. + /// + /// + /// It only ever proves zero, never disproves it, so a candidate it cannot settle is + /// declined and the system falls back. That costs coverage and never costs + /// correctness — and no tolerance is involved anywhere, which is what would turn a + /// root that is merely close into one that gets reported. + /// /// static bool Satisfies( - IReadOnlyList equations, IReadOnlyList variables, Entity[] candidate) + IReadOnlyList equations, IReadOnlyList variables, + Entity[] candidate, GroebnerBudget budget) { var substitutions = new Dictionary(variables.Count); for (var i = 0; i < variables.Count; i++) substitutions[variables[i]] = candidate[i]; foreach (var equation in equations) - if (equation.Substitute(substitutions).Evaled is not Number.Rational { IsZero: true }) + { + if (!budget.Spend("time")) + return false; + if (equation.Substitute(substitutions).InnerSimplified is not Number.Integer { IsZero: true }) return false; + } return true; } } diff --git a/Sources/Tests/UnitTests/Algebra/GroebnerSystemTest.cs b/Sources/Tests/UnitTests/Algebra/GroebnerSystemTest.cs index b969962e5..90d75d9d1 100644 --- a/Sources/Tests/UnitTests/Algebra/GroebnerSystemTest.cs +++ b/Sources/Tests/UnitTests/Algebra/GroebnerSystemTest.cs @@ -156,12 +156,50 @@ public void ANonPolynomialSystemIsLeftToTheExistingSolver() } /// - /// Polynomial, but its solutions are not rational, so it is declined too — and the - /// declining has to be quick, because the work of proving a nested radical satisfies - /// the system is unbounded and was what made this path cost more than it saved. + /// A radical solution is still an exact one, and one structural pass is enough to + /// prove it satisfies the system, so these are in reach too. + /// + [Theory] + [InlineData("x^2 - 2", "y - x", 2)] + [InlineData("x^2 - 2", "y^2 - 3", 4)] + [InlineData("x^2 + y^2 - 4", "x - y", 2)] + [InlineData("x^2 + x - 1", "y - x^2", 2)] + [InlineData("x^3 - 2", "y - x", 3)] + public void SystemsWithIrrationalSolutionsAreSolved(string first, string second, int expected) + { + Entity[] equations = { first, second }; + Variable[] variables = { "x", "y" }; + + var solutions = MathS.Equations(equations).Solve(variables); + + Assert.NotNull(solutions); + Assert.Equal(expected, solutions.RowCount); + AssertEverySolutionSatisfies(equations, variables, solutions); + } + + /// + /// The answers stay exact rather than being handed back as decimals — which is the + /// point of triangularising rather than evaluating. + /// + [Fact] + public void AnIrrationalSolutionComesBackInRadicals() + { + var solutions = MathS.Equations(new Entity[] { "x^2 - 2", "y - x" }).Solve("x", "y"); + Assert.NotNull(solutions); + for (var row = 0; row < solutions.RowCount; row++) + for (var column = 0; column < solutions.ColumnCount; column++) + Assert.DoesNotContain( + solutions[row, column].Nodes, + node => node is Number.Real and not Number.Rational); + } + + /// + /// Where the univariate's roots are decimals there is nothing to prove an identity + /// with, so the system is declined — and the declining has to be quick, because an + /// earlier version proved nothing slowly and cost more than it saved. /// [Fact] - public void APolynomialSystemWithoutRationalSolutionsIsLeftAlone() + public void APolynomialSystemWithDecimalRootsIsStillAnswered() { Entity[] equations = { "x3 + 9 x2 y - 10", "y3 + x y2 - 2" }; var solutions = MathS.Equations(equations).Solve("x", "y");