no special case for coeffs (#528)

* no special case for `coeffs`

* adjust to change in coeffs

* make coeffs type stable for ImmutablePolynomial; push instability to internal coeffs0

Author: jverzani

src/abstract-polynomial.jl

@@ -100,15 +100,15 @@ function coeffs(p::AbstractLaurentUnivariatePolynomial)
     [p[i] for i ∈ firstindex(p):lastindex(p)]
 end
 
-# a₀, …, aₙ or error even for Laurent
+# copy of a₀, …, aₙ or error even for Laurent
 function coeffs0(p::AbstractUnivariatePolynomial)
-    a,b = firstindex(p), lastindex(p)
+    a,b = firstindex(p), degree(p)
     a < 0 && throw(ArgumentError("Polynomial has negative index terms. Use `pairs` instead."))
     [p[i] for i ∈ 0:b]
 end
 function coeffs0(p::AbstractDenseUnivariatePolynomial)
-    a,b = firstindex(p), lastindex(p)
-    iszero(a) && return p.coeffs
+    a,b = firstindex(p), degree(p)
+    iszero(a) && return copy(p.coeffs)
     [p[i] for i ∈ 0:b]
 end
 

src/polynomial-container-types/immutable-dense-polynomial.jl

@@ -102,6 +102,9 @@ function trim_trailing_zeros!!(cs::Tuple)
     xs
 end
 
+coeffs0(p::ImmutableDensePolynomial) = trim_trailing_zeros!!(p.coeffs)
+
+
 ## chop. Also, not type stable
 function Base.chop(p::ImmutableDensePolynomial{B,T,X,N};
                    rtol::Real = Base.rtoldefault(real(T)),
@@ -208,10 +211,6 @@ function degree(p::ImmutableDensePolynomial{B,T,X,N}) where {B,T,X,N}
     return i - 1
 end
 
-function coeffs(p::P) where {P <: ImmutableDensePolynomial}
-    trim_trailing_zeros!!(p.coeffs)
-end
-
 # zero, one
 Base.zero(::Type{<:ImmutableDensePolynomial{B,T,X}}) where {B,T,X} =
     ImmutableDensePolynomial{B,T,X,0}(())

src/polynomials/multroot.jl

@@ -7,7 +7,7 @@ using ..Polynomials
 using LinearAlgebra
 
 
-import ..Polynomials: PnPolynomial, StandardBasisPolynomial
+import ..Polynomials: PnPolynomial, StandardBasisPolynomial, trim_trailing_zeros!!, chop!!, coeffs0
 
 """
     multroot(p; verbose=false, method=:direct, kwargs...)
@@ -102,10 +102,11 @@ function multroot(p::StandardBasisPolynomial{T}; verbose=false,
                   kwargs...) where {T}
 
     # degenerate case, constant
-    degree(p) == 0 && return (values=T[], multiplicities=Int[], κ=NaN, ϵ=NaN)
+    cs = coeffs0(p)
+    length(cs) <= 1 && return (values=T[], multiplicities=Int[], κ=NaN, ϵ=NaN)
 
     # degenerate case, all zeros
-    if (nz = findfirst(!iszero, coeffs(p))) == length(coeffs(p))
+    if (nz = findfirst(!iszero, cs)) == length(cs)
         return (values=zeros(T,1), multiplicities=[nz-1], κ=NaN, ϵ=NaN)
     end
 
@@ -146,7 +147,7 @@ function pejorative_manifold(
     )  where {T,X}
 
     S = float(T)
-    u = convert(PnPolynomial{S,X}, p)
+    u = convert(PnPolynomial{S,X}, coeffs0(p))
     nu₂ = norm(u, 2)
     θ2, ρ2 =  θ * nu₂, ρ * nu₂
     u, v, w, ρⱼ, κ = Polynomials.ngcd(
@@ -243,7 +244,7 @@ This follows Algorithm 1 of [Zeng](https://www.ams.org/journals/mcom/2005-74-250
 """
 function pejorative_root(p::StandardBasisPolynomial,
                          zs::Vector{S}, ls; kwargs...) where {S}
-    ps = reverse(coeffs(p))
+    ps = reverse(coeffs0(p))
     pejorative_root(ps, zs, ls; kwargs...)
 end
 
@@ -259,7 +260,6 @@ function pejorative_root(p, zs::Vector{S}, ls;
     ## using weights min(1/|aᵢ|), i ≠ 1
 
     m,n = sum(ls), length(zs)
-
     # storage
     a = p[2:end]./p[1]     # a ~ (p[n-1], p[n-2], ..., p[0])/p[n]
     W = Diagonal([min(1, 1/abs(aᵢ)) for aᵢ in a])
@@ -377,6 +377,7 @@ function backward_error(p, z̃s::Vector{S}, ls) where {S}
 end
 
 function stats(p::AbstractPolynomial, zs, ls)
+    p = chop(p; atol=0, rtol=0)
     cond_zl(p, zs, ls), backward_error(p, zs, ls)
 end
 
@@ -408,7 +409,7 @@ function pejorative_manifold(
     error("Does this get called?")
 
     S = float(T)
-    u = PnPolynomial{S,X}(S.(coeffs(p)))
+    u = PnPolynomial{S,X}(coeffs0(p))
 
     nu₂ = norm(u, 2)
 

test/StandardBasis.jl

@@ -984,7 +984,7 @@ end
 end
 
 @testset "multroot" begin
-    @testset for P in (Polynomial, ImmutablePolynomial, SparsePolynomial)
+    @testset for P in (Polynomial, ImmutablePolynomial, SparsePolynomial, LaurentPolynomial)
         rts = [1.0, sqrt(2), sqrt(3)]
         ls = [2, 3, 4]
         x = variable(P{Float64})
@@ -995,7 +995,7 @@ end
         @test out.ϵ <= sqrt(eps())
         @test out.κ * out.ϵ < sqrt(eps())  # small forward error
         # one for which the multiplicities are not correctly identified
-        n = 3 # was 4?
+        n = P == SparsePolynomial ? 2 : 4 # can be much higher (was 4), but SparsePolynomial loses accuracy.
         q = p^n
         out = Polynomials.Multroot.multroot(q)
         @test (out.multiplicities == n*ls) || (out.κ * out.ϵ > sqrt(eps()))  # large  forward error, l misidentified