From 96d95e754ff74f2d08332caa99dfda6e8b8a3454 Mon Sep 17 00:00:00 2001 From: jverzani Date: Fri, 15 Nov 2024 15:36:02 -0500 Subject: [PATCH 1/6] use leja order for factored polynomial, follow up on #568 --- Project.toml | 7 ++-- src/Polynomials.jl | 1 + src/polynomials/factored_polynomial.jl | 49 +++++++++++++++++--------- 3 files changed, 38 insertions(+), 19 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/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/polynomials/factored_polynomial.jl b/src/polynomials/factored_polynomial.jl index d9401bc7..b9f462d3 100644 --- a/src/polynomials/factored_polynomial.jl +++ b/src/polynomials/factored_polynomial.jl @@ -38,22 +38,22 @@ 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) From 835ff6b9c140fe0b9792040b3debddae356e6109 Mon Sep 17 00:00:00 2001 From: jverzani Date: Fri, 15 Nov 2024 17:28:20 -0500 Subject: [PATCH 2/6] doctest with fix --- docs/src/index.md | 113 +++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 101 insertions(+), 12 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index ab47ede5..2fed82f4 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -95,7 +95,13 @@ Polynomial(1 + 2*s + 3*s^2) julia> p + q ERROR: ArgumentError: Polynomials have different indeterminates -[...] +Stacktrace: + [1] assert_same_variable(X::Symbol, Y::Symbol) + @ Polynomials ~/julia/Polynomials/src/common.jl:522 + [2] +(p::Polynomial{Int64, :x}, q::Polynomial{Int64, :s}) + @ Polynomials ~/julia/Polynomials/src/common.jl:1114 + [3] top-level scope + @ none:1 ``` Except for operations involving constant polynomials. @@ -189,15 +195,70 @@ Polynomial(-6.0 + 11.0*x - 6.0*x^2 + 1.0*x^3) julia> roots(p) ERROR: MethodError: no method matching eigvals!(::Matrix{BigFloat}) -[...] +The function `eigvals!` exists, but no method is defined for this combination of argument types. + +Closest candidates are: + eigvals!(::AbstractMatrix{T}, !Matched::LinearAlgebra.Cholesky{T}; sortby) where T<:Number + @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/symmetriceigen.jl:260 + eigvals!(::StridedMatrix{T}, !Matched::LinearAlgebra.LU{T, <:StridedMatrix{T} where T}; sortby) where T + @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/symmetriceigen.jl:286 + eigvals!(!Matched::LinearAlgebra.SymTridiagonal{<:Union{Float32, Float64}, <:StridedVector{T} where T}) + @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/tridiag.jl:270 + ... + +Stacktrace: + [1] eigvals(A::Matrix{BigFloat}; kws::@Kwargs{}) + @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/eigen.jl:343 + [2] eigvals(A::Matrix{BigFloat}) + @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/eigen.jl:343 + [3] roots(p::Polynomial{BigFloat, :x}; kwargs::@Kwargs{}) + @ Polynomials ~/julia/Polynomials/src/polynomials/standard-basis/standard-basis.jl:500 + [4] roots(p::Polynomial{BigFloat, :x}) + @ Polynomials ~/julia/Polynomials/src/polynomials/standard-basis/standard-basis.jl:481 + [5] top-level scope + @ none:1 julia> using GenericLinearAlgebra +ERROR: ArgumentError: Package GenericLinearAlgebra not found in current path. +- Run `import Pkg; Pkg.add("GenericLinearAlgebra")` to install the GenericLinearAlgebra package. +Stacktrace: + [1] macro expansion + @ ./loading.jl:2223 [inlined] + [2] macro expansion + @ ./lock.jl:273 [inlined] + [3] __require(into::Module, mod::Symbol) + @ Base ./loading.jl:2198 + [4] #invoke_in_world#3 + @ ./essentials.jl:1089 [inlined] + [5] invoke_in_world + @ ./essentials.jl:1086 [inlined] + [6] require(into::Module, mod::Symbol) + @ Base ./loading.jl:2191 julia> roots(p) -3-element Vector{Complex{BigFloat}}: - 0.9999999999999999999999999999999999999999999999999999999999999999999999999999655 + 0.0im - 1.999999999999999999999999999999999999999999999999999999999999999999999999999931 - 0.0im - 2.999999999999999999999999999999999999999999999999999999999999999999999999999793 + 0.0im +ERROR: MethodError: no method matching eigvals!(::Matrix{BigFloat}) +The function `eigvals!` exists, but no method is defined for this combination of argument types. + +Closest candidates are: + eigvals!(::AbstractMatrix{T}, !Matched::LinearAlgebra.Cholesky{T}; sortby) where T<:Number + @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/symmetriceigen.jl:260 + eigvals!(::StridedMatrix{T}, !Matched::LinearAlgebra.LU{T, <:StridedMatrix{T} where T}; sortby) where T + @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/symmetriceigen.jl:286 + eigvals!(!Matched::LinearAlgebra.SymTridiagonal{<:Union{Float32, Float64}, <:StridedVector{T} where T}) + @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/tridiag.jl:270 + ... + +Stacktrace: + [1] eigvals(A::Matrix{BigFloat}; kws::@Kwargs{}) + @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/eigen.jl:343 + [2] eigvals(A::Matrix{BigFloat}) + @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/eigen.jl:343 + [3] roots(p::Polynomial{BigFloat, :x}; kwargs::@Kwargs{}) + @ Polynomials ~/julia/Polynomials/src/polynomials/standard-basis/standard-basis.jl:500 + [4] roots(p::Polynomial{BigFloat, :x}) + @ Polynomials ~/julia/Polynomials/src/polynomials/standard-basis/standard-basis.jl:481 + [5] top-level scope + @ none:1 ``` ### Comments on root finding @@ -576,7 +637,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)) @@ -651,7 +712,23 @@ However, as there would be an ambiguous outcome of the following ```jldoctest natural_inclusion julia> [one(p) one(q)] ERROR: ArgumentError: Polynomials have different indeterminates -[...] +Stacktrace: + [1] assert_same_variable(X::Symbol, Y::Symbol) + @ Polynomials ~/julia/Polynomials/src/common.jl:522 + [2] promote_rule(::Type{Polynomial{Int64, :x}}, ::Type{Polynomial{Int64, :y}}) + @ Polynomials ~/julia/Polynomials/src/promotions.jl:24 + [3] promote_type + @ ./promotion.jl:318 [inlined] + [4] promote_eltypeof(v1::Polynomial{Int64, :x}, v2::Polynomial{Int64, :y}) + @ Base ./abstractarray.jl:1618 + [5] _cat(::Val{2}, ::Polynomial{Int64, :x}, ::Polynomial{Int64, :y}) + @ Base ./abstractarray.jl:1820 + [6] cat + @ ./abstractarray.jl:2084 [inlined] + [7] hcat(::Polynomial{Int64, :x}, ::Polynomial{Int64, :y}) + @ Base ./abstractarray.jl:1968 + [8] top-level scope + @ none:1 ``` an error is thrown. @@ -683,7 +760,19 @@ but not this one: ```jldoctest natural_inclusion julia> [one(p), one(p) + one(q)] ERROR: ArgumentError: Polynomials have different indeterminates -[...] +Stacktrace: + [1] assert_same_variable(X::Symbol, Y::Symbol) + @ Polynomials ~/julia/Polynomials/src/common.jl:522 + [2] promote_rule(::Type{Polynomial{Int64, :x}}, ::Type{Polynomial{Int64, :y}}) + @ Polynomials ~/julia/Polynomials/src/promotions.jl:24 + [3] promote_type + @ ./promotion.jl:318 [inlined] + [4] promote_typeof(x::Polynomial{Int64, :x}, y::Polynomial{Int64, :y}) + @ Base ./promotion.jl:378 + [5] vect(::Polynomial{Int64, :x}, ::Vararg{Any}) + @ Base ./array.jl:184 + [6] top-level scope + @ none:1 ``` Also, mixing types can result in unspecific symbols, as this example shows: @@ -792,10 +881,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 +903,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. From 2d59f1f3c4352d72108dd3e7fb629c844ff7dcce Mon Sep 17 00:00:00 2001 From: jverzani Date: Fri, 15 Nov 2024 18:24:41 -0500 Subject: [PATCH 3/6] doctest --- docs/src/index.md | 42 ++----------------- src/common.jl | 12 +++++- src/polynomials/factored_polynomial.jl | 6 +-- .../standard-basis/immutable-polynomial.jl | 2 + .../standard-basis/laurent-polynomial.jl | 11 +++++ .../standard-basis/sparse-polynomial.jl | 2 + 6 files changed, 33 insertions(+), 42 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 2fed82f4..51c0dcf5 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -219,46 +219,12 @@ Stacktrace: @ none:1 julia> using GenericLinearAlgebra -ERROR: ArgumentError: Package GenericLinearAlgebra not found in current path. -- Run `import Pkg; Pkg.add("GenericLinearAlgebra")` to install the GenericLinearAlgebra package. -Stacktrace: - [1] macro expansion - @ ./loading.jl:2223 [inlined] - [2] macro expansion - @ ./lock.jl:273 [inlined] - [3] __require(into::Module, mod::Symbol) - @ Base ./loading.jl:2198 - [4] #invoke_in_world#3 - @ ./essentials.jl:1089 [inlined] - [5] invoke_in_world - @ ./essentials.jl:1086 [inlined] - [6] require(into::Module, mod::Symbol) - @ Base ./loading.jl:2191 julia> roots(p) -ERROR: MethodError: no method matching eigvals!(::Matrix{BigFloat}) -The function `eigvals!` exists, but no method is defined for this combination of argument types. - -Closest candidates are: - eigvals!(::AbstractMatrix{T}, !Matched::LinearAlgebra.Cholesky{T}; sortby) where T<:Number - @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/symmetriceigen.jl:260 - eigvals!(::StridedMatrix{T}, !Matched::LinearAlgebra.LU{T, <:StridedMatrix{T} where T}; sortby) where T - @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/symmetriceigen.jl:286 - eigvals!(!Matched::LinearAlgebra.SymTridiagonal{<:Union{Float32, Float64}, <:StridedVector{T} where T}) - @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/tridiag.jl:270 - ... - -Stacktrace: - [1] eigvals(A::Matrix{BigFloat}; kws::@Kwargs{}) - @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/eigen.jl:343 - [2] eigvals(A::Matrix{BigFloat}) - @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/eigen.jl:343 - [3] roots(p::Polynomial{BigFloat, :x}; kwargs::@Kwargs{}) - @ Polynomials ~/julia/Polynomials/src/polynomials/standard-basis/standard-basis.jl:500 - [4] roots(p::Polynomial{BigFloat, :x}) - @ Polynomials ~/julia/Polynomials/src/polynomials/standard-basis/standard-basis.jl:481 - [5] top-level scope - @ none:1 +3-element Vector{Complex{BigFloat}}: + 0.9999999999999999999999999999999999999999999999999999999999999999999999999999655 + 0.0im + 1.999999999999999999999999999999999999999999999999999999999999999999999999999931 - 0.0im + 2.999999999999999999999999999999999999999999999999999999999999999999999999999793 + 0.0im ``` ### Comments on root finding diff --git a/src/common.jl b/src/common.jl index d7c9eadd..5f350cc3 100644 --- a/src/common.jl +++ b/src/common.jl @@ -60,9 +60,11 @@ 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) +julia> Polynomials.fromroots(r) Polynomial(6 - 5*x + x^2) ``` """ @@ -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) ``` diff --git a/src/polynomials/factored_polynomial.jl b/src/polynomials/factored_polynomial.jl index b9f462d3..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,7 +31,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=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)) 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..1aea6e43 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,11 @@ 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 +Stacktrace: + [1] integrate(p::LaurentPolynomial{Int64, :x}) + @ Polynomials ~/julia/Polynomials/src/polynomials/standard-basis/standard-basis.jl:174 + [2] top-level scope + @ none:1 julia> integrate(P([1,1,1], -5)) LaurentPolynomial(-0.25*x⁻⁴ - 0.3333333333333333*x⁻³ - 0.5*x⁻²) @@ -184,6 +191,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 +232,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]) From 844ac2c187b8c1e1b2686b71f0ed907aa4900f6d Mon Sep 17 00:00:00 2001 From: jverzani Date: Fri, 15 Nov 2024 21:27:51 -0500 Subject: [PATCH 4/6] docfix --- docs/src/index.md | 63 ++----------------- .../standard-basis/laurent-polynomial.jl | 6 +- 2 files changed, 5 insertions(+), 64 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 51c0dcf5..6d27da12 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -95,13 +95,7 @@ Polynomial(1 + 2*s + 3*s^2) julia> p + q ERROR: ArgumentError: Polynomials have different indeterminates -Stacktrace: - [1] assert_same_variable(X::Symbol, Y::Symbol) - @ Polynomials ~/julia/Polynomials/src/common.jl:522 - [2] +(p::Polynomial{Int64, :x}, q::Polynomial{Int64, :s}) - @ Polynomials ~/julia/Polynomials/src/common.jl:1114 - [3] top-level scope - @ none:1 +[...] ``` Except for operations involving constant polynomials. @@ -195,28 +189,7 @@ Polynomial(-6.0 + 11.0*x - 6.0*x^2 + 1.0*x^3) julia> roots(p) ERROR: MethodError: no method matching eigvals!(::Matrix{BigFloat}) -The function `eigvals!` exists, but no method is defined for this combination of argument types. - -Closest candidates are: - eigvals!(::AbstractMatrix{T}, !Matched::LinearAlgebra.Cholesky{T}; sortby) where T<:Number - @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/symmetriceigen.jl:260 - eigvals!(::StridedMatrix{T}, !Matched::LinearAlgebra.LU{T, <:StridedMatrix{T} where T}; sortby) where T - @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/symmetriceigen.jl:286 - eigvals!(!Matched::LinearAlgebra.SymTridiagonal{<:Union{Float32, Float64}, <:StridedVector{T} where T}) - @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/tridiag.jl:270 - ... - -Stacktrace: - [1] eigvals(A::Matrix{BigFloat}; kws::@Kwargs{}) - @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/eigen.jl:343 - [2] eigvals(A::Matrix{BigFloat}) - @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/eigen.jl:343 - [3] roots(p::Polynomial{BigFloat, :x}; kwargs::@Kwargs{}) - @ Polynomials ~/julia/Polynomials/src/polynomials/standard-basis/standard-basis.jl:500 - [4] roots(p::Polynomial{BigFloat, :x}) - @ Polynomials ~/julia/Polynomials/src/polynomials/standard-basis/standard-basis.jl:481 - [5] top-level scope - @ none:1 +[...] julia> using GenericLinearAlgebra @@ -678,23 +651,7 @@ However, as there would be an ambiguous outcome of the following ```jldoctest natural_inclusion julia> [one(p) one(q)] ERROR: ArgumentError: Polynomials have different indeterminates -Stacktrace: - [1] assert_same_variable(X::Symbol, Y::Symbol) - @ Polynomials ~/julia/Polynomials/src/common.jl:522 - [2] promote_rule(::Type{Polynomial{Int64, :x}}, ::Type{Polynomial{Int64, :y}}) - @ Polynomials ~/julia/Polynomials/src/promotions.jl:24 - [3] promote_type - @ ./promotion.jl:318 [inlined] - [4] promote_eltypeof(v1::Polynomial{Int64, :x}, v2::Polynomial{Int64, :y}) - @ Base ./abstractarray.jl:1618 - [5] _cat(::Val{2}, ::Polynomial{Int64, :x}, ::Polynomial{Int64, :y}) - @ Base ./abstractarray.jl:1820 - [6] cat - @ ./abstractarray.jl:2084 [inlined] - [7] hcat(::Polynomial{Int64, :x}, ::Polynomial{Int64, :y}) - @ Base ./abstractarray.jl:1968 - [8] top-level scope - @ none:1 +[...] ``` an error is thrown. @@ -726,19 +683,7 @@ but not this one: ```jldoctest natural_inclusion julia> [one(p), one(p) + one(q)] ERROR: ArgumentError: Polynomials have different indeterminates -Stacktrace: - [1] assert_same_variable(X::Symbol, Y::Symbol) - @ Polynomials ~/julia/Polynomials/src/common.jl:522 - [2] promote_rule(::Type{Polynomial{Int64, :x}}, ::Type{Polynomial{Int64, :y}}) - @ Polynomials ~/julia/Polynomials/src/promotions.jl:24 - [3] promote_type - @ ./promotion.jl:318 [inlined] - [4] promote_typeof(x::Polynomial{Int64, :x}, y::Polynomial{Int64, :y}) - @ Base ./promotion.jl:378 - [5] vect(::Polynomial{Int64, :x}, ::Vararg{Any}) - @ Base ./array.jl:184 - [6] top-level scope - @ none:1 +[...] ``` Also, mixing types can result in unspecific symbols, as this example shows: diff --git a/src/polynomials/standard-basis/laurent-polynomial.jl b/src/polynomials/standard-basis/laurent-polynomial.jl index 1aea6e43..0779482b 100644 --- a/src/polynomials/standard-basis/laurent-polynomial.jl +++ b/src/polynomials/standard-basis/laurent-polynomial.jl @@ -50,11 +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 -Stacktrace: - [1] integrate(p::LaurentPolynomial{Int64, :x}) - @ Polynomials ~/julia/Polynomials/src/polynomials/standard-basis/standard-basis.jl:174 - [2] top-level scope - @ none:1 +[...] julia> integrate(P([1,1,1], -5)) LaurentPolynomial(-0.25*x⁻⁴ - 0.3333333333333333*x⁻³ - 0.5*x⁻²) From 625695ee117d730db4b25c0b9846d5eae039dace Mon Sep 17 00:00:00 2001 From: jverzani Date: Sat, 16 Nov 2024 07:30:43 -0500 Subject: [PATCH 5/6] allow test with FactoredPolynomial --- src/common.jl | 2 +- test/StandardBasis.jl | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/common.jl b/src/common.jl index 5f350cc3..1ab01959 100644 --- a/src/common.jl +++ b/src/common.jl @@ -1336,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/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 From a699d854ced60d558717d0380472edd200416c59 Mon Sep 17 00:00:00 2001 From: jverzani Date: Sat, 16 Nov 2024 07:47:32 -0500 Subject: [PATCH 6/6] remove cruft --- src/common.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/common.jl b/src/common.jl index 1ab01959..a7c2a8e5 100644 --- a/src/common.jl +++ b/src/common.jl @@ -64,7 +64,7 @@ julia> using Polynomials julia> r = [3, 2]; # (x - 3)(x - 2) -julia> Polynomials.fromroots(r) +julia> fromroots(r) Polynomial(6 - 5*x + x^2) ``` """