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