Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 1 commit into from
Aug 13, 2026
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 29 additions & 8 deletions Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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);
Expand All @@ -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;
Expand All @@ -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;
}
Expand Down Expand Up @@ -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;
Expand All @@ -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);
Expand All @@ -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;
Expand Down
18 changes: 9 additions & 9 deletions Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialGcd.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
}
Expand All @@ -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);
Expand All @@ -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;
}
Expand Down
31 changes: 20 additions & 11 deletions Sources/Tests/UnitTests/Algebra/Polynomials/MultivariateGcdTest.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading