fix(polys): fix multi factorisation over QQ<a>

sympy/polys/compatibility.py

@@ -101,6 +101,7 @@
 from sympy.polys.densetools import dup_mirror
 from sympy.polys.densetools import dup_scale
 from sympy.polys.densetools import dup_shift
+from sympy.polys.densetools import dmp_shift
 from sympy.polys.densetools import dup_transform
 from sympy.polys.densetools import dup_compose
 from sympy.polys.densetools import dmp_compose
@@ -515,6 +516,8 @@ def dup_scale(self, f, a):
         return self.from_dense(dup_scale(self.to_dense(f), a, self.domain))
     def dup_shift(self, f, a):
         return self.from_dense(dup_shift(self.to_dense(f), a, self.domain))
+    def dmp_shift(self, f, a):
+        return self.from_dense(dmp_shift(self.to_dense(f), a, self.ngens-1, self.domain))
     def dup_transform(self, f, p, q):
         return self.from_dense(dup_transform(self.to_dense(f), self.to_dense(p), self.to_dense(q), self.domain))
 

sympy/polys/densetools.py

@@ -894,6 +894,44 @@ def dup_shift(f, a, K):
     return f
 
 
+def dmp_shift(f, a, u, K):
+    """
+    Evaluate efficiently Taylor shift ``f(X + A)`` in ``K[X]``.
+
+    Examples
+    ========
+
+    >>> from sympy import symbols, ring, ZZ
+    >>> x, y = symbols('x y')
+    >>> R, _, _ = ring([x, y], ZZ)
+
+    >>> p = x**2*y + 2*x*y + 3*x + 4*y + 5
+
+    >>> R.dmp_shift(R(p), [ZZ(1), ZZ(2)])
+    x**2*y + 2*x**2 + 4*x*y + 11*x + 7*y + 22
+
+    >>> p.subs({x: x + 1, y: y + 2}).expand()
+    x**2*y + 2*x**2 + 4*x*y + 11*x + 7*y + 22
+    """
+    if not u:
+        return dup_shift(f, a[0], K)
+
+    if dmp_zero_p(f, u):
+        return f
+
+    a0, a1 = a[0], a[1:]
+
+    f = [ dmp_shift(c, a1, u-1, K) for c in f ]
+    n = len(f) - 1
+
+    for i in range(n, 0, -1):
+        for j in range(0, i):
+            afj = dmp_mul_ground(f[j], a0, u-1, K)
+            f[j + 1] = dmp_add(f[j + 1], afj, u-1, K)
+
+    return f
+
+
 def dup_transform(f, p, q, K):
     """
     Evaluate functional transformation ``q**n * f(p/q)`` in ``K[x]``.

sympy/polys/factortools.py

@@ -55,8 +55,7 @@
     dup_primitive, dmp_ground_primitive,
     dmp_eval_tail,
     dmp_eval_in, dmp_diff_eval_in,
-    dmp_compose,
-    dup_shift, dup_mirror)
+    dup_shift, dmp_shift, dup_mirror)
 
 from sympy.polys.euclidtools import (
     dmp_primitive,
@@ -125,6 +124,9 @@ def dmp_trial_division(f, factors, u, K):
             else:
                 break
 
+        if k == 0:
+            raise RuntimeError("trial division failed")
+
         result.append((factor, k))
 
     return _sort_factors(result)
@@ -1278,7 +1280,11 @@ def dup_ext_factor(f, K):
 
 
 def dmp_ext_factor(f, u, K):
-    """Factor multivariate polynomials over algebraic number fields. """
+    """Factor multivariate polynomials over algebraic number fields.
+
+    Algebraic Factoring and Rational Function Integration
+    Barry Trager 1976.
+    """
     if not u:
         return dup_ext_factor(f, K)
 
@@ -1296,12 +1302,11 @@ def dmp_ext_factor(f, u, K):
     if len(factors) == 1:
         factors = [f]
     else:
-        H = dmp_raise([K.one, s*K.unit], u, 0, K)
-
         for i, (factor, _) in enumerate(factors):
             h = dmp_convert(factor, u, K.dom, K)
             h, _, g = dmp_inner_gcd(h, g, u, K)
-            h = dmp_compose(h, H, u, K)
+            a = [s*K.unit]*(u+1)
+            h = dmp_shift(h, a, u, K)
             factors[i] = h
 
     return lc, dmp_trial_division(F, factors, u, K)

sympy/polys/polytools.py

@@ -6554,11 +6554,11 @@ def _try_factor(expr):
 
     try:
         return _generic_factor(f, gens, args, method='factor')
-    except PolynomialError as msg:
+    except PolynomialError:
         if not f.is_commutative:
             return factor_nc(f)
         else:
-            raise PolynomialError(msg)
+            raise
 
 
 @public

sympy/polys/sqfreetools.py

@@ -12,12 +12,12 @@
     dup_LC, dmp_ground_LC,
     dmp_zero_p,
     dmp_ground,
-    dup_degree, dmp_degree,
+    dup_degree, dmp_degree, dmp_degree_in,
     dmp_raise, dmp_inject,
     dup_convert)
 from sympy.polys.densetools import (
     dup_diff, dmp_diff, dmp_diff_in,
-    dup_shift, dmp_compose,
+    dup_shift, dmp_shift,
     dup_monic, dmp_ground_monic,
     dup_primitive, dmp_ground_primitive)
 from sympy.polys.euclidtools import (
@@ -70,8 +70,20 @@ def dmp_sqf_p(f, u, K):
     """
     if dmp_zero_p(f, u):
         return True
-    else:
-        return not dmp_degree(dmp_gcd(f, dmp_diff(f, 1, u, K), u, K), u)
+
+    for i in range(u+1):
+
+        fp = dmp_diff_in(f, 1, i, u, K)
+
+        if dmp_zero_p(fp, u):
+            continue
+
+        gcd = dmp_gcd(f, fp, u, K)
+
+        if dmp_degree_in(gcd, i, u) != 0:
+            return False
+
+    return True
 
 
 def dup_sqf_norm(f, K):
@@ -152,7 +164,7 @@ def dmp_sqf_norm(f, u, K):
         raise DomainError("ground domain must be algebraic")
 
     g = dmp_raise(K.mod.to_list(), u + 1, 0, K.dom)
-    F = dmp_raise([K.one, -K.unit], u, 0, K)
+    a = [-K.unit] * (u + 1)
 
     s = 0
 
@@ -163,7 +175,7 @@ def dmp_sqf_norm(f, u, K):
         if dmp_sqf_p(r, u, K.dom):
             break
         else:
-            f, s = dmp_compose(f, F, u, K), s + 1
+            f, s = dmp_shift(f, a, u, K), s + 1
 
     return s, f, r
 
@@ -381,6 +393,8 @@ def dmp_sqf_list(f, u, K, all=False):
     >>> R.dmp_sqf_list(f, all=True)
     (1, [(1, 1), (x + y, 2), (x, 3)])
 
+    Yun, David Y.Y. (1976). "On square-free decomposition algorithms".
+    doi:10.1145/800205.806320. ISBN 978-1-4503-7790-4.
     """
     if not u:
         return dup_sqf_list(f, K, all=all)