Repository: sympy/sympy
Pull request: fix(polys): fix multi factorisation over QQ<a>
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes gh-26497

#### Brief description of what is fixed or changed

There are various things here.

- Fix `dmp_sqf_p` to check properly that multivariate polynomials are square free over the ground domain. This is also fixes `Poly.is_sqf`.
- Change `dmp_sqf_norm` to use shifts in multiple variables.
- Add `_dmp_sqf_norm_shifts` to try to find efficient multivariate shifts.
- The `sqf_norm` functions now return a shift list rather than a single shift integer. This is to compute a shift like `f(x1-s1*a, x2-s2*a, ...)`.
- Add `dmp_shift` and `Poly.shift_list` for multivariate polynomial shifts.
- Add sanity checks that the degrees add up in various factorisation routines.
- Add `dmp_norm` to `PolyRing`.
- Add `norm` to `PolyElement`.
- Add tests for better general coverage of densetools.
- Add missing functions to polys docs and minor improvements to docstrings.
- Expand docstrings with explanation of some algorithms and references.

 The most important part is the fix for gh-26497 which is that `dmp_ext_factor` and hence `factor` over `ZZ_I` or `QQ_I` gave incorrect results for some multivariate polynomials. The bug is caused by `dmp_sqf_norm` choosing a polynomial that is not square-free which is in turn caused by `dmp_sqf_p` incorrectly reporting some polynomials as square-free.

The implementation of `dmp_sqf_p` on master treats a multivariate polynomial as if it were a univariate polynomial in its leading variable. This means that it only considers a polynomial to be square-free if it has squared factors of positive degree in the leading variable. Hence:
```python
In [1]: (y**2).as_poly(x, y).is_sqf
Out[1]: True
```
This is not the sense of "square-free" that is needed for `sqf_norm` i.e. a polynomial should only be square free if there are no squared factors having positive degree in *any* variable. Potentially there is some value in having two different notions of square-free but realistically the notion that we generally want is that of being square free in any variable so that is what I have changed it to.

Fixing `dmp_sqf_p` means that `dmp_sqf_norm` does not return an incorrect result with a norm that is not square-free. However then `dmp_sqf_norm` goes into an infinite loop. This is because it tries to shift the polynomial like `f(x,y,z,...) -> f(x-1,y,z,...) -> f(x-2,y,z,...)`. This is sequence of shifts is not guaranteed to result in a square-free norm for a multivariate polynomial. Instead we can use `f(x,y,z,...) -> f(x-1,y-1,z-1,...) -> f(x-2,y-2,z-2,...)`. Shifting all variables like this works and guarantees to find a square-free norm but can be very slow for some polynomials and caused some tests to hang.

A more efficient sequence of shifts requires that we use different shifts for each variable e.g. we try `f(x,y,z)`, `f(x-1,y,z)`, `f(x,y-1,z)` and so on before shifting all variables. This usually results in finding a shift such that the norm can be computed just as efficiently as before. The problem though is that this means that `s` returned from `dmp_sqf_norm` needs to be a list of integer shifts rather than a single integer.

Changing `dmp_sqf_norm` and also `Poly.sqf_norm` and `sqf_norm` to return a list of integer shifts is not backwards compatible. However I doubt that anyone is using these functions for anything. As far as I know their only real purpose is as an internal routine in Trager's algorithm for factorisation (`dup_ext_factor/dmp_ext_factor`). I don't see why anyone would be using these functions directly rather than e.g. `factor` so I think it should be okay to change the return type.

In turn I then added `dmp_shift` and `Poly.shift_list` so that the list of shifts can be used directly.

I then added various other things including tests, docs, etc as a general cleanup.

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* polys
    * A serious bug in factorisation of multivariate polynomials over extension fields was fixed. Previously, incorrect results were returned in some cases. The bug would most likely manifest when calling `factor` on an expression with more than one symbol and the complex unit `I`.
    * A bug in the Poly method `is_sqf` was fixed. Previously some multivariate polynomials were incorrectly reported as being square-free.
    * BREAKING CHANGE: The `sqf_norm` function and `Poly.sqf_norm` functions are now changed to return a list of integer shifts instead of a single integer shift. This is not a backwards compatible change but is needed because in general in the multivariate case different shifts are needed for each variable.
    * A new `shift_list` method is added to `Poly` and other polynomial types for computing shifts of multivariate polynomials.
    * Some documentation was added for internal routines in the polys module.
<!-- END RELEASE NOTES -->


Problem statement:
factor produces wrong output
`sympy.factor(expr)` produces the wrong output if an underscore is used in the name of one of the symbols in `expr`.
EDIT: The issue is more general. The correctness of the factorization output depends on the names of the variables. Using symbol names with more than one character also changes the output, for example.

**Observed behavior:**
For example,
```
a_n, m = sympy.symbols('a_n m', real=True, positive=True)
expr = (-a_n**4 *m**2 / (a_n**2 + 1) - 2 * a_n**4 *m / (a_n**2 + 1) - a_n**4 - a_n**4 / (a_n**2 + 1) + 2 * sympy.I * a_n**3 *m**2 / (a_n**2 + 1) + 4 * sympy.I * a_n**3 *m / (a_n**2 + 1) + 2 * sympy.I * a_n**3 + 2 * sympy.I * a_n**3 / (a_n**2 + 1) + 2 * sympy.I * a_n *m**2 / (a_n**2 + 1) + 4 * sympy.I * a_n *m / (a_n**2 + 1) + 2 * sympy.I * a_n + 2 * sympy.I * a_n / (a_n**2 + 1) +m**2 / (a_n**2 + 1) + 2 *m / (a_n**2 + 1) + 1 + 1 / (a_n**2 + 1))
```
gives different numerical results if I factor `expr` and replace `a_n` and `m` with numbers.
`expr.subs({a_n: 1, m: 1})` produces `12i`, whereas `expr.factor().subs({a_n: 1, m: 1})` produces `-6`.
A correct behavior would result in `12i`, or equivalent numerical expressions, in both cases.

On the other hand, defining
```
a_n, m = sympy.symbols('a m', real=True, positive=True)
```
correctly produces the equivalent outputs `12i` in both cases.

**Sympy version**
I observe this in the stable branch, but also in 1.12 and 1.10. Didn't check other versions.

**Other modifications that produce the same error**
Using `a_n, m = sympy.symbols('an m', real=True, positive=True)`, so a symbol with a name of more than one character, also reproduces the error.

Hints from issue discussion:
The bug seems to be related to factoring of multivariate expressions including `sympy.I`. A shorter reproduction example is

```py
import sympy
a, b = sympy.symbols("a b")
((b - sympy.I)**2 * (b + sympy.I) * (a + 1)).expand().factor()
```

which produces an incorrect expression $(a+1)(b^2 + 1)$. The originally reported role of variable naming seems to be related to the ordering of terms or factors: swapping $a$ and $b$ in the expression above makes the output correct.
This is a serious bug.

It has been around since SymPy started factorising such polynomials automatically: d3eab468768801dff2c810a5dc5d6d446d6db01c.

Before that you would need to use `extension=True` to factorise this. Now all mehods give the same:
```python
In [1]: import sympy
   ...: a, b = sympy.symbols("a b")
   ...: e = ((b - sympy.I)**2 * (b + sympy.I) * (a + 1))

In [6]: e.expand().factor()
Out[6]: 
        ⎛ 2    ⎞
(a + 1)⋅⎝b  + 1⎠

In [7]: e.expand().factor(extension=True)
Out[7]: 
        ⎛ 2    ⎞
(a + 1)⋅⎝b  + 1⎠

In [8]: e.expand().factor(domain=QQ_I)
Out[8]: 
        ⎛ 2    ⎞
(a + 1)⋅⎝b  + 1⎠

In [11]: e.expand().factor(domain=QQ.algebraic_field(I))
Out[11]: 
        ⎛ 2    ⎞
(a + 1)⋅⎝b  + 1⎠
```
There must be a bug in the code for factorisation over algebraic fields:
https://github.com/sympy/sympy/blob/d11d0b220e94dd04e0acaf339e784edb8e7d24fe/sympy/polys/factortools.py#L1280-L1307
The bug is in `dmp_sqf_list`:
```python
In [25]: sqf(expand(((y - I)**2 * (y + I) * (z + 1))))
Out[25]: 
       2        
(y - ⅈ) ⋅(y + ⅈ)

In [26]: sqf(expand(((y - I)**2 * (y + I) * (x + 1))))
Out[26]: x + 1
```
Both of those answers are wrong. The factors for one variable are missing in each case.

The bug is here:
https://github.com/sympy/sympy/blob/87bf9d4343202ee2d4d544a0ee8c6cdfe36feb6a/sympy/polys/sqfreetools.py#L409-L429

That seems to use the same logic as for the univariate case that you can compute the square free factorisation by differentiating and then computing the gcd. The difference though is that in the multivariate case there is more than one choice of variable to differentiate with respect to and this seems to make an arbitrary choice. We might then end up with a factor like `x + 1` and then differentiate wrt `y` and then exit the loop without collecting all factors.

What should probably happen is that once `p` gets down to degree 0 in the first variable it should recurse like it would at the start if the polynomial was already degree 0 in the first variable at the start:
https://github.com/sympy/sympy/blob/87bf9d4343202ee2d4d544a0ee8c6cdfe36feb6a/sympy/polys/sqfreetools.py#L404-L407

The implemented algorithm seems to be Yun's algorithm but applying the univariate version of the algorithm naively to multivariate polynomials. In Yun's paper it is stated that the algorithm requires the multivariate polynomials to be primitive in the given varaible that is used for differentiation:

> On square free decomposition algorithms, Yun, 1976
> https://dl.acm.org/doi/10.1145/800205.806320

