From dc6054b0551bfc12f70be62499a23252749e4c6e Mon Sep 17 00:00:00 2001 From: john verzani Date: Sat, 16 Nov 2024 10:19:24 -0500 Subject: [PATCH] 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 --- Project.toml | 7 ++- docs/src/index.md | 8 +-- src/Polynomials.jl | 1 + src/common.jl | 12 +++- src/polynomials/factored_polynomial.jl | 55 ++++++++++++------- .../standard-basis/immutable-polynomial.jl | 2 + .../standard-basis/laurent-polynomial.jl | 7 +++ .../standard-basis/sparse-polynomial.jl | 2 + test/StandardBasis.jl | 4 +- 9 files changed, 68 insertions(+), 30 deletions(-) diff --git a/Project.toml b/Project.toml index 160a89f8..94af7c06 100644 --- a/Project.toml +++ b/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"] diff --git a/docs/src/index.md b/docs/src/index.md index ab47ede5..6d27da12 100644 --- a/docs/src/index.md +++ b/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. diff --git a/src/Polynomials.jl b/src/Polynomials.jl index 821e511e..52fc1881 100644 --- a/src/Polynomials.jl +++ b/src/Polynomials.jl @@ -10,6 +10,7 @@ using LinearAlgebra import Base: evalpoly using Setfield using SparseArrays +using OrderedCollections include("abstract.jl") include("show.jl") diff --git a/src/common.jl b/src/common.jl index d7c9eadd..a7c2a8e5 100644 --- a/src/common.jl +++ b/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 diff --git a/src/polynomials/factored_polynomial.jl b/src/polynomials/factored_polynomial.jl index d9401bc7..f2fcdf14 100644 --- a/src/polynomials/factored_polynomial.jl +++ b/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,15 +106,23 @@ 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} @@ -116,6 +130,7 @@ function Base.convert(::Type{<:Polynomial}, p::FactoredPolynomial{T,X}) where {T 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) diff --git a/src/polynomials/standard-basis/immutable-polynomial.jl b/src/polynomials/standard-basis/immutable-polynomial.jl index abb63f61..aeadfbaf 100644 --- a/src/polynomials/standard-basis/immutable-polynomial.jl +++ b/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) diff --git a/src/polynomials/standard-basis/laurent-polynomial.jl b/src/polynomials/standard-basis/laurent-polynomial.jl index 62ee0ef5..0779482b 100644 --- a/src/polynomials/standard-basis/laurent-polynomial.jl +++ b/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 diff --git a/src/polynomials/standard-basis/sparse-polynomial.jl b/src/polynomials/standard-basis/sparse-polynomial.jl index 55d66dd3..0e6837f5 100644 --- a/src/polynomials/standard-basis/sparse-polynomial.jl +++ b/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]) diff --git a/test/StandardBasis.jl b/test/StandardBasis.jl index 7bb84446..d809a0da 100644 --- a/test/StandardBasis.jl +++ b/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