Add polynomial!, in-place version of polynomial

Author: blegat

src/differentiation.jl

@@ -37,7 +37,7 @@ differentiate(α::T, v::AbstractVariable) where T = zero(T)
 differentiate(v1::AbstractVariable, v2::AbstractVariable) = v1 == v2 ? 1 : 0
 differentiate(t::AbstractTermLike, v::AbstractVariable) = coefficient(t) * differentiate(monomial(t), v)
 # The polynomial function will take care of removing the zeros
-differentiate(p::APL, v::AbstractVariable) = polynomial(differentiate.(terms(p), v), SortedState())
+differentiate(p::APL, v::AbstractVariable) = polynomial!(differentiate.(terms(p), v), SortedState())
 differentiate(p::RationalPoly, v::AbstractVariable) = (differentiate(p.num, v) * p.den - p.num * differentiate(p.den, v)) / p.den^2
 
 const ARPL = Union{APL, RationalPoly}

src/operators.jl

@@ -182,15 +182,15 @@ function mul_to_terms!(ts::Vector{<:AbstractTerm}, p1::APL, p2::APL)
     return ts
 end
 function Base.:*(p::AbstractPolynomial, q::AbstractPolynomial)
-    polynomial(mul_to_terms!(MA.promote_operation(*, termtype(p), termtype(q))[], p, q))
+    polynomial!(mul_to_terms!(MA.promote_operation(*, termtype(p), termtype(q))[], p, q))
 end
 
 Base.isapprox(p::APL, α; kwargs...) = isapprox(promote(p, α)...; kwargs...)
 Base.isapprox(α, p::APL; kwargs...) = isapprox(promote(p, α)...; kwargs...)
 
 Base.:-(m::AbstractMonomialLike) = (-1) * m
 Base.:-(t::AbstractTermLike) = (-coefficient(t)) * monomial(t)
-Base.:-(p::APL) = polynomial((-).(terms(p)))
+Base.:-(p::APL) = polynomial!((-).(terms(p)))
 Base.:+(p::Union{APL, RationalPoly}) = p
 Base.:*(p::Union{APL, RationalPoly}) = p
 