I think this means that to apply this algorithm in the multivariate case you would first need to use `dmp_primitive` to factor out the content of the coefficients wrt x and then also compute the square free factorisation of the content recursively. Currently it uses `dmp_ground_primitive` which is not the same e.g.:
```python
In [12]: p = expand((y - I)**2 * (y + I) * (z + 1))

In [13]: p
Out[13]: 
 3      3      2        2                    
y ⋅z + y  - ⅈ⋅y ⋅z - ⅈ⋅y  + y⋅z + y - ⅈ⋅z - ⅈ

In [14]: p.as_poly().primitive()
Out[14]: (1, Poly(y**3*z + y**3 - I*y**2*z - I*y**2 + y*z + y - I*z - I, y, z, domain='ZZ_I'))

In [15]: p.as_poly(y).primitive()
Out[15]: (z + 1, Poly(y**3 - I*y**2 + y - I, y, domain='ZZ_I[z]'))
```
The bug for `sqf` does not require there to be any complex numbers i.e. it happens with ZZ as well as ZZ_I:
```python
In [1]: sqf(expand(((y - 2)**2 * (y + 2) * (x + 1))))
Out[1]: x + 1

In [2]: sqf(expand(((y - 2)**2 * (y + 2) * (z + 1))))
Out[2]: 
       2        
(y - 2) ⋅(y + 2)
```
I think that the difference in the case of `factor` is just that for `ZZ_I` it uses `dmp_sqf_list` whereas for `ZZ` it does not.

To be clear the expected output for sqf is like:
```python
In [3]: sqf(expand((x-1)*(x-2)*(x-3)**2*(x-4)**5))
Out[3]: 
       5        2 ⎛ 2          ⎞
(x - 4) ⋅(x - 3) ⋅⎝x  - 3⋅x + 2⎠
```
It is supposed to be a factorisation of the expression like with `factor` except that it is not a full factorisation into irreducibles.
This fixes the bug in `dmp_sqf_list` makes all examples above with `sqf` work as intended:
```diff
diff --git a/sympy/polys/sqfreetools.py b/sympy/polys/sqfreetools.py
index e11eebc7b7..2fd4f0b8ca 100644
--- a/sympy/polys/sqfreetools.py
+++ b/sympy/polys/sqfreetools.py
@@ -23,7 +23,7 @@
 from sympy.polys.euclidtools import (
     dup_inner_gcd, dmp_inner_gcd,
     dup_gcd, dmp_gcd,
-    dmp_resultant)
+    dmp_resultant, dmp_primitive)
 from sympy.polys.galoistools import (
     gf_sqf_list, gf_sqf_part)
 from sympy.polys.polyerrors import (
@@ -401,30 +401,36 @@ def dmp_sqf_list(f, u, K, all=False):
     deg = dmp_degree(f, u)
     if deg < 0:
         return coeff, []
-    elif deg == 0:
-        coeff2, factors = dmp_sqf_list(dmp_LC(f, u), u-1, K, all=all)
-        factors = [([fac], exp) for fac, exp in factors]
-        return coeff*coeff2, factors
 
-    result, i = [], 1
+    content, f = dmp_primitive(f, u, K)
 
-    h = dmp_diff(f, 1, u, K)
-    g, p, q = dmp_inner_gcd(f, h, u, K)
+    result = []
 
-    while True:
-        d = dmp_diff(p, 1, u, K)
-        h = dmp_sub(q, d, u, K)
+    if deg != 0:
 
-        if dmp_zero_p(h, u):
-            result.append((p, i))
-            break
+        h = dmp_diff(f, 1, u, K)
+        g, p, q = dmp_inner_gcd(f, h, u, K)
 
-        g, p, q = dmp_inner_gcd(p, h, u, K)
+        i = 1
 
-        if all or dmp_degree(g, u) > 0:
-            result.append((g, i))
+        while True:
+            d = dmp_diff(p, 1, u, K)
+            h = dmp_sub(q, d, u, K)
 
-        i += 1
+            if dmp_zero_p(h, u):
+                result.append((p, i))
+                break
+
+            g, p, q = dmp_inner_gcd(p, h, u, K)
+
+            if all or dmp_degree(g, u) > 0:
+                result.append((g, i))
+
+            i += 1
+
+    coeff_content, result_content = dmp_sqf_list(content, u-1, K, all=all)
+    coeff *= coeff_content
+    result += [([fac], exp) for fac, exp in result_content]
 
     return coeff, result
```
Unfortunately that does not fix the original problem. It turns out that `dmp_sqf_list` is not even being called but rather `dmp_sqf_part` which does seem to return correct output:
```python
In [3]: a, b = symbols("a b")

In [4]: e = ((b - I)**2 * (b + I) * (a + 1))

In [5]: e
Out[5]: 
               2        
(a + 1)⋅(b - ⅈ) ⋅(b + ⅈ)

In [6]: expand(e)
Out[6]: 
   3        2                3      2        
a⋅b  - ⅈ⋅a⋅b  + a⋅b - ⅈ⋅a + b  - ⅈ⋅b  + b - ⅈ

In [7]: sqf_part(expand(e))
Out[7]: 
   2        2    
a⋅b  + a + b  + 1

In [8]: sqf_part(expand(e)).factor()
Out[8]: 
        ⎛ 2    ⎞
(a + 1)⋅⎝b  + 1⎠

In [9]: expand(e).factor() # wrong
Out[9]: 
        ⎛ 2    ⎞
(a + 1)⋅⎝b  + 1⎠
```
Here the output from `sqf_part` is correct for `sqf_part` but the same output is not correct for `factor`. The problem then is that factor is returning the `sqf_part`.

The bug is somewhere in `dmp_ext_factor`:
https://github.com/sympy/sympy/blob/87bf9d4343202ee2d4d544a0ee8c6cdfe36feb6a/sympy/polys/factortools.py#L1280-L1307
I'm not familiar with the algorithm in use here. Again it looks very similar to the corresponding univariate algorithm. In the debugger we already have a wrong output from `dmp_factor_list_include` because we only get two factors when there should be three. The output of `dmp_factor_list` is correct for the given input so the problem needs to be earlier. The input `f` seems correct and `dmp_sqf_part` gets the right square free part.
Maybe the bug is triggered because the square free part has only real coefficients i.e. we're trying to factorise a polynomial in `QQ(I)[a,b]` but its square free part is in `QQ[a,b]`:
```python
In [11]: e
Out[11]: 
               2        
(a + 1)⋅(b - ⅈ) ⋅(b + ⅈ)

In [12]: sqf_part(e)
Out[12]: 
   2        2    
a⋅b  + a + b  + 1
```
It should be possible for `dmp_trial_division` to use the factors that are found to find all of the factors e.g.:
```python
In [20]: quo(e, sqf_part(expand(e)))
Out[20]: b - ⅈ
```
I'm not sure that's how the algorithm is supposed to work though. I don't actually know where this particular algorithm comes from.
The algorithm in `dup_ext_factor` seems to be algorithm 3.6.4 in Cohen's "A course in computational number theory":
http://tomlr.free.fr/Math%E9matiques/Math%20Complete/Number%20theory/A%20course%20in%20computational%20algebraic%20number%20theory%20-%20Cohen%20H..pdf

That algorithm is written for univariate polynomials though. I'm not sure if `dmp_ext_factor` is a valid generalisation of the algorithm to the multivariate case.
Possibly the `dmp_ext_factor` algorithm is from:

Algebraic factoring and rational function integration
Barry M. Trager 1976

https://dl.acm.org/doi/10.1145/800205.806338
The problem with `sqf` was fixed in gh-26514 but the problem with `factor` (`dmp_ext_factor`) remains.
The algorithm in `dmp_ext_factor` is the one in the Trager paper looks correct per the paper. It also calls `dmp_sqf_norm` which is also defined in the paper and also looks correct. However `dmp_sqf_norm` calls `dmp_sqf_p` to check if the polynomial is square free and that seems to have a similar bug to the one that `dmp_sqf_list` had:
https://github.com/sympy/sympy/blob/af0c8559a845df49a6d0f8855659cebd0ddc9bc8/sympy/polys/sqfreetools.py#L55-L74
It is not enough just to differentiate wrt one arbitrarily chosen variable in the multivariate case.

Hence:
```python
In [22]: p = x**2*y**4 + 2*x**2*y**2 + x**2 + 2*x*y**4 + 4*x*y**2 + 2*x + 2*y**4 + 4*y**2 + 2

In [23]: p
Out[23]: 
 2  4      2  2    2        4        2            4      2    
x ⋅y  + 2⋅x ⋅y  + x  + 2⋅x⋅y  + 4⋅x⋅y  + 2⋅x + 2⋅y  + 4⋅y  + 2

In [24]: p.factor()
Out[24]: 
        2               
⎛ 2    ⎞  ⎛ 2          ⎞
⎝y  + 1⎠ ⋅⎝x  + 2⋅x + 2⎠

In [25]: Poly(p).is_sqf  # wrong
Out[25]: True
```
I think that wrong result from `dmp_sqf_p` leads to the wrong result in `dmp_sqf_norm` and in turn `dmp_ext_factor`.

Commit message: polys: rename PolyElement.norm_algebraic -> norm

Patch to review:
diff --git a/doc/src/modules/polys/internals.rst b/doc/src/modules/polys/internals.rst
index 31c8659be326..25d24f6ac28c 100644
--- a/doc/src/modules/polys/internals.rst
+++ b/doc/src/modules/polys/internals.rst
@@ -545,7 +545,9 @@ Polynomial factorization in characteristic zero:
 
 .. currentmodule:: sympy.polys.factortools
 
+.. autofunction:: dup_trial_division
 .. autofunction:: dmp_trial_division
+.. autofunction:: dup_zz_mignotte_bound
 .. autofunction:: dmp_zz_mignotte_bound
 .. autofunction:: dup_zz_hensel_step
 .. autofunction:: dup_zz_hensel_lift
@@ -559,16 +561,49 @@ Polynomial factorization in characteristic zero:
 .. autofunction:: dmp_zz_wang_non_divisors
 .. autofunction:: dmp_zz_wang_test_points
 .. autofunction:: dmp_zz_wang_lead_coeffs
+.. autofunction:: dup_zz_diophantine
 .. autofunction:: dmp_zz_diophantine
 .. autofunction:: dmp_zz_wang_hensel_lifting
 .. autofunction:: dmp_zz_wang
 .. autofunction:: dmp_zz_factor
+.. autofunction:: dup_qq_i_factor
+.. autofunction:: dup_zz_i_factor
+.. autofunction:: dmp_qq_i_factor
+.. autofunction:: dmp_zz_i_factor
+.. autofunction:: dup_ext_factor
 .. autofunction:: dmp_ext_factor
 .. autofunction:: dup_gf_factor
