Ngcd bug (#333)

* multroot issue with zero as a root

Author: jverzani

Project.toml

@@ -2,14 +2,13 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "2.0.6"
+version = "2.0.7"
 
 [deps]
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 
-
 [compat]
 Intervals = "0.5, 1.0, 1.3"
 RecipesBase = "0.7, 0.8, 1"

src/polynomials/multroot.jl

@@ -79,12 +79,16 @@ is misidentified.
 function multroot(p::Polynomials.StandardBasisPolynomial{T}; verbose=false,
                   kwargs...) where {T}
 
+    # degenerate case, all zeros
+    if findfirst(!iszero, coeffs(p)) == length(coeffs(p))
+        return (values=zeros(T,1), multiplicities=nz, κ=NaN, ϵ=NaN)
+    end
+
     z, l = pejorative_manifold(p; kwargs...)
     z̃ = pejorative_root(p, z, l)
     κ, ϵ = stats(p, z̃, l)
 
     verbose && show_stats(κ, ϵ)
-
     (values = z̃, multiplicities = l, κ = κ, ϵ = ϵ)
 
 end

src/polynomials/ngcd.jl

@@ -11,7 +11,10 @@ function ngcd(p::P, q::Q,
               args...; kwargs...) where {T,X,P<:StandardBasisPolynomial{T,X},
                                          S,Y,Q<:StandardBasisPolynomial{S,Y}}
 
-    degree(q) > degree(p) && return ngcd(q,p,args...;kwargs...)
+    if (degree(q) > degree(p))
+        u,w,v,Θ,κ =  ngcd(q,p,args...;kwargs...)
+        return (u=u,v=v,w=w, Θ=Θ, κ=κ)
+    end
     if degree(p) > 5*(1+degree(q))
         a,b = divrem(p,q)
         return ngcd(q,b, args...; λ=100, kwargs...)
@@ -134,7 +137,7 @@ function ngcd(p::PnPolynomial{T,X},
               verbose=false,
               minⱼ = -1,
               λ = 1
-              ) where {T <: AbstractFloat, X}
+              ) where {T, X}
 
     m,n = length(p)-1, length(q)-1
     @assert m >= n
@@ -552,12 +555,12 @@ end
 
 
 function residual_error(p::P,q,uv,uw) where {T,X,P<:AbstractPolynomial{T,X}}
-    tot = zero(T)
+    tot = zero(real(T))
     for (pᵢ, uvᵢ) in zip(p,uv)
-        tot += (pᵢ-uvᵢ)^2
+        tot += norm(pᵢ-uvᵢ)^2
     end
     for (qᵢ, uwᵢ) in zip(q, uw)
-        tot += (qᵢ-uwᵢ)^2
+        tot += norm(qᵢ-uwᵢ)^2
     end
     sqrt(tot)
 end

test/StandardBasis.jl

@@ -1031,6 +1031,13 @@ end
         out3 = Polynomials.ngcd(p,q,3)
         @test out.Θ <= out1.Θ
         @test out.Θ <= out3.Θ
+
+        # check for correct output if degree p < degree q
+        x = variable(P{Float64})
+        p = -18.0 - 37.0*x - 54.0*x^2 - 36.0*x^3 - 16.0*x^4
+        q = 2.0 + 5.0*x + 8.0*x^2 + 7.0*x^3 + 4.0*x^4 + 1.0*x^5
+        out = Polynomials.ngcd(p,q)
+        @test out.u * out.v ≈ p
     end
 end