diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 0e3e67374..6cfc0e643 100644 --- a/BREAKING-CHANGES.md +++ b/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 diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 860684a90..471fe9c9c 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 + [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 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..71a85d3bc --- /dev/null +++ b/Sources/AngouriMath/Functions/Algebra/Groebner/GroebnerSystemSolver.cs @@ -0,0 +1,226 @@ +// +// 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 . + /// + /// + /// 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 + { + /// + /// 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, budget)) + 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) + { + 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. 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, 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 (!budget.Spend("time")) + return false; + if (equation.Substitute(substitutions).InnerSimplified is not Number.Integer { 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..90d75d9d1 --- /dev/null +++ b/Sources/Tests/UnitTests/Algebra/GroebnerSystemTest.cs @@ -0,0 +1,240 @@ +// +// 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); + } + + /// + /// 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 APolynomialSystemWithDecimalRootsIsStillAnswered() + { + 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); + } + } +}