+.. autofunction:: dmp_gf_factor
+.. autofunction:: dup_factor_list
+.. autofunction:: dup_factor_list_include
 .. autofunction:: dmp_factor_list
 .. autofunction:: dmp_factor_list_include
+.. autofunction:: dup_irreducible_p
 .. autofunction:: dmp_irreducible_p
 
+Square-free factorization:
+
+.. currentmodule:: sympy.polys.sqfreetools
+
+.. autofunction:: dup_sqf_p
+.. autofunction:: dmp_sqf_p
+.. autofunction:: dup_sqf_norm
+.. autofunction:: dmp_sqf_norm
+.. autofunction:: dmp_norm
+.. autofunction:: dup_gf_sqf_part
+.. autofunction:: dmp_gf_sqf_part
+.. autofunction:: dup_sqf_part
+.. autofunction:: dmp_sqf_part
+.. autofunction:: dup_gf_sqf_list
+.. autofunction:: dmp_gf_sqf_list
+.. autofunction:: dup_sqf_list
+.. autofunction:: dup_sqf_list_include
+.. autofunction:: dmp_sqf_list
+.. autofunction:: dmp_sqf_list_include
+.. autofunction:: dup_gff_list
+.. autofunction:: dmp_gff_list
+
+
 Groebner basis algorithms
 *************************
 
diff --git a/doc/src/modules/polys/literature.rst b/doc/src/modules/polys/literature.rst
index 47fa9453baa9..5cfc3336ed68 100644
--- a/doc/src/modules/polys/literature.rst
+++ b/doc/src/modules/polys/literature.rst
@@ -144,4 +144,16 @@ a theoretical foundation for implementing polynomials manipulation module.
     https://isc.tamu.edu/resources/preprints/1996/1996-02.pdf
 
 .. [Cohen93] Henri Cohen. "A Course in Computational Algebraic Number Theory",
-   Springer, 1993.
+    Springer, 1993.
+
+.. [Trager76] Barry M. Trager. "Algebraic factoriing and rational function
+    integration", Proceedings of SYMSAC 1976, pp. 219-226, ACM, 1976.
+    https://dl.acm.org/doi/abs/10.1145/800205.806338
+
+.. [Yun76] David Y.Y. Yun. "On square-free decomposition algorithms",
+    Proceedings of SYMSAC 1976, pp. 219-226, ACM, 1976.
+    https://dl.acm.org/doi/10.1145/800205.806320
+
+.. [Abbott13] John Abbott. "Bounds on factors in Z[x]".
+    Journal of Symbolic Computation 50 (2013), pp. 532-563
+    https://doi.org/10.1016/j.jsc.2012.09.004
diff --git a/sympy/polys/compatibility.py b/sympy/polys/compatibility.py
index 6635c61e569c..ee71d9efedae 100644
--- a/sympy/polys/compatibility.py
+++ b/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
@@ -209,9 +210,10 @@
 from sympy.polys.rootisolation import dup_isolate_all_roots
 
 from sympy.polys.sqfreetools import (
-    dup_sqf_p, dmp_sqf_p, dup_sqf_norm, dmp_sqf_norm, dup_gf_sqf_part, dmp_gf_sqf_part,
-    dup_sqf_part, dmp_sqf_part, dup_gf_sqf_list, dmp_gf_sqf_list, dup_sqf_list,
-    dup_sqf_list_include, dmp_sqf_list, dmp_sqf_list_include, dup_gff_list, dmp_gff_list)
+    dup_sqf_p, dmp_sqf_p, dmp_norm, dup_sqf_norm, dmp_sqf_norm,
+    dup_gf_sqf_part, dmp_gf_sqf_part, dup_sqf_part, dmp_sqf_part,
+    dup_gf_sqf_list, dmp_gf_sqf_list, dup_sqf_list, dup_sqf_list_include,
+    dmp_sqf_list, dmp_sqf_list_include, dup_gff_list, dmp_gff_list)
 
 from sympy.polys.galoistools import (
     gf_degree, gf_LC, gf_TC, gf_strip, gf_from_dict,
@@ -515,6 +517,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))
 
@@ -877,6 +881,10 @@ def dup_sqf_p(self, f):
     def dmp_sqf_p(self, f):
         return dmp_sqf_p(self.to_dense(f), self.ngens-1, self.domain)
 
+    def dmp_norm(self, f):
+        n = dmp_norm(self.to_dense(f), self.ngens-1, self.domain)
+        return self.to_ground().from_dense(n)
+
     def dup_sqf_norm(self, f):
         s, F, R = dup_sqf_norm(self.to_dense(f), self.domain)
         return (s, self.from_dense(F), self.to_ground().from_dense(R))
diff --git a/sympy/polys/densetools.py b/sympy/polys/densetools.py
index 30760c6559c6..0c9db99e9e68 100644
--- a/sympy/polys/densetools.py
+++ b/sympy/polys/densetools.py
@@ -782,7 +782,7 @@ def dmp_ground_extract(f, g, u, K):
 
 def dup_real_imag(f, K):
     """
-    Return bivariate polynomials ``f1`` and ``f2``, such that ``f = f1 + f2*I``.
+    Find ``f1`` and ``f2``, such that ``f(x+I*y) = f1(x,y) + f2(x,y)*I``.
 
     Examples
     ========
@@ -793,6 +793,11 @@ def dup_real_imag(f, K):
     >>> R.dup_real_imag(x**3 + x**2 + x + 1)
     (x**3 + x**2 - 3*x*y**2 + x - y**2 + 1, 3*x**2*y + 2*x*y - y**3 + y)
 
+    >>> from sympy.abc import x, y, z
+    >>> from sympy import I
+    >>> (z**3 + z**2 + z + 1).subs(z, x+I*y).expand().collect(I)
+    x**3 + x**2 - 3*x*y**2 + x - y**2 + I*(3*x**2*y + 2*x*y - y**3 + y) + 1
+
     """
     if not K.is_ZZ and not K.is_QQ:
         raise DomainError("computing real and imaginary parts is not supported over %s" % K)
@@ -894,6 +899,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]``.
diff --git a/sympy/polys/factortools.py b/sympy/polys/factortools.py
index c0316375ba1a..dc477f02e672 100644
--- a/sympy/polys/factortools.py
+++ b/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,
@@ -65,7 +64,9 @@
 from sympy.polys.sqfreetools import (
     dup_sqf_p,
     dup_sqf_norm, dmp_sqf_norm,
-    dup_sqf_part, dmp_sqf_part)
+    dup_sqf_part, dmp_sqf_part,
+    _dup_check_degrees, _dmp_check_degrees,
+    )
 
 from sympy.polys.polyutils import _sort_factors
 from sympy.polys.polyconfig import query
@@ -125,6 +126,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)
@@ -133,7 +137,7 @@ def dmp_trial_division(f, factors, u, K):
 def dup_zz_mignotte_bound(f, K):
     """
     The Knuth-Cohen variant of Mignotte bound for
-    univariate polynomials in `K[x]`.
+    univariate polynomials in ``K[x]``.
 
     Examples
     ========
@@ -145,17 +149,18 @@ def dup_zz_mignotte_bound(f, K):
     >>> R.dup_zz_mignotte_bound(f)
     152
 
-    By checking `factor(f)` we can see that max coeff is 8
+    By checking ``factor(f)`` we can see that max coeff is 8
 
-    Also consider a case that `f` is irreducible for example `f = 2*x**2 + 3*x + 4`
-    To avoid a bug for these cases, we return the bound plus the max coefficient of `f`
+    Also consider a case that ``f`` is irreducible for example
+    ``f = 2*x**2 + 3*x + 4``. To avoid a bug for these cases, we return the
+    bound plus the max coefficient of ``f``
 
     >>> f = 2*x**2 + 3*x + 4
     >>> R.dup_zz_mignotte_bound(f)
     6
 
-    Lastly,To see the difference between the new and the old Mignotte bound
-    consider the irreducible polynomial::
+    Lastly, to see the difference between the new and the old Mignotte bound
+    consider the irreducible polynomial:
 
     >>> f = 87*x**7 + 4*x**6 + 80*x**5 + 17*x**4 + 9*x**3 + 12*x**2 + 49*x + 26
     >>> R.dup_zz_mignotte_bound(f)
@@ -167,7 +172,7 @@ def dup_zz_mignotte_bound(f, K):
     References
     ==========
 
-    ..[1] [Abbott2013]_
+    ..[1] [Abbott13]_
 
     """
     from sympy.functions.combinatorial.factorials import binomial
@@ -702,6 +707,9 @@ def dup_zz_factor(f, K):
         H = dup_zz_zassenhaus(g, K)
 
     factors = dup_trial_division(f, H, K)
+
+    _dup_check_degrees(f, factors)
+
     return cont, factors
 
 
@@ -1177,6 +1185,8 @@ def dmp_zz_factor(f, u, K):
     for g, k in dmp_zz_factor(G, u - 1, K)[1]:
         factors.insert(0, ([g], k))
 
+    _dmp_check_degrees(f, u, factors)
+
     return cont, _sort_factors(factors)
 
 
@@ -1247,7 +1257,71 @@ def dmp_zz_i_factor(f, u, K0):
 
 
 def dup_ext_factor(f, K):
