Repository: diofant/diofant

Issue:
dmp_sqf_part(), dmp_sqf_p() and dmp_sqf_list() are wrong for multivariate polynomials
Simple example:
```
In [1]: R, x, y, z = ring('x y z', ZZ)                                                                                                  

In [2]: f = (x - y)*(z - 1)**2                                                                                                          

In [3]: f.is_squarefree   # same algorithm, as in dmp_sqf_p()                                                                                             
Out[3]: True

In [4]: R.dmp_sqf_part(f)                                                                                                               
Out[4]: x - y

In [5]: R.dmp_sqf_list(f)                                                                                                               
Out[5]: (1, [(x - y, 1)])
```

Apparently, functions are using algorithms, which really work only for univariate polynomials, e.g. [test for a square-free polynomial](https://github.com/diofant/diofant/blob/173ed6de7ab69a36ca209d6e9292d330c126b718/diofant/polys/sqfreetools.py#L27-L30) take formal derivative of poly ``f`` wrt the first argument.

Seems to be a long-standing issue for SymPy, e.g. it affects [this problem](https://github.com/sympy/sympy/issues/5786), though I can't find a dedicated bugreport.

@jksuom, you may want at least to track this problem for SymPy.

Discussion:

Patch:
diff --git a/diofant/domains/integerring.py b/diofant/domains/integerring.py
index f3fa96659..892e22d6f 100644
--- a/diofant/domains/integerring.py
+++ b/diofant/domains/integerring.py
@@ -1,5 +1,7 @@
 """Implementation of :class:`IntegerRing` class. """
 
+import abc
+
 from ..polys.polyerrors import CoercionFailed
 from .characteristiczero import CharacteristicZero
 from .groundtypes import (DiofantInteger, GMPYInteger, PythonInteger,
@@ -77,6 +79,11 @@ def _from_AlgebraicField(self, a, K0):
         if a.is_ground:
             return self.convert(a.LC(), K0.domain)
 
+    @abc.abstractmethod
+    def finite_field(self, p):
+        """Returns a finite field. """
+        raise NotImplementedError  # pragma: no cover
+
 
 class PythonIntegerRing(IntegerRing):
     """Integer ring based on Python's integers. """
@@ -105,6 +112,11 @@ def factorial(self, a):
         """Compute factorial of ``a``. """
         return python_factorial(a)
 
+    def finite_field(self, p):
+        """Returns a finite field. """
+        from .finitefield import PythonFiniteField
+        return PythonFiniteField(p, self)
+
 
 class GMPYIntegerRing(IntegerRing):
     """Integer ring based on GMPY's integers. """
@@ -134,6 +146,11 @@ def factorial(self, a):
         """Compute factorial of ``a``. """
         return gmpy_factorial(a)
 
+    def finite_field(self, p):
+        """Returns a finite field. """
+        from .finitefield import GMPYFiniteField
+        return GMPYFiniteField(p, self)
+
 
 ZZ_python = PythonIntegerRing()
 ZZ_gmpy = GMPYIntegerRing()
diff --git a/diofant/polys/compatibility.py b/diofant/polys/compatibility.py
index 846af6bb0..9e6001f81 100644
--- a/diofant/polys/compatibility.py
+++ b/diofant/polys/compatibility.py
@@ -48,7 +48,7 @@
                             dup_isolate_real_roots_sqf, dup_refine_real_root,
                             dup_root_upper_bound, dup_sturm)
 from .sqfreetools import (dmp_sqf_list, dmp_sqf_list_include, dmp_sqf_norm,
-                          dmp_sqf_part)
+                          dmp_sqf_p, dmp_sqf_part)
 
 
 __all__ = ('IPolys',)
@@ -525,6 +525,9 @@ def dup_sturm(self, f):
         seq = dup_sturm(self.to_dense(f), self.domain)
         return list(map(self.from_dense, seq))
 
+    def dmp_sqf_p(self, f):
+        return dmp_sqf_p(self.to_dense(f), self.ngens-1, self.domain)
+
     def dmp_sqf_norm(self, f):
         s, F, R = dmp_sqf_norm(self.to_dense(f), self.ngens-1, self.domain)
         return s, self.from_dense(F), self.to_ground().from_dense(R)
diff --git a/diofant/polys/factortools.py b/diofant/polys/factortools.py
index e40d2c1b9..9faf191ac 100644
--- a/diofant/polys/factortools.py
+++ b/diofant/polys/factortools.py
@@ -12,15 +12,16 @@
 from .densebasic import (dmp_convert, dmp_degree_in, dmp_degree_list,
                          dmp_eject, dmp_exclude, dmp_from_dict, dmp_ground_LC,
                          dmp_ground_p, dmp_include, dmp_inject, dmp_LC,
-                         dmp_nest, dmp_one, dmp_raise, dmp_strip, dmp_TC,
-                         dmp_terms_gcd, dmp_zero_p, dup_inflate)
+                         dmp_nest, dmp_normal, dmp_one, dmp_raise, dmp_strip,
+                         dmp_swap, dmp_TC, dmp_terms_gcd, dmp_zero_p,
+                         dup_inflate)
 from .densetools import (dmp_clear_denoms, dmp_compose, dmp_diff_eval_in,
                          dmp_eval_in, dmp_eval_tail, dmp_ground_content,
                          dmp_ground_monic, dmp_ground_primitive,
                          dmp_ground_trunc, dup_mirror, dup_trunc)
