Repository: sympy/sympy

Issue:
SymPy not computing eigenvalues for 5 x 5 matrix
Reported from https://stackoverflow.com/questions/65170764/is-it-possible-somehow-to-find-complex-eigenvalues-using-sympy

```
>>> A = Matrix([[5, -5, -3, 2, -7], [-2, -5, 0, 2, 1], [-2, -7, -5, -2, -6], [7, 10, 3, 9, -2], [4, -10, 3, -8, -4]])
>>> A.eigenvals()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "./sympy/matrices/matrices.py", line 372, in eigenvals
    return _eigenvals(self, error_when_incomplete=error_when_incomplete, **flags)
  File "./sympy/matrices/eigen.py", line 177, in _eigenvals
    **flags)
  File "./sympy/matrices/eigen.py", line 227, in _eigenvals_dict
    .format(block))
sympy.matrices.common.MatrixError: Could not compute eigenvalues for Matrix([[5, -5, 2, -3, -7], [-2, -5, 2, 0, 1], [7, 10, 9, 3, -2], [-2, -7, -2, -5, -6], [4, -10, -8, 3, -4]]). if you see this error, please report to SymPy issue tracker.
```

I believe the issue is the same as https://github.com/sympy/sympy/issues/20275. The characteristic polynomial is degree 5, so can't be solved explicitly. However, SymPy should just return a list of RootOf classes in this case, rather than an error (especially if the error is going to suggest that this is necessarily a bug). 

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v161). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* matrices
  - `eigenvals, eigenvects` without radical solution will be returned as `CRootOf` ([#20614](https://github.com/sympy/sympy/pull/20614) by [@sylee957](https://github.com/sylee957))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.8.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- 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 #20582

    #### Brief description of what is fixed or changed

    This uses the suggested quotient ring computation for eigenvector computation which can be applied for rational coefficients.
    I'd briefly do the integration test here, and if it succeeds, some previous issues can be identified and resolved as well.

    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below. 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 -->
    - matrices
      - `eigenvals, eigenvects` without radical solution will be returned as `CRootOf`
    <!-- END RELEASE NOTES -->

</details><p>

sympy-bot: ### 🟠

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v161). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at https://github.com/sympy/sympy-bot/issues/75.

The following commits **add new files**:
* e84493194552885a28be1101f054d3a08e22e32b:
  - `sympy/polys/eigentools.py`

If these files were added/deleted on purpose, you can ignore this message.

oscarbenjamin: We should think about where we want to put this kind of code. There is a tension between wanting it in matrices or in polys. I think this belongs in polys but higher level code built on top can be in mtraices. At least if we are going to have multiple files in polys for matrices then I guess there should be a `sympy.polys.matrices` package for them to go in though.
oscarbenjamin: Perhaps it's worth me explaining what I was thinking about in the design of the `ddm_*` functions and the `DDM` and `DomainMatrix` classes.

Firstly the `ddm_*` functions are designed to operate on a generic list-of-lists representation. It shouldn't be necessary to use any particular class with them. So if someone wants to make use of rref in another context they don't have to use the `DDM` or `DomainMatrix` classes e.g.:
```python
In [1]: from sympy.polys.domainmatrix import ddm_irref

In [2]: M = [[1, 2, 3], [4, 5, 6]]

In [3]: ddm_irref(M)
Out[3]: [0, 1]

In [4]: M
Out[4]: [[1.0, 0.0, -1.0], [0, 1.0, 2.0]]
```
I thought it was important to keep the core algorithms untied to the classes themselves so that they could be reused without needing to subclass `DDM` or `DomainMatrix`. This was based on seeing how `Matrix` had been abused by subclasses that put incompatible elements in to the matrix breaking the key internal invariants of the class. Had it been easy to reuse the core routines that wouldn't have been necessary. These functions work in place because they want to allow their callers to work with any list-of-lists style data structure: the caller has to provide the data structures.

Ideally the `ddm_*` functions don't need to know anything about domains and they don't do any checking of the types of the objects to make sure that they are as fast as possible. There are some cases where a domain argument is passed but I'd like to change that so that whatever is needed e.g. `K.zero` or `K.gcd` is passed but the functions should still be usable with objects that don't have a complete `Domain` object.

The `DDM` class is a list subclass that makes the list of list representation more palatable by adding basic operations like `+`, `*` etc. It also handles the boiler-plate involved in copying etc to make the inputs needed for each `ddm_*` function and ensuring that the operations don't need to work in place so e.g. `DDM.rref` returns a new matrix. The `DDM` class has an associated `domain` which it uses for checking that operations are compatible but it otherwise does not try to do any conversion. It expects that it's inputs are of the expected domain and it refuses mixed operations. No attempt is made to perform potentially slow conversions between domains. In this way `DDM` instances of a given domain represent a ring-like structure. Subclassing list and using a list of list representation makes operations like indexing with `M[i][j]` fast as well as iteration etc. We lose that benefit if we override `__getitem__` etc so we don't try to make e.g. slicing work with `DDM`. We could give it an `extract(rows, cols)` method though that can efficiently select rows and columns given as lists of indices.

The `DomainMatrix` class wraps `DDM` into something more palatable again that can handle conversions to and from `Expr` and unifying domains in mixed operations. The idea is that it holds a `rep` and has routines for converting that between different forms in order to make the operations work. At this level we don't want to have core routines indexing directly into `DomainMatrix` like `M[i][j]` so performance of indexing is not so important: we can add `__getitem__` to `DomainMatrix` and we can give it basic slicing (based on `DDM.extract`) but we should still try to keep this limited for speed.

Later it would be good to have `SDM` which represents a sparse matrix as a dict of dicts. The idea would then be to have `sdm_*` functions that work on dict of dict and `SDM` to wrap that up providing operations like `+`, `-` and `extract` for efficient slicing. A `DomainMatrix.rep` might be `DDM` or `SDM`. Operations on `DomainMatrix` can `unify` between dense and sparse representations if needed. We still have only one `DomainMatrix` class but it could have different underlying representations.
sylee957: > Ideally the `ddm_*` functions don't need to know anything about domains and they don't do any checking of the types of the objects to make sure that they are as fast as possible. There are some cases where a domain argument is passed but I'd like to change that so that whatever is needed e.g. `K.zero` or `K.gcd` is passed but the functions should still be usable with objects that don't have a complete `Domain` object.

I have thought in the other way that: 
When constants are needed, passing the context and operating with `K.zero, K.one` should be more robust with less conversions needed.
And when there are cases like `K.gcd` is needed, we would never achieve that consistency.
So I always thought in converse way that every DDM functions should take domain context as its arguments, but if you had different idea, the cases like `K.gcd` should be avoided anyway too.

Alternatively, we may make some necessary context function as its arguments, like `ddm_foo(M, K) -> ddm_foo(M, gcd)` or even like `ddm_foo(M, K) -> ddm_foo(M, zero, one, gcd)`,

The other problem I find is we may not always have overridden the arithmetics for domains. And working purely with how operators overrides should be fragile for that.
Of course we may fill in the missing places when `__bool__` is available or not, but comparisons like `>, <` may not be added for domains with complex numbers, but when they be used by matrix pivoting rules, we need to assume `>, <` exists anyway.
oscarbenjamin: > Alternatively, we may make some necessary context function as its arguments, like `ddm_foo(M, K) -> ddm_foo(M, gcd)` or even like `ddm_foo(M, K) -> ddm_foo(M, zero, one, gcd)`,

I think it makes more sense just to have separate functions for the different algorithms like:
```
ddm_foo_gcd(M, gcd)
ddm_foo_bareiss(M)
ddm_foo_field(M)
```
Probably the `gcd` is not needed if using Bareiss anyway.

> but when they be used by matrix pivoting rules, we need to assume `>, <` exists anyway.

In #20483 I added pivoting for one case like:
```
ddm_rref(M, pivot_score)
```
Where `pivot_score` can rank pivots. I'm not sure if that's best but anyway the most important thing is that we don't need pivoting to select "small" elements in the normal sense of small. Rather we should prefer a pivot like `10` over a pivot like `(1 + x + x**2)/(y-2)`. It does really need to be a domain specific scoring system but if the caller can pass in their own score function to the core routine then they don't need to use a Domain necessarily.
sylee957: For eigenvector computation, I think that it's not possible to break the tie with domain because I always need to make some algebraic extension and if I drop this information, it is impossible to work with.
I also think that you have added `ddm_irref(A, QQ)` in #20483 and this muddles the original idea that ddm methods should be agnostic about domains, but maybe this can be limited only for a few key places.
I could also have suggested what about manipulating the DDM directly in DomainMatrix methods, but if we need some abstraction layer to work with sparse classes, this can be dropped.
oscarbenjamin: > For eigenvector computation, I think that it's not possible to break the tie with domain because I always need to make some algebraic extension and if I drop this information, it is impossible to work with.

That's why I think it would be better to implement the eigenvalue calculation in terms of `DomainMatrix` rather than at the `ddm_*` level (I'm not sure why you're avoiding that).

Computing eigenvectors is an operation that can be expressed in high level terms: `nullspace(A - lambda*I)`. At the `ddm_*` level we should have a minimal set of computational kernels so I think `ddm_nullspace` makes sense but not `ddm_eigenvalues`. Keeping this set minimal is important because it needs to be duplicated for `sdm_*`. It's probably also worth considering that these `ddm_*` functions are the low-level primitives that would be reimplemented in something like symengine to make this faster.

As another example I now think that `ddm_charpoly` should probably be implemented at the `DomainMatrix` level so that it can be reused for dense and sparse. It is an algorithm that is naturally expressible in high-level operations but it has instead been implemented in very low-level operations. We should have a `dm_charpoly` function that takes a `DomainMatrix` and produces the list of coefficients. Possibly also `DomainMatrix.charpoly` should convert that to a `Poly`.

We need to be careful not repeat the mistake in polys that high-level algorithms are tied to low-level data structures. That's why it's so difficult to introduce a sparse representation for `Poly`: most of the complicated algorithms are tied to the dense representation. Hence e.g. `dup_factor_list` which I would consider analogous to `ddm_eigenvalues` in this respect. Of course we do need core routines for dense manipulation but there should be a minimal set of optimised primitives.

> I also think that you have added `ddm_irref(A, QQ)` in #20483

That PR is in an experimental state. I just pushed up the last two experimental commits so that you could see them. I think that the scorer function should be selected by the caller of `ddm_irref` and provided instead of the domain.
sylee957: What is the advantage that DomainMatrix should use Poly instead of list, or any other lower level polynomial object for chapoly?
It may be only one sympifiale object for DomainMatrix.
One of the reason that I can't easily break the tie of eigenvects and DDM is because I have to tie charpoly with `dup_factor` as charpoly barely returns a list.
oscarbenjamin: > What is the advantage that DomainMatrix should use Poly instead of list, or any other lower level polynomial object for chapoly?

I'm not sure. It might be better to have different methods like `charpoly_list`. These details can be worked out as more of the implementation is fleshed out.

> One of the reason that I can't easily break the tie of eigenvects and DDM is because I have to tie charpoly with `dup_factor` in the current representation.

It's not hard to convert the list to a Poly and vice versa e.g.:
```python
In [7]: from sympy.polys.domainmatrix import DomainMatrix

In [8]: dM = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)

In [9]: p = dM.charpoly()

In [10]: p
Out[10]: [mpz(1), mpz(-5), mpz(-2)]

In [11]: Poly(p, x, domain=ZZ)
Out[11]: Poly(x**2 - 5*x - 2, x, domain='ZZ')

In [12]: Poly(p, x, domain=ZZ).factor_list()
Out[12]: (1, [(Poly(x**2 - 5*x - 2, x, domain='ZZ'), 1)])

In [18]: Poly(p, x, domain=ZZ).factor_list()[1][0][0].rep.rep
Out[18]: [mpz(1), mpz(-5), mpz(-2)]
```
You can also just use `dup_factor_list` at the `DomainMatrix` level so I'm not sure what you mean that this is tied to `DDM`.
sylee957: I'd note that other major reason that I had tied to DDM was that DomainMatrix.rref had only accepted domains with is_field, but I think that I can drop this requirement arbitrarily and see if it works.
oscarbenjamin: > I'd note that other major reason that I had tied to DDM was that DomainMatrix.rref had only accepted domains with is_field, but I think that I can drop this requirement arbitrarily and see if it works.

Yes, for now I think it would be fine for `dm_eigenvals` to just do `DomainMatrix.from_ddm(ddm_irref(dM.rep.copy()))`.

Patch:
diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py
index eb1fca596c..04efdda461 100644
--- a/sympy/matrices/eigen.py
+++ b/sympy/matrices/eigen.py
@@ -10,7 +10,9 @@
 from sympy.core.numbers import Float
 from sympy.core.sympify import _sympify
 from sympy.functions.elementary.miscellaneous import sqrt
-from sympy.polys import roots
+from sympy.polys import roots, CRootOf, EX
+from sympy.polys.domainmatrix import (
+    DomainMatrix, dom_eigenvects, dom_eigenvects_to_sympy)
 from sympy.simplify import nsimplify, simplify as _simplify
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 
@@ -80,8 +82,7 @@ def _eigenvects_mpmath(M):
 def _eigenvals(
     M, error_when_incomplete=True, *, simplify=False, multiple=False,
     rational=False, **flags):
-    r"""Return eigenvalues using the Berkowitz algorithm to compute
-    the characteristic polynomial.
+    r"""Compute eigenvalues of the matrix.
 
     Parameters
     ==========
@@ -149,6 +150,14 @@ def _eigenvals(
 
     Eigenvalues of a matrix $A$ can be computed by solving a matrix
     equation $\det(A - \lambda I) = 0$
+
+    It's not always possible to return radical solutions for
+    eigenvalues for matrices larger than $4, 4$ shape due to
+    Abel-Ruffini theorem.
+
+    If there is no radical solution is found for the eigenvalue,
+    it may return eigenvalues in the form of
+    :class:`sympy.polys.rootoftools.ComplexRootOf`.
     """
     if not M:
         if multiple:
@@ -188,14 +197,20 @@ def _eigenvals_list(
             charpoly = block.charpoly(simplify=simplify)
         else:
             charpoly = block.charpoly()
+
         eigs = roots(charpoly, multiple=True, **flags)
 
-        if error_when_incomplete:
-            if len(eigs) != block.rows:
-                raise MatrixError(
-                    "Could not compute eigenvalues for {}. if you see this "
-                    "error, please report to SymPy issue tracker."
-                    .format(block))
+        if len(eigs) != block.rows:
+            degree = int(charpoly.degree())
+            f = charpoly.as_expr()
+            x = charpoly.gen
+            try:
+                eigs = [CRootOf(f, x, idx) for idx in range(degree)]
+            except NotImplementedError:
+                if error_when_incomplete:
+                    raise MatrixError
+                else:
+                    eigs = []
 
         all_eigs += eigs
 
@@ -217,14 +232,20 @@ def _eigenvals_dict(
             charpoly = block.charpoly(simplify=simplify)
         else:
             charpoly = block.charpoly()
+
         eigs = roots(charpoly, multiple=False, **flags)
 
-        if error_when_incomplete:
-            if sum(eigs.values()) != block.rows:
-                raise MatrixError(
-                    "Could not compute eigenvalues for {}. if you see this "
-                    "error, please report to SymPy issue tracker."
-                    .format(block))
+        if sum(eigs.values()) != block.rows:
+            degree = int(charpoly.degree())
+            f = charpoly.as_expr()
+            x = charpoly.gen
+            try:
+                eigs = {CRootOf(f, x, idx): 1 for idx in range(degree)}
+            except NotImplementedError:
+                if error_when_incomplete:
+                    raise MatrixError
+                else:
+                    eigs = {}
 
         for k, v in eigs.items():
             if k in all_eigs:
@@ -254,9 +275,39 @@ def _eigenspace(M, eigenval, iszerofunc=_iszero, simplify=False):
     return ret
 
 
+def _eigenvects_DOM(M, **kwargs):
+    DOM = DomainMatrix.from_Matrix(M, field=True, extension=True)
+    if DOM.domain != EX:
+        rational, algebraic = dom_eigenvects(DOM)
+        eigenvects = dom_eigenvects_to_sympy(
+            rational, algebraic, M.__class__, **kwargs)
+        eigenvects = sorted(eigenvects, key=lambda x: default_sort_key(x[0]))
+
+        return eigenvects
+    return None
+
+
+def _eigenvects_sympy(M, iszerofunc, simplify=True, **flags):
+    eigenvals = M.eigenvals(rational=False, **flags)
+
+    # Make sure that we have all roots in radical form
+    for x in eigenvals:
+        if x.has(CRootOf):
+            raise MatrixError(
+                "Eigenvector computation is not implemented if the matrix have "
+                "eigenvalues in CRootOf form")
+
+    eigenvals = sorted(eigenvals.items(), key=default_sort_key)
+    ret = []
+    for val, mult in eigenvals:
+        vects = _eigenspace(M, val, iszerofunc=iszerofunc, simplify=simplify)
+        ret.append((val, mult, vects))
+    return ret
+
+
 # This functions is a candidate for caching if it gets implemented for matrices.
 def _eigenvects(M, error_when_incomplete=True, iszerofunc=_iszero, *, chop=False, **flags):
-    """Return list of triples (eigenval, multiplicity, eigenspace).
+    """Compute eigenvectors of the matrix.
 
     Parameters
     ==========
@@ -291,6 +342,7 @@ def _eigenvects(M, error_when_incomplete=True, iszerofunc=_iszero, *, chop=False
 
     Returns
     =======
+
     ret : [(eigenval, multiplicity, eigenspace), ...]
         A ragged list containing tuples of data obtained by ``eigenvals``
         and ``nullspace``.
@@ -332,6 +384,7 @@ def _eigenvects(M, error_when_incomplete=True, iszerofunc=_iszero, *, chop=False
     """
     simplify = flags.get('simplify', True)
     primitive = flags.get('simplify', False)
+    flags.pop('simplify', None)  # remove this if it's there
     flags.pop('multiple', None)  # remove this if it's there
 
     if not isinstance(simplify, FunctionType):
@@ -343,15 +396,9 @@ def _eigenvects(M, error_when_incomplete=True, iszerofunc=_iszero, *, chop=False
             return _eigenvects_mpmath(M)
         M = M.applyfunc(lambda x: nsimplify(x, rational=True))
 
-    eigenvals = M.eigenvals(
-        rational=False, error_when_incomplete=error_when_incomplete,
-        **flags)
-
-    eigenvals = sorted(eigenvals.items(), key=default_sort_key)
-    ret = []
-    for val, mult in eigenvals:
-        vects = _eigenspace(M, val, iszerofunc=iszerofunc, simplify=simplify)
-        ret.append((val, mult, vects))
+    ret = _eigenvects_DOM(M)
+    if ret is None:
+        ret = _eigenvects_sympy(M, iszerofunc, simplify=simplify, **flags)
 
     if primitive:
         # if the primitive flag is set, get rid of any common
@@ -1120,10 +1167,12 @@ def pick_vec(small_basis, big_basis):
     # first calculate the jordan block structure
     eigs = mat.eigenvals()
 
-    # make sure that we found all the roots by counting
-    # the algebraic multiplicity
-    if sum(m for m in eigs.values()) != mat.cols:
-        raise MatrixError("Could not compute eigenvalues for {}".format(mat))
+    # Make sure that we have all roots in radical form
+    for x in eigs:
+        if x.has(CRootOf):
+            raise MatrixError(
+                "Jordan normal form is not implemented if the matrix have "
+                "eigenvalues in CRootOf form")
 
     # most matrices have distinct eigenvalues
     # and so are diagonalizable.  In this case, don't
diff --git a/sympy/polys/agca/extensions.py b/sympy/polys/agca/extensions.py
index f0be073f1b..0864da18f7 100644
--- a/sympy/polys/agca/extensions.py
+++ b/sympy/polys/agca/extensions.py
@@ -22,6 +22,9 @@ def __init__(self, rep, ext):
         self.rep = rep
         self.ext = ext
 
+    def __bool__(f):
+        return bool(f.rep)
+
     def __neg__(f):
         return ExtElem(-f.rep, f.ext)
 
@@ -229,6 +232,9 @@ def convert(self, f):
         rep = self.ring.convert(f)
         return ExtElem(rep % self.mod, self)
 
+    def to_sympy(self, f):
+        return self.ring.to_sympy(f.rep)
+
     __call__ = convert
 
     def __str__(self):
diff --git a/sympy/polys/domainmatrix.py b/sympy/polys/domainmatrix.py
index d5c3fac27f..dd97d0ba8b 100644
--- a/sympy/polys/domainmatrix.py
+++ b/sympy/polys/domainmatrix.py
@@ -1,9 +1,16 @@
 from operator import mul
 
+from sympy.core.symbol import Dummy
 from sympy.core.sympify import _sympify
+
 from sympy.matrices.common import (NonInvertibleMatrixError,
     NonSquareMatrixError, ShapeError)
+from sympy.polys import Poly
+from sympy.polys.agca.extensions import FiniteExtension
 from sympy.polys.constructor import construct_domain
+from sympy.polys.factortools import dup_factor_list
+from sympy.polys.polyroots import roots
+from sympy.polys.rootoftools import CRootOf
 
 
 class DDMError(Exception):
@@ -146,6 +153,22 @@ def rref(a):
         pivots = ddm_irref(b)
         return b, pivots
 
+    def nullspace(a):
+        rref, pivots = a.rref()
+        rows, cols = a.shape
+        domain = a.domain
+
+        basis = []
+        for i in range(cols):
+            if i in pivots:
+                continue
+            vec = [domain.one if i == j else domain.zero for j in range(cols)]
+            for ii, jj in enumerate(pivots):
+                vec[jj] -= rref[ii][i]
+            basis.append(vec)
+
+        return DDM(basis, (len(basis), cols), domain)
+
     def det(a):
         """Determinant of a"""
         m, n = a.shape
@@ -639,6 +662,9 @@ def rref(self):
         rref_ddm, pivots = self.rep.rref()
         return self.from_ddm(rref_ddm), tuple(pivots)
 
+    def nullspace(self):
+        return self.from_ddm(self.rep.nullspace())
+
     def inv(self):
         if not self.domain.is_Field:
             raise ValueError('Not a field')
@@ -683,3 +709,79 @@ def __eq__(A, B):
         if not isinstance(B, DomainMatrix):
             return NotImplemented
         return A.rep == B.rep
+
+
+def dom_eigenvects(A, l=Dummy('lambda')):
+    charpoly = A.charpoly()
+    rows, cols = A.shape
+    domain = A.domain
+    _, factors = dup_factor_list(charpoly, domain)
+
+    rational_eigenvects = []
+    algebraic_eigenvects = []
+    for base, exp in factors:
+        if len(base) == 2:
+            field = domain
+            eigenval = -base[1] / base[0]
+
+            EE_items = [
+                [eigenval if i == j else field.zero for j in range(cols)]
+                for i in range(rows)]
+            EE = DomainMatrix(EE_items, (rows, cols), field)
+
+            basis = (A - EE).nullspace()
+            rational_eigenvects.append((field, eigenval, exp, basis))
+        else:
+            minpoly = Poly.from_list(base, l, domain=domain)
+            field = FiniteExtension(minpoly)
+            eigenval = field(l)
+
+            AA_items = [
+                [Poly.from_list([item], l, domain=domain).rep for item in row]
+                for row in A.rep]
+            AA_items = [[field(item) for item in row] for row in AA_items]
+            AA = DomainMatrix(AA_items, (rows, cols), field)
+            EE_items = [
+                [eigenval if i == j else field.zero for j in range(cols)]
+                for i in range(rows)]
+            EE = DomainMatrix(EE_items, (rows, cols), field)
+
+            basis = (AA - EE).nullspace()
+            algebraic_eigenvects.append((field, minpoly, exp, basis))
+
+    return rational_eigenvects, algebraic_eigenvects
+
+
+def dom_eigenvects_to_sympy(
+    rational_eigenvects, algebraic_eigenvects,
+    Matrix, **kwargs
+):
+    result = []
+
+    for field, eigenvalue, multiplicity, eigenvects in rational_eigenvects:
+        eigenvects = eigenvects.rep
+        eigenvalue = field.to_sympy(eigenvalue)
+        new_eigenvects = [
+            Matrix([field.to_sympy(x) for x in vect])
+            for vect in eigenvects]
+        result.append((eigenvalue, multiplicity, new_eigenvects))
+
+    for field, minpoly, multiplicity, eigenvects in algebraic_eigenvects:
+        eigenvects = eigenvects.rep
+        l = minpoly.gens[0]
+
+        eigenvects = [[field.to_sympy(x) for x in vect] for vect in eigenvects]
+
+        degree = minpoly.degree()
+        minpoly = minpoly.as_expr()
+        eigenvals = roots(minpoly, l, **kwargs)
+        if len(eigenvals) != degree:
+            eigenvals = [CRootOf(minpoly, l, idx) for idx in range(degree)]
+
+        for eigenvalue in eigenvals:
+            new_eigenvects = [
+                Matrix([x.subs(l, eigenvalue) for x in vect])
+                for vect in eigenvects]
+            result.append((eigenvalue, multiplicity, new_eigenvects))
+
+    return result
diff --git a/sympy/polys/domains/quotientring.py b/sympy/polys/domains/quotientring.py
index cae78c3567..5a257f8966 100644
--- a/sympy/polys/domains/quotientring.py
+++ b/sympy/polys/domains/quotientring.py
@@ -30,6 +30,11 @@ def __str__(self):
         from sympy import sstr
         return sstr(self.data) + " + " + str(self.ring.base_ideal)
 
+    __repr__ = __str__
+
+    def __bool__(self):
+        return not self.ring.is_zero(self)
+
     def __add__(self, om):
         if not isinstance(om, self.__class__) or om.ring != self.ring:
             try:
@@ -71,7 +76,9 @@ def __truediv__(self, o):
         return self.ring.revert(o)*self
 
     def __pow__(self, oth):
-        return self.ring(self.data**oth)
+        if oth < 0:
+            return self.ring.revert(self) ** -oth
+        return self.ring(self.data ** oth)
 
     def __eq__(self, om):
         if not isinstance(om, self.__class__) or om.ring != self.ring:


Test patch:
diff --git a/sympy/matrices/tests/test_eigen.py b/sympy/matrices/tests/test_eigen.py
index 9e311cd422..b2a89e416c 100644
--- a/sympy/matrices/tests/test_eigen.py
+++ b/sympy/matrices/tests/test_eigen.py
@@ -1,10 +1,11 @@
 from sympy import (
     Rational, Symbol, N, I, Abs, sqrt, exp, Float, sin,
     cos, symbols)
-from sympy.matrices import eye, Matrix, dotprodsimp
+from sympy.matrices import eye, Matrix
 from sympy.core.singleton import S
 from sympy.testing.pytest import raises, XFAIL
 from sympy.matrices.matrices import NonSquareMatrixError, MatrixError
+from sympy.matrices.expressions.fourier import DFT
 from sympy.simplify.simplify import simplify
 from sympy.matrices.immutable import ImmutableMatrix
 from sympy.testing.pytest import slow
@@ -103,14 +104,11 @@ def NS(e, n):
     assert max(i.q for i in M._eigenvects[0][2][0]) > 1
     M._eigenvects = M.eigenvects(simplify=True)
     assert max(i.q for i in M._eigenvects[0][2][0]) == 1
+
     M = Matrix([[Rational(1, 4), 1], [1, 1]])
-    assert M.eigenvects(simplify=True) == [
+    assert M.eigenvects() == [
         (Rational(5, 8) - sqrt(73)/8, 1, [Matrix([[-sqrt(73)/8 - Rational(3, 8)], [1]])]),
         (Rational(5, 8) + sqrt(73)/8, 1, [Matrix([[Rational(-3, 8) + sqrt(73)/8], [1]])])]
-    with dotprodsimp(True):
-        assert M.eigenvects(simplify=False) == [
-            (Rational(5, 8) - sqrt(73)/8, 1, [Matrix([[-1/(-Rational(3, 8) + sqrt(73)/8)], [1]])]),
-            (Rational(5, 8) + sqrt(73)/8, 1, [Matrix([[8/(3 + sqrt(73))], [1]])])]
 
     # issue 10719
     assert Matrix([]).eigenvals() == {}
@@ -214,14 +212,16 @@ def test_eigenvals():
                 [1, 1, 1]])
     assert M.eigenvals() == {2*S.One: 1, -S.One: 1, S.Zero: 1}
 
-    # if we cannot factor the char poly, we raise an error
     m = Matrix([
         [3,  0,  0, 0, -3],
         [0, -3, -3, 0,  3],
         [0,  3,  0, 3,  0],
         [0,  0,  3, 0,  3],
         [3,  0,  0, 3,  0]])
-    raises(MatrixError, lambda: m.eigenvals())
+
+    # XXX Used dry-run test because arbitrary symbol that appears in
+    # CRootOf may not be unique.
+    assert m.eigenvals()
 
 
 def test_eigenvects():
@@ -453,7 +453,7 @@ def test_jordan_form():
     assert Matrix(1, 1, [1]).jordan_form() == (Matrix([1]), Matrix([1]))
     assert Matrix(1, 1, [1]).jordan_form(calc_transform=False) == Matrix([1])
 
-    # make sure if we cannot factor the characteristic polynomial, we raise an error
+    # If we have eigenvalues in CRootOf form, raise errors
     m = Matrix([[3, 0, 0, 0, -3], [0, -3, -3, 0, 3], [0, 3, 0, 3, 0], [0, 0, 3, 0, 3], [3, 0, 0, 3, 0]])
     raises(MatrixError, lambda: m.jordan_form())
 
@@ -620,3 +620,64 @@ def test_positive_semidefinite_cholesky():
     assert _is_positive_semidefinite_cholesky(m) == True
     m = Matrix([[1, 2, -1], [2, 5, 1], [-1, 1, 9]])
     assert _is_positive_semidefinite_cholesky(m) == False
+
+
+def test_issue_20582():
+    A = Matrix([
+        [5, -5, -3, 2, -7],
+        [-2, -5, 0, 2, 1],
+        [-2, -7, -5, -2, -6],
+        [7, 10, 3, 9, -2],
+        [4, -10, 3, -8, -4]
+    ])
+    # XXX Used dry-run test because arbitrary symbol that appears in
+    # CRootOf may not be unique.
+    assert A.eigenvects()
+
+
+def test_issue_20275():
+    # XXX We use complex expansions because complex exponentials are not
+    # recognized by polys.domains
+    A = DFT(3).as_explicit().expand(complex=True)
+    eigenvects = A.eigenvects()
+    assert eigenvects[0] == (
+        -1, 1,
+        [Matrix([[1 - sqrt(3)], [1], [1]])]
+    )
+    assert eigenvects[1] == (
+        1, 1,
+        [Matrix([[1 + sqrt(3)], [1], [1]])]
+    )
+    assert eigenvects[2] == (
+        -I, 1,
+        [Matrix([[0], [-1], [1]])]
+    )
+
+    A = DFT(4).as_explicit().expand(complex=True)
+    eigenvects = A.eigenvects()
+    assert eigenvects[0] == (
+        -1, 1,
+        [Matrix([[-1], [1], [1], [1]])]
+    )
+    assert eigenvects[1] == (
+        1, 2,
+        [Matrix([[1], [0], [1], [0]]), Matrix([[2], [1], [0], [1]])]
+    )
+    assert eigenvects[2] == (
+        -I, 1,
+        [Matrix([[0], [-1], [0], [1]])]
+    )
+
+    # XXX We skip test for some parts of eigenvectors which are very
+    # complicated and fragile under expression tree changes
+    A = DFT(5).as_explicit().expand(complex=True)
+    eigenvects = A.eigenvects()
+    assert eigenvects[0] == (
+        -1, 1,
+        [Matrix([[1 - sqrt(5)], [1], [1], [1], [1]])]
+    )
+    assert eigenvects[1] == (
+        1, 2,
+        [Matrix([[S(1)/2 + sqrt(5)/2], [0], [1], [1], [0]]),
+         Matrix([[S(1)/2 + sqrt(5)/2], [1], [0], [0], [1]])]
+    )
diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py
index 9a39d9f33d..0953ee7e9d 100644
--- a/sympy/matrices/tests/test_matrices.py
+++ b/sympy/matrices/tests/test_matrices.py
@@ -404,22 +404,27 @@ def test_issue_17247_expression_blowup_13():
         [1 - x, x + 1,     0, x + 1],
         [    0, 1 - x, x + 1, 1 - x],
         [    0,     0,     1 - x, 0]])