-    """Factor univariate polynomials over algebraic number fields. """
+    r"""Factor univariate polynomials over algebraic number fields.
+
+    The domain `K` must be an algebraic number field `k(a)` (see :ref:`QQ(a)`).
+
+    Examples
+    ========
+
+    First define the algebraic number field `K = \mathbb{Q}(\sqrt{2})`:
+
+    >>> from sympy import QQ, sqrt
+    >>> from sympy.polys.factortools import dup_ext_factor
+    >>> K = QQ.algebraic_field(sqrt(2))
+
+    We can now factorise the polynomial `x^2 - 2` over `K`:
+
+    >>> p = [K(1), K(0), K(-2)] # x^2 - 2
+    >>> p1 = [K(1), -K.unit]    # x - sqrt(2)
+    >>> p2 = [K(1), +K.unit]    # x + sqrt(2)
+    >>> dup_ext_factor(p, K) == (K.one, [(p1, 1), (p2, 1)])
+    True
+
+    Usually this would be done at a higher level:
+
+    >>> from sympy import factor
+    >>> from sympy.abc import x
+    >>> factor(x**2 - 2, extension=sqrt(2))
+    (x - sqrt(2))*(x + sqrt(2))
+
+    Explanation
+    ===========
+
+    Uses Trager's algorithm. In particular this function is algorithm
+    ``alg_factor`` from [Trager76]_.
+
+    If `f` is a polynomial in `k(a)[x]` then its norm `g(x)` is a polynomial in
+    `k[x]`. If `g(x)` is square-free and has irreducible factors `g_1(x)`,
+    `g_2(x)`, `\cdots` then the irreducible factors of `f` in `k(a)[x]` are
+    given by `f_i(x) = \gcd(f(x), g_i(x))` where the GCD is computed in
+    `k(a)[x]`.
+
+    The first step in Trager's algorithm is to find an integer shift `s` so
+    that `f(x-sa)` has square-free norm. Then the norm is factorized in `k[x]`
+    and the GCD of (shifted) `f` with each factor gives the shifted factors of
+    `f`. At the end the shift is undone to recover the unshifted factors of `f`
+    in `k(a)[x]`.
+
+    The algorithm reduces the problem of factorization in `k(a)[x]` to
+    factorization in `k[x]` with the main additional steps being to compute the
+    norm (a resultant calculation in `k[x,y]`) and some polynomial GCDs in
+    `k(a)[x]`.
+
+    In practice in SymPy the base field `k` will be the rationals :ref:`QQ` and
+    this function factorizes a polynomial with coefficients in an algebraic
+    number field  like `\mathbb{Q}(\sqrt{2})`.
+
+    See Also
+    ========
+
+    dmp_ext_factor:
+        Analogous function for multivariate polynomials over ``k(a)``.
+    dup_sqf_norm:
+        Subroutine ``sqfr_norm`` also from [Trager76]_.
+    sympy.polys.polytools.factor:
+        The high-level function that ultimately uses this function as needed.
+    """
     n, lc = dup_degree(f), dup_LC(f, K)
 
     f = dup_monic(f, K)
@@ -1274,11 +1348,59 @@ def dup_ext_factor(f, K):
         factors[i] = h
 
     factors = dup_trial_division(F, factors, K)
+
+    _dup_check_degrees(F, factors)
+
     return lc, factors
 
 
 def dmp_ext_factor(f, u, K):
-    """Factor multivariate polynomials over algebraic number fields. """
+    r"""Factor multivariate polynomials over algebraic number fields.
+
+    The domain `K` must be an algebraic number field `k(a)` (see :ref:`QQ(a)`).
+
+    Examples
+    ========
+
+    First define the algebraic number field `K = \mathbb{Q}(\sqrt{2})`:
+
+    >>> from sympy import QQ, sqrt
+    >>> from sympy.polys.factortools import dmp_ext_factor
+    >>> K = QQ.algebraic_field(sqrt(2))
+
+    We can now factorise the polynomial `x^2 y^2 - 2` over `K`:
+
+    >>> p = [[K(1),K(0),K(0)], [], [K(-2)]] # x**2*y**2 - 2
+    >>> p1 = [[K(1),K(0)], [-K.unit]]       # x*y - sqrt(2)
+    >>> p2 = [[K(1),K(0)], [+K.unit]]       # x*y + sqrt(2)
+    >>> dmp_ext_factor(p, 1, K) == (K.one, [(p1, 1), (p2, 1)])
+    True
+
+    Usually this would be done at a higher level:
+
+    >>> from sympy import factor
+    >>> from sympy.abc import x, y
+    >>> factor(x**2*y**2 - 2, extension=sqrt(2))
+    (x*y - sqrt(2))*(x*y + sqrt(2))
+
+    Explanation
+    ===========
+
+    This is Trager's algorithm for multivariate polynomials. In particular this
+    function is algorithm ``alg_factor`` from [Trager76]_.
+
+    See :func:`dup_ext_factor` for explanation.
+
+    See Also
+    ========
+
+    dmp_ext_factor:
+        Analogous function for multivariate polynomials over ``k(a)``.
+    dmp_sqf_norm:
+        Multivariate version of subroutine ``sqfr_norm`` also from [Trager76]_.
+    sympy.polys.polytools.factor:
+        The high-level function that ultimately uses this function as needed.
+    """
     if not u:
         return dup_ext_factor(f, K)
 
@@ -1296,15 +1418,18 @@ 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 = [si*K.unit for si in s]
+            h = dmp_shift(h, a, u, K)
             factors[i] = h
 
-    return lc, dmp_trial_division(F, factors, u, K)
+    result = dmp_trial_division(F, factors, u, K)
+
+    _dmp_check_degrees(F, u, result)
+
+    return lc, result
 
 
 def dup_gf_factor(f, K):
diff --git a/sympy/polys/numberfields/subfield.py b/sympy/polys/numberfields/subfield.py
index 4c6a1e147eb5..b959ddeb27a6 100644
--- a/sympy/polys/numberfields/subfield.py
+++ b/sympy/polys/numberfields/subfield.py
@@ -356,7 +356,7 @@ def primitive_element(extension, x=None, *, ex=False, polys=False):
                 continue
             _, factors = factor_list(g, extension=ext)
             g = _choose_factor(factors, x, gen)
-            s, _, g = g.sqf_norm()
+            [s], _, g = g.sqf_norm()
             gen += s*ext
             coeffs.append(s)
 
@@ -378,7 +378,7 @@ def primitive_element(extension, x=None, *, ex=False, polys=False):
         L = QQ.algebraic_field((p, ext))
         _, factors = factor_list(f, domain=L)
         f = _choose_factor(factors, x, gen)
-        s, g, f = f.sqf_norm()
+        [s], g, f = f.sqf_norm()
         gen += s*ext
         coeffs.append(s)
         K = QQ.algebraic_field((f, gen))
diff --git a/sympy/polys/polyclasses.py b/sympy/polys/polyclasses.py
index 058f32681e22..c19f72a227f8 100644
--- a/sympy/polys/polyclasses.py
+++ b/sympy/polys/polyclasses.py
@@ -81,6 +81,7 @@
     dmp_compose,
     dup_decompose,
     dup_shift,
+    dmp_shift,
     dup_transform,
     dmp_lift)
 
@@ -897,6 +898,11 @@ def shift(f, a):
 
         return f._shift(f.dom.convert(a))
 
+    def shift_list(f, a):
+        """Efficiently compute Taylor shift ``f(X + A)``. """
+        a = [f.dom.convert(ai) for ai in a]
+        return f._shift_list(a)
+
     def _shift(f, a):
         raise NotImplementedError
 
@@ -1574,6 +1580,10 @@ def _shift(f, a):
         """Efficiently compute Taylor shift ``f(x + a)``. """
         return f.per(dup_shift(f._rep, a, f.dom))
 
+    def _shift_list(f, a):
+        """Efficiently compute Taylor shift ``f(X + A)``. """
+        return f.per(dmp_shift(f._rep, a, f.lev, f.dom))
+
     def _transform(f, p, q):
         """Evaluate functional transformation ``q**n * f(p/q)``."""
         return f.per(dup_transform(f._rep, p._rep, q._rep, f.dom))
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index b17310651072..80f317d6798a 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -3102,13 +3102,32 @@ def shift(f, a):
         >>> Poly(x**2 - 2*x + 1, x).shift(2)
         Poly(x**2 + 2*x + 1, x, domain='ZZ')
 
+        See Also
+        ========
+
+        shift_list: Analogous method for multivariate polynomials.
         """
-        if hasattr(f.rep, 'shift'):
-            result = f.rep.shift(a)
-        else:  # pragma: no cover
-            raise OperationNotSupported(f, 'shift')
+        return f.per(f.rep.shift(a))
 
-        return f.per(result)
+    def shift_list(f, a):
+        """
+        Efficiently compute Taylor shift ``f(X + A)``.
+
+        Examples
+        ========
+
+        >>> from sympy import Poly
+        >>> from sympy.abc import x, y
+
+        >>> Poly(x*y, [x,y]).shift_list([1, 2]) == Poly((x+1)*(y+2), [x,y])
+        True
+
+        See Also
+        ========
+
+        shift: Analogous method for univariate polynomials.
+        """
+        return f.per(f.rep.shift_list(a))
 
     def transform(f, p, q):
         """
@@ -3240,7 +3259,7 @@ def sqf_norm(f):
         >>> s, f, r = Poly(x**2 + 1, x, extension=[sqrt(3)]).sqf_norm()
 
         >>> s
-        1
+        [1]
         >>> f
         Poly(x**2 - 2*sqrt(3)*x + 4, x, domain='QQ<sqrt(3)>')
         >>> r
@@ -6042,7 +6061,7 @@ def sqf_norm(f, *gens, **args):
     >>> from sympy.abc import x
 
     >>> sqf_norm(x**2 + 1, extension=[sqrt(3)])