-from .euclidtools import dmp_content, dmp_inner_gcd, dmp_primitive
+from .euclidtools import dmp_inner_gcd, dmp_primitive
 from .galoistools import (gf_add_mul, gf_div, gf_factor, gf_factor_sqf,
-                          gf_from_int_poly, gf_gcdex, gf_mul, gf_rem, gf_sqf_p,
+                          gf_from_int_poly, gf_gcdex, gf_mul, gf_rem,
                           gf_to_int_poly)
 from .polyconfig import query
 from .polyerrors import (CoercionFailed, DomainError, EvaluationFailed,
@@ -202,11 +203,12 @@ def dup_zz_zassenhaus(f, K):
             continue
 
         px = K.convert(px)
+        F = dmp_normal(f, 0, K.finite_field(px))
 
-        F = gf_from_int_poly(f, px)
-
-        if not gf_sqf_p(F, px, K):
+        if not dmp_sqf_p(F, 0, K.finite_field(px)):
             continue
+
+        F = gf_from_int_poly(f, px)
         fsqfx = gf_factor_sqf(F, px, K)[1]
         a.append((px, fsqfx))
         if len(fsqfx) < 15 or len(a) > 4:
@@ -1010,29 +1012,30 @@ def dmp_ext_factor(f, u, K):
     """Factor multivariate polynomials over algebraic number fields. """
     lc = dmp_ground_LC(f, u, K)
     f = dmp_ground_monic(f, u, K)
-    factors = []
 
-    if u:
-        for factor, k in dmp_ext_factor(dmp_content(f, u, K), u - 1, K)[1]:
-            factors.append(([factor], k))
+    if dmp_ground_p(f, None, u):
+        return lc, []
 
-    if dmp_degree_in(f, 0, u) > 0:
-        sqf = dmp_sqf_part(f, u, K)
-        s, g, r = dmp_sqf_norm(sqf, u, K)
+    f, F = dmp_sqf_part(f, u, K), f
+    s, g, r = dmp_sqf_norm(f, u, K)
 
-        _, sqf_factors = dmp_factor_list(r, u, K.domain)
+    _, factors = dmp_factor_list(r, u, K.domain)
 
+    if len(factors) == 1:
+        factors = [f]
+    else:
         H = dmp_raise([K.one, s*K.unit], u, 0, K)
 
-        for i, (factor, _) in enumerate(sqf_factors):
+        for i, (factor, _) in enumerate(factors):
             h = dmp_convert(factor, u, K.domain, K)
             h, _, g = dmp_inner_gcd(h, g, u, K)
-            h = dmp_compose(h, H, u, K)
-            sqf_factors[i] = h
-
-        factors.extend(dmp_trial_division(f, sqf_factors, u, K))
+            for j in range(u + 1):
+                h = dmp_swap(h, 0, j, u, K)
+                h = dmp_compose(h, H, u, K)
+                h = dmp_swap(h, 0, j, u, K)
+            factors[i] = h
 
-    return lc, factors
+    return lc, dmp_trial_division(F, factors, u, K)
 
 
 def dmp_gf_factor(f, u, K):
diff --git a/diofant/polys/galoistools.py b/diofant/polys/galoistools.py
index 61b6153fa..a107c0d96 100644
--- a/diofant/polys/galoistools.py
+++ b/diofant/polys/galoistools.py
@@ -1183,46 +1183,6 @@ def gf_irreducible_p(f, p, K):
     return _irred_methods[method](f, p, K)
 
 
-def gf_sqf_p(f, p, K):
-    """
-    Return ``True`` if ``f`` is square-free in ``GF(p)[x]``.
-
-    Examples
-    ========
-
-    >>> gf_sqf_p([3, 2, 4], 5, ZZ)
-    True
-    >>> gf_sqf_p([2, 4, 4, 2, 2, 1, 4], 5, ZZ)
-    False
-    """
-    _, f = gf_monic(f, p, K)
-
-    if not f:
-        return True
-    else:
-        return gf_gcd(f, gf_diff(f, p, K), p, K) == [K.one]
-
-
-def gf_sqf_part(f, p, K):
-    """
-    Return square-free part of a ``GF(p)[x]`` polynomial.
-
-    Examples
-    ========
-
-    >>> gf_sqf_part([1, 1, 3, 0, 1, 0, 2, 2, 1], 5, ZZ)
-    [1, 4, 3]
-    """
-    _, sqf = gf_sqf_list(f, p, K)
-
-    g = [K.one]
-
-    for f, _ in sqf:
-        g = gf_mul(g, f, p, K)
-
-    return g
-
-
 def gf_sqf_list(f, p, K):
     """
     Return the square-free decomposition of a ``GF(p)[x]`` polynomial.
diff --git a/diofant/polys/rings.py b/diofant/polys/rings.py
index 34769da81..c2942ebac 100644
--- a/diofant/polys/rings.py
+++ b/diofant/polys/rings.py
@@ -790,10 +790,7 @@ def is_quadratic(self):
 
     @property
     def is_squarefree(self):
-        if self.is_zero:
-            return True
-        else:
-            return not self.gcd(self.diff(0)).degree(0)
+        return self.ring.dmp_sqf_p(self)
 
     @property
     def is_irreducible(self):
diff --git a/diofant/polys/sqfreetools.py b/diofant/polys/sqfreetools.py
index e1c7e3bc9..619efbe9f 100644
--- a/diofant/polys/sqfreetools.py
+++ b/diofant/polys/sqfreetools.py
@@ -1,12 +1,13 @@
 """Square-free decomposition algorithms and related tools. """
 
-from .densearith import dmp_mul_ground, dmp_neg, dmp_quo, dmp_sub
-from .densebasic import (dmp_convert, dmp_degree_in, dmp_ground, dmp_ground_LC,
-                         dmp_inject, dmp_raise, dmp_zero_p)
-from .densetools import (dmp_compose, dmp_diff, dmp_ground_monic,
+from .densearith import dmp_mul, dmp_mul_ground, dmp_neg, dmp_quo, dmp_sub
+from .densebasic import (dmp_convert, dmp_ground, dmp_ground_LC, dmp_ground_p,
+                         dmp_inject, dmp_one_p, dmp_raise, dmp_swap,
+                         dmp_zero_p)
+from .densetools import (dmp_compose, dmp_diff_in, dmp_ground_monic,
                          dmp_ground_primitive)
-from .euclidtools import dmp_gcd, dmp_inner_gcd, dmp_resultant
-from .galoistools import gf_sqf_list, gf_sqf_part
+from .euclidtools import dmp_gcd, dmp_resultant
+from .galoistools import gf_sqf_list
 from .polyerrors import DomainError
 
 
@@ -17,17 +18,25 @@ def dmp_sqf_p(f, u, K):
     Examples
     ========
 
-    >>> dmp_sqf_p([[]], 1, ZZ)
+    >>> R, x, y = ring("x y", ZZ)
+
+    >>> R.dmp_sqf_p(R(0))
     True
-    >>> dmp_sqf_p([[1], [2, 0], [1, 0, 0]], 1, ZZ)
+    >>> R.dmp_sqf_p((x + y)**2)
     False
-    >>> dmp_sqf_p([[1], [], [1, 0, 0]], 1, ZZ)
+    >>> R.dmp_sqf_p(x**2 + y**2)
     True
     """
-    if dmp_zero_p(f, u):
+    if dmp_ground_p(f, None, u):
         return True
     else:
-        return not dmp_degree_in(dmp_gcd(f, dmp_diff(f, 1, u, K), u, K), 0, u)
+        g = f
+        for i in range(u + 1):
+            g = dmp_gcd(g, dmp_diff_in(f, 1, i, u, K), u, K)
+            if dmp_ground_p(g, None, u):
+                return True
+        else:
+            return False
 
 
 def dmp_sqf_norm(f, u, K):
@@ -42,17 +51,10 @@ def dmp_sqf_norm(f, u, K):
 
     >>> K = QQ.algebraic_field(I)
     >>> R, x, y = ring("x y", K)
-    >>> _, X, Y = ring("x y", QQ)
-
-    >>> s, f, r = R.dmp_sqf_norm(x*y + y**2)
-
-    >>> s == 1
-    True
-    >>> f == x*y + y**2 - I*y
-    True
-    >>> r == X**2*Y**2 + 2*X*Y**3 + Y**4 + Y**2
-    True
 
+    >>> R.dmp_sqf_norm(x*y + y**2)
+    (1, x*y - I*x + y**2 - 3*I*y - 2,
+     x**2*y**2 + x**2 + 2*x*y**3 + 2*x*y + y**4 + 5*y**2 + 4)
     """
     if not K.is_AlgebraicField:
         raise DomainError("ground domain must be algebraic")
@@ -67,21 +69,13 @@ def dmp_sqf_norm(f, u, K):
         r = dmp_resultant(g, h, u + 1, K.domain)
 
         if dmp_sqf_p(r, u, K.domain):
-            break
+            return s, f, r
         else:
-            f, s = dmp_compose(f, F, u, K), s + 1
-
-    return s, f, r
-
-
-def dmp_gf_sqf_part(f, u, K):
-    """Compute square-free part of ``f`` in ``GF(p)[X]``. """
-    if not u:
-        f = dmp_convert(f, u, K, K.domain)
-        g = gf_sqf_part(f, K.mod, K.domain)
-        return dmp_convert(g, u, K.domain, K)
-    else:  # pragma: no cover
-        raise NotImplementedError('multivariate polynomials over finite fields')
+            for j in range(u + 1):
+                f = dmp_swap(f, 0, j, u, K)
+                f = dmp_compose(f, F, u, K)
+                f = dmp_swap(f, 0, j, u, K)
+            s += 1
 
 
 def dmp_sqf_part(f, u, K):
@@ -98,7 +92,13 @@ def dmp_sqf_part(f, u, K):
 
     """
     if K.is_FiniteField:
-        return dmp_gf_sqf_part(f, u, K)
+        _, sqf = dmp_sqf_list(f, u, K)
+
+        g = [K.one]
+        for f, _ in sqf:
+            g = dmp_mul(g, f, u, K)
+
+        return g
 
     if dmp_zero_p(f, u):
         return f
@@ -106,7 +106,9 @@ def dmp_sqf_part(f, u, K):
     if K.is_negative(dmp_ground_LC(f, u, K)):
         f = dmp_neg(f, u, K)
 
-    gcd = dmp_gcd(f, dmp_diff(f, 1, u, K), u, K)
+    gcd = f
+    for i in range(u + 1):
+        gcd = dmp_gcd(gcd, dmp_diff_in(f, 1, i, u, K), u, K)
     sqf = dmp_quo(f, gcd, u, K)
 
     if K.is_Field:
@@ -131,6 +133,42 @@ def dmp_gf_sqf_list(f, u, K):
         raise NotImplementedError('multivariate polynomials over finite fields')
 
 
+def dmp_rr_yun0_sqf_list(f, u, K):
+    """Compute square-free decomposition of ``f`` in zero-characteristic ring ``K``.
+
+    References
+    ==========
+
+    * [LeeM13]_, page 8
+    """
+    if dmp_ground_p(f, None, u):
+        return []
+
+    result, count = [], 1
+    qs = [dmp_diff_in(f, 1, i, u, K) for i in range(u + 1)]
+
+    g = f
+    for q in qs:
+        g = dmp_gcd(g, q, u, K)
+
+    while not dmp_one_p(f, u, K):
+        for i in range(u + 1):
+            qs[i] = dmp_quo(qs[i], g, u, K)
+        f = dmp_quo(f, g, u, K)
+        for i in range(u + 1):
+            qs[i] = dmp_sub(qs[i], dmp_diff_in(f, 1, i, u, K), u, K)
+
+        g = f
+        for q in qs:
+            g = dmp_gcd(g, q, u, K)
+        if not dmp_one_p(g, u, K):
+            result.append((g, count))
+
+        count += 1
+
+    return result
+
+
 def dmp_sqf_list(f, u, K):
     """
     Return square-free decomposition of a polynomial in ``K[X]``.
@@ -140,9 +178,7 @@ def dmp_sqf_list(f, u, K):
 
     >>> R, x, y = ring("x y", ZZ)
 
-    >>> f = x**5 + 2*x**4*y + x**3*y**2
-
-    >>> R.dmp_sqf_list(f)
+    >>> R.dmp_sqf_list(x**5 + 2*x**4*y + x**3*y**2)
     (1, [(x + y, 2), (x, 3)])
     """
     if K.is_FiniteField:
@@ -158,30 +194,7 @@ def dmp_sqf_list(f, u, K):
             f = dmp_neg(f, u, K)
             coeff = -coeff
 
-    if dmp_degree_in(f, 0, u) <= 0:
-        return coeff, []
-
-    result, i = [], 1
-
-    h = dmp_diff(f, 1, u, K)
-    g, p, q = dmp_inner_gcd(f, h, u, K)
-
-    while True:
-        d = dmp_diff(p, 1, u, K)
-        h = dmp_sub(q, d, u, K)
-
-        if dmp_zero_p(h, u):
-            result.append((p, i))
-            break
-
-        g, p, q = dmp_inner_gcd(p, h, u, K)
-
-        if dmp_degree_in(g, 0, u) > 0:
-            result.append((g, i))
-
-        i += 1
-
-    return coeff, result
+    return coeff, dmp_rr_yun0_sqf_list(f, u, K)
 
 
 def dmp_sqf_list_include(f, u, K):
@@ -193,9 +206,7 @@ def dmp_sqf_list_include(f, u, K):
 
     >>> R, x, y = ring("x y", ZZ)
 
-    >>> f = x**5 + 2*x**4*y + x**3*y**2
-
-    >>> R.dmp_sqf_list_include(f)
+    >>> R.dmp_sqf_list_include(x**5 + 2*x**4*y + x**3*y**2)
     [(1, 1), (x + y, 2), (x, 3)]
     """
     coeff, factors = dmp_sqf_list(f, u, K)
diff --git a/docs/internals/polys.rst b/docs/internals/polys.rst
index eceab1cdd..29a08a908 100644
--- a/docs/internals/polys.rst
+++ b/docs/internals/polys.rst
@@ -84,6 +84,9 @@ may be slightly more efficient.)
 .. automodule:: diofant.polys.rootisolation
     :members:
 
+.. automodule:: diofant.polys.sqfreetools
+    :members:
+
 Manipulation of dense, univariate polynomials with finite field coefficients
 ****************************************************************************
 
diff --git a/docs/literature.rst b/docs/literature.rst
index dee9550c3..ee5f7d326 100644
--- a/docs/literature.rst
+++ b/docs/literature.rst
@@ -185,3 +185,6 @@ a theoretical foundation for implementing the Diofant.
     CA Systems, 1999, `<http://www.math.unm.edu/~wester/cas/book/Wester.pdf>`_
 
 .. [Gruntz96] Gruntz Thesis, http://www.cybertester.com/data/gruntz.pdf
+
+.. [LeeM13] M M-D Lee, Factorization of multivariate polynomials,
+    Ph.D. thesis, University of Kaiserslautern, 2013.
diff --git a/docs/release/notes-0.10.rst b/docs/release/notes-0.10.rst
index b4408b196..dbc8e7073 100644
--- a/docs/release/notes-0.10.rst
+++ b/docs/release/notes-0.10.rst
@@ -80,6 +80,7 @@ Minor changes
 * Support gmpy2.mpz ground type for numerator/denominator of :class:`~diofant.core.numbers.Rational`, see :pull:`694`.
 * Added ``FALLBACK_GCD_ZZ_METHOD`` configuration option to specify GCD algorithm for polynomials with integer coefficients if heuristic GCD was off or just unlucky, see :pull:`721`.
 * Added ``GCD_AA_METHOD`` configuration option to specify GCD algorithm for polynomials with algebraic coefficients, see :pull:`721`.
+* :meth:`~diofant.polys.polytools.Poly.sqf_part`, :meth:`~diofant.polys.polytools.Poly.sqf_norm`, :meth:`~diofant.polys.polytools.Poly.sqf_list` methods and  :attr:`~diofant.polys.polytools.Poly.is_squarefree` property use notion of being square-free w.r.t. to all polynomial variables, see :pull:`726`.
 
 Developer changes
 =================


Test patch:
diff --git a/diofant/domains/tests/test_gmpy.py b/diofant/domains/tests/test_gmpy.py
index 109d5699e..05d4ddddb 100644
--- a/diofant/domains/tests/test_gmpy.py
+++ b/diofant/domains/tests/test_gmpy.py
@@ -13,6 +13,9 @@
 
 @pytest.mark.skipif(gmpy is None, reason="no gmpy")
 def test_convert():
+    assert ZZ_gmpy.finite_field(2) == FF_gmpy(2)
+    assert ZZ_python.finite_field(2) == FF_python(2)
+
     F3 = FF(3)
     F3_gmpy = FF_gmpy(3)
     F3_python = FF_python(3)
diff --git a/diofant/polys/tests/test_galoistools.py b/diofant/polys/tests/test_galoistools.py
index e8f977830..9e8be4c40 100644
--- a/diofant/polys/tests/test_galoistools.py
+++ b/diofant/polys/tests/test_galoistools.py
@@ -19,8 +19,7 @@
                                        gf_irreducible_p, gf_lcm, gf_monic,
                                        gf_mul, gf_mul_ground, gf_neg, gf_pow,
                                        gf_pow_mod, gf_Qbasis, gf_Qmatrix,
-                                       gf_quo, gf_rem, gf_sqf_list, gf_sqf_p,
-                                       gf_sqf_part, gf_sqr, gf_sub,
+                                       gf_quo, gf_rem, gf_sqr, gf_sub,
                                        gf_sub_ground, gf_sub_mul, gf_to_dict,
                                        gf_to_int_poly, gf_trace_map, gf_value,
                                        linear_congruence)
@@ -408,34 +407,6 @@ def test_gf_irreducible_p():
     assert gf_irred_p_rabin(h, 17, ZZ) is False
 
 
-def test_gf_squarefree():
-    assert gf_sqf_list([], 11, ZZ) == (0, [])
-    assert gf_sqf_list([1], 11, ZZ) == (1, [])
-    assert gf_sqf_list([1, 1], 11, ZZ) == (1, [([1, 1], 1)])
-
-    assert gf_sqf_p([], 11, ZZ) is True
-    assert gf_sqf_p([1], 11, ZZ) is True
-    assert gf_sqf_p([1, 1], 11, ZZ) is True
-
-    f = gf_from_dict({11: 1, 0: 1}, 11, ZZ)
-
-    assert gf_sqf_p(f, 11, ZZ) is False
-
-    assert gf_sqf_list(f, 11, ZZ) == (1, [([1, 1], 11)])
-
-    f = [1, 5, 8, 4]
-
-    assert gf_sqf_p(f, 11, ZZ) is False
-
-    assert gf_sqf_list(f, 11, ZZ) == (1, [([1, 1], 1), ([1, 2], 2)])
-
-    assert gf_sqf_part(f, 11, ZZ) == [1, 3, 2]
-
-    f = [1, 0, 0, 2, 0, 0, 2, 0, 0, 1, 0]
-
-    assert gf_sqf_list(f, 3, ZZ) == (1, [([1, 0], 1), ([1, 1], 3), ([1, 2], 6)])
-
-
 def test_gf_frobenius_map():
     f = [2, 0, 1, 0, 2, 2, 0, 2, 2, 2]
     g = [1, 1, 0, 2, 0, 1, 0, 2, 0, 1]
@@ -698,8 +669,6 @@ def test_gf_factor():
     p = ZZ(nextprime(int((2**15 * pi).evalf())))
     f = gf_from_dict({15: 1, 1: 1, 0: 1}, p, ZZ)
 
-    assert gf_sqf_p(f, p, ZZ) is True
-
     g = (1, [([1, 22730, 68144], 1),
              ([1, 81553, 77449, 86810, 4724], 1),
              ([1, 86276, 56779, 14859, 31575], 1),
@@ -728,8 +697,6 @@ def test_gf_factor():
     p = ZZ(nextprime(int((2**4 * pi).evalf())))
     f = [1, 2, 5, 26, 41, 39, 38]
 
-    assert gf_sqf_p(f, p, ZZ) is True
-
     g = (1, [([1, 44, 26], 1),
              ([1, 11, 25, 18, 30], 1)])
 
diff --git a/diofant/polys/tests/test_sqfreetools.py b/diofant/polys/tests/test_sqfreetools.py
index 5a5ea0df2..6a27cdd10 100644
--- a/diofant/polys/tests/test_sqfreetools.py
+++ b/diofant/polys/tests/test_sqfreetools.py
@@ -15,58 +15,56 @@
 f_0, f_1, f_2, f_3, f_4, f_5, f_6 = f_polys()
 
 
-def test_dup_sqf():
+def test_dmp_sqf():
     R, x = ring("x", ZZ)
 
-    assert R.dmp_sqf_part(0) == 0
+    assert R(0).sqf_part() == 0
     assert R(0).is_squarefree is True
 
-    assert R.dmp_sqf_part(7) == 1
+    assert R(7).sqf_part() == 1
     assert R(7).is_squarefree is True
 
-    assert R.dmp_sqf_part(2*x + 2) == x + 1
+    assert (2*x + 2).sqf_part() == x + 1
     assert (2*x + 2).is_squarefree is True
 
-    assert R.dmp_sqf_part(x**3 + x + 1) == x**3 + x + 1
+    assert (x**3 + x + 1).sqf_part() == x**3 + x + 1
     assert (x**3 + x + 1).is_squarefree is True
 
-    assert R.dmp_sqf_part(-x**3 + x + 1) == x**3 - x - 1
+    assert (-x**3 + x + 1).sqf_part() == x**3 - x - 1
     assert (-x**3 + x + 1).is_squarefree is True
 
-    assert R.dmp_sqf_part(2*x**3 + 3*x**2) == 2*x**2 + 3*x
+    assert (2*x**3 + 3*x**2).sqf_part() == 2*x**2 + 3*x
     assert (2*x**3 + 3*x**2).is_squarefree is False
 
-    assert R.dmp_sqf_part(-2*x**3 + 3*x**2) == 2*x**2 - 3*x
+    assert (-2*x**3 + 3*x**2).sqf_part() == 2*x**2 - 3*x
     assert (-2*x**3 + 3*x**2).is_squarefree is False
 
-    assert R.dmp_sqf_part(x**3 - 3*x - 2) == x**2 - x - 2
+    assert (x**3 - 3*x - 2).sqf_part() == x**2 - x - 2
     assert (x**3 - 3*x - 2).is_squarefree is False
 
-    assert R.dmp_sqf_list(0) == (0, [])
-    assert R.dmp_sqf_list(1) == (1, [])
+    assert R(0).sqf_list() == (0, [])
+    assert R(1).sqf_list() == (1, [])
 
-    assert R.dmp_sqf_list(x) == (1, [(x, 1)])
-    assert R.dmp_sqf_list(2*x**2) == (2, [(x, 2)])
-    assert R.dmp_sqf_list(3*x**3) == (3, [(x, 3)])
+    assert x.sqf_list() == (1, [(x, 1)])
+    assert (2*x**2).sqf_list() == (2, [(x, 2)])
+    assert (3*x**3).sqf_list() == (3, [(x, 3)])
 
-    assert R.dmp_sqf_list(-x**5 + x**4 + x - 1) == \
-        (-1, [(x**3 + x**2 + x + 1, 1), (x - 1, 2)])
-    assert R.dmp_sqf_list(x**8 + 6*x**6 + 12*x**4 + 8*x**2) == \
-        ( 1, [(x, 2), (x**2 + 2, 3)])
+    assert (-x**5 + x**4 + x - 1).sqf_list() == (-1, [(x**3 + x**2 + x + 1, 1),
+                                                      (x - 1, 2)])
+    assert (x**8 + 6*x**6 + 12*x**4 + 8*x**2).sqf_list() == (1, [(x, 2),
+                                                                 (x**2 + 2, 3)])
 
-    assert R.dmp_sqf_list(2*x**2 + 4*x + 2) == (2, [(x + 1, 2)])
+    assert (2*x**2 + 4*x + 2).sqf_list() == (2, [(x + 1, 2)])
 
     R, x = ring("x", QQ)
-    assert R.dmp_sqf_list(2*x**2 + 4*x + 2) == (2, [(x + 1, 2)])
+    assert (2*x**2 + 4*x + 2).sqf_list() == (2, [(x + 1, 2)])
 
     R, x = ring("x", FF(2))
-    assert R.dmp_sqf_list(x**2 + 1) == (1, [(x + 1, 2)])
+    assert (x**2 + 1).sqf_list() == (1, [(x + 1, 2)])
 
     R, x = ring("x", FF(3))
-    assert R.dmp_sqf_list(x**10 + 2*x**7 + 2*x**4 + x) == \
-        (1, [(x, 1),
-             (x + 1, 3),
-             (x + 2, 6)])
+    assert (x**10 + 2*x**7 + 2*x**4 + x).sqf_list() == (1, [(x, 1), (x + 1, 3),
+                                                            (x + 2, 6)])
 
     R1, x = ring("x", ZZ)
     R2, y = ring("y", FF(3))
@@ -74,12 +72,70 @@ def test_dup_sqf():
     f = x**3 + 1
     g = y**3 + 1
 
-    assert R1.dmp_sqf_part(f) == f
-    assert R2.dmp_sqf_part(g) == y + 1
+    assert f.sqf_part() == f
+    assert g.sqf_part() == y + 1
 
     assert f.is_squarefree is True
     assert g.is_squarefree is False
 
+    R, x = ring("x", FF(5))
+
+    f = x**8 + x**7 + 3*x**6 + x**4 + 2*x**2 + 2*x + 1
+
+    assert f.sqf_part() == x**2 + 4*x + 3
+
+    f = 3*x**2 + 2*x + 4
+
+    assert f.is_squarefree is True
+
+    f = 2*x**6 + 4*x**5 + 4*x**4 + 2*x**3 + 2*x**2 + x + 4
+
+    assert f.is_squarefree is False
+
+    R, x = ring("x", FF(11))
+
+    f = x**3 + 5*x**2 + 8*x + 4
+
+    assert f.sqf_part() == x**2 + 3*x + 2
+
+    assert R(0).is_squarefree is True
+    assert R(0).sqf_list() == (0, [])
+
+    assert R(1).is_squarefree is True
+    assert R(1).sqf_list() == (1, [])
+
+    assert (x + 1).is_squarefree is True
+    assert (x + 1).sqf_list() == (1, [(x + 1, 1)])
+
+    f = x**11 + 1
+
+    assert f.is_squarefree is False
+    assert f.sqf_list() == (1, [(x + 1, 11)])
+
+    f = x**3 + 5*x**2 + 8*x + 4
+
+    assert f.is_squarefree is False
+    assert f.sqf_list() == (1, [(x + 1, 1), (x + 2, 2)])
+
+    R, x = ring('x', FF(3))
+
+    f = x**10 + 2*x**7 + 2*x**4 + x
+
+    assert f.is_squarefree is False
+    assert f.sqf_list() == (1, [(x, 1), (x + 1, 3), (x + 2, 6)])
+
+    R, x = ring('x', FF(53))
+
+    f = x**6 + 2*x**5 + 5*x**4 + 26*x**3 + 41*x**2 + 39*x + 38
+
+    assert f.is_squarefree is True
+
+    R, x = ring('x', FF(102953))
+
+    f = x**15 + x + 1
+
+    assert f.is_squarefree is True
+
     R, x, y = ring("x,y", ZZ)
 
     A = x**4 - 3*x**2 + 6
@@ -89,9 +145,9 @@ def test_dup_sqf():
     res = R.dmp_resultant(f, g)
     h = (4*y**2 + 1).drop(x)
 
-    assert R.drop(x).dmp_sqf_list(res) == (45796, [(h, 3)])
+    assert res.sqf_list() == (45796, [(h, 3)])
 
-    pytest.raises(DomainError, lambda: R.dmp_sqf_norm(x**2 - 1))
+    pytest.raises(DomainError, lambda: (x**2 - 1).sqf_norm())
 
     Rt, t = ring("t", ZZ)
     R, x = ring("x", Rt)
@@ -100,18 +156,17 @@ def test_dup_sqf():
     K = QQ.algebraic_field(sqrt(3))
     R, x = ring("x", K)
     _, X = ring("x", QQ)
-    assert R.dmp_sqf_norm(x**2 - 2) == (1, x**2 + K([QQ(-2), QQ(0)])*x + 1, X**4 - 10*X**2 + 1)
-
+    assert (x**2 - 2).sqf_norm() == (1, x**2 + K([QQ(-2), QQ(0)])*x + 1,
+                                     X**4 - 10*X**2 + 1)
 
-def test_dmp_sqf():
     R, x, y = ring("x,y", ZZ)
-    assert R.dmp_sqf_part(0) == 0
+    assert R(0).sqf_part() == 0
     assert R(0).is_squarefree is True
 
-    assert R.dmp_sqf_part(7) == 1
+    assert R(7).sqf_part() == 1
     assert R(7).is_squarefree is True
 
-    assert R.dmp_sqf_list(3) == (3, [])
+    assert R(3).sqf_list() == (3, [])
     assert R.dmp_sqf_list_include(3) == [(3, 1)]
 
     R, x, y, z = ring("x,y,z", ZZ)
@@ -127,41 +182,55 @@ def test_dmp_sqf():
     assert (f_5**2).is_squarefree is False
 
     assert f_4.is_squarefree is True
-    assert R.dmp_sqf_part(f_4) == -f_4
-
-    assert R.dmp_sqf_part(f_5) == x + y - z
+    assert f_4.sqf_part() == -f_4
+    assert f_5.sqf_part() == x + y - z
 
     R, x, y, z, t = ring("x,y,z,t", ZZ)
     assert f_6.is_squarefree is True
-    assert R.dmp_sqf_part(f_6) == f_6
+    assert f_6.sqf_part() == f_6
 
     R, x = ring("x", ZZ)
     f = -x**5 + x**4 + x - 1
 
-    assert R.dmp_sqf_list(f) == (-1, [(x**3 + x**2 + x + 1, 1), (x - 1, 2)])
+    assert f.sqf_list() == (-1, [(x**3 + x**2 + x + 1, 1), (x - 1, 2)])
     assert R.dmp_sqf_list_include(f) == [(-x**3 - x**2 - x - 1, 1), (x - 1, 2)]
 
     f = 2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16
 
-    assert R.dmp_sqf_list(f) == (2, [(x + 1, 2), (x + 2, 3)])
+    assert f.sqf_list() == (2, [(x + 1, 2), (x + 2, 3)])
     assert R.dmp_sqf_list_include(f) == [(2, 1), (x + 1, 2), (x + 2, 3)]
 
     R, x, y = ring("x,y", ZZ)
     f = -x**5 + x**4 + x - 1
 
-    assert R.dmp_sqf_list(f) == (-1, [(x**3 + x**2 + x + 1, 1), (x - 1, 2)])
+    assert f.sqf_list() == (-1, [(x**3 + x**2 + x + 1, 1), (x - 1, 2)])
     assert R.dmp_sqf_list_include(f) == [(-x**3 - x**2 - x - 1, 1), (x - 1, 2)]
 
-    pytest.raises(DomainError, lambda: R.dmp_sqf_norm(x**2 + y**2))
+    pytest.raises(DomainError, lambda: (x**2 + y**2).sqf_norm())
 
     f = -x**2 + 2*x - 1
     assert R.dmp_sqf_list_include(f) == [(-1, 1), (x - 1, 2)]
 
     R, x, y = ring("x,y", FF(2))
-    pytest.raises(NotImplementedError, lambda: R.dmp_sqf_list(y**2 + 1))
-    pytest.raises(NotImplementedError,
-                  lambda: R.dmp_sqf_part(x**3 + 2*x**2*y + x*y**2))
+    pytest.raises(NotImplementedError, lambda: (y**2 + 1).sqf_list())
+    pytest.raises(NotImplementedError, lambda: (x**3 + 2*x**2*y + x*y**2).sqf_part())
 
     R, x, y = ring("x,y", QQ.algebraic_field(I))
-    assert R.dmp_sqf_list(x**2 + 2*I*x - 1) == (R.one.to_dense()[0][0],
-                                                [(x + I, 2)])
+    assert (x**2 + 2*I*x - 1).sqf_list() == (1, [(x + I, 2)])
+
+
+def test_diofantissue_714():
+    R, x, y, z = ring('x y z', ZZ)
+
+    f = (x - y)*(z - 1)**2
+
+    assert f.is_squarefree is False
+    assert f.sqf_part() == (x - y)*(z - 1)
+    assert f.sqf_list() == (1, [(x - y, 1), (z - 1, 2)])
+
+    g = f
+    f = (x - y)*(z - 1)
+    assert (f*g).sqf_list() == (1, [(x - y, 2), (z - 1, 3)])
+
+    assert f.is_squarefree
+    assert ((x - y)*f).is_squarefree is False