Report index
StatProfilerHTML.jl report
Generated on Thu, 26 Mar 2020 19:09:01
File source code
Line Exclusive Inclusive Code
1 export polynomial, polynomial!, polynomialtype, terms, nterms, coefficients, monomials
2 export coefficienttype, monomialtype
3 export mindegree, maxdegree, extdegree, effective_variables
4 export leadingterm, leadingcoefficient, leadingmonomial
5 export removeleadingterm, removemonomials, monic
6
7 LinearAlgebra.norm(p::AbstractPolynomialLike, r::Int=2) = LinearAlgebra.norm(coefficients(p), r)
8
9 changecoefficienttype(::Type{TT}, ::Type{T}) where {TT<:AbstractTermLike, T} = termtype(TT, T)
10 changecoefficienttype(::Type{PT}, ::Type{T}) where {PT<:AbstractPolynomial, T} = polynomialtype(PT, T)
11
12 changecoefficienttype(p::PT, ::Type{T}) where {PT<:APL, T} = convert(changecoefficienttype(PT, T), p)
13
14 abstract type ListState end
15 abstract type UnsortedState <: ListState end
16 struct MessyState <: UnsortedState end
17 # No duplicates or zeros
18 struct UniqState <: UnsortedState end
19 sortstate(::MessyState) = SortedState()
20 sortstate(::UniqState) = SortedUniqState()
21 struct SortedState <: ListState end
22 struct SortedUniqState <: ListState end
23
24 """
25 polynomial(p::AbstractPolynomialLike)
26
27 Converts `p` to a value with polynomial type.
28
29 polynomial(p::AbstractPolynomialLike, ::Type{T}) where T
30
31 Converts `p` to a value with polynomial type with coefficient type `T`.
32
33 polynomial(a::AbstractVector, mv::AbstractVector{<:AbstractMonomialLike})
34
35 Creates a polynomial equal to `dot(a, mv)`.
36
37 polynomial(terms::AbstractVector{<:AbstractTerm}, s::ListState=MessyState())
38
39 Creates a polynomial equal to `sum(terms)` where `terms` are guaranteed to be in state `s`.
40
41 polynomial(f::Function, mv::AbstractVector{<:AbstractMonomialLike})
42
43 Creates a polynomial equal to `sum(f(i) * mv[i] for i in 1:length(mv))`.
44
45 ### Examples
46
47 Calling `polynomial([2, 4, 1], [x, x^2*y, x*y])` should return ``4x^2y + xy + 2x``.
48 """
49 polynomial(p::AbstractPolynomial) = p
50 polynomial(p::APL{T}, ::Type{T}) where T = polynomial(terms(p))
51 polynomial(p::APL{T}) where T = polynomial(p, T)
52 function polynomial(Q::AbstractMatrix, mv::AbstractVector)
53 LinearAlgebra.dot(mv, Q * mv)
54 end
55 function polynomial(Q::AbstractMatrix, mv::AbstractVector, ::Type{T}) where T
56 polynomial(polynomial(Q, mv), T)
57 end
58
59 polynomial(f::Function, mv::AbstractVector{<:AbstractMonomialLike}) = polynomial!([f(i) * mv[i] for i in 1:length(mv)])
60
61 polynomial(a::AbstractVector, x::AbstractVector, s::ListState=MessyState()) = polynomial([α * m for (α, m) in zip(a, x)], s)
62
63 polynomial(ts::AbstractVector, s::ListState=MessyState()) = sum(ts)
64 polynomial!(ts::AbstractVector, s::ListState=MessyState()) = sum(ts)
65
66 polynomial!(ts::AbstractVector{<:AbstractTerm}, s::SortedUniqState) = polynomial(coefficient.(ts), monomial.(ts), s)
67
68 function polynomial!(ts::AbstractVector{<:AbstractTerm}, s::SortedState)
69 92 (11 %)
92 (11 %) samples spent in polynomial!
92 (100 %) (incl.) when called from polynomial! line 72
13 (14 %) samples spent calling uniqterms!
79 (86 %) samples spent calling polynomial!
polynomial!(uniqterms!(ts), SortedUniqState())
70 end
71 function polynomial!(ts::AbstractVector{<:AbstractTerm}, s::UnsortedState=MessyState())
72 1388 (166 %)
1388 (166 %) samples spent in polynomial!
694 (50 %) (incl.) when called from polynomial! line 72
694 (50 %) (incl.) when called from * line 188
602 (43 %) samples spent calling sort!##kw
92 (7 %) samples spent calling polynomial!
694 (50 %) samples spent calling polynomial!
polynomial!(sort!(ts, lt=(>)), sortstate(s))
73 end
74
75 _collect(v::Vector) = v
76 _collect(v::AbstractVector) = collect(v)
77 polynomial(ts::AbstractVector{<:AbstractTerm}, args::Vararg{ListState, N}) where {N} = polynomial!(MA.mutable_copy(_collect(ts)), args...)
78
79 """
80 polynomialtype(p::AbstractPolynomialLike)
81
82 Returns the type that `p` would have if it was converted into a polynomial.
83
84 polynomialtype(::Type{PT}) where PT<:AbstractPolynomialLike
85
86 Returns the same as `polynomialtype(::PT)`.
87
88 polynomialtype(p::AbstractPolynomialLike, ::Type{T}) where T
89
90 Returns the type that `p` would have if it was converted into a polynomial of coefficient type `T`.
91
92 polynomialtype(::Type{PT}, ::Type{T}) where {PT<:AbstractPolynomialLike, T}
93
94 Returns the same as `polynomialtype(::PT, ::Type{T})`.
95 """
96 polynomialtype(::Union{P, Type{P}}) where P <: APL = Base.promote_op(polynomial, P)
97 polynomialtype(::Union{P, Type{P}}) where P <: AbstractPolynomial = P
98 polynomialtype(::Union{M, Type{M}}) where M<:AbstractMonomialLike = polynomialtype(termtype(M))
99 polynomialtype(::Union{M, Type{M}}, ::Type{T}) where {M<:AbstractMonomialLike, T} = polynomialtype(termtype(M, T))
100 polynomialtype(::Union{P, Type{P}}, ::Type{T}) where {P <: APL, T} = polynomialtype(polynomialtype(P), T)
101 polynomialtype(::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}) where PT <: APL = polynomialtype(PT)
102 polynomialtype(::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}, ::Type{T}) where {PT <: APL, T} = polynomialtype(PT, T)
103
104
13 (2 %) samples spent in uniqterms!
13 (100 %) (incl.) when called from polynomial! line 69
function uniqterms!(ts::AbstractVector{<: AbstractTerm})
105 i = firstindex(ts)
106 for j in Iterators.drop(eachindex(ts), 1)
107 4 (0 %)
4 (100 %) samples spent calling iszero
if !iszero(ts[j])
108 7 (1 %)
1 (14 %) samples spent calling ==
6 (86 %) samples spent calling getindex
if monomial(ts[i]) == monomial(ts[j])
109 1 (0 %)
1 (100 %) samples spent calling setindex!
ts[i] = MA.add!(coefficient(ts[i]), coefficient(ts[j])) * monomial(ts[i])
110 else
111 if !iszero(ts[i])
112 i += 1
113 end
114 1 (0 %)
1 (100 %) samples spent calling setindex!
ts[i] = MA.copy_if_mutable(ts[j])
115 end
116 end
117 end
118 if i < length(ts)
119 if iszero(ts[i])
120 i -= 1
121 end
122 resize!(ts, i)
123 end
124 ts
125 end
126
127 """
128 terms(p::AbstractPolynomialLike)
129
130 Returns an iterator over the nonzero terms of the polynomial `p` sorted in the decreasing monomial order.
131
132 ### Examples
133
134 Calling `terms` on ``4x^2y + xy + 2x`` should return an iterator of ``[4x^2y, xy, 2x]``.
135 """
136 terms(t::AbstractTermLike) = iszero(t) ? termtype(t)[] : [term(t)]
137 terms(p::AbstractPolynomialLike) = terms(polynomial(p))
138
139 """
140 nterms(p::AbstractPolynomialLike)
141
142 Returns the number of nonzero terms in `p`, i.e. `length(terms(p))`.
143
144 ### Examples
145
146 Calling `nterms` on ``4x^2y + xy + 2x`` should return 3.
147 """
148 function nterms end
149
150 nterms(t::AbstractTermLike) = iszero(t) ? 0 : 1
151 nterms(p::AbstractPolynomialLike) = length(terms(p))
152
153 """
154 coefficients(p::AbstractPolynomialLike)
155
156 Returns an iterator over the coefficients of `p` of the nonzero terms of the polynomial sorted in the decreasing monomial order.
157
158 coefficients(p::AbstractPolynomialLike, X::AbstractVector)
159
160 Returns an iterator over the coefficients of the monomials of `X` in `p` where `X` is a monomial vector not necessarily sorted but with no duplicate entry.
161
162 ### Examples
163
164 Calling `coefficients` on ``4x^2y + xy + 2x`` should return an iterator of ``[4, 1, 2]``.
165 Calling `coefficients(4x^2*y + x*y + 2x + 3, [x, 1, x*y, y])` should return an iterator of ``[2, 3, 1, 0]``.
166 """
167 coefficients(p::APL) = coefficient.(terms(p))
168 function coefficients(p::APL{T}, X::AbstractVector) where T
169 σ, mv = sortmonovec(X)
170 @assert length(mv) == length(X) # no duplicate in X
171 c = zeros(T, length(mv))
172 i = 1
173 for t in terms(p)
174 m = monomial(t)
175 while i <= length(mv) && mv[i] > m
176 c[σ[i]] = zero(T)
177 i += 1
178 end
179 if i <= length(mv) && mv[i] == m
180 c[σ[i]] = coefficient(t)
181 i += 1
182 end
183 end
184 c
185 end
186
187 """
188 monomials(p::AbstractPolynomialLike)
189
190 Returns an iterator over the monomials of `p` of the nonzero terms of the polynomial sorted in the decreasing order.
191
192 monomials(vars::Tuple, degs::AbstractVector{Int}, filter::Function = m -> true)
193
194 Builds the vector of all the monovec `m` with variables `vars` such that the degree `degree(m)` is in `degs` and `filter(m)` is `true`.
195
196 ### Examples
197
198 Calling `monomials` on ``4x^2y + xy + 2x`` should return an iterator of ``[x^2y, xy, x]``.
199
200 Calling `monomials((x, y), [1, 3], m -> degree(m, y) != 1)` should return `[x^3, x*y^2, y^3, x]` where `x^2*y` and `y` have been excluded by the filter.
201 """
202 monomials(p::APL) = monovec(monomial.(terms(p)))
203
204 #$(SIGNATURES)
205 """
206 mindegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})
207
208 Returns the minimal total degree of the monomials of `p`, i.e. `minimum(degree, terms(p))`.
209
210 mindegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)
211
212 Returns the minimal degree of the monomials of `p` in the variable `v`, i.e. `minimum(degree.(terms(p), v))`.
213
214 ### Examples
215 Calling `mindegree` on on ``4x^2y + xy + 2x`` should return 1, `mindegree(4x^2y + xy + 2x, x)` should return 1 and `mindegree(4x^2y + xy + 2x, y)` should return 0.
216 """
217 function mindegree(X::AbstractVector{<:AbstractTermLike}, args...)
218 isempty(X) ? 0 : minimum(t -> degree(t, args...), X)
219 end
220 function mindegree(p::AbstractPolynomialLike, args...)
221 mindegree(terms(p), args...)
222 end
223
224 #$(SIGNATURES)
225 """
226 maxdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})
227
228 Returns the maximal total degree of the monomials of `p`, i.e. `maximum(degree, terms(p))`.
229
230 maxdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)
231
232 Returns the maximal degree of the monomials of `p` in the variable `v`, i.e. `maximum(degree.(terms(p), v))`.
233
234 ### Examples
235 Calling `maxdegree` on on ``4x^2y + xy + 2x`` should return 3, `maxdegree(4x^2y + xy + 2x, x)` should return 2 and `maxdegree(4x^2y + xy + 2x, y)` should return 1.
236 """
237 function maxdegree(X::AbstractVector{<:AbstractTermLike}, args...)
238 isempty(X) ? 0 : maximum(t -> degree(t, args...), X)
239 end
240 function maxdegree(p::AbstractPolynomialLike, args...)
241 maxdegree(terms(p), args...)
242 end
243
244 #$(SIGNATURES)
245 """
246 extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})
247
248 Returns the extremal total degrees of the monomials of `p`, i.e. `(mindegree(p), maxdegree(p))`.
249
250 extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)
251
252 Returns the extremal degrees of the monomials of `p` in the variable `v`, i.e. `(mindegree(p, v), maxdegree(p, v))`.
253
254 ### Examples
255 Calling `extdegree` on on ``4x^2y + xy + 2x`` should return `(1, 3)`, `extdegree(4x^2y + xy + 2x, x)` should return `(1, 2)` and `maxdegree(4x^2y + xy + 2x, y)` should return `(0, 1)`.
256 """
257 function extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, args...)
258 (mindegree(p, args...), maxdegree(p, args...))
259 end
260
261 """
262 effective_variables(p::AbstractPolynomialLike)
263
264 Return a vector of `eltype` `variable_union_type(p)` (see [`variable_union_type`](@ref)),
265 containing all the variables that has nonzero degree in at least one term.
266 That is, return all the variables `v` such that `maxdegree(p, v)` is not zero.
267 """
268 function effective_variables(p::AbstractPolynomialLike)
269 VT = variable_union_type(p)
270 return VT[v for v in variables(p) if !iszero(maxdegree(p, v))]
271 end
272
273 """
274 leadingterm(p::AbstractPolynomialLike)
275
276 Returns the coefficient of the leading term, i.e. `first(terms(p))`.
277
278 ### Examples
279
280 Calling `leadingterm` on ``4x^2y + xy + 2x`` should return ``4x^2y``.
281 """
282 function leadingterm(p::AbstractPolynomialLike)
283 if iszero(p)
284 zeroterm(p)
285 else
286 first(terms(p))
287 end
288 end
289 leadingterm(t::AbstractTermLike) = term(t)
290
291 #$(SIGNATURES)
292 """
293 leadingcoefficient(p::AbstractPolynomialLike)
294
295 Returns the coefficient of the leading term of `p`, i.e. `coefficient(leadingterm(p))`.
296
297 ### Examples
298
299 Calling `leadingcoefficient` on ``4x^2y + xy + 2x`` should return ``4`` and calling it on ``0`` should return ``0``.
300 """
301 function leadingcoefficient(p::AbstractPolynomialLike)
302 coefficient(leadingterm(p))
303 end
304
305 #$(SIGNATURES)
306 """
307 leadingmonomial(p::AbstractPolynomialLike)
308
309 Returns the monomial of the leading term of `p`, i.e. `monomial(leadingterm(p))` or `first(monomials(p))`.
310
311 ### Examples
312
313 Calling `leadingmonomial` on ``4x^2y + xy + 2x`` should return ``x^2y``.
314 """
315 function leadingmonomial(p::AbstractPolynomialLike)
316 # first(monomials(p)) would be more efficient for DynamicPolynomials but
317 # monomial(leadingterm(p)) is more efficient for TypedPolynomials and is better if p is a term
318 monomial(leadingterm(p))
319 end
320
321 #$(SIGNATURES)
322 """
323 removeleadingterm(p::AbstractPolynomialLike)
324
325 Returns a polynomial with the leading term removed in the polynomial `p`.
326
327 ### Examples
328
329 Calling `removeleadingterm` on ``4x^2y + xy + 2x`` should return ``xy + 2x``.
330 """
331 function removeleadingterm(p::AbstractPolynomialLike)
332 # Iterators.drop returns an Interators.Drop which is not an AbstractVector
333 polynomial(terms(p)[2:end], SortedUniqState())
334 end
335
336 #$(SIGNATURES)
337 """
338
339 Returns a polynomial with the terms having their monomial in the monomial vector `mv` removed in the polynomial `p`.
340
341 ### Examples
342
343 Calling `removemonomials(4x^2*y + x*y + 2x, [x*y])` should return ``4x^2*y + 2x``.
344 """
345 function removemonomials(p::AbstractPolynomialLike, mv::AbstractVector{MT}) where {MT <: AbstractMonomialLike}
346 smv = monovec(mv) # Make sure it is sorted
347 i = 1
348 q = zero(p)
349 for t in terms(p)
350 m = monomial(t)
351 while i <= length(smv) && smv[i] > m
352 i += 1
353 end
354 if i > length(smv) || smv[i] != m
355 q += t
356 end
357 end
358 q
359 end
360
361 """
362 monic(p::AbstractPolynomialLike)
363
364 Returns `p / leadingcoefficient(p)` where the leading coefficient of the returned polynomials is made sure to be exactly one to avoid rounding error.
365 """
366 function monic(p::APL)
367 α = leadingcoefficient(p)
368 polynomial!(_divtoone.(terms(p), α))
369 end
370 monic(m::AbstractMonomialLike) = m
371 monic(t::AbstractTermLike{T}) where T = one(T) * monomial(t)
372
373 function _divtoone(t::AbstractTermLike{T}, α::S) where {T, S}
374 U = Base.promote_op(/, T, S)
375 β = coefficient(t)
376 if β == α
377 one(U) * monomial(t)
378 else
379 (β / α) * monomial(t)
380 end
381 end
382
383 """
384 mapcoefficientsnz(f::Function, p::AbstractPolynomialLike)
385
386 Returns the polynomial obtained by applying `f` to each coefficients where `f` is a function such that `f(x)` is nonzero if `x` is nonzero.
387
388 ### Examples
389
390 Calling `mapcoefficientsnz(α -> α^2, 2x*y + 3x + 1)` should return `4x*y + 9x + 1`.
391 """
392 function mapcoefficientsnz(f::Function, p::AbstractPolynomialLike)
393 # Invariant: p has only nonzero coefficient
394 # therefore f(α) will be nonzero for every coefficient α of p
395 # hence we can use Uniq
396 polynomial!(mapcoefficientsnz.(f, terms(p)), SortedUniqState())
397 end
398 mapcoefficientsnz(f::Function, t::AbstractTermLike) = f(coefficient(t)) * monomial(t)
399
400 Base.round(t::AbstractTermLike; args...) = round(coefficient(t); args...) * monomial(t)
401 function Base.round(p::AbstractPolynomialLike; args...)
402 # round(0.1) is zero so we cannot use SortedUniqState
403 polynomial!(round.(terms(p); args...), SortedState())
404 end
405
406 Base.ndims(::Type{<:AbstractPolynomialLike}) = 0
407 Base.broadcastable(p::AbstractPolynomialLike) = Ref(p)