-    (1, x**2 - 2*sqrt(3)*x + 4, x**4 - 4*x**2 + 16)
+    ([1], x**2 - 2*sqrt(3)*x + 4, x**4 - 4*x**2 + 16)
 
     """
     options.allowed_flags(args, ['polys'])
@@ -6054,10 +6073,12 @@ def sqf_norm(f, *gens, **args):
 
     s, g, r = F.sqf_norm()
 
+    s_expr = [Integer(si) for si in s]
+
     if not opt.polys:
-        return Integer(s), g.as_expr(), r.as_expr()
+        return s_expr, g.as_expr(), r.as_expr()
     else:
-        return Integer(s), g, r
+        return s_expr, g, r
 
 
 @public
@@ -6554,11 +6575,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
diff --git a/sympy/polys/rings.py b/sympy/polys/rings.py
index 57b0e92e6946..7e3db5d3997b 100644
--- a/sympy/polys/rings.py
+++ b/sympy/polys/rings.py
@@ -2983,7 +2983,10 @@ def shift(f, a):
         if f.ring.is_univariate:
             return f.ring.dup_shift(f, a)
         else:
-            raise MultivariatePolynomialError("polynomial shift")
+            raise MultivariatePolynomialError("shift: use shift_list instead")
+
+    def shift_list(f, a):
+        return f.ring.dmp_shift(f, a)
 
     def sturm(f):
         if f.ring.is_univariate:
@@ -2994,6 +2997,9 @@ def sturm(f):
     def gff_list(f):
         return f.ring.dmp_gff_list(f)
 
+    def norm(f):
+        return f.ring.dmp_norm(f)
+
     def sqf_norm(f):
         return f.ring.dmp_sqf_norm(f)
 
diff --git a/sympy/polys/sqfreetools.py b/sympy/polys/sqfreetools.py
index a27b6f573aa9..58dbec904ee7 100644
--- a/sympy/polys/sqfreetools.py
+++ b/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_degree_list,
     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 (
@@ -30,6 +30,22 @@
     MultivariatePolynomialError,
     DomainError)
 
+
+def _dup_check_degrees(f, result):
+    """Sanity check the degrees of a computed factorization in K[x]."""
+    deg = sum(k * dup_degree(fac) for (fac, k) in result)
+    assert deg == dup_degree(f)
+
+
+def _dmp_check_degrees(f, u, result):
+    """Sanity check the degrees of a computed factorization in K[X]."""
+    degs = [0] * (u + 1)
+    for fac, k in result:
+        degs_fac = dmp_degree_list(fac, u)
+        degs = [d1 + k * d2 for d1, d2 in zip(degs, degs_fac)]
+    assert tuple(degs) == dmp_degree_list(f, u)
+
+
 def dup_sqf_p(f, K):
     """
     Return ``True`` if ``f`` is a square-free polynomial in ``K[x]``.
@@ -70,20 +86,37 @@ 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):
-    """
-    Square-free norm of ``f`` in ``K[x]``, useful over algebraic domains.
+    r"""
+    Find a shift of `f` in `K[x]` that has square-free norm.
+
+    The domain `K` must be an algebraic number field `k(a)` (see :ref:`QQ(a)`).
 
-    Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and ``r(x) = Norm(g(x))``
-    is a square-free polynomial over K, where ``a`` is the algebraic extension of ``K``.
+    Returns `(s,g,r)`, such that `g(x)=f(x-sa)`, `r(x)=\text{Norm}(g(x))` and
+    `r` is a square-free polynomial over `k`.
 
     Examples
     ========
 
+    We first create the algebraic number field `K=k(a)=\mathbb{Q}(\sqrt{3})`
+    and rings `K[x]` and `k[x]`:
+
     >>> from sympy.polys import ring, QQ
     >>> from sympy import sqrt
 
@@ -91,15 +124,48 @@ def dup_sqf_norm(f, K):
     >>> R, x = ring("x", K)
     >>> _, X = ring("x", QQ)
 
-    >>> s, f, r = R.dup_sqf_norm(x**2 - 2)
+    We can now find a square free norm for a shift of `f`:
+
+    >>> f = x**2 - 1
+    >>> s, g, r = R.dup_sqf_norm(f)
+
+    The choice of shift `s` is arbitrary and the particular values returned for
+    `g` and `r` are determined by `s`.
 
     >>> s == 1
     True
-    >>> f == x**2 + K([QQ(-2), QQ(0)])*x + 1
+    >>> g == x**2 - 2*sqrt(3)*x + 2
     True
-    >>> r == X**4 - 10*X**2 + 1
+    >>> r == X**4 - 8*X**2 + 4
     True
 
+    The invariants are:
+
+    >>> g == f.shift(-s*K.unit)
+    True
+    >>> g.norm() == r
+    True
+    >>> r.is_squarefree
+    True
+
+    Explanation
+    ===========
+
+    This is part of Trager's algorithm for factorizing polynomials over
+    algebraic number fields. In particular this function is algorithm
+    ``sqfr_norm`` from [Trager76]_.
+
+    See Also
+    ========
+
+    dmp_sqf_norm:
+        Analogous function for multivariate polynomials over ``k(a)``.
+    dmp_norm:
+        Computes the norm of `f` directly without any shift.
+    dup_ext_factor:
+        Function implementing Trager's algorithm that uses this.
+    sympy.polys.polytools.sqf_norm:
+        High-level interface for using this function.
     """
     if not K.is_Algebraic:
         raise DomainError("ground domain must be algebraic")
@@ -118,16 +184,62 @@ def dup_sqf_norm(f, K):
     return s, f, r
 
 
+def _dmp_sqf_norm_shifts(f, u, K):
+    """Generate a sequence of candidate shifts for dmp_sqf_norm."""
+    #
+    # We want to find a minimal shift if possible because shifting high degree
+    # variables can be expensive e.g. x**10 -> (x + 1)**10. We try a few easy
+    # cases first before the final infinite loop that is guaranteed to give
+    # only finitely many bad shifts (see Trager76 for proof of this in the
+    # univariate case).
+    #
+
+    # First the trivial shift [0, 0, ...]
+    n = u + 1
+    s0 = [0] * n
+    yield s0, f
+
+    # Shift in multiples of the generator of the extension field K
+    a = K.unit
+
+    # Variables of degree > 0 ordered by increasing degree
+    d = dmp_degree_list(f, u)
+    var_indices = [i for di, i in sorted(zip(d, range(u+1))) if di > 0]
+
+    # Now try [1, 0, 0, ...], [0, 1, 0, ...]
+    for i in var_indices:
+        s1 = s0.copy()
+        s1[i] = 1
+        a1 = [-a*s1i for s1i in s1]
+        f1 = dmp_shift(f, a1, u, K)
+        yield s1, f1
+
+    # Now try [1, 1, 1, ...], [2, 2, 2, ...]
+    j = 0
+    while True:
+        j += 1
+        sj = [j] * n
+        aj = [-a*j] * n
+        fj = dmp_shift(f, aj, u, K)
+        yield sj, fj
+
+
 def dmp_sqf_norm(f, u, K):
-    """
-    Square-free norm of ``f`` in ``K[X]``, useful over algebraic domains.
+    r"""
+    Find a shift of ``f`` in ``K[X]`` that has square-free norm.
+
+    The domain `K` must be an algebraic number field `k(a)` (see :ref:`QQ(a)`).
 
-    Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and ``r(x) = Norm(g(x))``
-    is a square-free polynomial over K, where ``a`` is the algebraic extension of ``K``.
+    Returns `(s,g,r)`, such that `g(x_1,x_2,\cdots)=f(x_1-s_1 a, x_2 - s_2 a,
+    \cdots)`, `r(x)=\text{Norm}(g(x))` and `r` is a square-free polynomial over
+    `k`.
 
     Examples
     ========
 
+    We first create the algebraic number field `K=k(a)=\mathbb{Q}(i)` and rings
+    `K[x]` and `k[x]`:
+
     >>> from sympy.polys import ring, QQ
     >>> from sympy import I
 
@@ -135,42 +247,140 @@ def dmp_sqf_norm(f, u, K):
     >>> R, x, y = ring("x,y", K)
     >>> _, X, Y = ring("x,y", QQ)
 
-    >>> s, f, r = R.dmp_sqf_norm(x*y + y**2)
+    We can now find a square free norm for a shift of `f`:
 
-    >>> s == 1
+    >>> f = x*y + y**2
+    >>> s, g, r = R.dmp_sqf_norm(f)
+
+    The choice of shifts ``s`` is arbitrary and the particular values returned
+    for ``g`` and ``r`` are determined by ``s``.
+
+    >>> s
+    [0, 1]
+    >>> g == x*y - I*x + y**2 - 2*I*y - 1
+    True
+    >>> r == X**2*Y**2 + X**2 + 2*X*Y**3 + 2*X*Y + Y**4 + 2*Y**2 + 1
     True
-    >>> f == x*y + y**2 + K([QQ(-1), QQ(0)])*y
+
+    The required invariants are:
+
+    >>> g == f.shift_list([-si*K.unit for si in s])
     True
-    >>> r == X**2*Y**2 + 2*X*Y**3 + Y**4 + Y**2
+    >>> g.norm() == r
     True
+    >>> r.is_squarefree
+    True
+
+    Explanation
+    ===========
+
+    This is part of Trager's algorithm for factorizing polynomials over
+    algebraic number fields. In particular this function is a multivariate
+    generalization of algorithm ``sqfr_norm`` from [Trager76]_.
+
+    See Also
+    ========
 
+    dup_sqf_norm:
+        Analogous function for multivariate polynomials over ``k(a)``.
+    dmp_norm:
+        Computes the norm of `f` directly without any shift.
+    dmp_ext_factor:
+        Function implementing Trager's algorithm that uses this.
+    sympy.polys.polytools.sqf_norm:
+        High-level interface for using this function.
     """
     if not u:
-        return dup_sqf_norm(f, K)
+        s, g, r = dup_sqf_norm(f, K)
+        return [s], g, r
 
     if not K.is_Algebraic:
         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)
 
-    s = 0
+    for s, f in _dmp_sqf_norm_shifts(f, u, K):
 
-    while True:
         h, _ = dmp_inject(f, u, K, front=True)
         r = dmp_resultant(g, h, u + 1, K.dom)
 
         if dmp_sqf_p(r, u, K.dom):
             break
-        else:
-            f, s = dmp_compose(f, F, u, K), s + 1
 
     return s, f, r
 
 
 def dmp_norm(f, u, K):