@@ -249,8 +249,8 @@ Base.:*(m::AbstractMonomialLike, t::AbstractTermLike) = coefficient(t) * (m * mo
 Base.:*(t::AbstractTermLike, m::AbstractMonomialLike) = coefficient(t) * (monomial(t) * m)
 Base.:*(t1::AbstractTermLike, t2::AbstractTermLike) = (coefficient(t1) * coefficient(t2)) * (monomial(t1) * monomial(t2))
 
-Base.:*(t::AbstractTermLike, p::APL) = polynomial(map(te -> t * te, terms(p)))
-Base.:*(p::APL, t::AbstractTermLike) = polynomial(map(te -> te * t, terms(p)))
+Base.:*(t::AbstractTermLike, p::APL) = polynomial!(map(te -> t * te, terms(p)))
+Base.:*(p::APL, t::AbstractTermLike) = polynomial!(map(te -> te * t, terms(p)))
 Base.:*(p::APL, q::APL) = polynomial(p) * polynomial(q)
 
 # guaranteed that monomial(t1) > monomial(t2)
@@ -260,6 +260,7 @@ function _polynomial_2terms(t1::TT, t2::TT, ::Type{T}) where {TT<:AbstractTerm,
     elseif iszero(t2)
         polynomial(t1, T)
     else
+        # not `polynomial!` because we `t1` and `t2` cannot be modified
         polynomial(termtype(TT, T)[t1, t2], SortedUniqState())
     end
 end

src/polynomial.jl

@@ -1,4 +1,4 @@
-export polynomial, polynomialtype, terms, nterms, coefficients, monomials
+export polynomial, polynomial!, polynomialtype, terms, nterms, coefficients, monomials
 export coefficienttype, monomialtype
 export mindegree, maxdegree, extdegree
 export leadingterm, leadingcoefficient, leadingmonomial
@@ -49,17 +49,32 @@ Calling `polynomial([2, 4, 1], [x, x^2*y, x*y])` should return ``4x^2y + xy + 2x
 polynomial(p::AbstractPolynomial) = p
 polynomial(p::APL{T}, ::Type{T}) where T = polynomial(terms(p))
 polynomial(p::APL{T}) where T = polynomial(p, T)
-polynomial(ts::AbstractVector, s::ListState=MessyState()) = sum(ts)
-polynomial(ts::AbstractVector{<:AbstractTerm}, s::SortedUniqState) = polynomial(coefficient.(ts), monomial.(ts), s)
-polynomial(a::AbstractVector, x::AbstractVector, s::ListState=MessyState()) = polynomial([α * m for (α, m) in zip(a, x)], s)
-polynomial(f::Function, mv::AbstractVector{<:AbstractMonomialLike}) = polynomial([f(i) * mv[i] for i in 1:length(mv)])
 function polynomial(Q::AbstractMatrix, mv::AbstractVector)
     LinearAlgebra.dot(mv, Q * mv)
 end
 function polynomial(Q::AbstractMatrix, mv::AbstractVector, ::Type{T}) where T
     polynomial(polynomial(Q, mv), T)
 end
 
+polynomial(f::Function, mv::AbstractVector{<:AbstractMonomialLike}) = polynomial!([f(i) * mv[i] for i in 1:length(mv)])
+
+polynomial(a::AbstractVector, x::AbstractVector, s::ListState=MessyState()) = polynomial([α * m for (α, m) in zip(a, x)], s)
+
+polynomial(ts::AbstractVector, s::ListState=MessyState()) = sum(ts)
+
+polynomial!(ts::AbstractVector{<:AbstractTerm}, s::SortedUniqState) = polynomial(coefficient.(ts), monomial.(ts), s)
+
+function polynomial!(ts::AbstractVector{<:AbstractTerm}, s::SortedState)
+    polynomial!(uniqterms!(ts), SortedUniqState())
+end
+function polynomial!(ts::AbstractVector{<:AbstractTerm}, s::UnsortedState=MessyState())
+    polynomial!(sort!(ts, lt=(>)), sortstate(s))
+end
+
+_collect(v::Vector) = v
+_collect(v::AbstractVector) = collect(v)
+polynomial(ts::AbstractVector{<:AbstractTerm}, args::Vararg{ListState, N}) where {N} = polynomial!(MA.mutable_copy(_collect(ts)), args...)
+
 """
     polynomialtype(p::AbstractPolynomialLike)
 
@@ -85,25 +100,6 @@ polynomialtype(::Union{P, Type{P}}, ::Type{T}) where {P <: APL, T} = polynomialt
 polynomialtype(::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}) where PT <: APL = polynomialtype(PT)
 polynomialtype(::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}, ::Type{T}) where {PT <: APL, T} = polynomialtype(PT, T)
 
-function uniqterms(ts::AbstractVector{T}) where T <: AbstractTerm
-    result = T[]
-    sizehint!(result, length(ts))
-    for t in ts
-        if !iszero(t)
-            if isempty(result) || monomial(t) != monomial(last(result))
-                push!(result, t)
-            else
-                coef = coefficient(last(result)) + coefficient(t)
-                if iszero(coef)
-                    pop!(result)
-                else
-                    result[end] = coef * monomial(t)
-                end
-            end
-        end
-    end
-    result
-end
 function uniqterms!(ts::AbstractVector{<: AbstractTerm})
     i = firstindex(ts)
     for j in Iterators.drop(eachindex(ts), 1)
@@ -126,8 +122,6 @@ function uniqterms!(ts::AbstractVector{<: AbstractTerm})
     end
     ts
 end
-polynomial(ts::AbstractVector{<:AbstractTerm}, s::SortedState) = polynomial(uniqterms(ts), SortedUniqState())
-polynomial(ts::AbstractVector{<:AbstractTerm}, s::UnsortedState=MessyState()) = polynomial(sort(ts, lt=(>)), sortstate(s))
 
 """
     terms(p::AbstractPolynomialLike)
@@ -358,7 +352,7 @@ Returns `p / leadingcoefficient(p)` where the leading coefficient of the returne
 """
 function monic(p::APL)
     α = leadingcoefficient(p)
-    polynomial(_divtoone.(terms(p), α))
+    polynomial!(_divtoone.(terms(p), α))
 end
 monic(m::AbstractMonomialLike) = m
 monic(t::AbstractTermLike{T}) where T = one(T) * monomial(t)
@@ -386,14 +380,14 @@ function mapcoefficientsnz(f::Function, p::AbstractPolynomialLike)
     # Invariant: p has only nonzero coefficient
     # therefore f(α) will be nonzero for every coefficient α of p
     # hence we can use Uniq
-    polynomial(mapcoefficientsnz.(f, terms(p)), SortedUniqState())
+    polynomial!(mapcoefficientsnz.(f, terms(p)), SortedUniqState())
 end
 mapcoefficientsnz(f::Function, t::AbstractTermLike) = f(coefficient(t)) * monomial(t)
 
 Base.round(t::AbstractTermLike; args...) = round(coefficient(t); args...) * monomial(t)
 function Base.round(p::AbstractPolynomialLike; args...)
     # round(0.1) is zero so we cannot use SortedUniqState
-    polynomial(round.(terms(p); args...), SortedState())
+    polynomial!(round.(terms(p); args...), SortedState())
 end
 
 Base.ndims(::Type{<:AbstractPolynomialLike}) = 0