WIP: version 1 (#198)

This PR does:

* close #187; close #199; close #201

* require use of using Polynomials.PolyCompat for Poly and Pade types.

* Add ImmutablePolynomial type to take advantage of evalpoly function. Backed by an NTuple, not an array

* Add SparsePolynomial type backed by a dictionary, not an array
* Add evalpoly(x, p) method to enable broadcast for arrays of polynomials. Close #209
* refactor code from Polynomial.jl into standard-basis.jl to reuse amongst ImmutablePolynomial and Poly.jl
* add StandardBasisPolynomial abstract type
* modify + and * for immutable and sparse polynomials to close #206
* extend tests for Polynomial to cover ImmutablePolynomial, SparsePolynomial, rename test file
borrow constructorof from ConstructionBase.jl
* modify one, add oneunit
* remove restriction that T <: Number in AbstractPolynomial type. This is enforceable by subtypes
* Close #201 (issue with `round`)

Author: jverzani

Project.toml

@@ -2,7 +2,7 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "0.8.0"
+version = "1.0.0"
 
 [deps]
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"

README.md

@@ -11,7 +11,7 @@ Basic arithmetic, integration, differentiation, evaluation, and root finding ove
 ## Installation
 
 ```julia
-(v1.2) pkg> add Polynomials
+(v1.4) pkg> add Polynomials
 
 julia> using Polynomials
 ```
@@ -21,6 +21,8 @@ julia> using Polynomials
 #### Available Polynomials
 
 * `Polynomial` - Standard polynomials
+* `ImmutablePolynomial` - Standard polynomial backed by a tuple for faster evaluation of values
+* `SparsePolynomial` - Standard basis polynomial backed by a dictionary to hold  sparse high-degree  polynomials
 * `ChebyshevT` - Chebyshev polynomials of the first kind
 
 #### Construction and Evaluation
@@ -186,14 +188,24 @@ Polynomial objects also have other methods:
 
 ## Related Packages
 
+* [StaticUnivariatePolynomials.jl](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) Fixed-size univariate polynomials backed by a Tuple
+
 * [MultiPoly.jl](https://github.com/daviddelaat/MultiPoly.jl) for sparse multivariate polynomials
 
-* [MultivariatePolynomials.jl](https://github.com/blegat/MultivariatePolynomials.jl) for multivariate polynomials and moments of commutative or non-commutative variables
+* [DynamicPolynomals.jl](https://github.com/JuliaAlgebra/DynamicPolynomials.jl) Multivariate polynomials implementation of commutative and non-commutative variables
+
+* [MultivariatePolynomials.jl](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl) for multivariate polynomials and moments of commutative or non-commutative variables
+
+* [PolynomialRings](https://github.com/tkluck/PolynomialRings.jl) A library for arithmetic and algebra with multi-variable polynomials.
+
+* [AbstractAlgebra.jl](https://github.com/wbhart/AbstractAlgebra.jl) and [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series
+
+* [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) for a fast complex polynomial root finder. For larger degree problems, also [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) and [AMRVW](https://github.com/jverzani/AMRVW.jl).
 
-* [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series
 
-* [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) for a fast complex polynomial root finder
+## Legacy code
 
+As of v0.7, the internals of this package were greatly generalized and new types and method names were introduced. For compatability purposes, legacy code can be run after issuing `using Polynomials.PolyCompat`.
 
 ## Contributing
 

docs/src/index.md

@@ -120,7 +120,7 @@ resulting polynomial is one lower than the degree of `p`.
 
 ```jldoctest
 julia> derivative(Polynomial([1, 3, -1]))
-Polynomial(3.0 - 2.0*x)
+Polynomial(3 - 2*x)
 ```
 
 ### Root-finding
@@ -211,13 +211,21 @@ julia> as[3], p[2], collect(p)[3]
 
 ## Related Packages
 
+* [StaticUnivariatePolynomials.jl](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) Fixed-size univariate polynomials backed by a Tuple
+
 * [MultiPoly.jl](https://github.com/daviddelaat/MultiPoly.jl) for sparse multivariate polynomials
 
-* [MultivariatePolynomials.jl](https://github.com/blegat/MultivariatePolynomials.jl) for multivariate polynomials and moments of commutative or non-commutative variables
+* [DynamicPolynomals.jl](https://github.com/JuliaAlgebra/DynamicPolynomials.jl) Multivariate polynomials implementation of commutative and non-commutative variables
+
+* [MultivariatePolynomials.jl](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl) for multivariate polynomials and moments of commutative or non-commutative variables
+
+* [PolynomialRings](https://github.com/tkluck/PolynomialRings.jl) A library for arithmetic and algebra with multi-variable polynomials.
+
+* [AbstractAlgebra.jl](https://github.com/wbhart/AbstractAlgebra.jl) and [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series
+
+* [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) for a fast complex polynomial root finder. For larger degree problems, also [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) and [AMRVW](https://github.com/jverzani/AMRVW.jl).
 
-* [AbstractAlgeebra.jl](https://github.com/wbhart/AbstractAlgebra.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series.
 
-* [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) for a fast complex polynomial root finder
 
 ## Contributing
 

docs/src/polynomials/polynomial.md

@@ -11,7 +11,8 @@ Polynomial
 ```
 
 ```@docs
-Pade
-Pade(x)
+ImmutablePolynomial
+SparsePolynomial
 ```
 
+

src/Polynomials.jl

@@ -13,10 +13,14 @@ include("common.jl")
 
 
 # Polynomials
+include("polynomials/standard-basis.jl")
 include("polynomials/Polynomial.jl")
+include("polynomials/ImmutablePolynomial.jl")
+include("polynomials/SparsePolynomial.jl")
+
 include("polynomials/ChebyshevT.jl")
 
-# to be deprecated, then removed
-include("compat.jl")
+# compat; opt-in with `using Polynomials.PolyCompat`
+include("polynomials/Poly.jl")
 
 end # module

src/abstract.jl

@@ -35,17 +35,62 @@ macro register(name)
     poly = esc(name)
     quote
         Base.convert(::Type{P}, p::P) where {P<:$poly} = p
-        Base.convert(P::Type{<:$poly}, p::$poly) where {T} = P(coeffs(p), p.var)
+        Base.convert(P::Type{<:$poly}, p::$poly{T}) where {T} = P(coeffs(p), p.var)
         Base.promote_rule(::Type{$poly{T}}, ::Type{$poly{S}}) where {T,S} =
             $poly{promote_type(T, S)}
         Base.promote_rule(::Type{$poly{T}}, ::Type{S}) where {T,S<:Number} =
             $poly{promote_type(T, S)}
         $poly(coeffs::AbstractVector{T}, var::SymbolLike = :x) where {T} =
             $poly{T}(coeffs, Symbol(var))
-        $poly{T}(x::AbstractVector{S}, var = :x) where {T,S<:Number} = $poly(T.(x), var)
-        $poly{T}(n::Number, var = :x) where {T} = $poly(T(n), var)
+        $poly{T}(x::AbstractVector{S}, var = :x) where {T,S<:Number} =
+            $poly(T.(x), var)
+        $poly{T}(n::S, var = :x) where {T, S<:Number} =
+            $poly(T[n], var)
         $poly(n::Number, var = :x) = $poly([n], var)
-        $poly{T}(var=:x) where {T} = variable($poly{T}, var)
-        $poly(var=:x) = variable($poly, var)
+        $poly{T}(var::SymbolLike=:x) where {T} = variable($poly{T}, Symbol(var))
+        $poly(var::SymbolLike=:x) = variable($poly, Symbol(var))
     end
 end
+
+
+# Macros to register POLY{α, T} and POLY{α, β, T}
+macro register1(name)
+    poly = esc(name)
+    quote
+        Base.convert(::Type{P}, p::P) where {P<:$poly} = p
+        Base.promote_rule(::Type{$poly{α,T}}, ::Type{$poly{α,S}}) where {α,T,S} =
+            $poly{α,promote_type(T, S)}
+        Base.promote_rule(::Type{$poly{α,T}}, ::Type{S}) where {α,T,S<:Number} = 
+            $poly{α,promote_type(T,S)}
+        function $poly{α,T}(x::AbstractVector{S}, var::Polynomials.SymbolLike = :x) where {α,T,S}
+            $poly{α,T}(T.(x), Symbol(var))
+        end
+        $poly{α}(coeffs::AbstractVector{T}, var::Polynomials.SymbolLike=:x) where {α,T} =
+            $poly{α,T}(coeffs, Symbol(var))
+        $poly{α,T}(n::Number, var::Polynomials.SymbolLike = :x) where {α,T} = n*one($poly{α,T}, Symbol(var))
+        $poly{α}(n::Number, var::Polynomials.SymbolLike = :x) where {α} = n*one($poly{α}, Symbol(var))
+        $poly{α,T}(var::Polynomials.SymbolLike=:x) where {α, T} = variable($poly{α,T}, Symbol(var))
+        $poly{α}(var::Polynomials.SymbolLike=:x) where {α} = variable($poly{α}, Symbol(var))
+    end
+end
+
+
+# Macro to register POLY{α, β, T}
+macro register2(name)
+    poly = esc(name)
+    quote
+        Base.convert(::Type{P}, p::P) where {P<:$poly} = p
+        Base.promote_rule(::Type{$poly{α,β,T}}, ::Type{$poly{α,β,S}}) where {α,β,T,S} =
+            $poly{α,β,promote_type(T, S)}
+        Base.promote_rule(::Type{$poly{α,β,T}}, ::Type{S}) where {α,β,T,S<:Number} =
+            $poly{α,β,promote_type(T, S)}
+        $poly{α,β}(coeffs::AbstractVector{T}, var::Polynomials.SymbolLike = :x) where {α,β,T} =
+            $poly{α,β,T}(coeffs, Symbol(var))
+        $poly{α,β,T}(x::AbstractVector{S}, var = :x) where {α,β,T,S<:Number} = $poly{α,β,T}(T.(x), var)
+        $poly{α,β,T}(n::Number, var = :x) where {α,β,T} = n*one($poly{α,β,T}, var)
+        $poly{α,β}(n::Number, var = :x) where {α,β} = n*one($poly{α,β}, var)
+        $poly{α,β,T}(var=:x) where {α,β, T} = variable($poly{α,β,T}, var)
+        $poly{α,β}(var=:x) where {α,β} = variable($poly{α,β}, var)
+    end
+end
+

src/common.jl

@@ -110,7 +110,7 @@ Returns the roots of the given polynomial. This is calculated via the eigenvalue
 
 !!! note
 
-    The [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) package provides an alternative that is a bit faster and  abit more accurate; the [AMRVW.jl](https://github.com/jverzani/AMRVW.jl) package provides an alternative for high-degree polynomials.
+        The [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) package provides an alternative that is a bit faster and a bit more accurate; the [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) provides an interface to FORTRAN code implementing an algorithm that can handle very large polynomials (it is  `O(n^2)` not `O(n^3)`. the [AMRVW.jl](https://github.com/jverzani/AMRVW.jl) package implements the algorithm in Julia, allowing the use of other  number types.
 
 """
 function roots(q::AbstractPolynomial{T}; kwargs...) where {T <: Number}
@@ -364,8 +364,8 @@ Base.broadcastable(p::AbstractPolynomial) = Ref(p)
 
 # basis
 # return the kth basis polynomial for the given polynomial type, e.g. x^k for Polynomial{T}
-function basis(p::P, k::Int) where {P<:AbstractPolynomial}
-    basis(P, k)
+function basis(p::P, k::Int; var=:x) where {P<:AbstractPolynomial}
+    basis(P, k, var=var)
 end
 
 function basis(::Type{P}, k::Int; var=:x) where {P <: AbstractPolynomial}
@@ -384,14 +384,6 @@ end
 
 Base.collect(p::P) where {P <: AbstractPolynomial} = collect(P, p)
 
-function Base.getproperty(p::AbstractPolynomial, nm::Symbol)
-    if nm == :a
-        throw(ArgumentError("AbstractPolynomial.a is not supported, use coeffs(AbstractPolynomial) instead."))
-    end
-    return getfield(p, nm)
-end
-
-
 # getindex
 function Base.getindex(p::AbstractPolynomial{T}, idx::Int) where {T <: Number}
     idx < 0 && throw(BoundsError(p, idx))
@@ -434,6 +426,7 @@ Base.hash(p::AbstractPolynomial, h::UInt) = hash(p.var, hash(coeffs(p), h))
 
 Returns a representation of 0 as the given polynomial.
 """
+Base.zero(::Type{P}, var=:x) where {T, P <: AbstractPolynomial{T}} = P(zeros(T, 1), var)
 Base.zero(::Type{P}, var=:x) where {P <: AbstractPolynomial} = P(zeros(1), var)
 Base.zero(p::P) where {P <: AbstractPolynomial} = zero(P, p.var)
 """
@@ -442,9 +435,11 @@ Base.zero(p::P) where {P <: AbstractPolynomial} = zero(P, p.var)
 
 Returns a representation of 1 as the given polynomial.
 """
+Base.one(::Type{P}, var=:x) where {T, P <: AbstractPolynomial{T}} = P(ones(T, 1), var)
 Base.one(::Type{P}, var=:x) where {P <: AbstractPolynomial} = P(ones(1), var)
 Base.one(p::P) where {P <: AbstractPolynomial} = one(P, p.var)
-
+Base.oneunit(p::P, args...) where {P <: AbstractPolynomial} = one(p, args...)
+Base.oneunit(::Type{P}, args...) where {P <: AbstractPolynomial} = one(P, args...)
 #=
 arithmetic =#
 Base.:-(p::P) where {P <: AbstractPolynomial} = P(-coeffs(p), p.var)
@@ -457,10 +452,12 @@ function Base.:*(p::P, c::S) where {P <: AbstractPolynomial,S}
     T = promote_type(P, S)
     return T(coeffs(p) .* c, p.var)
 end
+
 function Base.:/(p::P, c::S) where {T,P <: AbstractPolynomial{T},S}
     R = promote_type(P, eltype(one(T) / one(S)))
     return R(coeffs(p) ./ c, p.var)
 end
+
 Base.:-(p1::AbstractPolynomial, p2::AbstractPolynomial) = +(p1, -p2)
 
 function Base.:+(p::P, n::Number) where {P <: AbstractPolynomial}
@@ -530,7 +527,7 @@ Comparisons =#
 Base.isequal(p1::P, p2::P) where {P <: AbstractPolynomial} = hash(p1) == hash(p2)
 Base.:(==)(p1::AbstractPolynomial, p2::AbstractPolynomial) =
     (p1.var == p2.var) && (coeffs(p1) == coeffs(p2))
-Base.:(==)(p::AbstractPolynomial, n::Number) = coeffs(p) == [n]
+Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && p[0] == n
 Base.:(==)(n::Number, p::AbstractPolynomial) = p == n
 
 function Base.isapprox(p1::AbstractPolynomial{T},
@@ -550,9 +547,9 @@ function Base.isapprox(p1::AbstractPolynomial{T},
 end
 
 function Base.isapprox(p1::AbstractPolynomial{T},
-    n::S;
-    rtol::Real = (Base.rtoldefault(T, S, 0)),
-    atol::Real = 0,) where {T,S}
+                       n::S;
+                       rtol::Real = (Base.rtoldefault(T, S, 0)),
+                       atol::Real = 0,) where {T,S}
     p1t = truncate(p1, rtol = rtol, atol = atol)
     if length(p1t) != 1
         return false

src/contrib.jl

@@ -31,11 +31,11 @@ end
 ## Code from Julia 1.4   (https://github.com/JuliaLang/julia/blob/master/base/math.jl#L101 on 4/8/20)
 ## cf. https://github.com/JuliaLang/julia/pull/32753
 ## Slight modification when `x` is a matrix
-## Remove once dependencies for Julia 1.0.0
+## Remove once dependencies for Julia 1.0.0 are dropped
 function evalpoly(x::S, p::Tuple) where {S}
     if @generated
         N = length(p.parameters)
-        ex = :(p[end])
+        ex = :(p[end]*_one(S))
         for i in N-1:-1:1
             ex = :(_muladd(x, $ex, p[$i]))
         end
@@ -49,7 +49,7 @@ evalpoly(x, p::AbstractVector) = _evalpoly(x, p)
 
 function _evalpoly(x::S, p) where {S}
     N = length(p)
-    ex = p[end]
+    ex = p[end]*_one(x)
     for i in N-1:-1:1
         ex = _muladd(x, ex, p[i])
     end
@@ -59,8 +59,8 @@ end
 function evalpoly(z::Complex, p::Tuple)
     if @generated
         N = length(p.parameters)
-        a = :(p[end])
-        b = :(p[end-1])
+        a = :(p[end]*_one(z))
+        b = :(p[end-1]*_one(z))
         as = []
         for i in N-2:-1:1
             ai = Symbol("a", i)
@@ -81,16 +81,16 @@ function evalpoly(z::Complex, p::Tuple)
         _evalpoly(z, p)
     end
 end
-evalpoly(z::Complex, p::Tuple{<:Any}) = p[1]
+evalpoly(z::Complex, p::Tuple{<:Any}) = p[1]*_one(z)
 
 
 evalpoly(z::Complex, p::AbstractVector) = _evalpoly(z, p)
 
 function _evalpoly(z::Complex, p)
-    length(p) == 1 && return p[1]
+    length(p) == 1 && return p[1]*_one(z)
     N = length(p)
-    a = p[end]
-    b = p[end-1]
+    a = p[end]*_one(z)
+    b = p[end-1]*_one(z)
 
     x = real(z)
     y = imag(z)
@@ -105,7 +105,25 @@ function _evalpoly(z::Complex, p)
     muladd(ai, z, b)
 end
 
-## modify muladd for matrices
+## 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
 
+# try to get y = P(c::T)(x::S) = P{T}(c)(x::S) to
+# have y == one(T)*one(S)*x
+_one(P::Type{<:Matrix}) = one(eltype(P))*I
+_one(x::Matrix) = one(eltype(x))*I
+_one(x) = one(x)
+
+## get type of parametric composite type without type parameters
+## this is needed when the underlying type changes, e.g. with integration
+## where T=Int might become T=Float64
+##
+## trick from [ConstructionBase.jl](https://github.com/JuliaObjects/ConstructionBase.jl/blob/b5686b755bd3bee29b181b3cb18fe2effa0f10a2/src/ConstructionBase.jl#L25)
+## as noted in https://discourse.julialang.org/t/get-new-type-with-different-parameter/37253/4
+##
+@generated function constructorof(::Type{T}) where T
+    getfield(parentmodule(T), nameof(T))
+end
+

src/pade.jl

@@ -27,7 +27,6 @@ struct Pade{T <: Number,S <: Number}
     q::Union{Poly{S}, Polynomial{S}}
     var::Symbol
     function Pade{T,S}(p::Union{Poly{T}, Polynomial{T}}, q::Union{Poly{S}, Polynomial{S}}) where {T,S}
-        Base.depwarn("Use of `Pade` from v1.0 forward will require `using Polynomials.PolyCompat`", :Pade)
         if p.var != q.var error("Polynomials must have same variable") end
         new{T,S}(p, q, p.var)
     end
@@ -89,7 +88,7 @@ Evaluate the Pade approximant at the given point.
 
 # Examples
 ```jldoctest pade
-julia> using SpecialFunctions, Polynomials
+julia> using Polynomials, Polynomials.PolyCompat, SpecialFunctions
 
 
 julia> p = Polynomial(@.(1 // BigInt(gamma(1:17))));