From 60ba67647ecbe9d9dced8733326f7360b03c724b Mon Sep 17 00:00:00 2001 From: john verzani Date: Wed, 7 Jul 2021 17:40:20 -0400 Subject: [PATCH] [WIP] Ring type (issue #349, #234) (#350) * drop T <: Number * add `conv` alternative for `fastconv` * restrict tests by VERSION * doc updates --- Project.toml | 2 +- docs/src/index.md | 97 ++++++++- src/abstract.jl | 46 ++++- src/common.jl | 100 ++++++---- src/contrib.jl | 10 +- src/polynomials/ChebyshevT.jl | 16 +- src/polynomials/ImmutablePolynomial.jl | 58 +++--- src/polynomials/LaurentPolynomial.jl | 39 ++-- src/polynomials/Polynomial.jl | 22 +-- src/polynomials/SparsePolynomial.jl | 61 +++--- src/polynomials/factored_polynomial.jl | 17 +- src/polynomials/standard-basis.jl | 92 +++++++-- src/show.jl | 5 +- test/Project.toml | 1 + test/StandardBasis.jl | 264 ++++++++++++++++++++++--- 15 files changed, 639 insertions(+), 191 deletions(-) diff --git a/Project.toml b/Project.toml index bb65eb73..7af3aee1 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "Polynomials" uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" license = "MIT" author = "JuliaMath" -version = "2.0.12" +version = "2.0.13" [deps] Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5" diff --git a/docs/src/index.md b/docs/src/index.md index 835e7c01..47616800 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -161,7 +161,7 @@ By design, this is not type-stable; the returned roots may be real or complex. The default `roots` function uses the eigenvalues of the [companion](https://en.wikipedia.org/wiki/Companion_matrix) matrix for -a polynomial. This is an `𝑶(n^3)` operation. +a polynomial. This is an `𝑶(n^3)` operation. For polynomials with `BigFloat` coefficients, the `GenericLinearAlgebra` package can be seamlessly used: @@ -204,7 +204,7 @@ julia> PolynomialRoots.roots(coeffs(p)) 1.0000000000000002 + 0.0im ``` -The roots are always returned as complex numbers. +The roots are always returned as complex numbers. * The @@ -226,7 +226,7 @@ julia> AMRVW.roots(float.(coeffs(p))) 2.9999999999999964 + 0.0im ``` -The roots are returned as complex numbers. +The roots are returned as complex numbers. Both `PolynomialRoots` and `AMRVW` are generic and work with `BigFloat` coefficients, for example. @@ -429,7 +429,7 @@ Fit a polynomial (of degree `deg`) to `x` and `y` using polynomial interpolation using Plots, Polynomials xs = range(0, 10, length=10) ys = @. exp(-xs) -f = fit(xs, ys) # degree = length(xs) - 1 +f = fit(xs, ys) # degree = length(xs) - 1 f2 = fit(xs, ys, 2) # degree = 2 scatter(xs, ys, markerstrokewidth=0, label="Data") @@ -583,7 +583,7 @@ Polynomial(4 + 2*x + 3*x^2) If `q` is non-constant, such as `variable(Polynomial, :y)`, then there would be an error due to the mismatched symbols. (The mathematical result would need a multivariate polynomial, not a univariate polynomial, as this package provides.) -The same conversion is done for polynomial multiplication: constant polynomials are treated as numbers; non-constant polynomials must have their symbols match. +The same conversion is done for polynomial multiplication: constant polynomials are treated as numbers; non-constant polynomials must have their symbols match. There is an oddity -- though the following two computations look the same, they are technically different: @@ -595,7 +595,7 @@ julia> one(Polynomial, :y) + one(Polynomial, :x) Polynomial(2.0) ``` -Both are constant polynomials over `Int`, but the first has the indeterminate `:y`, the second `:x`. +Both are constant polynomials over `Int`, but the first has the indeterminate `:y`, the second `:x`. This technical difference causes no issues with polynomial addition or multiplication, as there constant polynomials are treated as numbers, but can be an issue when constant polynomials are used as array elements. @@ -669,6 +669,91 @@ julia> [1 p; p 1] + [1 2one(q); 3 4] # array{P{T,:x}} + array{P{T,:y}} Though were a non-constant polynomial with indeterminate `y` replacing `2one(q)` above, that addition would throw an error. + +## Non-number types for `T` + +The coefficients of the polynomial may be non-number types, such as matrices or other polynomials, albeit not every operation is fully supported. + +For example, a polynomial with matrix coefficients, might be constructed with: + +```jldoctest non_number +julia> using Polynomials + +julia> a,b,c = [1 0;2 1], [1 0; 3 1], [1 0; 4 1] +([1 0; 2 1], [1 0; 3 1], [1 0; 4 1]) + +julia> p = Polynomial([a,b,c]) +Polynomial([1 0; 2 1] + [1 0; 3 1]*x + [1 0; 4 1]*x^2) + +julia> q = derivative(p) +Polynomial([1 0; 3 1] + [2 0; 8 2]*x) +``` + +Various operations are available, `derivative` was shown above, here are the vector-space operations: + +```jldoctest non_number +julia> 2p +Polynomial([2 0; 4 2] + [2 0; 6 2]*x + [2 0; 8 2]*x^2) + +julia> p + q +Polynomial([2 0; 5 2] + [3 0; 11 3]*x + [1 0; 4 1]*x^2) +``` + +polynomial multiplication: + +```jldoctest non_number +julia> p * q +Polynomial([1 0; 5 1] + [3 0; 18 3]*x + [3 0; 21 3]*x^2 + [2 0; 16 2]*x^3) +``` + +polynomial evaluation, here either with a scalar or a matrix: + +```jldoctest non_number +julia> p(2) +2×2 Matrix{Int64}: + 7 0 + 24 7 + +julia> p(b) +2×2 Matrix{Int64}: + 3 0 + 18 3 +``` + +But if the type `T` lacks support of some generic functions, such as `zero(T)` and `one(T)`, then there may be issues. For example, when `T <: AbstractMatrix` the output of `p-p` is an error, as the implementation assumes `zero(T)` is defined. For static arrays, this isn't an issue, as there is support for `zero(T)`. Other polynomial types, such as `SparsePolynomial` have less support, as some specialized methods assume more of the generic interface be implemented. + +Similarly, using polynomials for `T` is a possibility: + +```jldoctest non_number +julia> a,b,c = Polynomial([1],:y), Polynomial([0,1],:y), Polynomial([0,0,1],:y) +(Polynomial(1), Polynomial(y), Polynomial(y^2)) + +julia> p = Polynomial([a,b,c], :x) +Polynomial(Polynomial(1) + Polynomial(y)*x + Polynomial(y^2)*x^2) + +julia> q = derivative(p) +Polynomial(Polynomial(y) + Polynomial(2*y^2)*x) +``` + +Again, much works: + +```jldoctest non_number +julia> 2p +Polynomial(Polynomial(2) + Polynomial(2*y)*x + Polynomial(2*y^2)*x^2) + +julia> p + q +Polynomial(Polynomial(1 + y) + Polynomial(y + 2*y^2)*x + Polynomial(y^2)*x^2) + +julia> p(2) +Polynomial(1 + 2*y + 4*y^2) + +julia> p(b) +Polynomial(1 + y^2 + y^4) +``` + +But much doesn't. For example, implicit promotion can fail. For example, the scalar multiplication `p * b` will fail, as the methods assume this is the fallback polynomial multiplication and not the intended scalar multiplication. + + ## Rational functions The package provides support for rational functions -- fractions of polynomials (for most types). The construction of the basic type mirrors the construction of rational numbers. diff --git a/src/abstract.jl b/src/abstract.jl index 1cd2aec2..a32c8ce4 100644 --- a/src/abstract.jl +++ b/src/abstract.jl @@ -1,17 +1,50 @@ export AbstractPolynomial + + const SymbolLike = Union{AbstractString,Char,Symbol} """ AbstractPolynomial{T,X} -An abstract container for various polynomials. +An abstract type for various polynomials. + +A polynomial type holds an indeterminate `X`; coefficients of type `T`, stored in some container type; and an implicit basis, such as the standard basis. # Properties - `coeffs` - The coefficients of the polynomial + +# The type `T` + +`T` need not be `T <: Number`, at the moment it is not restricted + +Some `T`s will not be successful + +* scalar mult: `c::Number * p::Polynomial` should be defined +* scalar mult: `c::T * p:Polynomial{T}` An ambiguity when `T <: AbstractPolynomial` +* scalar mult: `p:Polynomial{T} * c::T` need not commute + +* scalar add/sub: `p::Polynomial{T} + q::Polynomial{T}` should be defined +* scalar sub: `p::Polynomial{T} - p::Polynomial{T}` generally needs `zeros(T,1)` defined for `zero(Polynomial{T})` + +* poly mult: `p::Polynomial{T} * q::Polynomial{T}` Needs "`T * T`" defined (e.g. `Base.promote_op(*, Vector{Int}, Vector{Int}))` needs to be something.) +* poly powers: `p::Polynomial{T}^2` needs "`T^2`" defined + +* implicit promotion: `p::Polynomial{T} + c::Number` needs `convert(T, c)` defined +* implicit promotion: `p::Polynomial{T} + c::T` ambiguity if `T <: AbstractPolynomial` + +* evaluation: `p::Polynomial{T}(s::Number)` +* evaluation `p::Polynomial{T}(c::T)` needs `T*T` defined + +* derivatives: `derivative(p::Polynomial{T})` +* integrals: `integrate(p::Polynomial{T})` + + """ abstract type AbstractPolynomial{T,X} end + +## ----- # We want ⟒(P{α…,T}) = P{α…}; this default # works for most cases ⟒(P::Type{<:AbstractPolynomial}) = constructorof(P) @@ -34,7 +67,7 @@ Polynomials.@register MyPolynomial ``` # Implementations -This will implement simple self-conversions like `convert(::Type{MyPoly}, p::MyPoly) = p` and creates two promote rules. The first allows promotion between two types (e.g. `promote(Polynomial, ChebyshevT)`) and the second allows promotion between parametrized types (e.g. `promote(Polynomial{T}, Polynomial{S})`). +This will implement simple self-conversions like `convert(::Type{MyPoly}, p::MyPoly) = p` and creates two promote rules. The first allows promotion between two types (e.g. `promote(Polynomial, ChebyshevT)`) and the second allows promotion between parametrized types (e.g. `promote(Polynomial{T}, Polynomial{S})`). For constructors, it implements the shortcut for `MyPoly(...) = MyPoly{T}(...)`, singleton constructor `MyPoly(x::Number, ...)`, conversion constructor `MyPoly{T}(n::S, ...)`, and `variable` alternative `MyPoly(var=:x)`. """ @@ -44,7 +77,7 @@ macro register(name) Base.convert(::Type{P}, p::P) where {P<:$poly} = p function Base.convert(P::Type{<:$poly}, p::$poly{T,X}) where {T,X} isconstant(p) && return constructorof(P){eltype(P),indeterminate(P)}(constantterm(p)) - constructorof(P){eltype(P), indeterminate(P,p)}(coeffs(p)) + constructorof(P){eltype(P), indeterminate(P,p)}(_coeffs(p)) end Base.promote(p::P, q::Q) where {X, T, P <:$poly{T,X}, Q <: $poly{T,X}} = p,q Base.promote_rule(::Type{<:$poly{T,X}}, ::Type{<:$poly{S,X}}) where {T,S,X} = $poly{promote_type(T, S),X} @@ -52,7 +85,7 @@ macro register(name) $poly{promote_type(T, S),X} $poly(coeffs::AbstractVector{T}, var::SymbolLike = :x) where {T} = $poly{T, Symbol(var)}(coeffs) - $poly{T}(x::AbstractVector{S}, var::SymbolLike = :x) where {T,S<:Number} = + $poly{T}(x::AbstractVector{S}, var::SymbolLike = :x) where {T,S} = $poly{T,Symbol(var)}(T.(x)) function $poly(coeffs::G, var::SymbolLike=:x) where {G} !Base.isiterable(G) && throw(ArgumentError("coeffs is not iterable")) @@ -76,11 +109,11 @@ macro registerN(name, params...) αs = tuple(esc.(params)...) quote Base.convert(::Type{P}, q::Q) where {$(αs...),T, P<:$poly{$(αs...),T}, Q <: $poly{$(αs...),T}} = q - Base.convert(::Type{$poly{$(αs...)}}, q::Q) where {$(αs...),T, Q <: $poly{$(αs...),T}} = q + Base.convert(::Type{$poly{$(αs...)}}, q::Q) where {$(αs...),T, Q <: $poly{$(αs...),T}} = q Base.promote(p::P, q::Q) where {$(αs...),T, X, P <:$poly{$(αs...),T,X}, Q <: $poly{$(αs...),T,X}} = p,q Base.promote_rule(::Type{<:$poly{$(αs...),T,X}}, ::Type{<:$poly{$(αs...),S,X}}) where {$(αs...),T,S,X} = $poly{$(αs...),promote_type(T, S),X} - Base.promote_rule(::Type{<:$poly{$(αs...),T,X}}, ::Type{S}) where {$(αs...),T,X,S<:Number} = + Base.promote_rule(::Type{<:$poly{$(αs...),T,X}}, ::Type{S}) where {$(αs...),T,X,S<:Number} = $poly{$(αs...),promote_type(T,S),X} function $poly{$(αs...),T}(x::AbstractVector{S}, var::SymbolLike = :x) where {$(αs...),T,S} @@ -98,4 +131,3 @@ macro registerN(name, params...) (p::$poly)(x) = evalpoly(x, p) end end - diff --git a/src/common.jl b/src/common.jl index db9d23bd..c473a88a 100644 --- a/src/common.jl +++ b/src/common.jl @@ -81,7 +81,7 @@ is [`Polynomial`](@ref). ## Weights -Weights may be assigned to the points by specifying a vector or matrix of weights. +Weights may be assigned to the points by specifying a vector or matrix of weights. When specified as a vector, `[w₁,…,wₙ]`, the weights should be non-negative as the minimization problem is `argmin_β Σᵢ wᵢ |yᵢ - Σⱼ @@ -89,9 +89,9 @@ Vᵢⱼ βⱼ|² = argmin_β || √(W)⋅(y - V(x)β)||²`, where, `W` the digon matrix formed from `[w₁,…,wₙ]`, is used for the solution, `V` being the Vandermonde matrix of `x` corresponding to the specified degree. This parameterization of the weights is different from that of -`numpy.polyfit`, where the weights would be specified through -`[ω₁,ω₂,…,ωₙ] = [√w₁, √w₂,…,√wₙ]` -with the answer solving +`numpy.polyfit`, where the weights would be specified through +`[ω₁,ω₂,…,ωₙ] = [√w₁, √w₂,…,√wₙ]` +with the answer solving `argminᵦ | (ωᵢ⋅yᵢ- ΣⱼVᵢⱼ(ω⋅x) βⱼ) |^2`. When specified as a matrix, `W`, the solution is through the normal @@ -131,8 +131,8 @@ fit′(P::Type{<:AbstractPolynomial}, x::AbstractVector{T}, y::AbstractVector{T}, args...; kwargs...) where {T} = fit(P,x,y,args...; kwargs...) - - + + fit(x::AbstractVector, y::AbstractVector, deg::Integer = length(x) - 1; @@ -168,7 +168,7 @@ _wlstsq(vand, y, W::AbstractMatrix) = qr(vand' * W * vand) \ (vand' * W * y) """ roots(::AbstractPolynomial; kwargs...) -Returns the roots, or zeros, of the given polynomial. +Returns the roots, or zeros, of the given polynomial. This is calculated via the eigenvalues of the companion matrix. The `kwargs` are passed to the `LinearAlgeebra.eigvals` call. @@ -181,7 +181,7 @@ This is calculated via the eigenvalues of the companion matrix. The `kwargs` are """ -function roots(q::AbstractPolynomial{T}; kwargs...) where {T <: Number} +function roots(q::AbstractPolynomial{T}; kwargs...) where {T} p = convert(Polynomial{T}, q) roots(p; kwargs...) @@ -231,7 +231,7 @@ end Compute the definite integral of the given polynomial from `a` to `b`. Will throw an error if either `a` or `b` are out of the polynomial's domain. """ -function integrate(p::AbstractPolynomial, a::Number, b::Number) +function integrate(p::AbstractPolynomial, a, b) P = integrate(p) return P(b) - P(a) end @@ -328,8 +328,8 @@ function chop!(ps::Vector{T}; tol = norm(ps) * rtol + atol for i = lastindex(ps):-1:1 val = ps[i] - if abs(val) > tol - resize!(ps, i); + if abs(val) > tol + resize!(ps, i); return nothing end end @@ -428,8 +428,12 @@ LinearAlgebra.transpose!(p::AbstractPolynomial) = p Conversions =# Base.convert(::Type{P}, p::P) where {P <: AbstractPolynomial} = p Base.convert(P::Type{<:AbstractPolynomial}, x) = P(x) -function Base.convert(::Type{S}, p::P) where {S <: Number,T, P<:Polynomials.AbstractPolynomial{T}} - Polynomials.isconstant(p) && return convert(S, Polynomials.constantterm(p)) +function Base.convert(::Type{S}, p::P) where {S <: Number,T, P<:AbstractPolynomial{T}} + isconstant(p) && return convert(S, constantterm(p)) + throw(ArgumentError("Can't convert a nonconstant polynomial to type $S")) +end +function Base.convert(::Type{T}, p::P) where {T, P<:AbstractPolynomial{T}} + isconstant(p) && return constantterm(p) throw(ArgumentError("Can't convert a nonconstant polynomial to type $S")) end @@ -537,8 +541,9 @@ Determine whether a polynomial is a monic polynomial, i.e., its leading coeffici ismonic(p::AbstractPolynomial) = isone(convert(Polynomial, p)[end]) "`hasnan(p::AbstractPolynomial)` are any coefficients `NaN`" -hasnan(p::AbstractPolynomial) = any(isnan, p) - +hasnan(p::AbstractPolynomial) = any(hasnan, p) +hasnan(p::AbstractArray) = any(hasnan.(p)) +hasnan(x) = isnan(x) """ isconstant(::AbstractPolynomial) @@ -554,6 +559,8 @@ Return the coefficient vector. For a standard basis polynomial these are `[a_0, """ coeffs(p::AbstractPolynomial) = p.coeffs +# hook in for offset coefficients of Laurent Polynomials +_coeffs(p::AbstractPolynomial) = coeffs(p) # specialize this to p[0] when basis vector is 1 @@ -570,7 +577,7 @@ constantterm(p::AbstractPolynomial{T}) where {T} = p(zero(T)) Return the degree of the polynomial, i.e. the highest exponent in the polynomial that has a nonzero coefficient. The degree of the zero polynomial is defined to be -1. The default method assumes the basis polynomial, `βₖ` has degree `k`. """ -degree(p::AbstractPolynomial) = iszero(p) ? -1 : lastindex(p) +degree(p::AbstractPolynomial) = iszero(p) ? -1 : lastindex(p) """ @@ -617,7 +624,7 @@ Base.eachindex(p::AbstractPolynomial) = 0:length(p) - 1 Base.broadcastable(p::AbstractPolynomial) = Ref(p) # getindex -function Base.getindex(p::AbstractPolynomial{T}, idx::Int) where {T <: Number} +function Base.getindex(p::AbstractPolynomial{T}, idx::Int) where {T} m,M = firstindex(p), lastindex(p) idx < m && throw(BoundsError(p, idx)) idx > M && return zero(T) @@ -628,7 +635,7 @@ Base.getindex(p::AbstractPolynomial, indices) = [p[i] for i in indices] Base.getindex(p::AbstractPolynomial, ::Colon) = coeffs(p) # setindex -function Base.setindex!(p::AbstractPolynomial, value::Number, idx::Int) +function Base.setindex!(p::AbstractPolynomial, value, idx::Int) n = length(coeffs(p)) if n ≤ idx resize!(p.coeffs, idx + 1) @@ -638,16 +645,18 @@ function Base.setindex!(p::AbstractPolynomial, value::Number, idx::Int) return p end -Base.setindex!(p::AbstractPolynomial, value::Number, idx::Number) = +Base.setindex!(p::AbstractPolynomial, value, idx::Number) = setindex!(p, value, convert(Int, idx)) -Base.setindex!(p::AbstractPolynomial, value::Number, indices) = - [setindex!(p, value, i) for i in indices] +#Base.setindex!(p::AbstractPolynomial, value::Number, indices) = +# [setindex!(p, value, i) for i in indices] Base.setindex!(p::AbstractPolynomial, values, indices) = - [setindex!(p, v, i) for (v, i) in zip(values, indices)] -Base.setindex!(p::AbstractPolynomial, value::Number, ::Colon) = - setindex!(p, value, eachindex(p)) + [setindex!(p, v, i) for (v, i) in tuple.(values, indices)] +# [setindex!(p, v, i) for (v, i) in zip(values, indices)] +#Base.setindex!(p::AbstractPolynomial, value, ::Colon) = +# setindex!(p, value, eachindex(p)) Base.setindex!(p::AbstractPolynomial, values, ::Colon) = - [setindex!(p, v, i) for (v, i) in zip(values, eachindex(p))] +# [setindex!(p, v, i) for (v, i) in zip(values, eachindex(p))] + [setindex!(p, v, i) for (v, i) in tuple.(values, eachindex(p))] #= Iteration =# @@ -826,11 +835,17 @@ basis(p::P, k::Int, _var=indeterminate(p); var=_var) where {P<:AbstractPolynomia #= arithmetic =# -Base.:-(p::P) where {P <: AbstractPolynomial} = _convert(p, -coeffs(p)) +Base.:-(p::P) where {P <: AbstractPolynomial} = _convert(p, -coeffs(p)) + +Base.:*(p::AbstractPolynomial, c::Number) = scalar_mult(p, c) +Base.:*(c::Number, p::AbstractPolynomial) = scalar_mult(c, p) +Base.:*(c::T, p::P) where {T, P <: AbstractPolynomial{T}} = scalar_mult(c, p) +Base.:*(p::P, c::T) where {T, P <: AbstractPolynomial{T}} = scalar_mult(p, c) + +# implicitly identify c::Number with a constant polynomials Base.:+(c::Number, p::AbstractPolynomial) = +(p, c) Base.:-(p::AbstractPolynomial, c::Number) = +(p, -c) Base.:-(c::Number, p::AbstractPolynomial) = +(-p, c) -Base.:*(c::Number, p::AbstractPolynomial) = *(p, c) # scalar operations # no generic p+c, as polynomial addition falls back to scalar ops @@ -849,7 +864,7 @@ Base.:-(p1::AbstractPolynomial, p2::AbstractPolynomial) = +(p1, -p2) ## +(p::P,c::Number) and +(p::P, q::Q) where {T,S,X,P<:SubtypePolynomial{T,X},Q<:SubtypePolynomial{S,X}} ## though the default for poly+poly isn't terrible -# polynomial + scalar +# polynomial + scalar; implicit identification of c with c*one(P) Base.:+(p::P, c::T) where {T,X, P<:AbstractPolynomial{T,X}} = p + c * one(P) function Base.:+(p::P, c::S) where {T,X, P<:AbstractPolynomial{T,X}, S} @@ -887,7 +902,7 @@ function ⊕(P::Type{<:AbstractPolynomial}, p1::Vector{T}, p2::Vector{S}) where for i in 1:n2 cs[i] += p2[i] end - + return cs end @@ -914,10 +929,10 @@ function ⊕(P::Type{<:AbstractPolynomial}, p1::Dict{Int,T}, p2::Dict{Int,S}) wh R = promote_type(T,S) p = Dict{Int, R}() - + # this allocates in the union -# for i in union(eachindex(p1), eachindex(p2)) +# for i in union(eachindex(p1), eachindex(p2)) # p[i] = p1[i] + p2[i] # end @@ -929,18 +944,28 @@ function ⊕(P::Type{<:AbstractPolynomial}, p1::Dict{Int,T}, p2::Dict{Int,S}) wh @inbounds p[i] = get(p1, i, zero(R)) + pi end end - + return p end ## -- multiplication -## scalar *, / -function Base.:*(p::P, c::S) where {P <: AbstractPolynomial,S} - _convert(p, coeffs(p) .* c) + +# this fall back not necessarily efficient (e.g., sparse) +function scalar_mult(p::P, c::S) where {S, T, X, P<:AbstractPolynomial{T,X}} + R = Base.promote_op(*, T, S) # typeof(one(T)*one(S))? + 𝐏 = ⟒(P){R,X} + 𝐏([pᵢ * c for pᵢ ∈ coeffs(p)]) end +function scalar_mult(c::S, p::P) where {S, T, X, P<:AbstractPolynomial{T, X}} + R = Base.promote_op(*, T, S) + 𝐏 = ⟒(P){R,X} + 𝐏([c * pᵢ for pᵢ ∈ coeffs(p)]) +end + + function Base.:/(p::P, c::S) where {P <: AbstractPolynomial,S} _convert(p, coeffs(p) ./ c) end @@ -955,7 +980,6 @@ function Base.:*(p1::P, p2::Q) where {T,X,P <: AbstractPolynomial{T,X},S,Y,Q <: return p1 * p2 end - Base.:^(p::AbstractPolynomial, n::Integer) = Base.power_by_squaring(p, n) function Base.divrem(num::P, den::O) where {P <: AbstractPolynomial,O <: AbstractPolynomial} @@ -1025,7 +1049,7 @@ Comparisons =# Base.isequal(p1::P, p2::P) where {P <: AbstractPolynomial} = hash(p1) == hash(p2) Base.:(==)(p1::AbstractPolynomial, p2::AbstractPolynomial) = check_same_variable(p1,p2) && (coeffs(p1) == coeffs(p2)) -Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && p[0] == n +Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && constantterm(p) == n Base.:(==)(n::Number, p::AbstractPolynomial) = p == n function Base.isapprox(p1::AbstractPolynomial{T,X}, @@ -1050,7 +1074,7 @@ function Base.isapprox(p1::AbstractPolynomial{T}, n::S; rtol::Real = (Base.rtoldefault(T, S, 0)), atol::Real = 0,) where {T,S} - return isapprox(p1, _convert(p1, [n])) + return isapprox(p1, _convert(p1, [n])) end Base.isapprox(n::S, diff --git a/src/contrib.jl b/src/contrib.jl index a5d20b7c..9897ff3a 100644 --- a/src/contrib.jl +++ b/src/contrib.jl @@ -2,6 +2,7 @@ using Base.Cartesian + # direct version (do not check if threshold is satisfied) @generated function fastconv(E::Array{T,N}, k::Array{T,N}) where {T,N} quote @@ -40,7 +41,7 @@ function evalpoly(x::S, p::Tuple) where {S} N = length(p.parameters) ex = :(p[end]*_one(S)) for i in N-1:-1:1 - ex = :(_muladd(x, $ex, p[$i])) + ex = :(_muladd($ex, x, p[$i])) end ex else @@ -54,7 +55,7 @@ function _evalpoly(x::S, p) where {S} N = length(p) ex = p[end]*_one(x) for i in N-1:-1:1 - ex = _muladd(x, ex, p[i]) + ex = _muladd(ex, x, p[i]) end ex end @@ -110,8 +111,8 @@ end ## modify muladd, as needed _muladd(a,b,c) = muladd(a,b,c) -_muladd(a::Vector, b, c) = a.*b .+ c -_muladd(a::Matrix, b, c) = a*(b*I) + c*I +_muladd(a, b::Vector, c) = a.*b .+ c +_muladd(a, b::Matrix, c) = (a*I)*b + c*I # try to get y = P(c::T)(x::S) = P{T}(c)(x::S) to # have y == one(T)*one(S)*x @@ -129,4 +130,3 @@ end @generated function constructorof(::Type{T}) where T getfield(parentmodule(T), nameof(T)) end - diff --git a/src/polynomials/ChebyshevT.jl b/src/polynomials/ChebyshevT.jl index 1ab386c9..b5a4e511 100644 --- a/src/polynomials/ChebyshevT.jl +++ b/src/polynomials/ChebyshevT.jl @@ -1,7 +1,7 @@ export ChebyshevT """ - ChebyshevT{T<:Number, X}(coeffs::AbstractVector) + ChebyshevT{T, X}(coeffs::AbstractVector) Chebyshev polynomial of the first kind. @@ -9,7 +9,7 @@ Construct a polynomial from its coefficients `coeffs`, lowest order first, optio terms of the given variable `var`, which can be a character, symbol, or string. !!! note - `ChebyshevT` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first + `ChebyshevT` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first index always corresponding to the coefficient of `T_0(x)`. # Examples @@ -39,9 +39,9 @@ The latter shows how to evaluate a `ChebyshevT` polynomial outside of its domain The Chebyshev polynomials are also implemented in `ApproxFun`, `ClassicalOrthogonalPolynomials.jl`, `FastTransforms.jl`, and `SpecialPolynomials.jl`. """ -struct ChebyshevT{T <: Number, X} <: AbstractPolynomial{T, X} +struct ChebyshevT{T, X} <: AbstractPolynomial{T, X} coeffs::Vector{T} - function ChebyshevT{T, X}(coeffs::AbstractVector{S}) where {T <: Number,X, S} + function ChebyshevT{T, X}(coeffs::AbstractVector{S}) where {T, X, S} if Base.has_offset_axes(coeffs) @warn "ignoring the axis offset of the coefficient vector" @@ -264,8 +264,12 @@ end function showterm(io::IO, ::Type{ChebyshevT{T,X}}, pj::T, var, j, first::Bool, mimetype) where {N, T,X} iszero(pj) && return false !first && print(io, " ") - print(io, hasneg(T) && isneg(pj) ? "- " : (!first ? "+ " : "")) - print(io, "$(abs(pj))⋅T_$j($var)") + if hasneg(T) + print(io, isneg(pj) ? "- " : (!first ? "+ " : "")) + print(io, "$(abs(pj))⋅T_$j($var)") + else + print(io, "+ ", "$(pj)⋅T_$j($var)") + end return true end diff --git a/src/polynomials/ImmutablePolynomial.jl b/src/polynomials/ImmutablePolynomial.jl index 5e81deec..1acfe344 100644 --- a/src/polynomials/ImmutablePolynomial.jl +++ b/src/polynomials/ImmutablePolynomial.jl @@ -1,9 +1,9 @@ export ImmutablePolynomial """ - ImmutablePolynomial{T<:Number, X, N}(coeffs::AbstractVector{T}) + ImmutablePolynomial{T, X, N}(coeffs::AbstractVector{T}) -Construct an immutable (static) polynomial from its coefficients +Construct an immutable (static) polynomial from its coefficients `a₀, a₁, …, aₙ`, lowest order first, optionally in terms of the given variable `x` where `x` can be a character, symbol, or string. @@ -29,7 +29,7 @@ immutable polynomials can not promote to a common type. As such, they are precluded from use in rational functions. !!! note - `ImmutablePolynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first + `ImmutablePolynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first index always corresponding to the constant term. # Examples @@ -51,9 +51,9 @@ ImmutablePolynomial(1.0) This was modeled after [StaticUnivariatePolynomials](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) by `@tkoolen`. """ -struct ImmutablePolynomial{T <: Number, X, N} <: StandardBasisPolynomial{T, X} +struct ImmutablePolynomial{T, X, N} <: StandardBasisPolynomial{T, X} coeffs::NTuple{N, T} - function ImmutablePolynomial{T,X,N}(coeffs::NTuple{N,T}) where {T <: Number, X, N} + function ImmutablePolynomial{T,X,N}(coeffs::NTuple{N,T}) where {T, X, N} N == 0 && return new{T,X, 0}(coeffs) iszero(coeffs[end]) && throw(ArgumentError("Leading term must be non-zero")) new{T,X,N}(coeffs) @@ -63,19 +63,29 @@ end @register ImmutablePolynomial ## Various interfaces - ## Abstract Vector coefficients +function ImmutablePolynomial{T,X, N}(coeffs::AbstractVector{T}) where {T,X, N} + cs = NTuple{N,T}(coeffs[i] for i ∈ firstindex(coeffs):N) + ImmutablePolynomial{T, X, N}(cs) + +end + +function ImmutablePolynomial{T,X, N}(coeffs::AbstractVector{S}) where {T,X, N, S} + cs = NTuple{N,T}(coeffs[i] for i ∈ firstindex(coeffs):N) + ImmutablePolynomial{T, X, N}(cs) + +end + function ImmutablePolynomial{T,X}(coeffs::AbstractVector{S}) where {T,X,S} R = promote_type(T,S) - + if Base.has_offset_axes(coeffs) @warn "ignoring the axis offset of the coefficient vector" end N = findlast(!iszero, coeffs) N == nothing && return ImmutablePolynomial{R,X,0}(()) N′ = N + 1 - firstindex(coeffs) - cs = NTuple{N′,T}(coeffs[i] for i ∈ firstindex(coeffs):N) - ImmutablePolynomial{T, X, N′}(cs) + ImmutablePolynomial{T, X, N′}([coeffs[i] for i ∈ firstindex(coeffs):N]) end ## -- Tuple arguments @@ -105,7 +115,7 @@ Base.copy(p::P) where {P <: ImmutablePolynomial} = P(coeffs(p)) degree(p::ImmutablePolynomial{T,X, N}) where {T,X,N} = N - 1 # no trailing zeros isconstant(p::ImmutablePolynomial{T,X,N}) where {T,X,N} = N <= 1 -Base.setindex!(p::ImmutablePolynomial, val::Number, idx::Int) = throw(ArgumentError("ImmutablePolynomials are immutable")) +Base.setindex!(p::ImmutablePolynomial, val, idx::Int) = throw(ArgumentError("ImmutablePolynomials are immutable")) for op in [:isequal, :(==)] @eval function Base.$op(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) where {T,N,S,M} @@ -116,7 +126,7 @@ for op in [:isequal, :(==)] n2 = findlast(!iszero, p2s) (n1 == nothing && n2 == nothing) && return true (n1 == nothing || n2 == nothing) && return false - $op(p1s[1:n1],p2s[1:n2]) && return true + $op(p1s[1:n1],p2s[1:n2]) && return true false end end @@ -154,7 +164,7 @@ function Base.:+(p::P, c::S) where {T, X, N, P <: ImmutablePolynomial{T,X,N}, S< q = ImmutablePolynomial{R,X,N}(cs) return q - + end Base.:-(p::ImmutablePolynomial{T,X,N}) where {T,X,N} = ImmutablePolynomial{T,X,N}(.-p.coeffs) @@ -179,15 +189,21 @@ end ## multiplication -function Base.:*(p::ImmutablePolynomial{T,X,N}, c::S) where {T, X,N, S <: Number} - R = eltype(one(T)*one(S)) - P = ImmutablePolynomial{R,X} - (N == 0 || iszero(c)) && return zero(P) +function scalar_mult(p::ImmutablePolynomial{T,X,N}, c::S) where {T, X,N, S <: Number} + R = eltype(p[0] * c * 0) + (N == 0 || iszero(c)) && return zero(ImmutablePolynomial{R,X}) cs = p.coeffs .* c - iszero(cs[end]) ? P(cs) : P{N}(cs) # more performant to specify when N is known + return ImmutablePolynomial(cs, X) + end + +function scalar_mult(c::S, p::ImmutablePolynomial{T,X,N}) where {T, X,N, S <: Number} + R = eltype(p[0] * c * 0) + (N == 0 || iszero(c)) && return zero(ImmutablePolynomial{R,X}) + cs = p.coeffs .* c + return ImmutablePolynomial(cs, X) end -function Base.:/(p::ImmutablePolynomial{T,X,N}, c::S) where {T,X,N,S <: Number} +function Base.:/(p::ImmutablePolynomial{T,X,N}, c::S) where {T,X,N,S} R = eltype(one(T)/one(S)) P = ImmutablePolynomial{R,X} (N == 0 || isinf(c)) && return zero(P) @@ -222,12 +238,12 @@ LinearAlgebra.norm(q::ImmutablePolynomial, p::Real) = _norm(q.coeffs, p) for j = 2:N expr = :($expr + abs2(a[$j])) end - + return quote $(Expr(:meta, :inline)) @inbounds return sqrt($expr) end - + end _norm_p0(x) = iszero(x) ? zero(x) : one(x) @@ -258,5 +274,3 @@ _norm_p0(x) = iszero(x) ? zero(x) : one(x) end end - - diff --git a/src/polynomials/LaurentPolynomial.jl b/src/polynomials/LaurentPolynomial.jl index 31403516..4527dc6a 100644 --- a/src/polynomials/LaurentPolynomial.jl +++ b/src/polynomials/LaurentPolynomial.jl @@ -5,7 +5,7 @@ export LaurentPolynomial A [Laurent](https://en.wikipedia.org/wiki/Laurent_polynomial) polynomial is of the form `a_{m}x^m + ... + a_{n}x^n` where `m,n` are integers (not necessarily positive) with ` m <= n`. -The `coeffs` specify `a_{m}, a_{m-1}, ..., a_{n}`. +The `coeffs` specify `a_{m}, a_{m-1}, ..., a_{n}`. The argument `m` represents the lowest exponent of the variable in the series, and is taken to be zero by default. Laurent polynomials and standard basis polynomials promote to Laurent polynomials. Laurent polynomials may be converted to a standard basis polynomial when `m >= 0` @@ -69,12 +69,12 @@ julia> x^degree(p) * p(x⁻¹) # reverses coefficients LaurentPolynomial(3.0 + 2.0*x + 1.0*x²) ``` """ -struct LaurentPolynomial{T <: Number, X} <: StandardBasisPolynomial{T, X} +struct LaurentPolynomial{T, X} <: StandardBasisPolynomial{T, X} coeffs::Vector{T} m::Base.RefValue{Int} n::Base.RefValue{Int} function LaurentPolynomial{T,X}(coeffs::AbstractVector{S}, - m::Union{Int, Nothing}=nothing) where {T <: Number, X, S} + m::Union{Int, Nothing}=nothing) where {T, X, S} fnz = findfirst(!iszero, coeffs) fnz == nothing && return new{T,X}(zeros(T,1), Ref(0), Ref(0)) @@ -84,7 +84,7 @@ struct LaurentPolynomial{T <: Number, X} <: StandardBasisPolynomial{T, X} cs = convert(Vector{T}, coeffs[fnz:lnz]) return new{T,X}(cs, Ref(fnz), Ref(lnz)) else - + c = convert(Vector{T}, coeffs[fnz:lnz]) m′ = fnz - 1 + (m == nothing ? 0 : m) @@ -100,16 +100,15 @@ end ## constructors function LaurentPolynomial{T}(coeffs::AbstractVector{S}, m::Int, var::SymbolLike=:x) where { - T <: Number, S <: Number} + T, S <: Number} LaurentPolynomial{T,Symbol(var)}(T.(coeffs), m) end -function LaurentPolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where { - T <: Number} +function LaurentPolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T} LaurentPolynomial{T, Symbol(var)}(coeffs, 0) end -function LaurentPolynomial(coeffs::AbstractVector{T}, m::Int, var::SymbolLike=:x) where {T <: Number} +function LaurentPolynomial(coeffs::AbstractVector{T}, m::Int, var::SymbolLike=:x) where {T} LaurentPolynomial{T, Symbol(var)}(coeffs, m) end @@ -133,21 +132,17 @@ function Base.convert(P::Type{<:Polynomial}, q::LaurentPolynomial) P([q[i] for i in 0:n], indeterminate(q)) end -# need to add p.m[], so abstract.jl method isn't sufficent -function Base.convert(::Type{P}, p::LaurentPolynomial) where {P<:LaurentPolynomial} - S′ = _eltype(P) - Y′ = _indeterminate(P) - S = S′ == nothing ? eltype(p) : S′ - Y = Y′ == nothing ? indeterminate(p) : Y′ - isconstant(p) && return LaurentPolynomial{S,Y}(constantterm(p)) - LaurentPolynomial{S,Y}(p.coeffs, p.m[]) + +# work around for non-applicable convert(::Type{<:P}, p::P{T,X}) in abstract.jl +struct OffsetCoeffs{V} + coeffs::V + m::Int end -# function Base.convert(::Type{P}, q::StandardBasisPolynomial{S}) where {T, P <:LaurentPolynomial{T},S} -# v′ = _indeterminate(P) -# X = v′ == nothing ? indeterminate(q) : v′ -# ⟒(P){T,X}([q[i] for i in 0:degree(q)], 0) -# end +_coeffs(p::LaurentPolynomial) = OffsetCoeffs(p.coeffs, p.m[]) +function LaurentPolynomial{T,X}(p::OffsetCoeffs) where {T, X} + LaurentPolynomial{T,X}(p.coeffs, p.m) +end function Base.convert(::Type{P}, q::StandardBasisPolynomial{S}) where {P <:LaurentPolynomial,S} @@ -528,7 +523,7 @@ function integrate(p::P) where {T, X, P<: LaurentPolynomial{T, X}} !iszero(p[-1]) && throw(ArgumentError("Can't integrate Laurent polynomial with `x⁻¹` term")) R = eltype(one(T)/1) Q = ⟒(P){R, X} - + if hasnan(p) return Q([NaN],0) end diff --git a/src/polynomials/Polynomial.jl b/src/polynomials/Polynomial.jl index 4a5ca6b8..5cb46223 100644 --- a/src/polynomials/Polynomial.jl +++ b/src/polynomials/Polynomial.jl @@ -1,7 +1,7 @@ export Polynomial """ - Polynomial{T<:Number, X}(coeffs::AbstractVector{T}, [var = :x]) + Polynomial{T, X}(coeffs::AbstractVector{T}, [var = :x]) Construct a polynomial from its coefficients `coeffs`, lowest order first, optionally in terms of the given variable `var` which may be a character, symbol, or a string. @@ -14,8 +14,8 @@ with combinations of polynomials and scalars. However, operations involving two polynomials of different variables causes an error except those involving a constant polynomial. !!! note - `Polynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first - index always corresponding to the constant term. In order to use the axis of `coeffs` as exponents, + `Polynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first + index always corresponding to the constant term. In order to use the axis of `coeffs` as exponents, consider using a [`LaurentPolynomial`](@ref) or possibly a [`SparsePolynomial`](@ref). # Examples @@ -32,9 +32,9 @@ julia> one(Polynomial) Polynomial(1.0) ``` """ -struct Polynomial{T <: Number, X} <: StandardBasisPolynomial{T, X} +struct Polynomial{T, X} <: StandardBasisPolynomial{T, X} coeffs::Vector{T} - function Polynomial{T, X}(coeffs::AbstractVector{S}) where {T <: Number, X, S} + function Polynomial{T, X}(coeffs::AbstractVector{S}) where {T, X, S} if Base.has_offset_axes(coeffs) @warn "ignoring the axis offset of the coefficient vector" end @@ -44,7 +44,7 @@ struct Polynomial{T <: Number, X} <: StandardBasisPolynomial{T, X} new{T,X}(cs) end # non-copying alternative assuming !iszero(coeffs[end]) - function Polynomial{T, X}(checked::Val{false}, coeffs::AbstractVector{T}) where {T <: Number, X} + function Polynomial{T, X}(checked::Val{false}, coeffs::AbstractVector{T}) where {T, X} new{T, X}(coeffs) end end @@ -56,7 +56,7 @@ end # scalar +,* faster than standard-basis/common versions as it avoids a copy function Base.:+(p::P, c::S) where {T, X, P <: Polynomial{T, X}, S<:Number} - R = promote_type(T, S) + R = Base.promote_op(+, T, S) Q = Polynomial{R,X} as = convert(Vector{R}, copy(coeffs(p))) as[1] += c @@ -64,7 +64,7 @@ function Base.:+(p::P, c::S) where {T, X, P <: Polynomial{T, X}, S<:Number} end function Base.:*(p::P, c::S) where {T, X, P <: Polynomial{T,X} , S <: Number} - R = promote_type(T,S) + R = Base.promote_op(*, T, S) #promote_type(T,S) Q = Polynomial{R, X} as = R[aᵢ * c for aᵢ ∈ coeffs(p)] iszero(as[end]) ? Q(as) : Q(Val(false), as) @@ -88,9 +88,9 @@ function Base.:+(p1::P1, p2::P2) where {T,X, P1<:Polynomial{T,X}, Q(Val(false), cs) end - -function Base.:*(p::P, q::P) where {T,X, P<:Polynomial{T,X}} + +# redundant, a bit faster +function Base.:*(p::P, q::P) where {T <: Number,X, P<:Polynomial{T,X}} c = fastconv(p.coeffs, q.coeffs) return iszero(c[end]) ? P(c) : P(Val(false), c) end - diff --git a/src/polynomials/SparsePolynomial.jl b/src/polynomials/SparsePolynomial.jl index a2e0a69f..aa7995e9 100644 --- a/src/polynomials/SparsePolynomial.jl +++ b/src/polynomials/SparsePolynomial.jl @@ -5,7 +5,7 @@ export SparsePolynomial Polynomials in the standard basis backed by a dictionary holding the non-zero coefficients. For polynomials of high degree, this might be -advantageous. +advantageous. # Examples: @@ -38,31 +38,31 @@ julia> p(1) ``` """ -struct SparsePolynomial{T <: Number, X} <: StandardBasisPolynomial{T, X} +struct SparsePolynomial{T, X} <: StandardBasisPolynomial{T, X} coeffs::Dict{Int, T} - function SparsePolynomial{T, X}(coeffs::AbstractDict{Int, S}) where {T <: Number, X, S} + function SparsePolynomial{T, X}(coeffs::AbstractDict{Int, S}) where {T, X, S} c = Dict{Int, T}(coeffs) for (k,v) in coeffs iszero(v) && pop!(c, k) end new{T, X}(c) end - function SparsePolynomial{T,X}(checked::Val{false}, coeffs::AbstractDict{Int, T}) where {T <: Number, X} + function SparsePolynomial{T,X}(checked::Val{false}, coeffs::AbstractDict{Int, T}) where {T, X} new{T,X}(coeffs) end end @register SparsePolynomial -function SparsePolynomial{T}(coeffs::AbstractDict{Int, S}, var::SymbolLike=:x) where {T <: Number, S} +function SparsePolynomial{T}(coeffs::AbstractDict{Int, S}, var::SymbolLike=:x) where {T, S} SparsePolynomial{T, Symbol(var)}(convert(Dict{Int,T}, coeffs)) end -function SparsePolynomial(coeffs::AbstractDict{Int, T}, var::SymbolLike=:x) where {T <: Number} +function SparsePolynomial(coeffs::AbstractDict{Int, T}, var::SymbolLike=:x) where {T} SparsePolynomial{T, Symbol(var)}(coeffs) end -function SparsePolynomial{T,X}(coeffs::AbstractVector{S}) where {T <: Number, X, S} +function SparsePolynomial{T,X}(coeffs::AbstractVector{S}) where {T, X, S} if Base.has_offset_axes(coeffs) @warn "ignoring the axis offset of the coefficient vector" @@ -109,11 +109,11 @@ function coeffs(p::SparsePolynomial{T}) where {T} cs[k+1]=v end cs - + end # get/set index -function Base.getindex(p::SparsePolynomial{T}, idx::Int) where {T <: Number} +function Base.getindex(p::SparsePolynomial{T}, idx::Int) where {T} get(p.coeffs, idx, zero(T)) end @@ -151,16 +151,16 @@ end ## ## ---- ## - + function evalpoly(x::S, p::SparsePolynomial{T}) where {T,S} - tot = zero(T) * EvalPoly._one(x) + tot = zero(T) * EvalPoly._one(x) for (k,v) in p.coeffs tot = EvalPoly._muladd(x^k, v, tot) end - + return tot - + end # map: over values -- not keys @@ -186,7 +186,7 @@ function Base.:+(p::SparsePolynomial{T,X}, c::S) where {T, X, S <: Number} iszero(D[0]) && pop!(D,0) return P(Val(false), D) - + end # Implement over fallback. A bit faster as it covers T != S @@ -195,28 +195,42 @@ function Base.:+(p1::P1, p2::P2) where {T,X, P1<:SparsePolynomial{T,X}, R = promote_type(T,S) Q = SparsePolynomial{R,X} - + d1, d2 = degree(p1), degree(p2) cs = d1 > d2 ? ⊕(P1, p1.coeffs, p2.coeffs) : ⊕(P1, p2.coeffs, p1.coeffs) return d1 != d2 ? Q(Val(false), cs) : Q(cs) - + end ## Multiplication -function Base.:*(p::P, c::S) where {T, X, P <: SparsePolynomial{T,X}, S <: Number} +function scalar_mult(p::P, c::S) where {T, X, P <: SparsePolynomial{T,X}, S<:Number} - R = promote_type(T,S) + R1 = promote_type(T,S) + R = typeof(zero(c)*zero(T)) Q = ⟒(P){R,X} - - q = zero(Q) + q = zero(Q) for (k,pₖ) ∈ pairs(p) q[k] = pₖ * c end - + return q end +function scalar_mult(c::S, p::P) where {T, X, P <: SparsePolynomial{T,X}, S<:Number} + + R1 = promote_type(T,S) + R = typeof(zero(c)*zero(T)) + Q = ⟒(P){R,X} + q = zero(Q) + for (k,pₖ) ∈ pairs(p) + q[k] = c * pₖ + end + + return q +end + + function Base.:*(p::P, q::Q) where {T,X,P<:SparsePolynomial{T,X}, S, Q<:SparsePolynomial{S,X}} R = promote_type(T,S) @@ -225,7 +239,7 @@ end function derivative(p::SparsePolynomial{T,X}, order::Integer = 1) where {T,X} - + order < 0 && throw(ArgumentError("Order of derivative must be non-negative")) order == 0 && return p @@ -247,7 +261,7 @@ end function integrate(p::P) where {T, X, P<:SparsePolynomial{T,X}} - + R = eltype(one(T)/1) Q = SparsePolynomial{R,X} @@ -263,4 +277,3 @@ function integrate(p::P) where {T, X, P<:SparsePolynomial{T,X}} return ∫p end - diff --git a/src/polynomials/factored_polynomial.jl b/src/polynomials/factored_polynomial.jl index 25dbe067..3f4c07db 100644 --- a/src/polynomials/factored_polynomial.jl +++ b/src/polynomials/factored_polynomial.jl @@ -142,7 +142,7 @@ function printpoly(io::IO, p::FactoredPolynomial{T,X}, mimetype=nothing) where { print(io, "(") print(io, x) if hasneg(T) - isneg(k) ? print(io, " + ", -k) : print(io, " - ", k) + isneg(k) ? print(io, " + ", -k) : print(io, " - ", k) else print(io, " - ", k) end @@ -191,13 +191,13 @@ function Base.isapprox(p1::FactoredPolynomial{T,X}, # # sorting roots below works only with real roots... # isapprox(p1.c, p2.c, rtol=rtol, atol=atol) || return false # k1,k2 = sort(collect(keys(p1.coeffs)),by = x -> (real(x), imag(x))), sort(collect(keys(p2.coeffs)),by = x -> (real(x), imag(x))) - + # length(k1) == length(k2) || return false # for (k₁, k₂) ∈ zip(k1, k2) # isapprox(k₁, k₂, atol=atol, rtol=rtol) || return false # p1.coeffs[k₁] == p2.coeffs[k₂] || return false # end - + # return true end @@ -247,7 +247,7 @@ roots(p::FactoredPolynomial{T}) where {T} = Base.typed_vcat(T,[repeat([k],v) for # unary subtraction Base.:-(p::P) where {T,X,P<:FactoredPolynomial{T,X}} = (-1)*p -# addition +# addition function Base.:+(p::P, q::P) where {T,X,P<:FactoredPolynomial{T,X}} 𝑷 = Polynomial{T,X} 𝒑,𝒒 = convert(𝑷, p), convert(𝑷, q) @@ -264,7 +264,7 @@ function Base.:*(p::P, q::P) where {T,X, P<:FactoredPolynomial{T,X}} end # scalar mult -function Base.:*(p::P, c::S) where {S<:Number, T, X, P <: FactoredPolynomial{T, X}} +function scalar_mult(p::P, c::S) where {S<:Number, T, X, P <: FactoredPolynomial{T, X}} R = promote_type(T,S) d = Dict{R, Int}() # wident copy!(d, p.coeffs) @@ -304,7 +304,7 @@ function uvw(p::P, q::P; kwargs...) where {T, X, P<:FactoredPolynomial{T,X}} du, dv, dw = Dict{T,Int}(), Dict{T,Int}(), Dict{T,Int}() dp,dq = p.coeffs, q.coeffs kp,kq = keys(dp), keys(dq) - + for k ∈ setdiff(kp, kq) dv[k] = dp[k] end @@ -315,7 +315,7 @@ function uvw(p::P, q::P; kwargs...) where {T, X, P<:FactoredPolynomial{T,X}} pₖ,qₖ = dp[k], dq[k] m = min(pₖ, qₖ) du[k] = m - dv[k] = pₖ - m; + dv[k] = pₖ - m; dw[k] = qₖ - m end P(du), P(dv, p.c), P(dw, q.c) @@ -332,7 +332,7 @@ function Base.divrem(p::P, q::P) where {T, X, P<:FactoredPolynomial{T,X}} return (d, u*r) -end +end ## ---- @@ -349,4 +349,3 @@ function derivative(p::P,n::Int) where {P <: FactoredPolynomial} 𝒑⁽ⁿ⁾ = derivative(𝒑, n) convert(P, 𝒑⁽ⁿ⁾) end - diff --git a/src/polynomials/standard-basis.jl b/src/polynomials/standard-basis.jl index 0023b91b..20e08e1c 100644 --- a/src/polynomials/standard-basis.jl +++ b/src/polynomials/standard-basis.jl @@ -64,6 +64,9 @@ function variable(::Type{P}) where {P<:StandardBasisPolynomial} ⟒(P){T,X}([zero(T),one(T)]) end +# can bypass c*one(P) +Base.:+(p::P, c::T) where {T, P<:StandardBasisPolynomial{T}} = p + ⟒(P)([c]) + ## multiplication algorithms for computing p * q. ## default multiplication between same type. ## subtypes might relax to match T,S to avoid one conversion @@ -71,13 +74,32 @@ function Base.:*(p::P, q::P) where {T,X, P<:StandardBasisPolynomial{T,X}} cs = ⊗(P, coeffs(p), coeffs(q)) P(cs) end - -## put here, not with type defintion, in case reuse is possible -function ⊗(P::Type{<:StandardBasisPolynomial}, p::Vector{T}, q::Vector{S}) where {T,S} + +function ⊗(P::Type{<:StandardBasisPolynomial}, p::Vector{T}, q::Vector{S}) where {T<:Number,S<:Number} R = promote_type(T,S) fastconv(convert(Vector{R}, p), convert(Vector{R},q)) end +## put here, not with type defintion, in case reuse is possible +## `conv` can be used with matrix entries, unlike `fastconv` +function conv(p::Vector{T}, q::Vector{S}) where {T,S} + as = [p[1]*q[1]] + z = 0 * as[1] + n,m = length(p)-1, length(q)-1 + for i ∈ 1:n+m + Σ = z + for j ∈ max(0, i-m):min(i,n) + Σ += p[1+j]*q[1 + i-j] + end + push!(as, Σ) + end + as +end + +function ⊗(P::Type{<:StandardBasisPolynomial}, p::Vector{T}, q::Vector{S}) where {T,S} + conv(p, q) +end + ## Static size of product makes generated functions a good choice ## from https://github.com/tkoolen/StaticUnivariatePolynomials.jl/blob/master/src/monomial_basis.jl ## convolution of two tuples @@ -97,7 +119,7 @@ end return quote Base.@_inline_meta - tuple($(exprs...)) + tuple($(exprs...)) end end @@ -129,6 +151,7 @@ function fromroots(P::Type{<:StandardBasisPolynomial}, r::AbstractVector{T}; var return P(reverse(c), var) end +## ---- function derivative(p::P, order::Integer = 1) where {T, X, P <: StandardBasisPolynomial{T, X}} order < 0 && throw(ArgumentError("Order of derivative must be non-negative")) @@ -136,9 +159,9 @@ function derivative(p::P, order::Integer = 1) where {T, X, P <: StandardBasisPol # we avoid usage like Base.promote_op(*, T, Int) here, say, as # Base.promote_op(*, Rational, Int) is Any, not Rational in analogy to # Base.promote_op(*, Complex, Int) - R = eltype(one(T)*1) + R = typeof(constantterm(p)*1) Q = ⟒(P){R,X} - + order == 0 && return p hasnan(p) && return Q(R[NaN]) order > length(p) && return zero(Q) @@ -153,15 +176,16 @@ function derivative(p::P, order::Integer = 1) where {T, X, P <: StandardBasisPol end function integrate(p::P) where {T, X, P <: StandardBasisPolynomial{T, X}} - R = eltype(one(T)/1) - Q = ⟒(P){R,X} + + R = typeof(constantterm(p)/1) + Q = ⟒(P){R,X} hasnan(p) && return Q([NaN]) iszero(p) && return zero(Q) n = length(p) as = Vector{R}(undef, n + 1) - as[1] = zero(R) + as[1] = 0*constantterm(p) for (i, pᵢ) ∈ pairs(p) i′ = i + 1 @inbounds as[i′+1] = pᵢ/i′ @@ -171,7 +195,7 @@ end function Base.divrem(num::P, den::Q) where {T, P <: StandardBasisPolynomial{T}, S, Q <: StandardBasisPolynomial{S}} - assert_same_variable(num, den) + assert_same_variable(num, den) X = indeterminate(num) @@ -363,6 +387,52 @@ function companion(p::P) where {T, P <: StandardBasisPolynomial{T}} return comp end +# block companion matrix +function companion(p::P) where {T, M <: AbstractMatrix{T}, P<:StandardBasisPolynomial{M}} + C₀, C₁ = companion_pencil(p) + C₀ * inv(C₁) # could be more efficient +end + +# the companion pencil is `C₀`, `C₁` where `λC₁ - C₀` is singular for +# eigen values of the companion matrix: `vᵀ(λC₁ - C₀) = 0` => `vᵀλ = vᵀ(C₀C₁⁻¹)`, where `C₀C₁⁻¹` +# is the companion matrix. +function companion_pencil(p::P) where {T, P<:StandardBasisPolynomial{T}} + n = degree(p) + C₁ = diagm(0 => ones(T, n)) + C₁[end,end] = p[end] + + C₀ = diagm(-1 => ones(T, n-1)) + for i ∈ 0:n-1 + C₀[1+i,end] = -p[i] + end + C₀, C₁ +end + +# block version +function companion_pencil(p::P) where {T, M <: AbstractMatrix{T}, P<:StandardBasisPolynomial{M}} + + m,m′ = size(p[0]) + @assert m == m′ # square matrix + + n = degree(p) + + C₀ = zeros(T, n*m, n*m) + C₁ = zeros(T, n*m, n*m) + + + Δ = 1:m + for i ∈ 1:n-1 + C₁[(i-1)*m .+ Δ, (i-1)*m .+ Δ] .= I(m) + C₀[i*m .+ Δ, (i-1)*m .+ Δ] .= I(m) + end + C₁[(n-1)*m .+ Δ, (n-1)*m .+ Δ] .= p[end] + for i ∈ 0:n-1 + C₀[i*m .+ Δ, (n-1)*m .+ Δ] .= -p[i] + end + + C₀, C₁ +end + function roots(p::P; kwargs...) where {T, P <: StandardBasisPolynomial{T}} R = eltype(one(T)/one(T)) @@ -381,7 +451,7 @@ function roots(p::P; kwargs...) where {T, P <: StandardBasisPolynomial{T}} k == K && return zeros(R, k-1) - # find eigenvalues of the companion matrix + # find eigenvalues of the companion matrix of the 0-deflated polynomial comp = companion(⟒(P)(as[k:K], indeterminate(p))) L = eigvals(comp; kwargs...) append!(L, zeros(eltype(L), k-1)) diff --git a/src/show.jl b/src/show.jl index 1992472b..7c716e0f 100644 --- a/src/show.jl +++ b/src/show.jl @@ -7,7 +7,10 @@ export printpoly ## which can be modified by users for other Ts _iszero(x::T) where {T} = (x == zero(T)) === true +_iszero(x::AbstractArray{T}) where {T} = all(isequal.(x, zero(T))) + _isone(x::T) where {T} = (x == one(T)) === true +_isone(x::AbstractArray{T}) where {T} = all(isequal.(x, one(T))) "`hasneg(::T)` attribute is true if: `pj < zero(T)` is defined." hasneg(::Type{T}) where {T} = false @@ -216,7 +219,7 @@ julia> Polynomial([Dual(1,2), Dual(3,4)]) Polynomial(1 + 2ɛ + 3 + 4ɛ*x) ``` -```jldoctest +```jldoctest julia> using DualNumbers, Polynomials julia> function Base.show_unquoted(io::IO, pj::Dual, indent::Int, prec::Int) diff --git a/test/Project.toml b/test/Project.toml index abc0a286..0f1df4db 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,4 +5,5 @@ OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/StandardBasis.jl b/test/StandardBasis.jl index 1c8c1778..a83b6c5f 100644 --- a/test/StandardBasis.jl +++ b/test/StandardBasis.jl @@ -60,7 +60,6 @@ isimmutable(::Type{<:ImmutablePolynomial}) = true end end - @testset "Mapdomain" begin for P in Ps @@ -141,11 +140,220 @@ Base.getindex(z::ZVector, I::Int) = parent(z)[I + z.offset] end end +if VERSION >= v"1.5" +@testset "Non-number type" begin + conv = Polynomials.conv + @testset "T=Polynomial{Int,:y}" begin + for P in (Polynomial,) + + T = P{Int, :y} + a,b,c = T([1]), T([1,2]), T([1,2,3]) + p = P([a,b,c]) + q = P([a,b]) + s = 2 + d = c + + # scalar product + @test s*p == P([s*cᵢ for cᵢ ∈ [a,b,c]]) + @test p*s == P([cᵢ*s for cᵢ ∈ [a,b,c]]) + @test_throws ArgumentError d*p == P([d*cᵢ for cᵢ ∈ [a,b,c]]) # can't fix + @test_throws ArgumentError p*d == P([cᵢ*d for cᵢ ∈ [a,b,c]]) # can't fix + + # poly add + @test p + q == P([a+a,b+b,c]) + @test p - q == P([a-a,b-b,c]) + @test p - p == P([0*a]) + + # poly mult + @test p * q == P(conv([a,b,c], [a,b])) + @test q * p == P(conv([a,b], [a,b, c])) + + # poly powers + @test p^2 == p * p + + # evalution + @test p(s) == a + b * s + c * s * s + @test p(c) == a + b * c + c * c * c + + # ∂, ∫ + @test derivative(p) == P([b, 2c]) + @test integrate(p) == P([0*a, a, b/2, c/3]) + + # matrix element + pq = [p q] + @test pq[1] == p + @test pq[2] == q + + # implicit promotion + @test p + s == P([a+s, b, c]) + @test_throws Union{ArgumentError, MethodError} p + d == P([a+d, b, c]) # can't fix + @test p + P([d]) == P([a+d,b,c]) + + ps = [p s] + @test ps[1] == p + @test ps[2] == s + end + end + + @testset "T=Matrix (2x2)" begin + for P ∈ (Polynomial, ImmutablePolynomial) + a,b,c = [1 0; 1 1], [1 0; 2 1], [1 0; 3 1] + p = P([a,b,c]) + q = P([a,b]) + s = 2 + d = [4 1; 1 0] + + # scalar product + @test s*p == P([s*cᵢ for cᵢ ∈ [a,b,c]]) + @test p*s == P([cᵢ*s for cᵢ ∈ [a,b,c]]) + @test d*p == P([d*cᵢ for cᵢ ∈ [a,b,c]]) + @test p*d == P([cᵢ*d for cᵢ ∈ [a,b,c]]) + + # poly add + @test p + q == P([a+a,b+b,c]) + @test p - q == P([a-a,b-b,c]) + @test_throws MethodError p - p == P([0*a]) # no zeros to make zero polynomial + + # poly mult + @test p * q == P(conv([a,b,c], [a,b])) + @test q * p == P(conv([a,b], [a,b, c])) + + # poly powers + @test p^2 == p * p + + # evalution + @test p(s) == a + b * s + c * s * s + @test p(c) == a + b * c + c * c * c + + # ∂, ∫ + @test derivative(p) == P([b, 2c]) + @test integrate(p) == P([0*a, a, b/2, c/3]) + + # matrix element + @test [p q][1] == p + @test [p q][2] == q + + # implicit promotion + @test_throws MethodError p + s == P([a+s, b, c]) # OK, no a + s + @test p + d == P([a+d, b, c]) + @test p + P([d]) == P([a+d,b,c]) + + @test_throws MethodError [p s][1] == p # no promotion T(s) + @test_throws MethodError [p s][2] == s + end + end + + + @testset "T=Vector{Int}" begin + for P ∈ (Polynomial, ImmutablePolynomial) + a,b,c = [1,0,0], [1,1,0], [1,1,1] + p = P([a,b,c]) + q = P([a,b]) + s = 2 + d = [1,2,3] + + # scalar product + @test s*p == P([s*cᵢ for cᵢ ∈ [a,b,c]]) + @test p*s == P([cᵢ*s for cᵢ ∈ [a,b,c]]) + @test_throws MethodError d*p == P([d*cᵢ for cᵢ ∈ [a,b,c]]) # Ok, no * for T + @test_throws MethodError p*d == P([cᵢ*d for cᵢ ∈ [a,b,c]]) # Ok, no * for T + + # poly add + @test p + q == P([a+a,b+b,c]) + @test p - q == P([a-a,b-b,c]) + @test_throws MethodError p - p == P([0*a]) # no zero(T) to make zero polynomial + + # poly mult + @test_throws MethodError p * q == P(conv([a,b,c], [a,b])) # Ok, no * for T + @test_throws MethodError q * p == P(conv([a,b], [a,b, c])) # Ok, no * for T + + # poly powers + @test_throws MethodError p^2 == p * p # Ok, no * for T + + # evalution + @test p(s) == a + b * s + c * s * s + @test_throws MethodError p(c) == a + b * c + c * c * c # OK, no b * c + + # ∂, ∫ + @test derivative(p) == P([b, 2c]) + @test integrate(p) == P([0*a, a, b/2, c/3]) + + + # matrix element + @test [p q][1] == p + @test [p q][2] == q + + # implicit promotion + @test_throws MethodError p + s == P([a+s, b, c]) # OK, no a + s + @test p + d == P([a+d, b, c]) + @test p + P([d]) == P([a+d,b,c]) + + @test_throws MethodError [p s][1] == p # no promotion T(s) + @test_throws MethodError [p s][2] == s + end + end + + + if VERSION >= v"1.5" + eval(quote + using StaticArrays + end) + @testset "T=SA" begin + for P ∈ (Polynomial, ImmutablePolynomial ) + a,b,c = SA[1 0; 1 1], SA[1 0; 2 1], SA[1 0; 3 1] + p = P([a,b,c]) + q = P([a,b]) + s = 2 + d = SA[4 1; 1 0] + + # scalar product + @test s*p == P([s*cᵢ for cᵢ ∈ [a,b,c]]) + @test p*s == P([cᵢ*s for cᵢ ∈ [a,b,c]]) + @test d*p == P([d*cᵢ for cᵢ ∈ [a,b,c]]) + @test p*d == P([cᵢ*d for cᵢ ∈ [a,b,c]]) + + # poly add + @test p + q == P([a+a,b+b,c]) + @test p - p == P([0*a]) + + # poly mult + @test p * q == P(conv([a,b,c], [a,b])) + @test q * p == P(conv([a,b], [a,b, c])) + + # poly powers + @test p^2 == p * p + + + # evalution + @test p(s) == a + b * s + c * s * s + @test p(c) == a + b * c + c * c * c + + # ∂, ∫ + @test derivative(p) == P([b, 2c]) + @test integrate(p) == P([0*a, a, b/2, c/3]) + + # matrix element + @test [p q][1] == p + @test [p q][2] == q + + # implicit promotion + # @test_broken p + s == P([a .+ s, b, c]) # should error, doesn't + @test p + d == P([a + d, b, c]) + @test p + P([d]) == P([a + d,b,c]) + + @test_throws MethodError [p s][1] == p # no promotion T(s) + @test_throws MethodError [p s][2] == s # + end + end + end +end +end + @testset "OffsetVector" begin as = ones(3:4) bs = parent(as) - - + + for P in Ps # LaurentPolynomial accepts OffsetArrays; others throw warning if P == LaurentPolynomial @@ -154,11 +362,11 @@ end @test P(as) == P(bs) @test P{eltype(as)}(as) == P{eltype(as)}(bs) # (Or throw an error?) - # @test_throws ArgumentError P(as) + # @test_throws ArgumentError P(as) # @test P{eltype(as)}(as) == P{eltype(as)}(bs) end end - + a = [1,1] b = OffsetVector(a, axes(a)) c = ZVector(a) @@ -167,7 +375,7 @@ end if P == LaurentPolynomial && continue @test P(a) == P(b) == P(c) == P(d) end - + end end @@ -266,7 +474,7 @@ end @test pX != pS1 @test pS1 == pS2 @test pS1 == pS3 - + @test indeterminate(pS1 + pS1) == indeterminate(pS1) @test indeterminate(pS1 - pS1) == indeterminate(pS1) @test indeterminate(pS1 * pS1) == indeterminate(pS1) @@ -446,8 +654,8 @@ end @test Polynomials.evalpoly.(1/2, ps) ≈ [p(1/2) for p in ps] end - - + + # constant polynomials and type Ts = (Int, Float32, Float64, Complex{Int}, Complex{Float64}) for P in (Polynomial, ImmutablePolynomial, SparsePolynomial) @@ -761,21 +969,21 @@ end end meths = (Base.vect, Base.vcat, Base.hcat) for P in (Polynomial, ImmutablePolynomial, SparsePolynomial, LaurentPolynomial) - + p,q = P([1,2], :x), P([1,2], :y) P′′ = P == LaurentPolynomial ? P : P′ # different promotion rule - + # * should promote to Polynomial type if mixed (save Laurent Polynomial) @testset "promote mixed polys" begin for m ∈ meths - @test _test(m(p,p), P, :x) + @test _test(m(p,p), P, :x) @test _test(m(p,r), P′′, :x) end - + @test _test(Base.hvcat((2,1), p, r,[p r]), P′′, :x) - + end - + # * numeric constants should promote to a polynomial, when mixed @testset "promote numbers to poly" begin for m ∈ meths @@ -784,20 +992,20 @@ end @test _test(m(1,1,p), P, :x) @test _test(m(p,1,1), P, :x) end - - @test _test(Base.hvcat((3,1), 1, p, r,[1 p r]), P′′, :x) + + @test _test(Base.hvcat((3,1), 1, p, r,[1 p r]), P′′, :x) end - + # * non-constant polynomials must share the same indeterminate @testset "non constant polys share same X" begin for m ∈ meths @test_throws ArgumentError m(p,q) @test_throws ArgumentError m(p,s) end - + @test_throws ArgumentError Base.hvcat((2,1), p, q,[p q]) end - + # * constant polynomials are treated as `P{T,X}`, not elements of `T` @testset "constant polys" begin @@ -807,10 +1015,10 @@ end @test _test(m(1,1,one(p)), P, :x) @test _test(m(one(p),1,1), P, :x) end - - @test _test(Base.hvcat((3,1), 1, p, r,[1 p r]), P′′, :x) + + @test _test(Base.hvcat((3,1), 1, p, r,[1 p r]), P′′, :x) end - + # * Promotion can be forced to mix constant-polynomials @testset "Use typed constructor to mix constant polynomals" begin 𝑷,𝑸 = P{Int,:x}, P{Int,:y} # not typeof(p),... as Immutable carries N @@ -827,10 +1035,10 @@ end @test eltype(𝑷[1 one(p); one(p) one(q)]) == 𝑷 @test eltype(𝑸[1 one(p); one(p) one(q)]) == 𝑸 end - + @testset "hvcat" begin p,q = P([1,2],:x), P([1,2],:y) - + q1 = [q 1] q11 = [one(q) 1] @@ -843,10 +1051,10 @@ end @test eltype(𝑷[1 p; q11]) == 𝑷 @test eltype(Base.typed_hvcat(𝑷, (2, 1), 1, p, q11)) == 𝑷 end - + end - + end @testset "Linear Algebra" begin @@ -1143,7 +1351,7 @@ end p = Polynomial([1.234567890, 2.34567890]) io=IOBuffer(); printpoly(io, p, compact=true); @test String(take!(io)) == "1.23457 + 2.34568*x" io=IOBuffer(); printpoly(io, p, compact=true, mulsymbol=""); @test String(take!(io)) == "1.23457 + 2.34568x" - + ## issue 278 with complex @test printpoly_to_string(Polynomial([1 + im, 1, 2, im, 2im, 1+im, 1-im])) == "1 + im + x + 2*x^2 + im*x^3 + 2im*x^4 + (1 + im)x^5 + (1 - im)x^6"