error with NaN; allow method to pass through (#594)

Address issues #592 and #593 


* close #592 with an error message
* add option for `multroot` method to calling `poles` function. This is a partial fix. The underlying algorithm is designed with Float64 in mind, and the example is using BigFloats. Might make more sense to have specific choices based on type of element, but that needs some thinking through.

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"

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]

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))

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
 

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