one part of #602

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,30 @@ 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)
+
+    @show ps, qs
+
     ## 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

@@ -459,6 +459,7 @@ function roots(pq::AbstractRationalFunction{T};
                kwargs...) where {T}
     pq′ = lowest_terms(pq; method=method, kwargs...)
     den = numerator(pq′)
+    (hasnan(den) || any(isinf, den)) && throw(ArgumentError("Reduced expression has numeric issues"))
     mmethod = something(multroot_method, default_multroot_method(T))
     mr = Multroot.multroot(den; method=mmethod)
     (zs=mr.values, multiplicities = mr.multiplicities)

test/rational-functions.jl

@@ -124,6 +124,20 @@ 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)//(s)
+    num = numerator(Polynomials.lowest_terms(r))
+    @test length(roots(num)) == 14
+    @test_throws DimensionMismatch Polynomials.Multroot.multroot(num) # where to fix
+
+    r = (15s^14 + 1e-16s^15)//(1s + 14s^14)
+    num = numerator(Polynomials.lowest_terms(r))
+    @test length(roots(num)) == 14
+    @test_throws DimensionMismatch Polynomials.Multroot.multroot(num) # where to fix
+
+
 end
 
 @testset "As matrix elements" begin