-    """
-    Norm of ``f`` in ``K[X1, ..., Xn]``, often not square-free.
+    r"""
+    Norm of ``f`` in ``K[X]``, often not square-free.
+
+    The domain `K` must be an algebraic number field `k(a)` (see :ref:`QQ(a)`).
+
+    Examples
+    ========
+
+    We first define the algebraic number field `K = k(a) = \mathbb{Q}(\sqrt{2})`:
+
+    >>> from sympy import QQ, sqrt
+    >>> from sympy.polys.sqfreetools import dmp_norm
+    >>> k = QQ
+    >>> K = k.algebraic_field(sqrt(2))
+
+    We can now compute the norm of a polynomial `p` in `K[x,y]`:
+
+    >>> p = [[K(1)], [K(1),K.unit]]                  # x + y + sqrt(2)
+    >>> N = [[k(1)], [k(2),k(0)], [k(1),k(0),k(-2)]] # x**2 + 2*x*y + y**2 - 2
+    >>> dmp_norm(p, 1, K) == N
+    True
+
+    In higher level functions that is:
+
+    >>> from sympy import expand, roots, minpoly
+    >>> from sympy.abc import x, y
+    >>> from math import prod
+    >>> a = sqrt(2)
+    >>> e = (x + y + a)
+    >>> e.as_poly([x, y], extension=a).norm()
+    Poly(x**2 + 2*x*y + y**2 - 2, x, y, domain='QQ')
+
+    This is equal to the product of the expressions `x + y + a_i` where the
+    `a_i` are the conjugates of `a`:
+
+    >>> pa = minpoly(a)
+    >>> pa
+    _x**2 - 2
+    >>> rs = roots(pa, multiple=True)
+    >>> rs
+    [sqrt(2), -sqrt(2)]
+    >>> n = prod(e.subs(a, r) for r in rs)
+    >>> n
+    (x + y - sqrt(2))*(x + y + sqrt(2))
+    >>> expand(n)
+    x**2 + 2*x*y + y**2 - 2
+
+    Explanation
+    ===========
+
+    Given an algebraic number field `K = k(a)` any element `b` of `K` can be
+    represented as polynomial function `b=g(a)` where `g` is in `k[x]`. If the
+    minimal polynomial of `a` over `k` is `p_a` then the roots `a_1`, `a_2`,
+    `\cdots` of `p_a(x)` are the conjugates of `a`. The norm of `b` is the
+    product `g(a1) \times g(a2) \times \cdots` and is an element of `k`.
+
+    As in [Trager76]_ we extend this norm to multivariate polynomials over `K`.
+    If `b(x)` is a polynomial in `k(a)[X]` then we can think of `b` as being
+    alternately a function `g_X(a)` where `g_X` is an element of `k[X][y]` i.e.
+    a polynomial function with coefficients that are elements of `k[X]`. Then
+    the norm of `b` is the product `g_X(a1) \times g_X(a2) \times \cdots` and
+    will be an element of `k[X]`.
+
+    See Also
+    ========
+
+    dmp_sqf_norm:
+        Compute a shift of `f` so that the `\text{Norm}(f)` is square-free.
+    sympy.polys.polytools.Poly.norm:
+        Higher-level function that calls this.
     """
     if not K.is_Algebraic:
         raise DomainError("ground domain must be algebraic")
@@ -206,6 +416,10 @@ def dup_sqf_part(f, K):
     >>> R.dup_sqf_part(x**3 - 3*x - 2)
     x**2 - x - 2
 
+    See Also
+    ========
+
+    sympy.polys.polytools.Poly.sqf_part
     """
     if K.is_FiniteField:
         return dup_gf_sqf_part(f, K)
@@ -264,6 +478,8 @@ def dmp_sqf_part(f, u, K):
 
 def dup_gf_sqf_list(f, K, all=False):
     """Compute square-free decomposition of ``f`` in ``GF(p)[x]``. """
+    f_orig = f
+
     f = dup_convert(f, K, K.dom)
 
     coeff, factors = gf_sqf_list(f, K.mod, K.dom, all=all)
@@ -271,6 +487,8 @@ def dup_gf_sqf_list(f, K, all=False):
     for i, (f, k) in enumerate(factors):
         factors[i] = (dup_convert(f, K.dom, K), k)
 
+    _dup_check_degrees(f_orig, factors)
+
     return K.convert(coeff, K.dom), factors
 
 
@@ -283,6 +501,8 @@ def dup_sqf_list(f, K, all=False):
     """
     Return square-free decomposition of a polynomial in ``K[x]``.
 
+    Uses Yun's algorithm from [Yun76]_.
+
     Examples
     ========
 
@@ -296,10 +516,26 @@ def dup_sqf_list(f, K, all=False):
     >>> R.dup_sqf_list(f, all=True)
     (2, [(1, 1), (x + 1, 2), (x + 2, 3)])
 
+    See Also
+    ========
+
+    dmp_sqf_list:
+        Corresponding function for multivariate polynomials.
+    sympy.polys.polytools.sqf_list:
+        High-level function for square-free factorization of expressions.
+    sympy.polys.polytools.Poly.sqf_list:
+        Analogous method on :class:`~.Poly`.
+
+    References
+    ==========
+
+    [Yun76]_
     """
     if K.is_FiniteField:
         return dup_gf_sqf_list(f, K, all=all)
 
+    f_orig = f
+
     if K.is_Field:
         coeff = dup_LC(f, K)
         f = dup_monic(f, K)
@@ -333,6 +569,8 @@ def dup_sqf_list(f, K, all=False):
 
         i += 1
 
+    _dup_check_degrees(f_orig, result)
+
     return coeff, result
 
 
@@ -366,7 +604,7 @@ def dup_sqf_list_include(f, K, all=False):
 
 def dmp_sqf_list(f, u, K, all=False):
     """
-    Return square-free decomposition of a polynomial in ``K[X]``.
+    Return square-free decomposition of a polynomial in `K[X]`.
 
     Examples
     ========
@@ -381,6 +619,26 @@ def dmp_sqf_list(f, u, K, all=False):
     >>> R.dmp_sqf_list(f, all=True)
     (1, [(1, 1), (x + y, 2), (x, 3)])
 
+    Explanation
+    ===========
+
+    Uses Yun's algorithm for univariate polynomials from [Yun76]_ recrusively.
+    The multivariate polynomial is treated as a univariate polynomial in its
+    leading variable. Then Yun's algorithm computes the square-free
+    factorization of the primitive and the content is factored recursively.
+
+    It would be better to use a dedicated algorithm for multivariate
+    polynomials instead.
+
+    See Also
+    ========
+
+    dup_sqf_list:
+        Corresponding function for univariate polynomials.
+    sympy.polys.polytools.sqf_list:
+        High-level function for square-free factorization of expressions.
+    sympy.polys.polytools.Poly.sqf_list:
+        Analogous method on :class:`~.Poly`.
     """
     if not u:
         return dup_sqf_list(f, K, all=all)
@@ -388,6 +646,8 @@ def dmp_sqf_list(f, u, K, all=False):
     if K.is_FiniteField:
         return dmp_gf_sqf_list(f, u, K, all=all)
 
+    f_orig = f
+
     if K.is_Field:
         coeff = dmp_ground_LC(f, u, K)
         f = dmp_ground_monic(f, u, K)
@@ -445,6 +705,8 @@ def dmp_sqf_list(f, u, K, all=False):
 
     result = [(result[i], i) for i in sorted(result)]
 
+    _dmp_check_degrees(f_orig, u, result)
+
     return coeff, result
 
 
diff --git a/sympy/polys/tests/test_densetools.py b/sympy/polys/tests/test_densetools.py
index d9c9cf6e56b1..43dae691f52d 100644
--- a/sympy/polys/tests/test_densetools.py
+++ b/sympy/polys/tests/test_densetools.py
@@ -20,7 +20,7 @@
     dup_primitive, dmp_ground_primitive,
     dup_extract, dmp_ground_extract,
     dup_real_imag,
-    dup_mirror, dup_scale, dup_shift,
+    dup_mirror, dup_scale, dup_shift, dmp_shift,
     dup_transform,
     dup_compose, dmp_compose,
     dup_decompose,
@@ -39,7 +39,7 @@
 
 from sympy.polys.specialpolys import f_polys
 
-from sympy.polys.domains import FF, ZZ, QQ, EX, RR
+from sympy.polys.domains import FF, ZZ, QQ, ZZ_I, QQ_I, EX, RR
 from sympy.polys.rings import ring
 
 from sympy.core.numbers import I
@@ -75,6 +75,8 @@ def test_dup_integrate():
 
 
 def test_dmp_integrate():
+    assert dmp_integrate([QQ(1)], 2, 0, QQ) == [QQ(1, 2), QQ(0), QQ(0)]
+
     assert dmp_integrate([[[]]], 1, 2, QQ) == [[[]]]
     assert dmp_integrate([[[]]], 2, 2, QQ) == [[[]]]
 
@@ -107,6 +109,9 @@ def test_dmp_integrate_in():
         dmp_swap(
             dmp_integrate(dmp_swap(f, 0, 2, 3, QQ), 3, 3, QQ), 0, 2, 3, QQ)
 
+    raises(IndexError, lambda: dmp_integrate_in(f, 1, -1, 3, QQ))
+    raises(IndexError, lambda: dmp_integrate_in(f, 1, 4, 3, QQ))
+
 
 def test_dup_diff():
     assert dup_diff([], 1, ZZ) == []
@@ -180,6 +185,8 @@ def test_dmp_diff_in():
     assert dmp_diff_in(f_6, 3, 2, 3, ZZ) == \
         dmp_swap(dmp_diff(dmp_swap(f_6, 0, 2, 3, ZZ), 3, 3, ZZ), 0, 2, 3, ZZ)
 
+    raises(IndexError, lambda: dmp_diff_in(f_6, 1, -1, 3, ZZ))
+    raises(IndexError, lambda: dmp_diff_in(f_6, 1, 4, 3, ZZ))
 
 def test_dup_eval():
     assert dup_eval([], 7, ZZ) == 0
@@ -217,6 +224,9 @@ def test_dmp_eval_in():
     assert dmp_eval_in(f, -2, 2, 2, ZZ) == \
         [[45], [], [], [-9, -1, 0, -44]]
 
+    raises(IndexError, lambda: dmp_eval_in(f_6, ZZ(1), -1, 3, ZZ))
+    raises(IndexError, lambda: dmp_eval_in(f_6, ZZ(1), 4, 3, ZZ))
+
 
 def test_dmp_eval_tail():
     assert dmp_eval_tail([[]], [1], 1, ZZ) == []
@@ -248,6 +258,11 @@ def test_dmp_diff_eval_in():
     assert dmp_diff_eval_in(f_6, 2, 7, 1, 3, ZZ) == \
         dmp_eval(dmp_diff(dmp_swap(f_6, 0, 1, 3, ZZ), 2, 3, ZZ), 7, 3, ZZ)
 
