Merge pull request #241 from jverzani/issue_240

add tolerance to gcd; adjust chop; close issue #240

Author: jverzani

Project.toml

@@ -3,7 +3,7 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "1.1.2"
+version = "1.1.3"
 
 [deps]
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"

src/common.jl

@@ -209,9 +209,10 @@ function chop!(p::AbstractPolynomial{T};
     rtol::Real = Base.rtoldefault(real(T)),
                atol::Real = 0,) where {T}
     isempty(coeffs(p)) && return p
+    tol = norm(coeffs(p)) * rtol + atol
     for i = lastindex(p):-1:0
         val = p[i]
-        if !isapprox(val, zero(T); rtol = rtol, atol = atol)
+        if abs(val) > tol #!isapprox(val, zero(T); rtol = rtol, atol = atol)
             resize!(p.coeffs, i + 1); 
             return p
         end
@@ -520,7 +521,7 @@ function Base.divrem(num::P, den::O) where {P <: AbstractPolynomial,O <: Abstrac
 end
 
 """
-    gcd(a::AbstractPolynomial, b::AbstractPolynomial)
+    gcd(a::AbstractPolynomial, b::AbstractPolynomial; atol::Real=0, rtol::Real=Base.rtoldefault)
 
 Find the greatest common denominator of two polynomials recursively using
 [Euclid's algorithm](http://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Euclid.27s_algorithm).
@@ -535,15 +536,18 @@ Polynomial(4.0 - 6.0*x + 2.0*x^2)
 
 ```
 """
-function Base.gcd(p1::AbstractPolynomial{T}, p2::AbstractPolynomial{S}) where {T,S}
+function Base.gcd(p1::AbstractPolynomial{T}, p2::AbstractPolynomial{S};
+                  atol::Real=zero(real(promote_type(T,S))),
+                  rtol::Real=Base.rtoldefault(real(promote_type(T,S)))
+                  ) where {T,S}
     r₀, r₁ = promote(p1, p2)
     iter = 1
     itermax = length(r₁)
 
     while !iszero(r₁) && iter ≤ itermax
         _, rtemp = divrem(r₀, r₁)
         r₀ = r₁
-        r₁ = truncate(rtemp)  
+        r₁ = truncate(rtemp; atol=atol, rtol=rtol)  
         iter += 1
     end
     return r₀

test/StandardBasis.jl

@@ -673,6 +673,13 @@ end
         @test degree(gcd(p1, P(eps(0.)))) == 0          # ditto
         @test degree(gcd(p1, P(0))) == degree(p1) # P(0) has the roots of p1
         @test degree(gcd(p1 + p2 * 170.10734737144486, p2)) == 0          # see, c.f., #122
+
+        ## Issue #240; add tolerance to gcd
+        a = P([0.8457170323029561, 0.47175077674705257,  0.9775441940117577]);
+        b = P([0.5410010714904849, 0.533604905984294]);
+        d = P([0.5490673726445683, 0.15991109487875477]);
+        @test Polynomials.isconstant(gcd(a*d,b*d))
+        @test !Polynomials.isconstant(gcd(a*d, b*d, atol=sqrt(eps())))
 
         p1 = fromroots(P, [1.,2.,3.])
         p2 = fromroots(P, [1.,2.,6.])