diff --git a/benchmarks/polys.py b/benchmarks/polys.py
index 60318c2..c2953a2 100644
--- a/benchmarks/polys.py
+++ b/benchmarks/polys.py
@@ -1,4 +1,4 @@
-from sympy import symbols, prod, prem, rem, degree, LC
+from sympy import symbols, prod, prem, rem, degree, LC, pquo, quo
 from sympy.polys import QQ, Poly
 
 
@@ -50,7 +50,6 @@ def as_poly(self):
     def as_ring(self):
         return (self.to_ring(self.f), self.to_ring(self.g), self.to_ring(self.d), self.ring)
 
-
 class _LinearDenseQuadraticGCD(_GCDExample):
     """A pair of linearly dense quartic inputs with quadratic GCDs"""
 
@@ -71,7 +70,6 @@ def make_poly(self, n):
         g = d * (2 + x ** (n + 1) + sum([y[i] ** (n + 1) for i in range(n)]))
         return f, g, d, syms
 
-
 class _QuadraticNonMonicGCD(_GCDExample):
     """A pair of quadratic polynomials with a non-monic GCD."""
 
@@ -82,7 +80,6 @@ def make_poly(self, n):
         g = d * (2 + x * y[0] + sum(y[1:n])) ** 2
         return f, g, d, syms
 
-
 class _SparseNonMonicQuadratic(_GCDExample):
     """A pair of sparse non-monic quadratic polynomials with linear GCDs."""
 
@@ -127,8 +124,8 @@ def time_op(self, n, impl):
     def teardown(self, n, impl):
         assert self.expected_result == self.returned_result
 
-
 class _TimePREM(_TimeOP):
+    """Benchmarks for pseudo-remainder method"""
 
     def expected(self, f, g, d, syms):
         x = syms[0]
@@ -162,5 +159,44 @@ class TimePREM_QuadraticNonMonicGCD(_TimePREM):
 
 
 class TimePREM_SparseNonMonicQuadratic(_TimePREM):
+    GCDExampleCLS = _SparseNonMonicQuadratic
+    params = [(1, 3, 5, 8), ('expr', 'dense', 'sparse')]
+
+
+class _TimePQUO(_TimeOP):
+    """Benchmarks for pseudo-quotient method"""
+
+    def expected(self, f, g, d, syms):
+        x = syms[0]
+        pquo_f_g_x = quo(f * LC(g, x) ** (degree(f, x) - degree(g, x) + 1), g, x)
+        return pquo_f_g_x
+
+    def get_func_expr(self, f, g, d, syms):
+        x = syms[0]
+        return lambda: pquo(f, g, x)
+
+    def get_func_poly(self, f, g, d):
+        return lambda: f.pquo(g)
+
+    def get_func_sparse(self, f, g, d, ring):
+        return lambda: f.pquo(g)
+
+
+class TimePQUO_LinearDenseQuadraticGCD(_TimePQUO):
+    GCDExampleCLS = _LinearDenseQuadraticGCD
+    params = [(1, 3, 5), ('expr', 'dense', 'sparse')] # This case is slow for n=8.
+
+
+class TimePQUO_SparseGCDHighDegree(_TimePQUO):
+    GCDExampleCLS = _SparseGCDHighDegree
+    params = [(1, 3, 5, 8), ('expr', 'dense', 'sparse')]
+
+
+class TimePQUO_QuadraticNonMonicGCD(_TimePQUO):
+    GCDExampleCLS = _QuadraticNonMonicGCD
+    params = [(1, 3, 5), ('expr', 'dense', 'sparse')] # This case is slow for n=8.
+
+
+class TimePQUO_SparseNonMonicQuadratic(_TimePQUO):
     GCDExampleCLS = _SparseNonMonicQuadratic
     params = [(1, 3, 5, 8), ('expr', 'dense', 'sparse')]
\ No newline at end of file