diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialResultant.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialResultant.cs
index c221abd22..640c629a4 100644
--- a/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialResultant.cs
+++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialResultant.cs
@@ -54,13 +54,51 @@ namespace AngouriMath.Functions
internal static class PolynomialResultant
{
///
- /// The elimination is cubic in the size of the Sylvester matrix, with a
- /// multiplication of two multivariate polynomials at every step, so the size is what
- /// decides whether this finishes. The bound is on deg f + deg g; 24 leaves a
- /// 13824-step elimination as the worst case admitted, and admits the discriminant of
- /// anything up to degree 12, since deg f + deg f' is 2 deg f - 1.
+ /// The bound on deg f + deg g. It exists to stop the matrix being built at
+ /// all; what the elimination then costs is bounded by
+ /// instead, and the two are not the same question.
///
- private const int MaxSylvesterSize = 24;
+ ///
+ /// Size is the cheap axis, which is why this is generous. With scalar entries the
+ /// elimination does size^3 / 3 units of the work counted below, and 40 measured
+ /// 31ms and 35MB of allocation — where the same size with three-term entries costs 22
+ /// seconds and 79GB. So a ceiling on size alone can only be right for one shape of
+ /// input at a time, and this one is set for the shape it can actually decide: 40
+ /// admits the discriminant of anything up to degree 20, since deg f + deg f' is
+ /// 2 deg f - 1. https://github.com/asc-community/AngouriMath/issues/921
+ ///
+ private const int MaxSylvesterSize = 40;
+
+ ///
+ /// The budget on the elimination itself, in units of one coefficient-pair
+ /// multiplication: each step adds the product of the term counts of the two operands
+ /// of each of its two products, which is what a sparse multiplication costs.
+ ///
+ ///
+ ///
+ /// Measured rather than reasoned, across the two axes that matter — Sylvester size,
+ /// and terms per entry. Over the whole range the unit predicts both resources
+ /// linearly and tightly: 1.4KB of allocation and 0.4µs per unit. So this budget
+ /// is about a second, and about three gigabytes of garbage, and moving it moves both
+ /// together.
+ ///
+ ///
+ /// A bound on the input shape cannot do this job, which is the measurement's real
+ /// finding. The cost is mild in size and violent in terms per entry — between the
+ /// fifth and the eighth power of it over the range — so no product of the two is the
+ /// right law. Worse, the most expensive input is not the widest: entries wide
+ /// enough trip early and the elimination
+ /// declines cheaply, while three-term entries grow just slowly enough to spend 6
+ /// seconds and 21GB before declining at size 24. The shape that costs the most is the
+ /// one just inside the term ceiling, and no bound read off the degrees can see it.
+ ///
+ ///
+ /// Set to admit every case that answered at all in the sweep bar one — two-term
+ /// entries at size 40, which answers in 2.9 seconds having allocated 8.4GB — and to
+ /// refuse the expensive refusals, which is where it earns its keep.
+ ///
+ ///
+ private const long MaxEliminationWork = 2_500_000;
[ConstantField] private static readonly ERational MinusOne = ERational.One.Negate();
@@ -217,6 +255,7 @@ private static bool TryDivideOut(
matrix[rightDegree + row][row + rightDegree - power] = coefficient;
var negated = false;
+ var work = 0L;
var previous = MultivariatePolynomial.One(variableCount);
for (var pivot = 0; pivot + 1 < size; pivot++)
{
@@ -243,6 +282,13 @@ private static bool TryDivideOut(
var leading = matrix[row][pivot];
for (var column = pivot + 1; column < size; column++)
{
+ // Charged before the multiplication rather than after it, so that the
+ // budget cannot be overshot by the one step that was going to be the
+ // most expensive of them.
+ work += (long)head.TermCount * matrix[row][column].TermCount
+ + (long)leading.TermCount * matrix[pivot][column].TermCount;
+ if (work > MaxEliminationWork)
+ return null;
if (head.Multiply(matrix[row][column]) is not { } kept
|| leading.Multiply(matrix[pivot][column]) is not { } removed
|| kept.Subtract(removed).DivideExact(previous) is not { } reduced)
diff --git a/Sources/Tests/UnitTests/Algebra/Polynomials/PolynomialResultantTest.cs b/Sources/Tests/UnitTests/Algebra/Polynomials/PolynomialResultantTest.cs
index b11d9dfcd..c8c76e339 100644
--- a/Sources/Tests/UnitTests/Algebra/Polynomials/PolynomialResultantTest.cs
+++ b/Sources/Tests/UnitTests/Algebra/Polynomials/PolynomialResultantTest.cs
@@ -520,7 +520,8 @@ public void TheDiscriminantOfAQuadraticInOneVariableSeesTheOther()
[Fact]
public void InputTooLargeForTheEliminationIsRefusedRatherThanAttempted()
{
- var coefficients = new ERational[21];
+ // deg f + deg g of 42, one past the ceiling on size.
+ var coefficients = new ERational[22];
for (var power = 0; power < coefficients.Length; power++)
coefficients[power] = Rational(power % 5 + 1);
var poly = Univariate(coefficients);
@@ -533,8 +534,8 @@ public void TheLargestAdmittedEliminationStillAnswers()
{
// deg f + deg g of exactly the ceiling, so that the refusal above reads as a
// bound rather than as a description of everything past a handful of terms.
- var leftRoots = new[] { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 };
- var rightRoots = new[] { 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24 };
+ var leftRoots = new[] { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 };
+ var rightRoots = new[] { 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40 };
AssertConstant(
ResultantFromRoots(1, leftRoots, 1, rightRoots),
PolynomialResultant.Resultant(
@@ -542,5 +543,43 @@ public void TheLargestAdmittedEliminationStillAnswers()
0, NoOtherVariables),
"Res(f, g) at the ceiling");
}
+
+ ///
+ /// Three variables, every coefficient in the main one carrying two monomials in the
+ /// other two, at a Sylvester size the ceiling admits. The elimination is refused all
+ /// the same, on the budget rather than on the size — which is the whole point of
+ /// there being two bounds: the cost of an elimination is not readable off its degrees.
+ /// https://github.com/asc-community/AngouriMath/issues/921
+ ///
+ [Fact]
+ public void AnEliminationPastTheWorkBudgetIsRefusedThoughItsSizeIsAdmitted()
+ {
+ var left = Wide(20, 2, seed: 1);
+ var right = Wide(20, 2, seed: 2);
+ Assert.Equal(20, left.DegreeIn(0));
+ Assert.Equal(20, right.DegreeIn(0));
+ Assert.Null(PolynomialResultant.Resultant(left, right, 0, new[] { 1, 2 }));
+ }
+
+ ///
+ /// A polynomial of the given degree in variable 0, each of whose coefficients carries
+ /// distinct monomials in variables 1 and 2.
+ ///
+ private static MultivariatePolynomial Wide(int degree, int width, int seed)
+ {
+ var shape = new[] { (A: 0, B: 0), (A: 1, B: 0), (A: 0, B: 1), (A: 1, B: 1) };
+ var result = MultivariatePolynomial.Zero(3);
+ for (var power = 0; power <= degree; power++)
+ for (var term = 0; term < width; term++)
+ {
+ var value = Rational(1 + (seed * 7 + power * 3 + term * 5) % 11);
+ var monomial = MultivariatePolynomial.Constant(3, value).ShiftedBy(0, power)
+ ?.ShiftedBy(1, shape[term % shape.Length].A)
+ ?.ShiftedBy(2, shape[term % shape.Length].B);
+ Assert.NotNull(monomial);
+ result = result.Add(monomial!);
+ }
+ return result;
+ }
}
}