|
| 1 | +struct Polynomial{CoeffType,T<:AbstractTerm{CoeffType},V<:AbstractVector{T}} <: AbstractPolynomial{CoeffType} |
| 2 | + terms::V |
| 3 | + |
| 4 | + Polynomial{C, T, V}(terms::AbstractVector{T}) where {C, T, V} = new{C, T, V}(terms) |
| 5 | +end |
| 6 | + |
| 7 | +function coefficients(p::Polynomial) |
| 8 | + return LazyMap{coefficienttype(p)}(coefficient, terms(p)) |
| 9 | +end |
| 10 | +function monomials(p::Polynomial) |
| 11 | + return LazyMap{monomialtype(p)}(monomial, terms(p)) |
| 12 | +end |
| 13 | + |
| 14 | +const VectorPolynomial{C,T} = Polynomial{C,T,Vector{T}} |
| 15 | + |
| 16 | +termtype(::Type{<:Polynomial{C,T}}) where {C,T} = T |
| 17 | +terms(p::Polynomial) = p.terms |
| 18 | +constantmonomial(::Union{Polynomial{C,TT},Type{Polynomial{C,TT}}}) where {C,TT} = constantmonomial(TT) |
| 19 | +Base.convert(::Type{Polynomial{C,TT,V}}, p::Polynomial{C,TT,V}) where {C,TT,V} = p |
| 20 | +function Base.convert(PT::Type{Polynomial{C,TT,V}}, p::AbstractPolynomialLike) where {C,TT,V} |
| 21 | + return PT(convert(V, terms(p))) |
| 22 | +end |
| 23 | +function Base.convert(::Type{Polynomial{C}}, p::AbstractPolynomialLike) where {C} |
| 24 | + TT = termtype(p, C) |
| 25 | + return convert(Polynomial{C,TT,Vector{TT}}, p) |
| 26 | +end |
| 27 | +function Base.convert(::Type{Polynomial}, p::AbstractPolynomialLike) |
| 28 | + return convert(Polynomial{coefficienttype(p)}, p) |
| 29 | +end |
| 30 | + |
| 31 | + |
| 32 | +_change_eltype(::Type{<:Vector}, ::Type{T}) where T = Vector{T} |
| 33 | +function polynomialtype(::Type{Polynomial{C, T, V}}, ::Type{NewC}) where {C, T, V, NewC} |
| 34 | + NewT = termtype(T, NewC) |
| 35 | + Polynomial{NewC, NewT, _change_eltype(V, NewT)} |
| 36 | +end |
| 37 | + |
| 38 | +function Base.promote_rule(::Type{Polynomial}, ::Type{<:AbstractPolynomialLike}) |
| 39 | + return Polynomial |
| 40 | +end |
| 41 | +function Base.promote_rule(::Type{<:AbstractPolynomialLike}, ::Type{Polynomial}) |
| 42 | + return Polynomial |
| 43 | +end |
| 44 | +function promote_rule_constant(::Type{T}, ::Type{Polynomial}) where T |
| 45 | + return Any |
| 46 | + # `convert(Polynomial{T}, ::T)` cannot be implemented as we don't know the type of the term |
| 47 | + #return Polynomial{T} |
| 48 | +end |
| 49 | + |
| 50 | +(p::Polynomial)(s...) = substitute(Eval(), p, s) |
| 51 | + |
| 52 | +function Base.one(::Type{P}) where {C,TT,P<:Polynomial{C,TT}} |
| 53 | + return convert(P, one(TT)) |
| 54 | +end |
| 55 | +Base.one(p::Polynomial) = one(typeof(p)) |
| 56 | + |
| 57 | +Base.zero(::Type{Polynomial{C,T,A}}) where {C,T,A} = Polynomial{C,T,A}(A()) |
| 58 | +Base.zero(t::Polynomial) = zero(typeof(t)) |
| 59 | + |
| 60 | +jointerms(terms1::AbstractArray{<:Term}, terms2::AbstractArray{<:Term}) = Sequences.mergesorted(terms1, terms2, compare, combine) |
| 61 | +function jointerms!(output::AbstractArray{<:Term}, terms1::AbstractArray{<:Term}, terms2::AbstractArray{<:Term}) |
| 62 | + resize!(output, length(terms1) + length(terms2)) |
| 63 | + Sequences.mergesorted!(output, terms1, terms2, compare, combine) |
| 64 | +end |
| 65 | + |
| 66 | +Base.:(+)(p1::Polynomial, p2::Polynomial) = polynomial!(jointerms(terms(p1), terms(p2)), SortedUniqState()) |
| 67 | +Base.:(+)(p1::Polynomial, p2::Polynomial{<:LinearAlgebra.UniformScaling}) = p1 + mapcoefficientsnz(J -> J.λ, p2) |
| 68 | +Base.:(+)(p1::Polynomial{<:LinearAlgebra.UniformScaling}, p2::Polynomial) = mapcoefficientsnz(J -> J.λ, p1) + p2 |
| 69 | +Base.:(+)(p1::Polynomial{<:LinearAlgebra.UniformScaling}, p2::Polynomial{<:LinearAlgebra.UniformScaling}) = mapcoefficientsnz(J -> J.λ, p1) + p2 |
| 70 | +function MA.operate_to!(result::Polynomial, ::typeof(+), p1::Polynomial, p2::Polynomial) |
| 71 | + if result === p1 || result === p2 |
| 72 | + error("Cannot call `operate_to!(output, +, p, q)` with `output` equal to `p` or `q`, call `operate!` instead.") |
| 73 | + end |
| 74 | + jointerms!(result.terms, terms(p1), terms(p2)) |
| 75 | + result |
| 76 | +end |
| 77 | +function MA.operate_to!(result::Polynomial, ::typeof(*), p::Polynomial, t::AbstractTermLike) |
| 78 | + if iszero(t) |
| 79 | + MA.operate!(zero, result) |
| 80 | + else |
| 81 | + resize!(result.terms, nterms(p)) |
| 82 | + for i in eachindex(p.terms) |
| 83 | + # TODO could use MA.mul_to!! for indices that were presents in `result` before the `resize!`. |
| 84 | + result.terms[i] = p.terms[i] * t |
| 85 | + end |
| 86 | + return result |
| 87 | + end |
| 88 | +end |
| 89 | +function MA.operate_to!(result::Polynomial, ::typeof(*), t::AbstractTermLike, p::Polynomial) |
| 90 | + if iszero(t) |
| 91 | + MA.operate!(zero, result) |
| 92 | + else |
| 93 | + resize!(result.terms, nterms(p)) |
| 94 | + for i in eachindex(p.terms) |
| 95 | + # TODO could use MA.mul_to!! for indices that were presents in `result` before the `resize!`. |
| 96 | + result.terms[i] = t * p.terms[i] |
| 97 | + end |
| 98 | + return result |
| 99 | + end |
| 100 | +end |
| 101 | +Base.:(-)(p1::Polynomial, p2::Polynomial) = polynomial!(jointerms(terms(p1), (-).(terms(p2)))) |
| 102 | +Base.:(-)(p1::Polynomial, p2::Polynomial{<:LinearAlgebra.UniformScaling}) = p1 - mapcoefficientsnz(J -> J.λ, p2) |
| 103 | +Base.:(-)(p1::Polynomial{<:LinearAlgebra.UniformScaling}, p2::Polynomial) = mapcoefficientsnz(J -> J.λ, p1) - p2 |
| 104 | +Base.:(-)(p1::Polynomial{<:LinearAlgebra.UniformScaling}, p2::Polynomial{<:LinearAlgebra.UniformScaling}) = mapcoefficientsnz(J -> J.λ, p1) - p2 |
| 105 | + |
| 106 | +LinearAlgebra.adjoint(x::Polynomial) = polynomial!(adjoint.(terms(x))) |
| 107 | + |
| 108 | +function mapcoefficients(f::Function, p::Polynomial; nonzero = false) |
| 109 | + terms = map(p.terms) do term |
| 110 | + mapcoefficients(f, term) |
| 111 | + end |
| 112 | + if !nonzero |
| 113 | + filter!(!iszero, terms) |
| 114 | + end |
| 115 | + return polynomial!(terms) |
| 116 | +end |
| 117 | +function mapcoefficients!(f::Function, p::Polynomial; nonzero = false) |
| 118 | + for i in eachindex(p.terms) |
| 119 | + t = p.terms[i] |
| 120 | + p.terms[i] = Term(f(coefficient(t)), monomial(t)) |
| 121 | + end |
| 122 | + if !nonzero |
| 123 | + filter!(!iszero, p.terms) |
| 124 | + end |
| 125 | + return p |
| 126 | +end |
| 127 | + |
| 128 | +function mapcoefficients_to!(output::Polynomial, f::Function, p::Polynomial; nonzero = false) |
| 129 | + resize!(output.terms, nterms(p)) |
| 130 | + for i in eachindex(p.terms) |
| 131 | + t = p.terms[i] |
| 132 | + output.terms[i] = Term(f(coefficient(t)), monomial(t)) |
| 133 | + end |
| 134 | + if !nonzero |
| 135 | + filter!(!iszero, output.terms) |
| 136 | + end |
| 137 | + return output |
| 138 | +end |
| 139 | +function mapcoefficients_to!(output::Polynomial, f::Function, p::AbstractPolynomialLike; nonzero = false) |
| 140 | + return mapcoefficients_to!(output, f, polynomial(p); nonzero = false) |
| 141 | +end |
| 142 | + |
| 143 | +# The polynomials can be mutated. |
| 144 | +MA.mutability(::Type{<:VectorPolynomial}) = MA.IsMutable() |
| 145 | + |
| 146 | +# Terms are not mutable under the MutableArithmetics API |
| 147 | +function MA.mutable_copy(p::VectorPolynomial{C,TT}) where {C,TT} |
| 148 | + return VectorPolynomial{C,TT}([TT(MA.copy_if_mutable(coefficient(term)), monomial(term)) for term in terms(p)]) |
| 149 | +end |
| 150 | +Base.copy(p::VectorPolynomial) = MA.mutable_copy(p) |
| 151 | + |
| 152 | +function grlex end |
| 153 | +function MA.operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{T,TT}, q::Polynomial) where {T,TT} |
| 154 | + get1(i) = p.terms[i] |
| 155 | + function get2(i) |
| 156 | + t = q.terms[i] |
| 157 | + TT(MA.scaling_convert(T, MA.operate(op, coefficient(t))), monomial(t)) |
| 158 | + end |
| 159 | + set(i, t::TT) = p.terms[i] = t |
| 160 | + push(t::TT) = push!(p.terms, t) |
| 161 | + compare_monomials(t::TT, j::Int) = grlex(monomial(q.terms[j]), monomial(t)) |
| 162 | + compare_monomials(i::Int, j::Int) = compare_monomials(get1(i), j) |
| 163 | + combine(i::Int, j::Int) = p.terms[i] = Term(MA.operate!!(op, coefficient(p.terms[i]), coefficient(q.terms[j])), monomial(p.terms[i])) |
| 164 | + combine(t::TT, j::Int) = TT(MA.operate!!(op, coefficient(t), coefficient(q.terms[j])), monomial(t)) |
| 165 | + resize(n) = resize!(p.terms, n) |
| 166 | + # We can modify the coefficient since it's the result of `combine`. |
| 167 | + keep(t::Term) = !MA.iszero!!(coefficient(t)) |
| 168 | + keep(i::Int) = !MA.iszero!!(coefficient(p.terms[i])) |
| 169 | + polynomial_merge!( |
| 170 | + nterms(p), nterms(q), get1, get2, set, push, |
| 171 | + compare_monomials, combine, keep, resize |
| 172 | + ) |
| 173 | + return p |
| 174 | +end |
| 175 | +function MA.operate_to!(output::Polynomial, ::typeof(*), p::Polynomial, q::Polynomial) |
| 176 | + empty!(output.terms) |
| 177 | + mul_to_terms!(output.terms, p, q) |
| 178 | + sort!(output.terms, lt=(>)) |
| 179 | + uniqterms!(output.terms) |
| 180 | + return output |
| 181 | +end |
| 182 | +function MA.operate!(::typeof(*), p::Polynomial, q::Polynomial) |
| 183 | + return MA.operate_to!(p, *, MA.mutable_copy(p), q) |
| 184 | +end |
| 185 | + |
| 186 | +function MA.operate!(::typeof(zero), p::Polynomial) |
| 187 | + empty!(p.terms) |
| 188 | + return p |
| 189 | +end |
| 190 | +function MA.operate!(::typeof(one), p::Polynomial{T}) where T |
| 191 | + if isempty(p.terms) |
| 192 | + push!(p.terms, constantterm(one(T), p)) |
| 193 | + else |
| 194 | + t = p.terms[1] |
| 195 | + p.terms[1] = Term(MA.one!!(coefficient(t)), constantmonomial(t)) |
| 196 | + resize!(p.terms, 1) |
| 197 | + end |
| 198 | + return p |
| 199 | +end |
| 200 | + |
| 201 | +function MA.operate!(::typeof(removeleadingterm), p::Polynomial) |
| 202 | + deleteat!(p.terms, 1) |
| 203 | + return p |
| 204 | +end |
0 commit comments