From 1887de5762fa6cd4344f6371db596b3f16944100 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Thu, 13 Aug 2026 19:11:58 +0000 Subject: [PATCH] Bound a gcd's intermediates separately from its inputs (#920) PolynomialGcd.Gcd declined on inputs well inside every documented limit, because an intermediate of the subresultant remainder sequence went past MultivariatePolynomial.MaxTerms while neither input nor answer came close: left (b + c + 1) * (a + b) * (a + b + c + d) 19 terms right (a^2 + b*c + d) * (a + b + c + d) * (a + b + c + d) 29 terms gcd a + b + c + d 4 terms was declined A multivariate pseudo-remainder multiplies through by a leading coefficient that is itself a polynomial. The subresultant divisions bound how large the coefficients get -- which is what the class remark is about, and it holds -- but nothing bounds the monomial count, and three steps in a 25-term lead times a 159-term remainder passes 512. Raising MaxTerms is the one-line version and it is the wrong shape, which the measurement shows rather than the reasoning: at 2048 three documented refusals become answers, including TermCountsPastTheCeilingAreRefused, whose subject is a direct product of two 495-term inputs, and RefusedRatherThanAnsweredWrongly, whose input has 1001 terms and is meant to be stopped at the door. The bound on what may be *asked* was never too small. What was too small was the bound on what may be passed through on the way to an answer that is itself small. So MaxIntermediateTerms is a second constant, 4096, threaded as an optional argument through Multiply, Power and DivideExact and passed only by the gcd's remainder sequence. Everything reached any other way keeps MaxTerms by default, so the entry bound that protects the hot path is untouched and all three of those refusals still refuse. 3000 drawn triples 7 declined -> 0 declined the pair above declined -> a + b + c + d, checked as a greatest common divisor The sweep now asserts the count is zero rather than at most seven, so a refusal coming back fails there instead of passing quietly under a ceiling. Performance, which is the part #920 asked not be guessed at: TryCancel runs on every quotient the simplifier builds, so simpsweep is the hot path. Two runs at 4096 took 3m26.8s and 3m40.7s, one at 512 took 3m45.3s. The spread between repeats of the same setting is as large as the gap between the settings, so this supports no regression and does not support a speedup. The gcd test class is unchanged at 1 s. Verified: 6944 C# tests and 130 F# tests pass, 0 fail. propcheck 1340 checks 0 failures, simpsweep 10463/10463 agree 0 disagreements, rootcheck 596/596 clean, casbench 117/119 with 0 wrong, 0 error and 0 timeout. Co-authored-by: Claude Opus 5 --- .../Polynomials/MultivariatePolynomial.cs | 37 +++++++++++++++---- .../Algebra/Polynomials/PolynomialGcd.cs | 18 ++++----- .../Polynomials/MultivariateGcdTest.cs | 31 ++++++++++------ 3 files changed, 58 insertions(+), 28 deletions(-) diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs index 851fba84a..c41749132 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs @@ -63,6 +63,27 @@ internal sealed partial class MultivariatePolynomial /// internal const int MaxTerms = 512; + /// + /// The ceiling for an intermediate of a calculation whose input and answer are both + /// already inside . + /// + /// + /// A multivariate pseudo-remainder multiplies through by a leading coefficient that is + /// itself a polynomial, so its intermediates grow in monomial count even when neither + /// side is anywhere near the bound — the subresultant divisions bound how big the + /// coefficients get, not how many terms there are. Holding those intermediates to the + /// input bound refused a gcd of a 19-term and a 29-term pair. + /// https://github.com/asc-community/AngouriMath/issues/920 + /// + /// It is deliberately a second constant rather than a larger . + /// The input bound is what protects the hot path — TryCancel runs on every + /// quotient the simplifier builds — and raising it turned three documented refusals + /// into answers, including a direct product of two 495-term inputs. What was too small + /// was never the bound on what may be asked, only the bound on what may be passed + /// through on the way to an answer that is itself small. + /// + internal const int MaxIntermediateTerms = 4096; + private const int BitsPerVariable = 8; private const ulong PowerMask = 0xFF; @@ -163,7 +184,7 @@ private static void Accumulate(Dictionary into, ulong monomial into[monomial] = sum; } - internal MultivariatePolynomial? Multiply(MultivariatePolynomial other) + internal MultivariatePolynomial? Multiply(MultivariatePolynomial other, int maxTerms = MaxTerms) { if (IsZero || other.IsZero) return Zero(VariableCount); @@ -174,13 +195,13 @@ private static void Accumulate(Dictionary into, ulong monomial if (!TryMultiplyMonomials(left.Key, right.Key, VariableCount, out var monomial)) return null; Accumulate(result, monomial, left.Value.Multiply(right.Value).ToLowestTerms()); - if (result.Count > MaxTerms) + if (result.Count > maxTerms) return null; } return new(VariableCount, result); } - internal MultivariatePolynomial? Power(int exponent) + internal MultivariatePolynomial? Power(int exponent, int maxTerms = MaxTerms) { if (exponent < 0 || exponent > MaxDegree) return null; @@ -190,14 +211,14 @@ private static void Accumulate(Dictionary into, ulong monomial { if ((exponent & 1) == 1) { - if (result.Multiply(square) is not { } multiplied) + if (result.Multiply(square, maxTerms) is not { } multiplied) return null; result = multiplied; } exponent >>= 1; if (exponent == 0) break; - if (square.Multiply(square) is not { } squared) + if (square.Multiply(square, maxTerms) is not { } squared) return null; square = squared; } @@ -292,7 +313,7 @@ internal MultivariatePolynomial LeadingCoefficientIn(int variable) /// That is the check the caller relies on: nothing is cancelled that has not been /// divided out and seen to leave nothing behind. /// - internal MultivariatePolynomial? DivideExact(MultivariatePolynomial divisor) + internal MultivariatePolynomial? DivideExact(MultivariatePolynomial divisor, int maxTerms = MaxTerms) { if (divisor.IsZero) return null; @@ -305,7 +326,7 @@ internal MultivariatePolynomial LeadingCoefficientIn(int variable) var divisorValue = divisor.terms[divisorLead]; var quotient = new Dictionary(); var rest = this; - for (var step = 0; step <= MaxTerms; step++) + for (var step = 0; step <= maxTerms; step++) { if (rest.IsZero) return new(VariableCount, quotient); @@ -317,7 +338,7 @@ internal MultivariatePolynomial LeadingCoefficientIn(int variable) if (divisor.MultiplyByTerm(monomial, value) is not { } product) return null; rest = rest.Subtract(product); - if (rest.TermCount > MaxTerms) + if (rest.TermCount > maxTerms) return null; } return null; diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialGcd.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialGcd.cs index 7865fb64b..ee3cc164a 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialGcd.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialGcd.cs @@ -240,9 +240,9 @@ internal static bool TryCancel(Entity numerator, Entity denominator, // The division below is the whole point of the subresultant sequence: what // it leaves is a subresultant, so it comes out exact, and the coefficients // stay the size of the subresultants instead of compounding. - if (scale.Power(delta) is not { } scalePower - || previousLead.Multiply(scalePower) is not { } factor - || remainder.DivideExact(factor) is not { } next) + if (scale.Power(delta, MultivariatePolynomial.MaxIntermediateTerms) is not { } scalePower + || previousLead.Multiply(scalePower, MultivariatePolynomial.MaxIntermediateTerms) is not { } factor + || remainder.DivideExact(factor, MultivariatePolynomial.MaxIntermediateTerms) is not { } next) return null; left = right; @@ -252,9 +252,9 @@ internal static bool TryCancel(Entity numerator, Entity denominator, scale = previousLead; else if (delta > 1) { - if (previousLead.Power(delta) is not { } raised - || scale.Power(delta - 1) is not { } divisor - || raised.DivideExact(divisor) is not { } updated) + if (previousLead.Power(delta, MultivariatePolynomial.MaxIntermediateTerms) is not { } raised + || scale.Power(delta - 1, MultivariatePolynomial.MaxIntermediateTerms) is not { } divisor + || raised.DivideExact(divisor, MultivariatePolynomial.MaxIntermediateTerms) is not { } updated) return null; scale = updated; } @@ -281,8 +281,8 @@ internal static bool TryCancel(Entity numerator, Entity denominator, break; MultithreadingFunctional.ExitIfCancelled(); var shift = remainder.DegreeIn(main) - divisorDegree; - if (divisorLead.Multiply(remainder) is not { } scaled - || remainder.LeadingCoefficientIn(main).Multiply(divisor) is not { } cancelling + if (divisorLead.Multiply(remainder, MultivariatePolynomial.MaxIntermediateTerms) is not { } scaled + || remainder.LeadingCoefficientIn(main).Multiply(divisor, MultivariatePolynomial.MaxIntermediateTerms) is not { } cancelling || cancelling.ShiftedBy(main, shift) is not { } shifted) return null; remainder = scaled.Subtract(shifted); @@ -292,7 +292,7 @@ internal static bool TryCancel(Entity numerator, Entity denominator, return remainder; for (var i = 0; i < outstanding; i++) { - if (divisorLead.Multiply(remainder) is not { } scaled) + if (divisorLead.Multiply(remainder, MultivariatePolynomial.MaxIntermediateTerms) is not { } scaled) return null; remainder = scaled; } diff --git a/Sources/Tests/UnitTests/Algebra/Polynomials/MultivariateGcdTest.cs b/Sources/Tests/UnitTests/Algebra/Polynomials/MultivariateGcdTest.cs index e77710520..f902c85b8 100644 --- a/Sources/Tests/UnitTests/Algebra/Polynomials/MultivariateGcdTest.cs +++ b/Sources/Tests/UnitTests/Algebra/Polynomials/MultivariateGcdTest.cs @@ -295,10 +295,13 @@ static int[] Draw(Random random) /// clock-seeded generator out. /// /// - /// Seven of these three thousand are declined rather than answered, and the count is - /// asserted rather than ignored: - /// is what reaches the ceiling and why the inputs that do it are so small. The bound - /// is an upper one, so finding fewer never fails. + /// Seven of these three thousand used to be declined rather than answered, all of them + /// for the reason in — + /// an intermediate past the input bound on the way to a small answer. Since + /// that is none of them, and + /// the count is asserted at zero rather than bounded, so a refusal reappearing fails + /// here instead of passing quietly under a ceiling. + /// https://github.com/asc-community/AngouriMath/issues/920 /// [Fact] public void GcdIsMultiplicativeOverManyDrawnTriples() @@ -322,7 +325,7 @@ public void GcdIsMultiplicativeOverManyDrawnTriples() Assert.True(divisor.SameAs(expected), $"trial {trial}: the divisor is not the expected one"); AssertIsGreatestCommonDivisor(first, second, divisor, SweepVariables.Length); } - Assert.True(declined <= 7, $"{declined} of 3000 triples were declined"); + Assert.True(declined == 0, $"{declined} of 3000 triples were declined"); } #endregion @@ -501,12 +504,14 @@ public void TermCountsPastTheCeilingAreRefused() /// wrong answer is that the step declines. /// /// - /// Pinned as a refusal, which is a legitimate answer. Should the ceiling be raised, or - /// the sequence learn to keep its intermediates primitive, this becomes an answer and - /// the test should be changed to assert that answer deliberately. + /// This was pinned as a refusal, and is now the answer. The intermediate is held to + /// rather than to the input + /// bound, which is the distinction the refusal was missing: nothing about the question + /// asked here is large, only something passed through on the way to it. + /// https://github.com/asc-community/AngouriMath/issues/920 /// [Fact] - public void TheTermCeilingIsReachedByAnIntermediate() + public void AnIntermediatePastTheInputCeilingIsStillAnswered() { var variables = new[] { "a", "b", "c", "d" }; var left = Polynomial("(b + c + 1) * (a + b) * (a + b + c + d)", variables); @@ -514,11 +519,15 @@ public void TheTermCeilingIsReachedByAnIntermediate() Assert.Equal(19, left.TermCount); Assert.Equal(29, right.TermCount); - // a + b + c + d divides both, and is not found. + // a + b + c + d divides both, and is now found rather than declined. var common = Polynomial("a + b + c + d", variables); Assert.NotNull(left.DivideExact(common)); Assert.NotNull(right.DivideExact(common)); - Assert.Null(Gcd(left, right, variables.Length)); + + var divisor = Gcd(left, right, variables.Length); + Assert.NotNull(divisor); + Assert.True(divisor.SameAs(common), "the divisor found is not a + b + c + d"); + AssertIsGreatestCommonDivisor(left, right, divisor, variables.Length); } #endregion