Error with `Polynomial(Polynomial)`

[ParadaCarleton]

Iterated polynomial application seems to be broken, but only for Polynomials with different types. MWE:

julia> f = Polynomial((1, -1, -1.1))
Polynomial(1 - x - 1.1*x^2)

julia> f(f)
ERROR: MethodError: no method matching one(::Type{Polynomial{Float64, 3}}, ::Int64)

Closest candidates are:
  one(::P, ::Any) where P<:AbstractPolynomial
   @ Polynomials ~/.julia/packages/Polynomials/tVgyD/src/common.jl:781
  one(::Type{P}) where P<:Polynomials.StandardBasisPolynomial
   @ Polynomials ~/.julia/packages/Polynomials/tVgyD/src/polynomials/standard-basis.jl:65
  one(::Type{P}) where P<:AbstractPolynomial
   @ Polynomials ~/.julia/packages/Polynomials/tVgyD/src/common.jl:779
  ...

Stacktrace:
 [1] one(p::Polynomial{Float64, 3}, var::Int64)
   @ Polynomials ~/.julia/packages/Polynomials/tVgyD/src/common.jl:781
 [2] one(p::Polynomial{Float64, 3})
   @ Polynomials ~/.julia/packages/Polynomials/tVgyD/src/common.jl:781
 [3] _one(x::Polynomial{Float64, 3})
   @ Polynomials.EvalPoly ~/.julia/packages/Polynomials/tVgyD/src/contrib.jl:126
 [4] _evalpoly(x::Polynomial{Float64, 3}, p::Vector{Float64})
   @ Polynomials.EvalPoly ~/.julia/packages/Polynomials/tVgyD/src/contrib.jl:56
 [5] evalpoly(x::Polynomial{Float64, 3}, p::Vector{Float64})
   @ Polynomials.EvalPoly ~/.julia/packages/Polynomials/tVgyD/src/contrib.jl:50
 [6] evalpoly(x::Polynomial{Float64, 3}, p::Polynomial{Float64, 3})
   @ Polynomials ~/.julia/packages/Polynomials/tVgyD/src/polynomials/standard-basis.jl:43
 [7] (::Polynomial{Float64, 3})(x::Polynomial{Float64, 3})
   @ Polynomials ~/.julia/packages/Polynomials/tVgyD/src/abstract.jl:122
 [8] top-level scope
   @ REPL[15]:1

Also relevant is that this seems to be type-unstable:

julia> typeof(f(f))
Polynomial{Any, :x}

This could probably be fixed by promoting all the parameters to a common type.

[jverzani]

I can't replicate, though I haven't considered doing this before, so it may be broken. But this seems to be off:

julia> f = Polynomial((1, -1, -1.1))
Polynomial(1.0 - 1.0*3 - 1.1*3^2)

(As in where does 3 come from?)

For me, I see this:

julia> f = Polynomial((1, -1, -1.1))
Polynomial(1 - x - 1.1*x^2)

julia> f(f)
Polynomial(-1.1 + 3.2*x + 2.420000000000001*x^2 - 2.4200000000000004*x^3 - 1.3310000000000004*x^4)

[ParadaCarleton]

I honestly have no idea where the 3 came from :sweat_smile:

Retried like this, which makes it more clear that the problem is with a type instability:

julia> f = Polynomial((-1.1, 0, 1))
Polynomial(-1.1 + x^2)

julia> f(f)
Polynomial(0.1100000000000001 - 2.2*x^2 + 1.0*x^4)

julia> roots(f(f))
ERROR: MethodError: no method matching one(::Type{Any}

[jverzani]

Yeah, this shows the issue. Thanks. Let me look to see what needs to be done.

[jverzani]

Yeah, the issue is the use of a tuple to define the coefficients. Had you used (1.1, 0.0, 1.0) this wouldn't have popped up. Using a tuple to define coefficients isn't discouraged, but I'm not sure a special code path to accommodate that is warranted.

Curious though, what is the natural question where iteration of polynomials is part of the solution?

[ParadaCarleton]

Curious though, what is the natural question where iteration of polynomials is part of the solution?

Dynamical systems; iteration is a common way to get chaotic behavior. In this case, I happened to be making a bifurcation plot for the quadratic map.

[jverzani]

Yes, nice plot. There isn’t really savings using a tuple here, as they are converted to a vector in construction. But that shouldn’t be a vector of Real, rather Float64. Let me look.

[ParadaCarleton]

Yeah, the issue is the use of a tuple to define the coefficients. Had you used (1.1, 0.0, 1.0) this wouldn't have popped up. Using a tuple to define coefficients isn't discouraged, but I'm not sure a special code path to accommodate that is warranted.

I think adding an extra convenience method here to convert is probably the best approach, given it's a one-liner. But either that or throwing an error for Tuples is fine--the problem is the silent type instability.

[jverzani]

The problem sits with collect which promotes this tuple to Real in this case, not the tightest type possible. I am going to push a fix that calls promote on the tuple argument prior to collection. Just waiting to see what that breaks, if anything. Maybe I need to be educated about a base function that identifies the minimal eltype of a collection.