andreasnoack / FastPolynomialRoots.jl

Fast and backward stable computation of roots of polynomials in Julia
MIT License
15 stars 4 forks source link

Missing working methods after `Polynomials.roots()` is overwritten #19

Closed johanbluecreek closed 2 years ago

johanbluecreek commented 4 years ago

With the original Polynomials.roots() there are several working methods supporting polynomials with float and integer coefficients:

julia> using GenericLinearAlgebra

julia> using Polynomials

julia> Polynomial(Float64[1,10,100,1000]) |> roots
3-element Array{Complex{Float64},1}:
  -0.09999999999999999 + 0.0im
 -2.42861286636753e-17 - 0.1im
 -2.42861286636753e-17 + 0.1im

julia> Polynomial(BigFloat[1,10,100,1000]) |> roots |> x -> Complex{Float64}.(x)
3-element Array{Complex{Float64},1}:
                   -0.1 + 0.0im
 4.3180842775472223e-78 - 0.1im
 4.3180842775472223e-78 + 0.1im

julia> Polynomial(Int64[1,10,100,1000]) |> roots
3-element Array{Complex{Float64},1}:
  -0.09999999999999999 + 0.0im
 -2.42861286636753e-17 - 0.1im
 -2.42861286636753e-17 + 0.1im

julia> Polynomial(BigInt[1,10,100,1000]) |> roots |> x -> Complex{Float64}.(x)
3-element Array{Complex{Float64},1}:
                   -0.1 + 0.0im
 4.3180842775472223e-78 - 0.1im
 4.3180842775472223e-78 + 0.1im

when roots is overwritten by FastPolynomialRoots, the methods for integers break:

julia> using FastPolynomialRoots

julia> Polynomial(Float64[1,10,100,1000]) |> roots
3-element Array{Complex{Float64},1}:
 1.6219664500383146e-16 + 0.1000000000000001im
 1.6219664500383146e-16 - 0.1000000000000001im
   -0.09999999999999977 + 0.0im

julia> Polynomial(BigFloat[1,10,100,1000]) |> roots |> x -> Complex{Float64}.(x)
3-element Array{Complex{Float64},1}:
                   -0.1 + 0.0im
 4.3180842775472223e-78 - 0.1im
 4.3180842775472223e-78 + 0.1im

julia> Polynomial(Int64[1,10,100,1000]) |> roots
ERROR: MethodError: no method matching rootsFastPolynomialRoots(::Polynomial{Float64})
Closest candidates are:
  rootsFastPolynomialRoots(::Array{Complex{Float64},1}) at /home/user/.julia/packages/FastPolynomialRoots/Meom5/src/FastPolynomialRoots.jl:29
  rootsFastPolynomialRoots(::Array{Float64,1}) at /home/user/.julia/packages/FastPolynomialRoots/Meom5/src/FastPolynomialRoots.jl:10
Stacktrace:
 [1] roots(::Polynomial{Int64}) at /home/user/.julia/packages/FastPolynomialRoots/Meom5/src/FastPolynomialRoots.jl:6
 [2] |>(::Polynomial{Int64}, ::typeof(roots)) at ./operators.jl:823
 [3] top-level scope at REPL[10]:1

julia> Polynomial(BigInt[1,10,100,1000]) |> roots |> x -> Complex{Float64}.(x)
ERROR: MethodError: no method matching rootsFastPolynomialRoots(::Polynomial{Float64})
Closest candidates are:
  rootsFastPolynomialRoots(::Array{Complex{Float64},1}) at /home/user/.julia/packages/FastPolynomialRoots/Meom5/src/FastPolynomialRoots.jl:29
  rootsFastPolynomialRoots(::Array{Float64,1}) at /home/user/.julia/packages/FastPolynomialRoots/Meom5/src/FastPolynomialRoots.jl:10
Stacktrace:
 [1] roots(::Polynomial{BigInt}) at /home/user/.julia/packages/FastPolynomialRoots/Meom5/src/FastPolynomialRoots.jl:6
 [2] |>(::Polynomial{BigInt}, ::typeof(roots)) at ./operators.jl:823
 [3] top-level scope at REPL[11]:1

The expected behavior I would say is that the methods that were working before to still be working after FastPolynomialRoots has been used.

johanbluecreek commented 4 years ago

This seems to do it

diff --git a/src/FastPolynomialRoots.jl b/src/FastPolynomialRoots.jl
index 87f7baf..2077096 100644
--- a/src/FastPolynomialRoots.jl
+++ b/src/FastPolynomialRoots.jl
@@ -3,7 +3,7 @@ module FastPolynomialRoots
 using LibAMVW_jll, Polynomials

 Polynomials.roots(p::Union{Polynomial{Float64},Polynomial{Complex{Float64}}}) = rootsFastPolynomialRoots(coeffs(p))
-Polynomials.roots(p::Polynomial{<:Integer}) = rootsFastPolynomialRoots(convert(Polynomial{Float64}, p))
+Polynomials.roots(p::Polynomial{<:Integer}) = rootsFastPolynomialRoots(convert(Vector{Float64}, coeffs(p)))

 function rootsFastPolynomialRoots(a::Vector{Float64})

since then (continued session from above)

julia> Polynomials.roots(p::Polynomial{<:Integer}) = FastPolynomialRoots.rootsFastPolynomialRoots(convert(Vector{Float64}, coeffs(p)))

julia> Polynomial(Int64[1,10,100,1000]) |> roots
3-element Array{Complex{Float64},1}:
 1.6219664500383146e-16 + 0.1000000000000001im
 1.6219664500383146e-16 - 0.1000000000000001im
   -0.09999999999999977 + 0.0im

julia> Polynomial(BigInt[1,10,100,1000]) |> roots
3-element Array{Complex{Float64},1}:
 1.6219664500383146e-16 + 0.1000000000000001im
 1.6219664500383146e-16 - 0.1000000000000001im
   -0.09999999999999977 + 0.0im

or maybe BigInt should be expected to be using the original Polynomials.roots with GenericLinearAlgebra's eigvals?

andreasnoack commented 4 years ago

Usually, we convert BigInts to BigFloats. I think we should only overwrite methods for which this package has implementations so I suspect that the Polynomials.roots(p::Polynomial{<:Integer}) method should be narrowed to just Int.

Do you know how Roots solves the BigFloat case right now?

johanbluecreek commented 4 years ago

I think we should only overwrite methods for which this package has implementations so I suspect that the Polynomials.roots(p::Polynomial{<:Integer}) method should be narrowed to just Int.

I'd agree.

Do you know how Roots solves the BigFloat case right now?

I'm quite sure that the original Polynomials.roots (if that is what you mean by "Roots") solves the problem by deriving the companion matrix and sending it to eigvals. Since if I do not have GenericLinearAlgebra loaded the call with a Polynomial{BigInt} or BigFloat just gives me e.g. ERROR: MethodError: no method matching eigvals!(::Array{BigFloat,2}).