Fix Euclid's algorithm in gcd

Fixes #122.

Author: Arda Aytekin

src/Polynomials.jl

@@ -665,11 +665,19 @@ Poly(4.0 - 6.0⋅x + 2.0⋅x^2)
 ```
 """
 function gcd{T, S}(a::Poly{T}, b::Poly{S})
-  R = typeof(one(T)/one(S))
-  degree(b) == 0 && b ≈ zero(b) && return convert(Poly{R}, a)
+  U       = typeof(one(T)/one(S))
+  r₀, r₁  = convert(Poly{U}, a), truncate(convert(Poly{U}, b))
+  iter    = 1
+  itermax = degree(b) + 1
+
+  while r₁ ≉ zero(r₁) && iter ≤ itermax   # just to avoid unnecessary recursion
+    _, rtemp  = divrem(r₀, r₁)
+    r₀        = r₁
+    r₁        = truncate(rtemp)
+    iter      += 1
+  end
 
-  s, r = divrem(a, b)
-  return gcd(b, r)
+  return r₀
 end
 
 

test/runtests.jl

@@ -345,3 +345,17 @@ end
 @test length(collect(p1)) == degree(p1)+1
 
 @test [p1[idx] for idx in eachindex(p1)] == [1,2,0,3]
+
+p1 = Poly([2.,5.,1.])
+p2 = Poly([1.,2.,3.])
+
+@test degree(gcd(p1, p2))                       == 0          # no common roots
+@test degree(gcd(p1, Poly(5)))                  == 0          # ditto
+@test degree(gcd(p1, Poly(eps(0.))))            == 0          # ditto
+@test degree(gcd(p1, Poly(0)))                  == degree(p1) # Poly(0) has the roots of p1
+@test degree(gcd(p1+p2*170.10734737144486, p2)) == 0          # see, c.f., #122
+
+p1 = poly([1.,2.,3.])
+p2 = poly([1.,2.,6.])
+
+@test (res = roots(gcd(p1, p2)); 1. ∈ res && 2. ∈ res)