Ngcd cleanup

src/polynomials/multroot.jl

@@ -133,7 +133,7 @@ function pejorative_manifold(p::Polynomials.StandardBasisPolynomial{T};
 
         for z in rts
             tmp, ind = findmin(abs.(zs .- z))
-            ls[ind] = ls[ind] + 1
+            ls[ind] += 1
         end
 
     end
@@ -212,7 +212,7 @@ function pejorative_root(p, zs::Vector{S}, ls;
         @info ("""
 The multiplicity count may be in error: the initial guess for the roots failed
 to converge to a pejorative root.
-""")
+        """)
         return(zₖs)
     end
 
@@ -268,8 +268,8 @@ function cond_zl(p, zs::Vector{S}, ls) where {S}
     J = zeros(S, sum(ls), length(zs))
     W = diagm(0 => [min(1, 1/abs(aᵢ)) for aᵢ in p[2:end]])
     evalJ!(J, zs, ls)
-    F = qr(W*J)
-    σ = Polynomials.NGCD.smallest_singular_value(F.R)
+    _,R = qr(W*J)
+    σ = Polynomials.NGCD.smallest_singular_value(UpperTriangular(R))
     1 / σ
 end
 

src/polynomials/ngcd.jl

@@ -204,8 +204,6 @@ function ngcd(p::PnPolynomial{T,X},
               minⱼ = -1
               ) where {T, X}
 
-    verbose = false
-
     m,n = length(p)-1, length(q)-1
     (m == 1 || n == 0) && return trivial_gcd(p, q)
 
@@ -254,20 +252,20 @@ function ngcd(p::PnPolynomial{T,X},
         xx = view(x, 1:nc)
         λ = norm(Q,Inf)*norm(R,Inf)
         σ₋₁ = smallest_singular_value!(xx, V, ρ *  sqrt(m - j + 1), λ*ϵₘ)
-        verbose && @show j, σ₋₁, ρ *  sqrt(m - j + 1)
+        # @show j, σ₋₁, ρ *  sqrt(m - j + 1)
 
         # Lemma 7.1: If (p,q) is w/in ϵ of P^k_{mn} then σ₋₁ < ϵ√(m-j+1)
         if σ₋₁ ≤ ρ *  sqrt(m - j + 1)
             # candidate for degree; refine u₀, vₒ, w₀ to see if ρ < ϵ
             if iszero(σ₋₁)
-                @info "Determinant is zero, which shouldn't be the case. Treat results with scrutiny"
+                # determinant is 0
                 u, v, w = initial_uvw(Val(:iszero), j, p, q, xx)
             else
                 u, v, w = initial_uvw(Val(:ispossible), j, p, q, xx)
             end
             ϵₖ, κ = refine_uvw!(u, v, w, p, q, uv, uw)
             ϵ = max(atol, npq₂ * κ * rtol)
-            verbose && @show ϵₖ, ϵ
+            # @show ϵₖ, ϵ
             if ϵₖ ≤ ϵ
                 return (u=u, v=v, w=w, Θ=ϵₖ, κ=κ)
             end
@@ -430,28 +428,29 @@ function initial_uvw(::Val{:ispossible}, j, p::P, q::Q, x) where {T,X,
 
 end
 
-# find u₎, v₀. w₀ when R is singular
+# find u₀, v₀. w₀ when R is singular.
 function initial_uvw(::Val{:iszero}, j, p::P, q::Q, x) where {T,X,
                                                               P<:PnPolynomial{T,X},
                                                               Q<:PnPolynomial{T,X}}
 
     m,n = length(p)-1, length(q)-1
-    S = [convmtx(p, n-j+1) convmtx(q, m-j+1)]
-
+    S = SylvesterMatrix(p,q,j)
     F = qr(S)
     R = UpperTriangular(F.R)
 
     if iszero(det(R))
         x .= eigvals(R)[:,1]
     else
-        x .= ones(T, size(R,2))
-        ldiv!(R', x)
-        x ./= norm(x,2)
-        ldiv!(R, x)
-        x ./= norm(x)
+        rand!(x)
+        smallest_singular_value!(x, R)
+#        x .= ones(T, size(R,2))
+#        ldiv!(R', x)
+#        x ./= norm(x,2)
+#        ldiv!(R, x)
+#        x ./= norm(x)
     end
 
-    w = P(x[1:n-j+1])
+    w = P(x[1:n-j+1]) # ordering of S is not interlaced
     v = P(-x[(n-j+2):end])
 
     u = solve_u(v,w,p,q,j)
@@ -735,7 +734,7 @@ end
 
 function SylvesterMatrix(p, q, j)
     m, n = length(p)-1, length(q) - 1
-    Sₓ = hcat(convmtx(p, n-j + 1 ),  convmtx(q, m-j + 1))
+    Sₓ = hcat(convmtx(p, n - j + 1 ),  convmtx(q, m - j + 1))
 end
 
 ## ----

test/StandardBasis.jl

@@ -1292,7 +1292,7 @@ end
     @test degree(gcd(a*d, b*d, atol=sqrt(eps()))) > 0
     @test  degree(gcd(a*d,b*d, method=:noda_sasaki)) == degree(d)
     @test_skip degree(gcd(a*d,b*d, method=:numerical)) == degree(d) # issues on some architectures
-    l,m,n = (5,5,5) # realiable, though for larger l,m,n only **usually** correct
+    l,m,n = (2,2,2) # realiable, though for larger l,m,n only **usually** correct.
     u,v,w = fromroots.(rand.((l,m,n)))
     @test degree(gcd(u*v, u*w, method=:numerical)) == degree(u)