diff --git a/numth/types/polynomial.py b/numth/types/polynomial.py
index 6c96f3a..c625ea5 100644
--- a/numth/types/polynomial.py
+++ b/numth/types/polynomial.py
@@ -255,6 +255,38 @@ def derivative(self, order=None):
 
         return deriv
 
+    #-------------------------
+
+    def mod_gcd(self, other, modulus):
+        if self == other == Polynomial({0: 0}):
+            raise ValueError('gcd(0, 0) is undefined')
+
+        if self.degree == other.degree == 1:
+            s_m = self.coeffs[0] % modulus
+            o_m = other.coeffs[0] % modulus
+            return gcd(s_m, o_m) % modulus
+
+        a, b = self % modulus, other % modulus
+        while b.degree >= 0:
+            a, b = b, (a % b) % modulus
+
+        return (mod_inverse(a.leading_coeff, modulus) * a) % modulus
+    
+    #-------------------------
+    
+    def gcd(self, other):
+        if self == other == polyn('0'):
+            raise ValueError('gcd(0, 0) is undefined')
+        
+        if self.degree == other.degree == 1:
+            return gcd(self.coeffs[0], other.coeffs[0])
+
+        a, b = self.to_integer_polyn(), other.to_integer_polyn()
+        while b.degree >= 0:
+            a, b = b, a % b
+
+        return a.canonical() / gcd(*a.coeffs.values())
+
 #===========================================================
 
 def _term_pattern():
diff --git a/tests/test_types_polynomial.py b/tests/test_types_polynomial.py
index 84d138f..a31e61f 100644
--- a/tests/test_types_polynomial.py
+++ b/tests/test_types_polynomial.py
@@ -388,6 +388,19 @@ def test_mod_eval(e1, c1, e2, c2, val, mod):
     p = make(e1, c1, e2, c2)
     assume( p.mod_eval(val, mod) == p.eval(val) % mod )
 
+#-----------------------------
+
+    @given(*coords(4, coeff_min=1))
+    def test_mod_gcd(e1, c1, e2, c2, e3, c3, e4, c4):
+        modulus = choice([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47])
+        assume( are_distinct(e1, e2) and are_distinct(e3, e4) )
+        p1 = make(e1, c1, e2, c2) % modulus
+        p2 = make(e3, c3, e4, c4) % modulus
+        d = p1.mod_gcd(p2, modulus)
+        assert( d == p2.mod_gcd(p1, modulus) )
+        assert( p1 % d % modulus == p2 % d % modulus == polyn('0') )
+        assert( (p1 // d).mod_gcd(p2 // d, modulus) == polyn('1') )
+
 #=============================
 
 @given(