one part of #602 (#603)

* one part of #602

* add deflation to ngcd

* move where deflation happens

* deflate early

* bump compat

* bump compat

* slight cleanup

Author: jverzani

docs/Project.toml

@@ -8,4 +8,3 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 Documenter = "1"
-GR_jll = "< 0.58"

src/polynomials/multroot.jl

@@ -99,9 +99,15 @@ function multroot(p::StandardBasisPolynomial{T}; verbose=false,
     cs = coeffs0(p)
     length(cs) <= 1 && return (values=T[], multiplicities=Int[], κ=NaN, ϵ=NaN)
 
-    # degenerate case, all zeros
-    if (nz = findfirst(!iszero, cs)) == length(cs)
+    # deflate leading zeros?
+    i = findfirst(!Polynomials.iscoeffzero,cs)
+    if isnothing(i)
+        # degenerate case, all zeros
+        throw(ArgumentError("Zero polynomial has no roots defined"))
+    elseif i == length(cs)
         return (values=zeros(T,1), multiplicities=[nz-1], κ=NaN, ϵ=NaN)
+    elseif i > 1
+        p = Polynomial(cs[i:end])
     end
 
     # Basic algorithm is two steps
@@ -110,6 +116,11 @@ function multroot(p::StandardBasisPolynomial{T}; verbose=false,
     κ, ϵ = stats(p, z̃, l)
 
     verbose && show_stats(κ, ϵ)
+    if i > 1
+        push!(z̃, zero(T))
+        push!(l, i-1)
+    end
+
 
     (values = z̃, multiplicities = l, κ = κ, ϵ = ϵ)
 

src/polynomials/ngcd.jl

@@ -13,26 +13,18 @@ function ngcd(p::P, q::Q,
               kwargs...) where {T,X,P<:StandardBasisPolynomial{T,X},
                                 S,Y,Q<:StandardBasisPolynomial{S,Y}}
 
+
     # easy cases
     degree(p) < 0  && return (u=q,      v=p, w=one(q),  θ=NaN, κ=NaN)
     degree(p) == 0 && return (u=one(q), v=p, w=q,       θ=NaN, κ=NaN)
     degree(q) < 0  && return (u=one(q), v=p, w=zero(q), θ=NaN, κ=NaN)
     degree(q) == 0 && return (u=one(p), v=p, w=q,       Θ=NaN, κ=NaN)
 
+
     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...)
-    end
-
-    # other easy cases
-    p ≈ q          && return (u=p,v=one(p),  w=one(p),  θ=NaN, κ=NaN)
-    Polynomials.assert_same_variable(p,q)
-
-
 
     R = promote_type(float(T))
     𝑷 = Polynomials.constructorof(P){R,X}
@@ -45,14 +37,28 @@ function ngcd(p::P, q::Q,
     if nz == length(qs)
         x = variable(p)
         u = x^(nz-1)
-        v,w = 𝑷(ps[nz:end]), 𝑷(qs[nz:end])
+        ps,qs = ps[nz:end], qs[nz:end]
+        v,w = 𝑷(ps), 𝑷(qs)
         return (u=u, v=v, w=w, Θ=NaN, κ=NaN)
+    elseif 1 <= nz < length(qs)
+        ps,qs = ps[nz:end], qs[nz:end]
+        p,q = P(ps), Q(qs)
+    end
+
+    if degree(p) > 5*(1+degree(q))
+        a,b = divrem(p,q)
+        return ngcd(q, b, args...; λ=100,  kwargs...)
     end
 
+    # other easy cases
+    p ≈ q          && return (u=p,v=one(p),  w=one(p),  θ=NaN, κ=NaN)
+    Polynomials.assert_same_variable(p,q)
+
     ## call ngcd
     P′ = PnPolynomial
-    p′ = P′{R,X}(ps[nz:end])
-    q′ = P′{R,X}(qs[nz:end])
+    p′ = P′{R,X}(ps)
+    q′ = P′{R,X}(qs)
+
     out = NGCD.ngcd(p′, q′, args...; kwargs...)
     ## convert to original polynomial type
     𝑷 = Polynomials.constructorof(P){R,X}

src/rational-functions/common.jl

@@ -458,9 +458,10 @@ function roots(pq::AbstractRationalFunction{T};
                multroot_method=nothing, # :direct or:iterative
                kwargs...) where {T}
     pq′ = lowest_terms(pq; method=method, kwargs...)
-    den = numerator(pq′)
+    num = numerator(pq′)
+    (hasnan(num) || any(isinf, num)) && throw(ArgumentError("Reduced expression has numeric issues"))
     mmethod = something(multroot_method, default_multroot_method(T))
-    mr = Multroot.multroot(den; method=mmethod)
+    mr = Multroot.multroot(num; method=mmethod)
     (zs=mr.values, multiplicities = mr.multiplicities)
 end
 default_multroot_method(::Type{T}) where {T<:Union{BigFloat, Complex{BigFloat}, BigInt, Complex{BigInt}}} = :iterative

test/StandardBasis.jl

@@ -1046,6 +1046,7 @@ end
     multroot_method = :direct # fails if direct
     @test_throws ArgumentError poles(r; multroot_method)
     @test_throws ArgumentError roots(r; multroot_method)
+
 end
 
 @testset "critical points" begin

test/rational-functions.jl

@@ -124,6 +124,19 @@ end
     end
     @test norm(numerator(lowest_terms(d - pq)), Inf) <= sqrt(eps())
 
+        ## Issue #602
+    s = Polynomial([0,1],:s)
+    r = (15s^14 + 1e-16s^15)//(1s + 14s^14)
+    out = roots(r)
+    @test 0 ∈ out.zs
+    @test out.multiplicities == [1,13]
+
+    r = (15s^14 + 1e-16s^15)//(s)
+    out = roots(r)
+    @test 0 ∈ out.zs
+    @test out.multiplicities == [1,13]
+
+
     ## issue #605
     s = x//1
     X = 1//(s^4*(s+2))
@@ -145,7 +158,6 @@ end
         end
     end
     @test X ≈ Y
-
 end
 
 @testset "As matrix elements" begin