use leja order for factored polynomial, follow up on #568 (#587)

* use leja order for factored polynomial, follow up on #568

* doctest with fix

* doctest

* docfix

* allow test with FactoredPolynomial

* remove cruft

Author: jverzani

Project.toml

@@ -2,10 +2,11 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "4.0.11"
+version = "4.0.12"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
@@ -33,6 +34,7 @@ LinearAlgebra = "<0.0.1, 1"
 MakieCore = "0.6,0.7, 0.8"
 MutableArithmetics = "1"
 OffsetArrays = "1"
+OrderedCollections = "1.6.3"
 RecipesBase = "0.7, 0.8, 1"
 Requires = "1.1"
 Setfield = "1"
@@ -56,5 +58,4 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ChainRulesCore", "DualNumbers", "FFTW", "LinearAlgebra",
-"MutableArithmetics", "SparseArrays", "OffsetArrays", "SpecialFunctions", "Test"]
+test = ["Aqua", "ChainRulesCore", "DualNumbers", "FFTW", "LinearAlgebra", "MutableArithmetics", "SparseArrays", "OffsetArrays", "SpecialFunctions", "Test"]

docs/src/index.md

@@ -576,7 +576,7 @@ julia> p = Polynomial([24, -50, 35, -10, 1])
 Polynomial(24 - 50*x + 35*x^2 - 10*x^3 + x^4)
 
 julia> q = convert(FactoredPolynomial, p) # noisy form of `factor`:
-FactoredPolynomial((x - 4.0000000000000036) * (x - 2.9999999999999942) * (x - 1.0000000000000002) * (x - 2.0000000000000018))
+FactoredPolynomial((x - 4.0000000000000036) * (x - 1.0000000000000002) * (x - 2.9999999999999942) * (x - 2.0000000000000018))
 
 julia> map(x -> round(x, digits=10), q)
 FactoredPolynomial((x - 4.0) * (x - 2.0) * (x - 3.0) * (x - 1.0))
@@ -792,10 +792,10 @@ julia> P = FactoredPolynomial
 FactoredPolynomial
 
 julia> p,q = fromroots(P, [1,2,3,4]), fromroots(P, [2,2,3,5])
-(FactoredPolynomial((x - 4) * (x - 2) * (x - 3) * (x - 1)), FactoredPolynomial((x - 5) * (x - 2)² * (x - 3)))
+(FactoredPolynomial((x - 4) * (x - 1) * (x - 3) * (x - 2)), FactoredPolynomial((x - 5) * (x - 2)² * (x - 3)))
 
 julia> pq = p // q
-((x - 4) * (x - 2) * (x - 3) * (x - 1)) // ((x - 5) * (x - 2)² * (x - 3))
+((x - 4) * (x - 1) * (x - 3) * (x - 2)) // ((x - 5) * (x - 2)² * (x - 3))
 
 julia> lowest_terms(pq)
 ((x - 4.0) * (x - 1.0)) // ((x - 5.0) * (x - 2.0))
@@ -814,7 +814,7 @@ julia> for (λ, rs) ∈ r # reconstruct p/q from output of `residues`
        end
 
 julia> d
-((x - 4.0) * (x - 1.0000000000000002)) // ((x - 5.0) * (x - 2.0))
+((x - 1.0000000000000002) * (x - 4.0)) // ((x - 5.0) * (x - 2.0))
 ```
 
 A basic plot recipe is provided.

src/Polynomials.jl

@@ -10,6 +10,7 @@ using LinearAlgebra
 import Base: evalpoly
 using Setfield
 using SparseArrays
+using OrderedCollections
 
 include("abstract.jl")
 include("show.jl")

src/common.jl

@@ -60,6 +60,8 @@ Construct a polynomial of the given type given the roots. If no type is given, d
 
 # Examples
 ```jldoctest
+julia> using Polynomials
+
 julia> r = [3, 2]; # (x - 3)(x - 2)
 
 julia> fromroots(r)