-    with dotprodsimp(True):
-        ev = M.eigenvects()
-        assert ev[0][:2] == (0, 2)
-        assert ev[0][2][0] == Matrix([[0],[-1],[0],[1]])
-        assert ev[1][:2] == (x - sqrt(2)*(x - 1) + 1, 1)
-        assert (ev[1][2][0] - Matrix([
-            [-(-17*x**4 + 12*sqrt(2)*x**4 - 4*sqrt(2)*x**3 + 6*x**3 - 6*x - 4*sqrt(2)*x + 12*sqrt(2) + 17)/(-7*x**4 + 5*sqrt(2)*x**4 - 6*sqrt(2)*x**3 + 8*x**3 - 2*x**2 + 8*x + 6*sqrt(2)*x - 5*sqrt(2) - 7)],
-            [                      (-7*x**3 + 5*sqrt(2)*x**3 - x**2 + sqrt(2)*x**2 - sqrt(2)*x - x - 5*sqrt(2) - 7)/(-3*x**3 + 2*sqrt(2)*x**3 - 2*sqrt(2)*x**2 + 3*x**2 + 2*sqrt(2)*x + 3*x - 3 - 2*sqrt(2))],
-            [                                                                                           -(-3*x**2 + 2*sqrt(2)*x**2 + 2*x - 3 - 2*sqrt(2))/(-x**2 + sqrt(2)*x**2 - 2*sqrt(2)*x + 1 + sqrt(2))],
-            [                                                                                                                                                                                              1]])).expand() == Matrix([[0],[0],[0],[0]])
-        assert ev[2][:2] == (x + sqrt(2)*(x - 1) + 1, 1)
-        assert (ev[2][2][0] - Matrix([
-            [-(12*sqrt(2)*x**4 + 17*x**4 - 6*x**3 - 4*sqrt(2)*x**3 - 4*sqrt(2)*x + 6*x - 17 + 12*sqrt(2))/(7*x**4 + 5*sqrt(2)*x**4 - 6*sqrt(2)*x**3 - 8*x**3 + 2*x**2 - 8*x + 6*sqrt(2)*x - 5*sqrt(2) + 7)],
-            [                      (7*x**3 + 5*sqrt(2)*x**3 + x**2 + sqrt(2)*x**2 - sqrt(2)*x + x - 5*sqrt(2) + 7)/(2*sqrt(2)*x**3 + 3*x**3 - 3*x**2 - 2*sqrt(2)*x**2 - 3*x + 2*sqrt(2)*x - 2*sqrt(2) + 3)],
-            [                                                                                           -(2*sqrt(2)*x**2 + 3*x**2 - 2*x - 2*sqrt(2) + 3)/(x**2 + sqrt(2)*x**2 - 2*sqrt(2)*x - 1 + sqrt(2))],
-            [                                                                                                                                                                                            1]])).expand() == Matrix([[0],[0],[0],[0]])
+
+    ev = M.eigenvects()
+    assert ev[0] == (0, 2, [Matrix([0, -1, 0, 1])])
+    assert ev[1][0] == x - sqrt(2)*(x - 1) + 1
+    assert ev[1][1] == 1
+    assert ev[1][2][0].expand(deep=False, numer=True) == Matrix([
+        [(-x + sqrt(2)*(x - 1) - 1)/(x - 1)],
+        [-4*x/(x**2 - 2*x + 1) + (x + 1)*(x - sqrt(2)*(x - 1) + 1)/(x**2 - 2*x + 1)],
+        [(-x + sqrt(2)*(x - 1) - 1)/(x - 1)],
+        [1]
+    ])
+
+    assert ev[2][0] == x + sqrt(2)*(x - 1) + 1
+    assert ev[2][1] == 1
+    assert ev[2][2][0].expand(deep=False, numer=True) == Matrix([
+        [(-x - sqrt(2)*(x - 1) - 1)/(x - 1)],
+        [-4*x/(x**2 - 2*x + 1) + (x + 1)*(x + sqrt(2)*(x - 1) + 1)/(x**2 - 2*x + 1)],
+        [(-x - sqrt(2)*(x - 1) - 1)/(x - 1)],
+        [1]
+    ])
+
 
 def test_issue_17247_expression_blowup_14():
     M = Matrix(8, 8, ([1+x, 1-x]*4 + [1-x, 1+x]*4)*4)
@@ -2591,19 +2596,19 @@ def test_pinv():
 def test_pinv_rank_deficient_when_diagonalization_fails():
     # Test the four properties of the pseudoinverse for matrices when
     # diagonalization of A.H*A fails.
-    As = [Matrix([
-        [61, 89, 55, 20, 71, 0],
-        [62, 96, 85, 85, 16, 0],
-        [69, 56, 17,  4, 54, 0],
-        [10, 54, 91, 41, 71, 0],
-        [ 7, 30, 10, 48, 90, 0],
-        [0,0,0,0,0,0]])]
+    As = [
+        Matrix([
+            [61, 89, 55, 20, 71, 0],
+            [62, 96, 85, 85, 16, 0],
+            [69, 56, 17,  4, 54, 0],
+            [10, 54, 91, 41, 71, 0],
+            [ 7, 30, 10, 48, 90, 0],
+            [0, 0, 0, 0, 0, 0]])
+    ]
     for A in As:
         A_pinv = A.pinv(method="ED")
         AAp = A * A_pinv
         ApA = A_pinv * A
-        assert simplify(AAp * A) == A
-        assert simplify(ApA * A_pinv) == A_pinv
         assert AAp.H == AAp
         assert ApA.H == ApA