tweak defaults with ngcd with one euclidean step

Author: jverzani

src/polynomials/ngcd.jl

@@ -14,19 +14,37 @@ function ngcd(p::P, q::Q, args...;kwargs...) where {T, S, P<:StandardBasisPolyno
 
     p′,q′ = promote(p,q)
     u,v,w,Θ,κ = NGCD.ngcd(coeffs(p′), coeffs(q′), args...; kwargs...)
+
     PP = ⟒(typeof(p′))
     (u=PP(u, p.var), v=PP(v, p.var), w = PP(w, p.var), Θ=Θ, κ=κ)
+    
 end
 
 """
     ngcd′(p,q)
 
 When degree(p) ≫ degree(q), this uses a early call to `divrem`.
 """
-function ngcd′(p::P, q::P; kwargs...) where {P <: StandardBasisPolynomial}
+function ngcd′(p::P, q::P;
+               atol = eps(real(float(T))),
+               # rtol=eps(real(float(T))),
+               rtol = Base.rtoldefault(real(float(T))), # relax this
+               satol= atol,
+               srtol= rtol,
+               kwargs...
+               ) where {T, P <: StandardBasisPolynomial{T}}
+
+
     a, b = divrem(p,q)
-    u,v,w,Θ, κ = ngcd(q, b; kwargs...)
-    ngcd(p, q, degree(u))
+
+    # check if a=u (p,q) ≈ (aq,q)
+    if isapprox(p, a*q, atol=atol, rtol=rtol)
+        return a
+    else
+        u,v,w,Θ, κ = ngcd(q, b; atol=100atol, rtol=100rtol, kwargs...)
+        u
+        #ngcd(p, q, degree(u))
+    end
 end
 
 
@@ -111,10 +129,10 @@ Note: Based on work by Andreas Varga
 function ngcd(ps::Vector{T},
               qs::Vector{T};
               scale::Bool=true,
-              atol=(length(ps)+length(qs))*eps(real(T)),  
-              rtol=eps(real(T)), #Base.rtoldefault(real(T)),
-              satol=atol,
-              srtol=eps(real(T)),
+              atol = eps(real(T)),  
+              rtol = eps(real(T)),
+              satol = atol,
+              srtol = rtol,
               verbose=false
               ) where {T <: AbstractFloat}
 

test/StandardBasis.jl

@@ -714,7 +714,7 @@ end
     r1, r2 = 1/2, 3/2
     U(n) = prod( (x-r1*alpha(j,n))^2 + r1^2*beta(j,n)^2 for j in 1:n)
     V(n) = prod( (x-r2*alpha(j,n))^2 + r2^2*beta(j,n)^2 for j in 1:n)
-W(n) = prod( (x-r1*alpha(j,n))^2 + r1^2*beta(j,n)^2 for j in (n+1):2n)
+    W(n) = prod( (x-r1*alpha(j,n))^2 + r1^2*beta(j,n)^2 for j in (n+1):2n)
     for n in 2:2:20
         p = U(n) * V(n); q = U(n) * W(n)
         @test degree(gcd(p,q, method=:numerical)) == degree(U(n))
@@ -729,7 +729,16 @@ W(n) = prod( (x-r1*alpha(j,n))^2 + r1^2*beta(j,n)^2 for j in (n+1):2n)
         @test degree(gcd(p,dp, method=:numerical)) == sum(max.(ms .- 1, 0))
     end
 
-
+    #
+    x =  variable(P{Float64})
+    for n in (10,20,50, 100)
+        p = (x-1)^n * (x-2)^n * (x-3)
+        q = (x-1) * (x-2) * (x-4)
+        @test degree(gcd(p,q, method=:numerical)) != 2
+        @test degree(Polynomials.ngcd′(p,q)) == 2
+    end
+    
+    
 
 end