+    assert dmp_diff_eval_in(f_6, 2, 7, 0, 3, ZZ) == \
+        dmp_eval(dmp_diff(f_6, 2, 3, ZZ), 7, 3, ZZ)
+
+    raises(IndexError, lambda: dmp_diff_eval_in(f_6, 1, ZZ(1), 4, 3, ZZ))
+
 
 def test_dup_revert():
     f = [-QQ(1, 720), QQ(0), QQ(1, 24), QQ(0), -QQ(1, 2), QQ(0), QQ(1)]
@@ -271,6 +286,12 @@ def test_dup_trunc():
     assert dup_trunc([1, 2, 3, 4, 5, 6], ZZ(3), ZZ) == [1, -1, 0, 1, -1, 0]
     assert dup_trunc([6, 5, 4, 3, 2, 1], ZZ(3), ZZ) == [-1, 1, 0, -1, 1]
 
+    R = ZZ_I
+    assert dup_trunc([R(3), R(4), R(5)], R(3), R) == [R(1), R(-1)]
+
+    K = FF(5)
+    assert dup_trunc([K(3), K(4), K(5)], K(3), K) == [K(1), K(0)]
+
 
 def test_dmp_trunc():
     assert dmp_trunc([[]], [1, 2], 2, ZZ) == [[]]
@@ -294,6 +315,8 @@ def test_dup_monic():
 
 
 def test_dmp_ground_monic():
