Bound a gcd's intermediates separately from its inputs (#920) (#923)

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 <noreply@anthropic.com>

Author: Rafael-SOWNet

Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs

@@ -63,6 +63,27 @@ internal sealed partial class MultivariatePolynomial
         /// </summary>
         internal const int MaxTerms = 512;
 
+        /// <summary>
+        /// The ceiling for an intermediate of a calculation whose input and answer are both
+        /// already inside <see cref="MaxTerms"/>.
+        /// </summary>
+        /// <remarks>
+        /// 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 <see cref="MaxTerms"/>.
+        /// The input bound is what protects the hot path — <c>TryCancel</c> 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.
+        /// </remarks>
+        internal const int MaxIntermediateTerms = 4096;
+
         private const int BitsPerVariable = 8;
         private const ulong PowerMask = 0xFF;
 
@@ -163,7 +184,7 @@ private static void Accumulate(Dictionary<ulong, ERational> 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<ulong, ERational> 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<ulong, ERational> 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.
         /// </remarks>
-        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<ulong, ERational>();
             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;

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;
             }

Sources/Tests/UnitTests/Algebra/Polynomials/MultivariateGcdTest.cs

@@ -295,10 +295,13 @@ static int[] Draw(Random random)
         /// clock-seeded generator out.
         /// </summary>
         /// <remarks>
-        /// Seven of these three thousand are declined rather than answered, and the count is
-        /// asserted rather than ignored: <see cref="TheTermCeilingIsReachedByAnIntermediate"/>
-        /// 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 <see cref="AnIntermediatePastTheInputCeilingIsStillAnswered"/> —
+        /// an intermediate past the input bound on the way to a small answer. Since
+        /// <see cref="MultivariatePolynomial.MaxIntermediateTerms"/> 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
         /// </remarks>
         [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,24 +504,30 @@ public void TermCountsPastTheCeilingAreRefused()
         /// wrong answer is that the step declines.
         /// </summary>
         /// <remarks>
-        /// 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
+        /// <see cref="MultivariatePolynomial.MaxIntermediateTerms"/> 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
         /// </remarks>
         [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);
             var right = Polynomial("(a ^ 2 + b * c + d) * (a + b + c + d) * (a + b + c + d)", variables);
             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