From 2d00591ffc5dbe671ae7583d39dc8fbe915285bb Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 15 Aug 2025 13:17:42 +0530 Subject: [PATCH 1/4] feat: add fast in-place `Polynomial + constant` implementation --- src/operators.jl | 41 +++++++++++++++++++++++++++++++++++++ test/mutable_arithmetics.jl | 39 +++++++++++++++++++++++++++++++++++ 2 files changed, 80 insertions(+) diff --git a/src/operators.jl b/src/operators.jl index e5d4e0c..bb52f0b 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -92,6 +92,47 @@ function MA.operate_to!( return output end +""" + _constant_term_idx(p::Polynomial) + +Return the index of the constant term in `p` according to its monomial ordering. +""" +_constant_term_idx(p::Polynomial{V, M, T}) where {V, M <: Reverse, T} = lastindex(p.x) +_constant_term_idx(p::Polynomial) = firstindex(p.x) + +""" + _insert_constant_term!(p::Polynomial) + +Insert a constant (degree 0) term into polynomial `p` at the appropriate position for the +monomial ordering of `p`. Does not check if a constant term already exists. +""" +function _insert_constant_term!(p::Polynomial{V, M, T}) where {V, M <: Reverse, T} + push!(MP.coefficients(p), zero(T)) + push!(MP.monomials(p).Z, zeros(Int, length(MP.variables(p)))) + return p +end + +function _insert_constant_term!(p::Polynomial{V, M, T}) where {V, M, T} + insert!(MP.coefficients(p), 1, zero(T)) + insert!(MP.monomials(p).Z, 1, zeros(Int, length(MP.variables(p)))) + return p +end + +function MA.operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{V, M, T}, x::T) where {V, M, T} + c_idx = _constant_term_idx(p) + if MP.nterms(p) == 0 || !MP.isconstant(MP.terms(p)[c_idx]) + _insert_constant_term!(p) + c_idx = _constant_term_idx(p) + end + coeffs = MP.coefficients(p) + coeffs[c_idx] = op(coeffs[c_idx], x) + if iszero(coeffs[c_idx]) + deleteat!(coeffs, c_idx) + deleteat!(MP.monomials(p), c_idx) + end + return p +end + function MA.operate!( op::Union{typeof(+),typeof(-)}, p::Polynomial, diff --git a/test/mutable_arithmetics.jl b/test/mutable_arithmetics.jl index d2a908c..f65aed1 100644 --- a/test/mutable_arithmetics.jl +++ b/test/mutable_arithmetics.jl @@ -35,3 +35,42 @@ using DynamicPolynomials @test q == x + y + 1 end end + +@testset "Fast path cases: $ord" for ord in [ + InverseLexOrder, + LexOrder, + Graded{InverseLexOrder}, + Graded{LexOrder}, + Reverse{InverseLexOrder}, + Reverse{LexOrder}, + Graded{Reverse{InverseLexOrder}}, + Graded{Reverse{LexOrder}}, + Reverse{Graded{InverseLexOrder}}, + Reverse{Graded{LexOrder}}, + ] + @polyvar x y z monomial_order=ord + + # allocation tests vary between Julia versions, so they're upper bounds + @testset "Polynomial + constant" begin + poly = 2 * x^2 + 3 * x * y + z * y^2 + poly2 = copy(poly) + result = poly + 2 + MA.operate!(+, poly, 2) + @test isequal(poly, result) + # down from 576 using the generic method + @test (@allocated MA.operate!(+, poly2, 2)) <= 272 + # subsequent additions don't allocate + @test (@allocated MA.operate!(+, poly2, 2)) == 0 + + # also test `-` + poly = 2 * x^2 + 3 * x * y + z * y^2 + poly2 = copy(poly) + result = poly - 2 + MA.operate!(-, poly, 2) + @test isequal(poly, result) + # down from 576 using the generic method + @test (@allocated MA.operate!(-, poly2, 2)) <= 272 + # subsequent additions don't allocate + @test (@allocated MA.operate!(-, poly2, 2)) == 0 + end +end From fd8154fb50ec4aec2553317961760533c0cd284c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 15 Aug 2025 18:23:23 +0530 Subject: [PATCH 2/4] feat: add fast in-place `Polynomial + Variable` implementation --- src/operators.jl | 27 +++++++++++++++++++++++++++ test/mutable_arithmetics.jl | 26 ++++++++++++++++++++++++++ 2 files changed, 53 insertions(+) diff --git a/src/operators.jl b/src/operators.jl index bb52f0b..5daac6c 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -133,6 +133,33 @@ function MA.operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{V, M, T}, x: return p end +function MA.operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{V, M, T}, x::Variable{V, M}) where {V, M, T} + vars = MP.variables(p) + idx = searchsortedfirst(vars, x; rev = true) + monos = MP.monomials(p) + if idx > length(vars) || !isequal(vars[idx], x) + for mono in monos + insert!(MP.exponents(mono), idx, 0) + end + insert!(vars, idx, x) + end + mono = Monomial{V, M}(vars, zeros(Int, length(vars))) + mono.z[idx] = 1 + idx = searchsortedfirst(monos, mono) + coeffs = MP.coefficients(p) + N = MP.nterms(p) + if idx > N || !isequal(monos[idx], mono) + insert!(monos.Z, idx, MP.exponents(mono)) + insert!(coeffs, idx, zero(T)) + end + coeffs[idx] = op(coeffs[idx], one(T)) + if iszero(coeffs[idx]) + deleteat!(coeffs, idx) + deleteat!(monos, idx) + end + return p +end + function MA.operate!( op::Union{typeof(+),typeof(-)}, p::Polynomial, diff --git a/test/mutable_arithmetics.jl b/test/mutable_arithmetics.jl index f65aed1..5529013 100644 --- a/test/mutable_arithmetics.jl +++ b/test/mutable_arithmetics.jl @@ -73,4 +73,30 @@ end # subsequent additions don't allocate @test (@allocated MA.operate!(-, poly2, 2)) == 0 end + + @testset "Polynomial + Variable" begin + poly = 2 * x^2 + 3 * x * y + z * y^2 + poly2 = copy(poly) + result = poly + x + MA.operate!(+, poly, x) + @test isequal(poly, result) + # down from 18752 using the generic method + # 368 or 304 depending on ordering, more for different version + # pre is especially bad + @test (@allocated MA.operate!(+, poly2, x)) <= 400 + # down from 1904 using the generic method + @test (@allocated MA.operate!(+, poly2, x)) <= 144 + + # also test `-` + poly = 2 * x^2 + 3 * x * y + z * y^2 + poly2 = copy(poly) + result = poly - x + MA.operate!(-, poly, x) + @test isequal(poly, result) + # down from 18752 using the generic method + # 368 or 304 depending on ordering + @test (@allocated MA.operate!(-, poly2, x)) <= 400 + # down from 1904 using the generic method + @test (@allocated MA.operate!(-, poly2, x)) <= 144 + end end From 7e99f398f0b6ea8e2bc6653b04e392690436e61c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sun, 17 Aug 2025 10:40:09 +0530 Subject: [PATCH 3/4] Update src/operators.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Benoît Legat --- src/operators.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/operators.jl b/src/operators.jl index 5daac6c..cbacb4f 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -104,7 +104,7 @@ _constant_term_idx(p::Polynomial) = firstindex(p.x) _insert_constant_term!(p::Polynomial) Insert a constant (degree 0) term into polynomial `p` at the appropriate position for the -monomial ordering of `p`. Does not check if a constant term already exists. +monomial ordering of `p`. Assume that a constant term does not already exists. """ function _insert_constant_term!(p::Polynomial{V, M, T}) where {V, M <: Reverse, T} push!(MP.coefficients(p), zero(T)) From e6f14583e5ab70ace28e7276cc9884998d5b9941 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sun, 17 Aug 2025 10:41:43 +0530 Subject: [PATCH 4/4] refactor: rename `_constant_term_idx` to `_lowest_term_idx` --- src/operators.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/operators.jl b/src/operators.jl index cbacb4f..dd27226 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -93,12 +93,12 @@ function MA.operate_to!( end """ - _constant_term_idx(p::Polynomial) + _lowest_term_idx(p::Polynomial) -Return the index of the constant term in `p` according to its monomial ordering. +Return the index of the lowest term in `p` according to its monomial ordering. """ -_constant_term_idx(p::Polynomial{V, M, T}) where {V, M <: Reverse, T} = lastindex(p.x) -_constant_term_idx(p::Polynomial) = firstindex(p.x) +_lowest_term_idx(p::Polynomial{V, M, T}) where {V, M <: Reverse, T} = lastindex(p.x) +_lowest_term_idx(p::Polynomial) = firstindex(p.x) """ _insert_constant_term!(p::Polynomial) @@ -119,10 +119,10 @@ function _insert_constant_term!(p::Polynomial{V, M, T}) where {V, M, T} end function MA.operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{V, M, T}, x::T) where {V, M, T} - c_idx = _constant_term_idx(p) + c_idx = _lowest_term_idx(p) if MP.nterms(p) == 0 || !MP.isconstant(MP.terms(p)[c_idx]) _insert_constant_term!(p) - c_idx = _constant_term_idx(p) + c_idx = _lowest_term_idx(p) end coeffs = MP.coefficients(p) coeffs[c_idx] = op(coeffs[c_idx], x)