diff --git a/Project.toml b/Project.toml
index 2493809f..e7551c07 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,7 +2,7 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "4.0.14"
+version = "4.0.15"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/polynomials/multroot.jl b/src/polynomials/multroot.jl
index 260e6d1d..aad8884c 100644
--- a/src/polynomials/multroot.jl
+++ b/src/polynomials/multroot.jl
@@ -145,7 +145,6 @@ function pejorative_manifold(
     ϕ = 100,      # residual growth factor
     kwargs...
     )  where {T,X}
-
     S = float(T)
     u = convert(PnPolynomial{S,X}, coeffs0(p))
     nu₂ = norm(u, 2)
@@ -156,9 +155,9 @@ function pejorative_manifold(
         atol  = ρ2, rtol  = zero(ρ2))
     ρⱼ /= nu₂
 
+    hasnan(v) && throw(ArgumentError("NaN in reduced polynomial"))
     # root approximations
     zs = roots(v)
-
     # recover multiplicities
     ls = pejorative_manifold_multiplicities(Val(method),
                                             u, v, w,
@@ -206,7 +205,6 @@ function pejorative_manifold_multiplicities(
         end
 
     end
-
     ls
 
 end
@@ -225,6 +223,7 @@ function pejorative_manifold_multiplicities(
 
     dv = derivative(v)
     ls = w.(zs) ./ dv.(zs)
+
     ls = round.(Int, real.(ls))
 
     return ls
@@ -258,7 +257,6 @@ function pejorative_root(p, zs::Vector{S}, ls;
 
     ## Solve WJ Δz = W(Gl(z) - a)
     ## 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]
diff --git a/src/polynomials/ngcd.jl b/src/polynomials/ngcd.jl
index 514b2736..b9e971a8 100644
--- a/src/polynomials/ngcd.jl
+++ b/src/polynomials/ngcd.jl
@@ -49,7 +49,6 @@ function ngcd(p::P, q::Q,
     p′ = P′{R,X}(ps[nz:end])
     q′ = P′{R,X}(qs[nz:end])
     out = NGCD.ngcd(p′, q′, args...; kwargs...)
-
     ## convert to original polynomial type
     𝑷 = Polynomials.constructorof(P){R,X}
     u,v,w = convert.(𝑷, (out.u,out.v,out.w))
diff --git a/src/rational-functions/common.jl b/src/rational-functions/common.jl
index dae8bb5a..a6e85c50 100644
--- a/src/rational-functions/common.jl
+++ b/src/rational-functions/common.jl
@@ -427,15 +427,22 @@ end
 
 ## ---- zeros, poles, ...
 """
-    poles(pq::AbstractRationalFunction; method=:numerical, kwargs...)
+    poles(pq::AbstractRationalFunction;
+        method=:numerical, multroot_method=:direct, kwargs...)
 
 For a rational function `p/q`, first reduces to normal form, then finds the roots and multiplicities of the resulting denominator.
 
+* `method` is used to pass to `lowest_terms`
+* `multroot_method` is passed to the method argument of `multroot`, which can be `:direct` (the faster default) or `:iterative` (the slower, and possibly more robust alternate)
+
 """
-function poles(pq::AbstractRationalFunction; method=:numerical,  kwargs...)
+function poles(pq::AbstractRationalFunction;
+               method=:numerical, # for lowest_terms
+               multroot_method=:direct, # or :iterative
+               kwargs...)
     pq′ = lowest_terms(pq; method=method, kwargs...)
     den = denominator(pq′)
-    mr = Multroot.multroot(den)
+    mr = Multroot.multroot(den; method=multroot_method)
     (zs=mr.values, multiplicities = mr.multiplicities)
 end
 
diff --git a/test/StandardBasis.jl b/test/StandardBasis.jl
index 2c69bc73..cf9b5c22 100644
--- a/test/StandardBasis.jl
+++ b/test/StandardBasis.jl
@@ -1027,6 +1027,20 @@ end
     out = Polynomials.Multroot.multroot(Polynomials.Polynomial(pb))
     @test out.values ≈ [1.0] && out.multiplicities == [2]
 
+    ## Issue #592
+    p = Polynomial([NaN,NaN,NaN,NaN,NaN,NaN,NaN])
+    @test_throws ArgumentError roots(p//p)
+
+    ## issue #593
+    numcoeffs = Complex{BigFloat}[-0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, -6.4095e+257 - 3.314e+243im, -5.6118e+244 - 3.9514e+230im, -4.0102e+236 - 2.0735e+222im, -3.0725e+223 - 2.1631e+209im, -1.0975e+215 - 5.6751e+200im, -7.2072e+201 - 5.0751e+187im, -1.7167e+193 - 8.877e+178im, -9.3936e+179 - 6.6134e+165im, -1.678e+171 - 8.6782e+156im, -7.3448e+157 - 5.1693e+143im, -1.0499e+149 - 5.43e+134im, -3.4468e+135 - 2.4259e+121im, -4.1052e+126 - 2.1231e+112im, -8.9846e+112 - 6.3225e+98im, -9.1725e+103 - 4.742e+89im, -1.0036e+90 - 7.0656e+75im, -8.9646e+80 - 4.6354e+66im]
+    dencoeffs = Complex{BigFloat}[0.0 + 0.0im, 0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, -3.0837e+258 + 0.0im, -3.5997e+245 - 9.309e+230im, -1.9292e+237 - 7.0217e+217im, -1.9709e+224 - 5.0951e+209im, -5.2803e+215 - 3.2932e+196im, -4.6214e+202 - 1.1952e+188im, -8.2587e+193 - 6.4383e+174im, -6.0237e+180 - 1.5577e+166im, -8.0763e+171 - 6.7125e+152im, -4.711e+158 - 1.2183e+144im, -5.0507e+149 - 3.937e+130im, -2.2106e+136 - 5.7173e+121im, -1.9752e+127 - 1.2316e+108im, -5.7628e+113 - 1.4904e+99im, -4.4151e+104 - 1.6057e+85im, -6.4393e+90 - 1.6651e+76im, -4.3167e+81 + 0.0im]
+    p = Polynomial(numcoeffs)
+    q = Polynomial(dencoeffs)
+    r = p//q
+
+    @test_throws ArgumentError poles(r)
+    out = poles(r; multroot_method=:iterative)
+    @test out.multiplicities == [3]
 end
 
 @testset "critical points" begin