+    assert dmp_ground_monic([3, 6, 9], 0, ZZ) == [1, 2, 3]
+
     assert dmp_ground_monic([[3], [6], [9]], 1, ZZ) == [[1], [2], [3]]
 
     raises(
@@ -386,6 +409,8 @@ def test_dup_primitive():
 
 
 def test_dmp_ground_primitive():
+    assert dmp_ground_primitive([ZZ(1)], 0, ZZ) == (ZZ(1), [ZZ(1)])
+
     assert dmp_ground_primitive([[]], 1, ZZ) == (ZZ(0), [[]])
 
     assert dmp_ground_primitive(f_0, 2, ZZ) == (ZZ(1), f_0)
@@ -456,9 +481,15 @@ def test_dup_real_imag():
     assert dup_real_imag(
         [1, 2, 3], ZZ) == ([[1], [2], [-1, 0, 3]], [[2, 0], [2, 0]])
 
+    assert dup_real_imag([ZZ(1), ZZ(0), ZZ(1), ZZ(3)], ZZ) == (
+        [[ZZ(1)], [], [ZZ(-3), ZZ(0), ZZ(1)], [ZZ(3)]],
+        [[ZZ(3), ZZ(0)], [], [ZZ(-1), ZZ(0), ZZ(1), ZZ(0)]]
+    )
+
     raises(DomainError, lambda: dup_real_imag([EX(1), EX(2)], EX))
 
 
+
 def test_dup_mirror():
     assert dup_mirror([], ZZ) == []
     assert dup_mirror([1], ZZ) == [1]
@@ -483,6 +514,16 @@ def test_dup_shift():
     assert dup_shift([1, 2, 3, 4, 5], 7, ZZ) == [1, 30, 339, 1712, 3267]
 
 
+def test_dmp_shift():
+    assert dmp_shift([ZZ(1), ZZ(2)], [ZZ(1)], 0, ZZ) == [ZZ(1), ZZ(3)]
+
+    assert dmp_shift([[]], [ZZ(1), ZZ(2)], 1, ZZ) == [[]]
+
+    xy = [[ZZ(1), ZZ(0)], []]               # x*y
+    x1y2 = [[ZZ(1), ZZ(2)], [ZZ(1), ZZ(2)]] # (x+1)*(y+2)
+    assert dmp_shift(xy, [ZZ(1), ZZ(2)], 1, ZZ) == x1y2
+
+
 def test_dup_transform():
     assert dup_transform([], [], [1, 1], ZZ) == []
     assert dup_transform([], [1], [1, 1], ZZ) == []
@@ -570,12 +611,17 @@ def test_dup_decompose():
 def test_dmp_lift():
     q = [QQ(1, 1), QQ(0, 1), QQ(1, 1)]
 
-    f = [ANP([QQ(1, 1)], q, QQ), ANP([], q, QQ), ANP([], q, QQ),
+    f_a = [ANP([QQ(1, 1)], q, QQ), ANP([], q, QQ), ANP([], q, QQ),
          ANP([QQ(1, 1), QQ(0, 1)], q, QQ), ANP([QQ(17, 1), QQ(0, 1)], q, QQ)]
 
-    assert dmp_lift(f, 0, QQ.algebraic_field(I)) == \
-        [QQ(1), QQ(0), QQ(0), QQ(0), QQ(0), QQ(0), QQ(2), QQ(0), QQ(578),
-         QQ(0), QQ(0), QQ(0), QQ(1), QQ(0), QQ(-578), QQ(0), QQ(83521)]
+    f_lift = [QQ(1), QQ(0), QQ(0), QQ(0), QQ(0), QQ(0), QQ(2), QQ(0), QQ(578),
+              QQ(0), QQ(0), QQ(0), QQ(1), QQ(0), QQ(-578), QQ(0), QQ(83521)]
+
+    assert dmp_lift(f_a, 0, QQ.algebraic_field(I)) == f_lift
+
+    f_g = [QQ_I(1), QQ_I(0), QQ_I(0), QQ_I(0, 1), QQ_I(0, 17)]
+
+    assert dmp_lift(f_g, 0, QQ_I) == f_lift
 
     raises(DomainError, lambda: dmp_lift([EX(1), EX(2)], 0, EX))
 
diff --git a/sympy/polys/tests/test_factortools.py b/sympy/polys/tests/test_factortools.py
index 84133d4137e4..7f99097c71e9 100644
--- a/sympy/polys/tests/test_factortools.py
+++ b/sympy/polys/tests/test_factortools.py
@@ -562,7 +562,10 @@ def anp(element):
 
 
 def test_dmp_ext_factor():
-    R, x,y = ring("x,y", QQ.algebraic_field(sqrt(2)))
+    K = QQ.algebraic_field(sqrt(2))
+    R, x,y = ring("x,y", K)
+    sqrt2 = K.unit
+
     def anp(x):
         return ANP(x, [QQ(1), QQ(0), QQ(-2)], QQ)
 
@@ -588,6 +591,12 @@ def anp(x):
         (anp([QQ(2)]), [(anp([QQ(1)])*x + anp([QQ(-1), QQ(0)])*y, 1),
                         (anp([QQ(1)])*x + anp([QQ( 1), QQ(0)])*y, 1)])
 
+    f1 = y + 1
+    f2 = y + sqrt2
+    f3 = x**2 + x + 2 + 3*sqrt2
+    f = f1**2 * f2**2 * f3**2
+    assert R.dmp_ext_factor(f) == (K.one, [(f1, 2), (f2, 2), (f3, 2)])
+
 
 def test_dup_factor_list():
     R, x = ring("x", ZZ)
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index 2bab24433e43..0cce69d0ffb1 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -2313,6 +2313,11 @@ def test_compose():
 def test_shift():
     assert Poly(x**2 - 2*x + 1, x).shift(2) == Poly(x**2 + 2*x + 1, x)
 
+
+def test_shift_list():
+    assert Poly(x*y, [x,y]).shift_list([1,2]) == Poly((x+1)*(y+2), [x,y])
+
+
 def test_transform():
     # Also test that 3-way unification is done correctly
     assert Poly(x**2 - 2*x + 1, x).transform(Poly(x + 1), Poly(x - 1)) == \
@@ -2397,17 +2402,17 @@ def test_norm():
 
 def test_sqf_norm():
     assert sqf_norm(x**2 - 2, extension=sqrt(3)) == \
-        (1, x**2 - 2*sqrt(3)*x + 1, x**4 - 10*x**2 + 1)
+        ([1], x**2 - 2*sqrt(3)*x + 1, x**4 - 10*x**2 + 1)
     assert sqf_norm(x**2 - 3, extension=sqrt(2)) == \
-        (1, x**2 - 2*sqrt(2)*x - 1, x**4 - 10*x**2 + 1)
+        ([1], x**2 - 2*sqrt(2)*x - 1, x**4 - 10*x**2 + 1)
 
     assert Poly(x**2 - 2, extension=sqrt(3)).sqf_norm() == \
-        (1, Poly(x**2 - 2*sqrt(3)*x + 1, x, extension=sqrt(3)),
-            Poly(x**4 - 10*x**2 + 1, x, domain='QQ'))
+        ([1], Poly(x**2 - 2*sqrt(3)*x + 1, x, extension=sqrt(3)),
+              Poly(x**4 - 10*x**2 + 1, x, domain='QQ'))
 
     assert Poly(x**2 - 3, extension=sqrt(2)).sqf_norm() == \
-        (1, Poly(x**2 - 2*sqrt(2)*x - 1, x, extension=sqrt(2)),
-            Poly(x**4 - 10*x**2 + 1, x, domain='QQ'))
+        ([1], Poly(x**2 - 2*sqrt(2)*x - 1, x, extension=sqrt(2)),
+              Poly(x**4 - 10*x**2 + 1, x, domain='QQ'))
 
 
 def test_sqf():
@@ -2693,6 +2698,24 @@ def test_factor():
     assert factor_list((x - sqrt(2)*pi)*(x + sqrt(2)*pi), x) == (
         1, [(x - sqrt(2)*pi, 1), (x + sqrt(2)*pi, 1)])
 
+    # https://github.com/sympy/sympy/issues/26497
+    p = ((y - I)**2 * (y + I) * (x + 1))
+    assert factor(expand(p)) == p
+
+    p = ((x - I)**2 * (x + I) * (y + 1))
+    assert factor(expand(p)) == p
+
+    p = (y + 1)**2*(y + sqrt(2))**2*(x**2 + x + 2 + 3*sqrt(2))**2
+    assert factor(expand(p), extension=True) == p
+
+    e = (
+        -x**2*y**4/(y**2 + 1) + 2*I*x**2*y**3/(y**2 + 1) + 2*I*x**2*y/(y**2 + 1) +
+        x**2/(y**2 + 1) - 2*x*y**4/(y**2 + 1) + 4*I*x*y**3/(y**2 + 1) +
+        4*I*x*y/(y**2 + 1) + 2*x/(y**2 + 1) - y**4 - y**4/(y**2 + 1) + 2*I*y**3
+        + 2*I*y**3/(y**2 + 1) + 2*I*y + 2*I*y/(y**2 + 1) + 1 + 1/(y**2 + 1)
+    )
+    assert factor(e) == -(y - I)**3*(y + I)*(x**2 + 2*x + y**2 + 2)/(y**2 + 1)
+
 
 def test_factor_large():
     f = (x**2 + 4*x + 4)**10000000*(x**2 + 1)*(x**2 + 2*x + 1)**1234567
diff --git a/sympy/polys/tests/test_rings.py b/sympy/polys/tests/test_rings.py
index 2f560922ea81..d2cc34da5a83 100644
--- a/sympy/polys/tests/test_rings.py
+++ b/sympy/polys/tests/test_rings.py
@@ -1490,6 +1490,12 @@ def test_PolyElement_decompose():
 def test_PolyElement_shift():
     _, x = ring("x", ZZ)
     assert (x**2 - 2*x + 1).shift(2) == x**2 + 2*x + 1
+    assert (x**2 - 2*x + 1).shift_list([2]) == x**2 + 2*x + 1
+
+    R, x, y = ring("x, y", ZZ)
+    assert (x*y).shift_list([1, 2]) == (x+1)*(y+2)
+
+    raises(MultivariatePolynomialError, lambda: (x*y).shift(1))
 
 def test_PolyElement_sturm():
     F, t = field("t", ZZ)
@@ -1517,12 +1523,12 @@ def test_PolyElement_sqf_norm():
     R, x = ring("x", QQ.algebraic_field(sqrt(3)))
     X = R.to_ground().x
 
-    assert (x**2 - 2).sqf_norm() == (1, x**2 - 2*sqrt(3)*x + 1, X**4 - 10*X**2 + 1)
+    assert (x**2 - 2).sqf_norm() == ([1], x**2 - 2*sqrt(3)*x + 1, X**4 - 10*X**2 + 1)
 
     R, x = ring("x", QQ.algebraic_field(sqrt(2)))
     X = R.to_ground().x
 
-    assert (x**2 - 3).sqf_norm() == (1, x**2 - 2*sqrt(2)*x - 1, X**4 - 10*X**2 + 1)
+    assert (x**2 - 3).sqf_norm() == ([1], x**2 - 2*sqrt(2)*x - 1, X**4 - 10*X**2 + 1)
 
 def test_PolyElement_sqf_list():
     _, x = ring("x", ZZ)
diff --git a/sympy/polys/tests/test_sqfreetools.py b/sympy/polys/tests/test_sqfreetools.py
index 976b035163ed..b772a05a50e2 100644
--- a/sympy/polys/tests/test_sqfreetools.py
+++ b/sympy/polys/tests/test_sqfreetools.py
@@ -134,6 +134,10 @@ def test_dmp_sqf():
     f = -x**2 + 2*x - 1
     assert R.dmp_sqf_list_include(f) == [(-1, 1), (x - 1, 2)]
 
+    f = (y**2 + 1)**2*(x**2 + 2*x + 2)
+    assert R.dmp_sqf_p(f) is False
+    assert R.dmp_sqf_list(f) == (1, [(x**2 + 2*x + 2, 1), (y**2 + 1, 2)])
+
     R, x, y = ring("x,y", FF(2))
     raises(NotImplementedError, lambda: R.dmp_sqf_list(y**2 + 1))
 


Code review comments:
diff_hunk: @@ -118,59 +184,203 @@ def dup_sqf_norm(f, K):
     return s, f, r
 
 
+def _dmp_sqf_norm_shifts(f, u, K):
+    """Generate a sequence of candidate shifts for dmp_sqf_norm."""
+    #
+    # We want to find a minimal shift if possible because shifting high degree
+    # variables can be expensive e.g. x**10 -> (x + 1)**10. We try a few easy
+    # cases first before the final infinite loop that is guaranteed to give
+    # only finitely many bad shifts (see Trager76 for proof of this in the
+    # univariate case).
+    #
+
+    # First the trivial shift [0, 0, ...]
+    n = u + 1
+    s0 = [0] * n
+    yield s0, f
+
+    # Shift in multiples of the generator of the extension field K
+    a = K.unit
+
+    # Variables of degree > 0 ordered by increasing degree
+    d = dmp_degree_list(f, u)
+    var_indices = [i for di, i in sorted(zip(d, range(u+1))) if di > 0]
+
+    # Now try [1, 0, 0, ...], [0, 1, 0, ...]
+    for i in var_indices:
+        s1 = s0.copy()
+        s1[i] = 1
+        a1 = [-a*s1i for s1i in s1]
+        f1 = dmp_shift(f, a1, u, K)
+        yield s1, f1
+
+    # Now try [1, 1, 1, ...], [2, 2, 2, ...]
+    j = 0
+    while True:
+        j += 1
+        sj = [j] * n
+        aj = [-a*j] * n
+        fj = dmp_shift(f, aj, u, K)
+        yield sj, fj
+
+
 def dmp_sqf_norm(f, u, K):
-    """
-    Square-free norm of ``f`` in ``K[X]``, useful over algebraic domains.
+    r"""
+    Find a shift of ``f`` in ``K[X]`` that has square-free norm.
+
+    The domain `K` must be an algebraic number field `k(a)` (see :ref:`QQ(a)`).
 
-    Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and ``r(x) = Norm(g(x))``
-    is a square-free polynomial over K, where ``a`` is the algebraic extension of ``K``.
+    Returns `(s,g,r)`, such that `g(x_1,x_2,\cdots)=f(x_1-s_1 a, x_2 - s_2 a,
+    \cdots)`, `r(x)=\text{Norm}(g(x))` and `r` is a square-free polynomial over
+    `k`.
 
     Examples
     ========
 
+    We first create the algebraic number field `K=k(a)=\mathbb{Q}(i)` and rings
+    `K[x]` and `k[x]`:
path: sympy/polys/sqfreetools.py
text: @user1:
``K[x, y]`` and ``k[x, y]``

@author:
Fixed.

diff_hunk: @@ -118,59 +184,203 @@ def dup_sqf_norm(f, K):
     return s, f, r
 
 
+def _dmp_sqf_norm_shifts(f, u, K):
+    """Generate a sequence of candidate shifts for dmp_sqf_norm."""
+    #
+    # We want to find a minimal shift if possible because shifting high degree
+    # variables can be expensive e.g. x**10 -> (x + 1)**10. We try a few easy
+    # cases first before the final infinite loop that is guaranteed to give
+    # only finitely many bad shifts (see Trager76 for proof of this in the
+    # univariate case).
+    #
+
+    # First the trivial shift [0, 0, ...]
+    n = u + 1
+    s0 = [0] * n
+    yield s0, f
+
+    # Shift in multiples of the generator of the extension field K
+    a = K.unit
+
+    # Variables of degree > 0 ordered by increasing degree
+    d = dmp_degree_list(f, u)
+    var_indices = [i for di, i in sorted(zip(d, range(u+1))) if di > 0]
+
+    # Now try [1, 0, 0, ...], [0, 1, 0, ...]
+    for i in var_indices:
+        s1 = s0.copy()
+        s1[i] = 1
+        a1 = [-a*s1i for s1i in s1]
+        f1 = dmp_shift(f, a1, u, K)
+        yield s1, f1
+
+    # Now try [1, 1, 1, ...], [2, 2, 2, ...]
+    j = 0
+    while True:
+        j += 1
+        sj = [j] * n
+        aj = [-a*j] * n
+        fj = dmp_shift(f, aj, u, K)
+        yield sj, fj
+
+
 def dmp_sqf_norm(f, u, K):
-    """
-    Square-free norm of ``f`` in ``K[X]``, useful over algebraic domains.
+    r"""
+    Find a shift of ``f`` in ``K[X]`` that has square-free norm.
+
+    The domain `K` must be an algebraic number field `k(a)` (see :ref:`QQ(a)`).
 
-    Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and ``r(x) = Norm(g(x))``
-    is a square-free polynomial over K, where ``a`` is the algebraic extension of ``K``.
+    Returns `(s,g,r)`, such that `g(x_1,x_2,\cdots)=f(x_1-s_1 a, x_2 - s_2 a,
+    \cdots)`, `r(x)=\text{Norm}(g(x))` and `r` is a square-free polynomial over
+    `k`.
 
     Examples
     ========
 
+    We first create the algebraic number field `K=k(a)=\mathbb{Q}(i)` and rings
+    `K[x]` and `k[x]`:
+
     >>> from sympy.polys import ring, QQ
     >>> from sympy import I
 
     >>> 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)
+    We can now find a square free norm for a shift of `f`:
 
-    >>> s == 1
+    >>> f = x*y + y**2
+    >>> s, g, r = R.dmp_sqf_norm(f)
+
+    The choice of shifts ``s`` is arbitrary and the particular values returned
+    for ``g`` and ``r`` are determined by ``s``.
+
+    >>> s
+    [0, 1]
+    >>> g == x*y - I*x + y**2 - 2*I*y - 1
+    True
+    >>> r == X**2*Y**2 + X**2 + 2*X*Y**3 + 2*X*Y + Y**4 + 2*Y**2 + 1
     True
-    >>> f == x*y + y**2 + K([QQ(-1), QQ(0)])*y
+
+    The required invariants are:
+
+    >>> g == f.shift_list([-si*K.unit for si in s])
     True
-    >>> r == X**2*Y**2 + 2*X*Y**3 + Y**4 + Y**2
+    >>> g.norm() == r
     True
+    >>> r.is_squarefree
+    True
+
+    Explanation
+    ===========
+
+    This is part of Trager's algorithm for factorizing polynomials over
+    algebraic number fields. In particular this function is a multivariate
+    generalization of algorithm ``sqfr_norm`` from [Trager76]_.
+
+    See Also
+    ========
 
+    dup_sqf_norm:
+        Analogous function for multivariate polynomials over ``k(a)``.
path: sympy/polys/sqfreetools.py
text: @user1:
"multivariate" --> "univariate"

@author:
Fixed.