From 56afa5499682bfa05f0a00a6543d0527c95314ad Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 14 Aug 2026 02:07:04 +0000 Subject: [PATCH] Bound the resultant's elimination by measured work rather than by reasoned size (#921) MaxSylvesterSize = 24 was arithmetic about a step count, not a measurement of anything, and #921 asked for the two axes that actually matter to be measured: Sylvester size, and terms per entry. They were, across 5 widths and 9 sizes, timing and allocation both, and the answer is that a bound on the input shape cannot do the job. width size ms allocated outcome (before this change; ceiling lifted to sweep) 1 40 31 35 MB answered 2 24 422 689 MB answered 2 40 2880 8.4 GB answered 3 16 838 2.9 GB answered 3 24 6248 21 GB declined 3 40 21980 79 GB declined 8 24 2233 8.4 GB declined Three things follow, none of them visible from the step count. **Size is the cheap axis and 24 was far too mean for it.** With scalar entries the elimination does size^3/3 units of work; size 40 is 31ms and 35MB. The ceiling is now 40, which admits the discriminant of anything up to degree 20 rather than 12. The largest admitted elimination -- two degree-20 polynomials with roots 1..20 and 21..40, checked against the product over the differences of the roots -- takes 35ms. **Cost is violent in terms per entry**, between the fifth and the eighth power of it over the range, so no product of the two axes is the right law either. **And the most expensive input is not the widest**, which is what makes any bound read off the degrees wrong in principle. Entries wide enough trip MultivariatePolynomial.MaxTerms early and the elimination declines cheaply; 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. So the second bound is on the work done rather than on the work predicted. Each step charges the product of the term counts of the operands of each of its two products -- what a sparse multiplication costs -- and the elimination stops past MaxEliminationWork = 2,500,000, charged before the multiplication so the budget cannot be overshot by the step that was going to be the worst. Over the whole sweep that unit predicts both resources linearly and tightly, 1.4KB of allocation and 0.4us per unit, so the budget is about a second and about three gigabytes and moving it moves both together. every refusal in the sweep up to 22s and 79 GB -> 1.0-1.2s and 3.5-4.2 GB whole sweep 2m 50s -> 37s Nothing that answered within the old ceiling stops answering, and the cases that declined there decline up to six times faster. Above it the raise adds answers: scalar entries to size 40, two-term entries to size 32, the discriminant of a degree-16 polynomial. One case is refused that the sweep showed answering -- two-term entries at size 40, 2.9 seconds and 8.4GB -- and the old ceiling refused that one too. The cost of the raise is that a large wide input now takes up to a second to be refused where a ceiling of 24 refused it for nothing; that is the price of no longer refusing the cheap ones with it. No BREAKING-CHANGES entry: PolynomialResultant is internal and nothing in the library calls it yet, which is also why #921 said this blocks nothing and is worth settling before something does. Measured: suite 6961 passed / 0 failed; F# wrapper 130/130; casbench 116/119 with 0 wrong, 0 error, 0 timeout; propcheck 1340 checks / 0 failures; rootcheck 596/596 clean; simpsweep 10463/10463 agree. https://github.com/asc-community/AngouriMath/issues/921 Co-authored-by: Claude Opus 5 --- .../Polynomials/PolynomialResultant.cs | 58 +++++++++++++++++-- .../Polynomials/PolynomialResultantTest.cs | 45 +++++++++++++- 2 files changed, 94 insertions(+), 9 deletions(-) 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; + } } }