Skip to content

Commit

Permalink
use leja order for factored polynomial, follow up on #568 (#587)
Browse files Browse the repository at this point in the history 
* use leja order for factored polynomial, follow up on #568

* doctest with fix

* doctest

* docfix

* allow test with FactoredPolynomial

* remove cruft
  • Loading branch information
@jverzani
jverzani authored Nov 16, 2024
1 parent d5a472a commit dc6054b
Show file tree
Hide file tree
Showing 9 changed files with 68 additions and 30 deletions.
7 changes: 4 additions & 3 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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"
Expand All @@ -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"]
8 changes: 4 additions & 4 deletions docs/src/index.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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))
Expand All @@ -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.
Expand Down
1 change: 1 addition & 0 deletions src/Polynomials.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ using LinearAlgebra
import Base: evalpoly
using Setfield
using SparseArrays
using OrderedCollections

include("abstract.jl")
include("show.jl")
Expand Down
12 changes: 11 additions & 1 deletion src/common.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
```
Expand Down Expand Up @@ -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
Expand Down
55 changes: 36 additions & 19 deletions src/polynomials/factored_polynomial.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)²⁰⁰⁰)
Expand All @@ -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
Expand All @@ -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}
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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])
Expand All @@ -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)

Expand Down
2 changes: 2 additions & 0 deletions src/polynomials/standard-basis/immutable-polynomial.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
7 changes: 7 additions & 0 deletions src/polynomials/standard-basis/laurent-polynomial.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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⁻²)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
2 changes: 2 additions & 0 deletions src/polynomials/standard-basis/sparse-polynomial.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand Down
4 changes: 1 addition & 3 deletions test/StandardBasis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down

2 comments on commit dc6054b

@jverzani
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/119579

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v4.0.12 -m "<description of version>" dc6054b0551bfc12f70be62499a23252749e4c6e
git push origin v4.0.12

Please sign in to comment.