@@ -83,6 +85,8 @@ Construct a polynomial of the given type using the eigenvalues of the given matr
 
 # Examples
 ```jldoctest
+julia> using Polynomials
+
 julia> A = [1 2; 3 4]; # (x - 5.37228)(x + 0.37228)
 
 julia> fromroots(A)
@@ -754,6 +758,8 @@ Given values of x that are assumed to be unbounded (-∞, ∞), return values re
 
 # Examples
 ```jldoctest
+julia> using Polynomials
+
 julia> x = -10:10
 -10:10
 
@@ -976,6 +982,8 @@ Return the monomial `x` in the indicated polynomial basis.  If no type is give,
 
 # Examples
 ```jldoctest
+julia> using Polynomials
+
 julia> x = variable()
 Polynomial(x)
 
@@ -1237,6 +1245,8 @@ Find the greatest common denominator of two polynomials recursively using
 # Examples
 
 ```jldoctest
+julia> using Polynomials
+
 julia> gcd(fromroots([1, 1, 2]), fromroots([1, 2, 3]))
 Polynomial(4.0 - 6.0*x + 2.0*x^2)
 ```
@@ -1326,7 +1336,7 @@ function Base.isapprox(p1::AbstractPolynomial{T,X},
                        atol::Real = 0,) where {T,S,X}
     (hasnan(p1) || hasnan(p2)) && return false  # NaN poisons comparisons
     # copy over from abstractarray.jl
-    Δ  = norm(p1-p2)
+    Δ  = norm(p1-p2) # p1 - p2 does conversion to common type
     if isfinite(Δ)
         return Δ <= max(atol, rtol*max(norm(p1), norm(p2)))
     else

src/polynomials/factored_polynomial.jl

@@ -16,10 +16,10 @@ julia> p = FactoredPolynomial(Dict([0=>1, 1=>2, 3=>4]))
 FactoredPolynomial(x * (x - 3)⁴ * (x - 1)²)
 
 julia> q = fromroots(FactoredPolynomial, [0,1,2,3])
-FactoredPolynomial(x * (x - 2) * (x - 3) * (x - 1))
+FactoredPolynomial((x - 3) * x * (x - 2) * (x - 1))
 
 julia> p*q
-FactoredPolynomial(x² * (x - 2) * (x - 3)⁵ * (x - 1)³)
+FactoredPolynomial(x² * (x - 3)⁵ * (x - 1)³ * (x - 2))
 
 julia> p^1000
 FactoredPolynomial(x¹⁰⁰⁰ * (x - 3)⁴⁰⁰⁰ * (x - 1)²⁰⁰⁰)
@@ -31,29 +31,29 @@ julia> p = Polynomial([24, -50, 35, -10, 1])
 Polynomial(24 - 50*x + 35*x^2 - 10*x^3 + x^4)
 
 julia> q = convert(FactoredPolynomial, p) # noisy form of `factor`:
-FactoredPolynomial((x - 4.0000000000000036) * (x - 2.9999999999999942) * (x - 1.0000000000000002) * (x - 2.0000000000000018))
+FactoredPolynomial((x - 4.0000000000000036) * (x - 1.0000000000000002) * (x - 2.9999999999999942) * (x - 2.0000000000000018))
 
 julia> map(x->round(x, digits=12), q) # map works over factors and leading coefficient -- not coefficients in the standard basis
 FactoredPolynomial((x - 4.0) * (x - 2.0) * (x - 3.0) * (x - 1.0))
 ```
 """
 struct FactoredPolynomial{T <: Number, X} <: AbstractPolynomial{T, X}
-    coeffs::Dict{T,Int}
+    coeffs::OrderedDict{T,Int}
     c::T
-    function FactoredPolynomial{T, X}(checked::Val{false}, cs::Dict{T,Int}, c::T) where {T, X}
-        new{T,X}(cs,T(c))
+    function FactoredPolynomial{T, X}(checked::Val{false}, cs::AbstractDict{T,Int}, c::T) where {T, X}
+        new{T,X}(convert(OrderedDict,cs),T(c))
     end
-    function FactoredPolynomial{T, X}(cs::Dict{T,Int}, c=one(T)) where {T, X}
-        D = Dict{T,Int}()
+    function FactoredPolynomial{T, X}(cs::AbstractDict{T,Int}, c=one(T)) where {T, X}
+        D = OrderedDict{T,Int}()
         for (k,v) ∈ cs
             v > 0 && (D[k] = v)
         end
         FactoredPolynomial{T,X}(Val(false), D,T(c))
     end
-    function FactoredPolynomial(cs::Dict{T,Int}, c::S=1, var::SymbolLike=:x) where {T,S}
+    function FactoredPolynomial(cs::AbstractDict{T,Int}, c::S=1, var::SymbolLike=:x) where {T,S}
         X = Symbol(var)
         R = promote_type(T,S)
-        D = convert(Dict{R,Int}, cs)
+        D = convert(OrderedDict{R,Int}, cs)
         FactoredPolynomial{R,X}(D, R(c))
     end
 end
@@ -72,8 +72,14 @@ function FactoredPolynomial{T,X}(coeffs::AbstractVector{S}) where {T,S,X}
 
     zs = Multroot.multroot(p)
     c = p[end]
-    D = Dict(zip(zs.values, zs.multiplicities))
-    FactoredPolynomial(D, c, X)
+
+    D = LittleDict(k=>v for (k,v) ∈ zip(zs.values, zs.multiplicities))
+    D′ = OrderedDict{eltype(zs.values), Int}()
+    for z ∈ lejaorder(zs.values)
+        D′[z] = D[z]
+    end
+
+    FactoredPolynomial(D′, c, X)
 end
 
 function FactoredPolynomial{T}(coeffs::AbstractVector{S}, var::SymbolLike=:x) where {T,S}
@@ -100,22 +106,31 @@ function Base.convert(P::Type{<:FactoredPolynomial}, p::FactoredPolynomial{T,X})
     copy!(d, p.coeffs)
     FactoredPolynomial{𝑻,𝑿}(d, p.c)
 end
+
 Base.promote(p::P,q::Q) where {X,T,P<:FactoredPolynomial{T,X},Q<:FactoredPolynomial{T,X}} = p,q
+
 Base.promote_rule(::Type{<:FactoredPolynomial{T,X}}, ::Type{<:FactoredPolynomial{S,X}}) where {T,S,X} =
     FactoredPolynomial{promote_type(T,S), X}
+
 Base.promote_rule(::Type{<:FactoredPolynomial{T,X}}, ::Type{S}) where {T,S<:Number,X} =
     FactoredPolynomial{promote_type(T,S), X}
+
 FactoredPolynomial{T,X}(n::S) where {T,X,S<:Number} = T(n) * one(FactoredPolynomial{T,X})
+
 FactoredPolynomial{T}(n::S, var::SymbolLike=:x) where {T,S<:Number} = T(n) * one(FactoredPolynomial{T,Symbol(var)})
+
 FactoredPolynomial(n::S, var::SymbolLike=:x) where {S<:Number} = n * one(FactoredPolynomial{S,Symbol(var)})
+
 FactoredPolynomial(var::SymbolLike=:x) = variable(FactoredPolynomial, Symbol(var))
+
 (p::FactoredPolynomial)(x) = evalpoly(x, p)
 
 function Base.convert(::Type{<:Polynomial}, p::FactoredPolynomial{T,X}) where {T,X}
     x = variable(Polynomial{T,X})
     isconstant(p) && return Polynomial{T,X}(p.c)
     p(x)
 end
+
 function Base.convert(P::Type{<:FactoredPolynomial}, p::Polynomial{T,X}) where {T,X}
     isconstant(p) && return ⟒(P)(constantterm(p), X)
     ⟒(P)(coeffs(p), X)
@@ -239,8 +254,8 @@ end
 
 function fromroots(::Type{P}, r::AbstractVector{T}; var::SymbolLike=:x) where {T <: Number, P<:FactoredPolynomial}
     X = Symbol(var)
-    d = Dict{T,Int}()
-    for rᵢ ∈ r
+    d = OrderedDict{T,Int}()
+    for rᵢ ∈ lejaorder(r)
         d[rᵢ] = get(d, rᵢ, 0) + 1
     end
     FactoredPolynomial{T, X}(d)
@@ -291,8 +306,10 @@ end
 # scalar mult
 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)
+    d = OrderedDict{R, Int}() # wident
+    for (k,v) ∈ p.coeffs
+        d[k] = v
+    end
     FactoredPolynomial{R,X}(d, c * p.c)
 end
 scalar_mult(c::S, p::P) where {S<:Number, T, X, P <: FactoredPolynomial{T, X}} = scalar_mult(p, c) # assume commutative, as we have S <: Number
@@ -304,7 +321,7 @@ end
 
 function Base.:^(p::P, n::Integer) where {T,X, P<:FactoredPolynomial{T,X}}
     n >= 0 || throw(ArgumentError("n must be non-negative"))
-    d = Dict{T,Int}()
+    d = OrderedDict{T,Int}()
     for (k,v) ∈ p.coeffs
         d[k] = v*n
     end
@@ -316,7 +333,7 @@ end
 function Base.gcd(p::P, q::P) where {T, X, P<:FactoredPolynomial{T,X}}
     iszero(p) && return q
     iszero(q) && return p
-    d = Dict{T,Int}()
+    d = OrderedDict{T,Int}()
 
     for k ∈ intersect(keys(p.coeffs), keys(q.coeffs))
         d[k] = min(p.coeffs[k], q.coeffs[k])
@@ -327,7 +344,7 @@ end
 
 # return u,v,w with p = u*v , q = u*w
 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}()
+    du, dv, dw = OrderedDict{T,Int}(), OrderedDict{T,Int}(), OrderedDict{T,Int}()
     dp,dq = p.coeffs, q.coeffs
     kp,kq = keys(dp), keys(dq)
 

src/polynomials/standard-basis/immutable-polynomial.jl

@@ -35,6 +35,8 @@ are precluded from use in rational functions.
 # Examples
 
 ```jldoctest
+julia> using Polynomials
+
 julia> ImmutablePolynomial((1, 0, 3, 4))
 ImmutablePolynomial(1 + 3*x^2 + 4*x^3)
 

src/polynomials/standard-basis/laurent-polynomial.jl

@@ -17,6 +17,8 @@ Integration will fail if there is a `x⁻¹` term in the polynomial.
 
 # Examples:
 ```jldoctest
+julia> using Polynomials
+
 julia> P = LaurentPolynomial;
 
 julia> p = P([1,1,1],  -1)
@@ -48,6 +50,7 @@ LaurentPolynomial(1.0*x + 0.5*x² + 0.3333333333333333*x³)
 
 julia> integrate(p)  # x⁻¹  term is an issue
 ERROR: ArgumentError: Can't integrate Laurent polynomial with  `x⁻¹` term
+[...]
 
 julia> integrate(P([1,1,1], -5))
 LaurentPolynomial(-0.25*x⁻⁴ - 0.3333333333333333*x⁻³ - 0.5*x⁻²)
@@ -184,6 +187,8 @@ Call `p̂ = paraconj(p)` and `p̄` = conj(p)`, then this satisfies
 Examples:
 
 ```jldoctest
+julia> using Polynomials
+
 julia> z = variable(LaurentPolynomial, :z)
 LaurentPolynomial(z)
 
@@ -223,6 +228,8 @@ This satisfies for *imaginary* `s`: `conj(p(s)) = p̃(s) = (conj ∘ p)(s) = cco
 
 Examples:
 ```jldoctest
+julia> using Polynomials
+
 julia> s = 2im
 0 + 2im
 

src/polynomials/standard-basis/sparse-polynomial.jl

@@ -10,6 +10,8 @@ advantageous.
 # Examples:
 
 ```jldoctest
+julia> using Polynomials
+
 julia> P = SparsePolynomial;
 
 julia> p, q = P([1,2,3]), P([4,3,2,1])

test/StandardBasis.jl

@@ -985,9 +985,7 @@ end
             (b ≈ a) & isreal(coeffs(b))    # the coeff should be real
         end
         p = Polynomial([1; zeros(99); -1])
-        if P !== FactoredPolynomial
-            @test fromroots(P, roots(p)) * p[end] ≈ p
-        end
+        @test fromroots(P, roots(p)) * p[end] ≈ p
     end
 end