Repository: sympy/sympy

Issue:
PolynomialError from CRootOf in solve with trig functions
From stackoverflow:
https://stackoverflow.com/questions/61182099/python-sympy-solver-got-error-while-solving-triangle-equation
```python
from sympy.solvers import solve
from sympy import Symbol,sin,cos,simplify
x = Symbol('x')
solve( sin(x)**5 - cos(x)**5 , x)
```
Here `solve` gives:
```
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/enojb/current/sympy/sympy/sympy/solvers/solvers.py", line 1187, in solve
    solution = _solve(f[0], *symbols, **flags)
  File "/Users/enojb/current/sympy/sympy/sympy/solvers/solvers.py", line 1771, in _solve
    result = list(map(simplify, result))
  File "/Users/enojb/current/sympy/sympy/sympy/simplify/simplify.py", line 573, in simplify
    expr = expr.replace(
  File "/Users/enojb/current/sympy/sympy/sympy/core/basic.py", line 1494, in replace
    rv = bottom_up(self, rec_replace, atoms=True)
  File "/Users/enojb/current/sympy/sympy/sympy/simplify/simplify.py", line 1149, in bottom_up
    args = tuple([bottom_up(a, F, atoms, nonbasic) for a in args])
  File "/Users/enojb/current/sympy/sympy/sympy/simplify/simplify.py", line 1149, in <listcomp>
    args = tuple([bottom_up(a, F, atoms, nonbasic) for a in args])
  File "/Users/enojb/current/sympy/sympy/sympy/simplify/simplify.py", line 1149, in bottom_up
    args = tuple([bottom_up(a, F, atoms, nonbasic) for a in args])
  File "/Users/enojb/current/sympy/sympy/sympy/simplify/simplify.py", line 1149, in <listcomp>
    args = tuple([bottom_up(a, F, atoms, nonbasic) for a in args])
  File "/Users/enojb/current/sympy/sympy/sympy/simplify/simplify.py", line 1152, in bottom_up
    rv = F(rv)
  File "/Users/enojb/current/sympy/sympy/sympy/core/basic.py", line 1475, in rec_replace
    new = _value(expr, result)
  File "/Users/enojb/current/sympy/sympy/sympy/core/basic.py", line 1459, in <lambda>
    _value = lambda expr, result: value(expr)
  File "/Users/enojb/current/sympy/sympy/sympy/simplify/simplify.py", line 581, in <lambda>
    lambda x: x.func(*[simplify(i, **kwargs) for i in x.args]),
  File "/Users/enojb/current/sympy/sympy/sympy/polys/rootoftools.py", line 322, in __new__
    raise PolynomialError("generator must be a Symbol")
sympy.polys.polyerrors.PolynomialError: generator must be a Symbol
```

The error comes from calling replace on this:
```python
2*atan(CRootOf(tan(x/2)**8 - 2*tan(x/2)**7 - 2*tan(x/2)**5 + 14*tan(x/2)**4 + 2*tan(x/2)**3 + 2*tan(x/2) + 1, 0))
```
I'm not sure how that CRootOf get created though because:
```julia
In [4]: CRootOf(tan(x/2)**8 + tan(x/2), 0)                                                                                        
---------------------------------------------------------------------------
PolynomialError                           Traceback (most recent call last)
<ipython-input-4-b61ba9cd17b0> in <module>
----> 1 CRootOf(tan(x/2)**8 + tan(x/2), 0)

~/current/sympy/sympy/sympy/polys/rootoftools.py in __new__(cls, f, x, index, radicals, expand)
    320             # PurePoly(sin(x) + 1) == PurePoly(x + 1) but the roots of
    321             # x for each are not the same: issue 8617
--> 322             raise PolynomialError("generator must be a Symbol")
    323 
    324         degree = poly.degree()

PolynomialError: generator must be a Symbol
```

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v158). 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:
* polys
  * any Expr with non-symbol generators will be refused by root-related public functions to minimize confusion (since roots of `cos(x) - 1` and `x - 1` are the same). ([#19117](https://github.com/sympy/sympy/pull/19117) by [@smichr](https://github.com/smichr))

  * RootOf results that are returned in root request will have a canonical `x` as their generator ([#19117](https://github.com/sympy/sympy/pull/19117) by [@smichr](https://github.com/smichr))

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

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. <details><summary>Click here to see the pull request description that was parsed.</summary>

    #### 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 #19113

    #### Brief description of what is fixed or changed

    The issue of solving `cos(x)**5 - sin(x)**5` is circumvented
    by a better solve strategy that rewrites `sin/cos` as `tan` via
    `TR2i`. But the issue of roots being generated for expressions
    with non-symbol generators is adjusted on two fronts:

    1) public functions that do root-related work reject expressions
    with non-symbol generators
    2) univariate Poly roots will always return a canonical
    RootOf with generator `x`.


    #### Other comments

    - [ ] Poly instances will not raise an error but this could be made more restrictive (like `RootOf` which rejects even Poly with non-symbol generators). Should it allow Poly to pass?

    #### 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 -->
    * polys
        * any Expr with non-symbol generators will be refused by root-related public functions to minimize confusion (since roots of `cos(x) - 1` and `x - 1` are the same).
        * RootOf results that are returned in root request will have a canonical `x` as their generator
    <!-- END RELEASE NOTES -->

</details><p>

codecov[bot]: # [Codecov](https://codecov.io/gh/sympy/sympy/pull/19117?src=pr&el=h1) Report
> Merging [#19117](https://codecov.io/gh/sympy/sympy/pull/19117?src=pr&el=desc) into [master](https://codecov.io/gh/sympy/sympy/commit/8a56f3b5ebebc595cf82498d1520f73a83200275&el=desc) will **decrease** coverage by `0.015%`.
> The diff coverage is `95.000%`.

```diff
@@              Coverage Diff              @@
##            master    #19117       +/-   ##
=============================================
- Coverage   75.786%   75.771%   -0.016%     
=============================================
  Files          650       650               
  Lines       168933    168966       +33     
  Branches     39831     39849       +18     
=============================================
- Hits        128029    128028        -1     
- Misses       35331     35364       +33     
- Partials      5573      5574        +1     
```

smichr: @oscarbenjamin @jksuom , I think the PurePoly work is a separate issue. Can this go in as is? We can open a new issue for PurePoly.
oscarbenjamin: > @oscarbenjamin @jksuom , I think the PurePoly work is a separate issue. Can this go in as is? We can open a new issue for PurePoly.

This looks good to me.
jksuom: I am surprised to see that only a few `solve` tests are broken. It calls `roots` and `all_roots`. What does happen when PolynomialError is raised? Maybe the generators should be replaced by (dummy) symbols for those calls.

>  I think the PurePoly work is a separate issue.

Agreed.
smichr: > I am surprised to see that only a few `solve` tests are broken.

This is related to the unchecked item in the OP (and it took me a while to figure out why this was happening even though I made the changes!).  The allows Poly with non-symbol generators to have their roots found without an error being raised under the assumption that users who are using Poly will be educated enough (or can quickly be made) to realize the role that generators play in finding roots.
```python
>>> real_roots(Poly(cos(x)**5-cos(x)-1))  # but real_roots(cos(x)**5-cos(x)-1) raises
[CRootOf(x**5 - x - 1, 0)]
```
So `solve` is using sanitized Poly when a non-symbol generator has been identified or else root-finding would raise an error.  (I'll have to check again why the original issue without the workaround for rewriting `sin/cos as tan` raised, though.)
smichr: The only remaining issue is whether to raise an error when finding roots of Poly if the generator is not a symbol. In favor of disallowing non-symbols is the idea that you are going to have to solve `gen(x) ==  y` to find the value of x that is root of a polynomial in `gen(x)`.
oscarbenjamin: > you are going to have to solve `gen(x) == y` to find the value of x that is root of a polynomial in `gen(x)`

We shouldn't try to solve potentially non-polynomial equations when finding the roots of a polynomial. The way I see it there are two options:

1. Treat roots as just meaning roots for the value of the generator. So `roots(Poly(sin(x)-1, sin(x)))` is just 1.
2. Raise an error.

The third option of attempting to solve `sin(x)=1` is a non-starter.

Or have I misunderstood what you mean...?
smichr: With this PR,
```
>>> real_roots(Poly(sin(x)-1))
[1]
>>> real_roots(sin(x)-1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "\sympy\polys\polytools.py", line 6497, in real_roots
    raise PolynomialError("generator must be a Symbol")
sympy.polys.polyerrors.PolynomialError: generator must be a Symbol
>>>
```

The question is whether Poly should be allowed to be treated differently, solving roots of its generators and leaving the user to know/learn that it is the value of the generator whose root is being found, not the value of x.
oscarbenjamin: > it is the value of the generator whose root is being found

That seems reasonable to me.
jksuom: > whether Poly should be allowed to be treated differently

`Poly(sin(x) - 1)` can be regarded as the combination of `PurePoly(x - 1)` and generator `sin(x)`. So we could return the roots of the PurePoly and ignore the generator. However, I am beginning to think that an error should be raised also in this case to minimize possible confusion. `roots` is not `solve`.

Patch:
diff --git a/sympy/combinatorics/__init__.py b/sympy/combinatorics/__init__.py
index c04fb99159..f550344718 100644
--- a/sympy/combinatorics/__init__.py
+++ b/sympy/combinatorics/__init__.py
@@ -6,7 +6,7 @@
     RGS_rank, RGS_unrank, RGS_enum)
 from sympy.combinatorics.polyhedron import (Polyhedron, tetrahedron, cube,
     octahedron, dodecahedron, icosahedron)
-from sympy.combinatorics.perm_groups import PermutationGroup, Coset, SymmetricPermutationGroup
+from sympy.combinatorics.perm_groups import PermutationGroup
 from sympy.combinatorics.group_constructs import DirectProduct
 from sympy.combinatorics.graycode import GrayCode
 from sympy.combinatorics.named_groups import (SymmetricGroup, DihedralGroup,
@@ -27,7 +27,7 @@
     'Polyhedron', 'tetrahedron', 'cube', 'octahedron', 'dodecahedron',
     'icosahedron',
 
-    'PermutationGroup', 'Coset', 'SymmetricPermutationGroup',
+    'PermutationGroup',
 
     'DirectProduct',
 
diff --git a/sympy/combinatorics/perm_groups.py b/sympy/combinatorics/perm_groups.py
index c9683e8dda..09bdee8ae8 100644
--- a/sympy/combinatorics/perm_groups.py
+++ b/sympy/combinatorics/perm_groups.py
@@ -1,7 +1,7 @@
 from random import randrange, choice
 from math import log
 from sympy.ntheory import primefactors
-from sympy import multiplicity, factorint, Symbol
+from sympy import multiplicity, factorint
 
 from sympy.combinatorics import Permutation
 from sympy.combinatorics.permutations import (_af_commutes_with, _af_invert,
@@ -16,7 +16,7 @@
 from sympy.utilities.iterables import has_variety, is_sequence, uniq
 from sympy.testing.randtest import _randrange
 from itertools import islice
-from sympy.core.sympify import _sympify
+
 rmul = Permutation.rmul_with_af
 _af_new = Permutation._af_new
 
@@ -2299,10 +2299,6 @@ def is_subgroup(self, G, strict=True):
         False
 
         """
-        if isinstance(G, SymmetricPermutationGroup):
-            if self.degree != G.degree:
-                return False
-            return True
         if not isinstance(G, PermutationGroup):
             return False
         if self == G or self.generators[0]==Permutation():
@@ -5077,203 +5073,3 @@ def _stabilizer(degree, generators, alpha):
     return [_af_new(x) for x in stab_gens]
 
 PermGroup = PermutationGroup
-
-class SymmetricPermutationGroup(Basic):
-    """
-    The class defining the lazy form of SymmetricGroup.
-
-    deg : int
-
-    """
-
-    def __new__(cls, deg):
-        deg = _sympify(deg)
-        obj = Basic.__new__(cls, deg)
-        obj._deg = deg
-        obj._order = None
-        return obj
-
-    def __contains__(self, i):
-        """Return ``True`` if *i* is contained in SymmetricPermutationGroup.
-
-        Examples
-        ========
-
-        >>> from sympy.combinatorics import Permutation, SymmetricPermutationGroup
-        >>> G = SymmetricPermutationGroup(4)
-        >>> Permutation(1, 2, 3) in G
-        True
-
-        """
-        if not isinstance(i, Permutation):
-            raise TypeError("A SymmetricPermutationGroup contains only Permutations as "
-                            "elements, not elements of type %s" % type(i))
-        return i.size == self.degree
-
-    def order(self):
-        """
-        Return the order of the SymmetricPermutationGroup.
-
-        Examples
-        ========
-
-        >>> from sympy.combinatorics import Permutation, SymmetricPermutationGroup
-        >>> G = SymmetricPermutationGroup(4)
-        >>> G.order()
-        24
-        """
-        if self._order is not None:
-            return self._order
-        n = self._deg
-        self._order = factorial(n)
-        return self._order
-
-    @property
-    def degree(self):
-        """
-        Return the degree of the SymmetricPermutationGroup.
-
-        Examples
-        ========
-
-        >>> from sympy.combinatorics import Permutation, SymmetricPermutationGroup
-        >>> G = SymmetricPermutationGroup(4)
-        >>> G.degree
-        4
-
-        """
-        return self._deg
-
-    @property
-    def identity(self):
-        '''
-        Return the identity element of the SymmetricPermutationGroup.
-
-        Examples
-        ========
-
-        >>> from sympy.combinatorics import Permutation, SymmetricPermutationGroup
-        >>> G = SymmetricPermutationGroup(4)
-        >>> G.identity()
-        (3)
-
-        '''
-        return _af_new(list(range(self._deg)))
-
-
-class Coset(Basic):
-    """A left coset of a permutation group with respect to an element.
-
-    Parameters
-    ==========
-
-    g : Permutation
-
-    H : PermutationGroup
-
-    dir : "+" or "-", If not specified by default it will be "+"
-        here ``dir`` specified the type of coset "+" represent the
-        right coset and "-" represent the left coset.
-
-    G : PermutationGroup, optional
-        The group which contains *H* as its subgroup and *g* as its
-        element.
-
-        If not specified, it would automatically become a symmetric
-        group ``SymmetricPermutationGroup(g.size)`` and
-        ``SymmetricPermutationGroup(H.degree)`` if ``g.size`` and ``H.degree``
-        are matching.``SymmetricPermutationGroup`` is a lazy form of SymmetricGroup
-        used for representation purpose.
-
-    """
-
-    def __new__(cls, g, H, G=None, dir="+"):
-        g = _sympify(g)
-        if not isinstance(g, Permutation):
-            raise NotImplementedError
-
-        H = _sympify(H)
-        if not isinstance(H, PermutationGroup):
-            raise NotImplementedError
-
-        if G is not None:
-            G = _sympify(G)
-            if not isinstance(G, PermutationGroup) and not isinstance(G, SymmetricPermutationGroup):
-                raise NotImplementedError
-            if not H.is_subgroup(G):
-                raise ValueError("{} must be a subgroup of {}.".format(H, G))
-            if g not in G:
-                raise ValueError("{} must be an element of {}.".format(g, G))
-        else:
-            g_size = g.size
-            h_degree = H.degree
-            if g_size != h_degree:
-                raise ValueError(
-                    "The size of the permutation {} and the degree of "
-                    "the permutation group {} should be matching "
-                    .format(g, H))
-            G = SymmetricPermutationGroup(g.size)
-
-        if isinstance(dir, str):
-            dir = Symbol(dir)
-        elif not isinstance(dir, Symbol):
-            raise TypeError("dir must be of type basestring or "
-                    "Symbol, not %s" % type(dir))
-        if str(dir) not in ('+', '-'):
-            raise ValueError("dir must be one of '+' or '-' not %s" % dir)
-        obj = Basic.__new__(cls, g, H, G, dir)
-        obj._dir = dir
-        return obj
-
-    @property
-    def is_left_coset(self):
-        """
-        Check if the coset is left coset that is ``gH``.
-
-        Examples
-        ========
-
-        >>> from sympy.combinatorics import Permutation, PermutationGroup, Coset
-        >>> a = Permutation(1, 2)
-        >>> b = Permutation(0, 1)
-        >>> G = PermutationGroup([a, b])
-        >>> cst = Coset(a, G, dir="-")
-        >>> cst.is_left_coset
-        True
-
-        """
-        return str(self._dir) == '-'
-
-    @property
-    def is_right_coset(self):
-        """
-        Check if the coset is right coset that is ``Hg``.
-
-        Examples
-        ========
-
-        >>> from sympy.combinatorics import Permutation, PermutationGroup, Coset
-        >>> a = Permutation(1, 2)
-        >>> b = Permutation(0, 1)
-        >>> G = PermutationGroup([a, b])
-        >>> cst = Coset(a, G, dir="+")
-        >>> cst.is_right_coset
-        True
-
-        """
-        return str(self._dir) == '+'
-
-    def as_list(self):
-        """
-        Return all the elements of coset in the form of list.
-        """
-        g = self.args[0]
-        H = self.args[1]
-        cst = []
-        if str(self._dir) == '+':
-            for h in H.elements:
-                cst.append(h*g)
-        else:
-            for h in H.elements:
-                cst.append(g*h)
-        return cst
diff --git a/sympy/matrices/expressions/inverse.py b/sympy/matrices/expressions/inverse.py
index 37757830d6..97d7c9f7a4 100644
--- a/sympy/matrices/expressions/inverse.py
+++ b/sympy/matrices/expressions/inverse.py
@@ -67,8 +67,11 @@ def doit(self, **hints):
         if hints.get('deep', True):
             arg = arg.doit(**hints)
 
+        if isinstance(arg, MatPow):
+            return MatPow(arg.base, self.exp * arg.exp)
         return arg.inverse()
 
+
     def _eval_derivative_matrix_lines(self, x):
         arg = self.args[0]
         lines = arg._eval_derivative_matrix_lines(x)
diff --git a/sympy/matrices/expressions/kronecker.py b/sympy/matrices/expressions/kronecker.py
index d9c4b49601..7ede71fd13 100644
--- a/sympy/matrices/expressions/kronecker.py
+++ b/sympy/matrices/expressions/kronecker.py
@@ -391,7 +391,7 @@ def kronecker_mat_mul(expr):
 
 
 def kronecker_mat_pow(expr):
-    if isinstance(expr.base, KroneckerProduct) and all(a.is_square for a in expr.base.args):
+    if isinstance(expr.base, KroneckerProduct):
         return KroneckerProduct(*[MatPow(a, expr.exp) for a in expr.base.args])
     else:
         return expr
@@ -415,10 +415,8 @@ def combine_kronecker(expr):
     KroneckerProduct(A*B, B*A)
     >>> combine_kronecker(KroneckerProduct(A, B)+KroneckerProduct(B.T, A.T))
     KroneckerProduct(A + B.T, B + A.T)
-    >>> C = MatrixSymbol('C', n, n)
-    >>> D = MatrixSymbol('D', m, m)
-    >>> combine_kronecker(KroneckerProduct(C, D)**m)
-    KroneckerProduct(C**m, D**m)
+    >>> combine_kronecker(KroneckerProduct(A, B)**m)
+    KroneckerProduct(A**m, B**m)
     """
     def haskron(expr):
         return isinstance(expr, MatrixExpr) and expr.has(KroneckerProduct)
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 624e6e973d..14b4facef5 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -15,7 +15,6 @@
 from sympy.matrices.common import NonSquareMatrixError, NonInvertibleMatrixError
 from sympy.simplify import simplify
 from sympy.utilities.misc import filldedent
-from sympy.assumptions.ask import ask, Q
 
 
 def _sympifyit(arg, retval=None):
@@ -129,7 +128,15 @@ def __rmatmul__(self, other):
     @_sympifyit('other', NotImplemented)
     @call_highest_priority('__rpow__')
     def __pow__(self, other):
-        return MatPow(self, other).doit()
+        if not self.is_square:
+            raise NonSquareMatrixError("Power of non-square matrix %s" % self)
+        elif self.is_Identity:
+            return self
+        elif other == S.Zero:
+            return Identity(self.rows)
+        elif other == S.One:
+            return self
+        return MatPow(self, other).doit(deep=False)
 
     @_sympifyit('other', NotImplemented)
     @call_highest_priority('__pow__')
@@ -181,10 +188,6 @@ def _eval_transpose(self):
         return Transpose(self)
 
     def _eval_power(self, exp):
-        """
-        Override this in sub-classes to implement simplification of powers.  The cases where the exponent
-        is -1, 0, 1 are already covered in MatPow.doit(), so implementations can exclude these cases.
-        """
         return MatPow(self, exp)
 
     def _eval_simplify(self, **kwargs):
@@ -901,9 +904,6 @@ def _entry(self, i, j, **kwargs):
     def _eval_determinant(self):
         return S.One
 
-    def _eval_power(self, exp):
-        return self
-
 
 class GenericIdentity(Identity):
     """
@@ -967,9 +967,14 @@ def __new__(cls, m, n):
     def shape(self):
         return (self.args[0], self.args[1])
 
-    def _eval_power(self, exp):
-        # exp = -1, 0, 1 are already handled at this stage
-        if (exp < 0) == True:
+    @_sympifyit('other', NotImplemented)
+    @call_highest_priority('__rpow__')
+    def __pow__(self, other):
+        if other != 1 and not self.is_square:
+            raise NonSquareMatrixError("Power of non-square matrix %s" % self)
+        if other == 0:
+            return Identity(self.rows)
+        if other < 1:
             raise NonInvertibleMatrixError("Matrix det == 0; not invertible")
         return self
 
@@ -1047,10 +1052,6 @@ def __new__(cls, m, n, evaluate=False):
     def shape(self):
         return self._args
 
-    @property
-    def is_Identity(self):
-        return self._is_1x1() == True
-
     def as_explicit(self):
         from sympy import ImmutableDenseMatrix
         return ImmutableDenseMatrix.ones(*self.shape)
@@ -1061,16 +1062,6 @@ def doit(self, **hints):
             args = [a.doit(**hints) for a in args]
         return self.func(*args, evaluate=True)
 
-    def _eval_power(self, exp):
-        # exp = -1, 0, 1 are already handled at this stage
-        if self._is_1x1() == True:
-            return Identity(1)
-        if (exp < 0) == True:
-            raise NonInvertibleMatrixError("Matrix det == 0; not invertible")
-        if ask(Q.integer(exp)):
-            return self.shape[0] ** (exp - 1) * OneMatrix(*self.shape)
-        return super(OneMatrix, self)._eval_power(exp)
-
     def _eval_transpose(self):
         return OneMatrix(self.cols, self.rows)
 
diff --git a/sympy/matrices/expressions/matmul.py b/sympy/matrices/expressions/matmul.py
index 8f8472004d..adeed1d7dc 100644
--- a/sympy/matrices/expressions/matmul.py
+++ b/sympy/matrices/expressions/matmul.py
@@ -5,7 +5,7 @@
 from sympy.functions import adjoint
 from sympy.strategies import (rm_id, unpack, typed, flatten, exhaust,
         do_one, new)
-from sympy.matrices.common import ShapeError, NonInvertibleMatrixError
+from sympy.matrices.common import ShapeError
 from sympy.matrices.matrices import MatrixBase
 
 from .inverse import Inverse
@@ -333,17 +333,12 @@ def combine_powers(mul):
         if A_base == B_base:
             new_exp = A_exp + B_exp
             new_args[-1] = MatPow(A_base, new_exp).doit(deep=False)
-            continue
-        elif not isinstance(B_base, MatrixBase):
-            try:
-                B_base_inv = B_base.inverse()
-            except NonInvertibleMatrixError:
-                B_base_inv = None
-            if B_base_inv is not None and A_base == B_base_inv:
-                new_exp = A_exp - B_exp
-                new_args[-1] = MatPow(A_base, new_exp).doit(deep=False)
-                continue
-        new_args.append(B)
+        elif not isinstance(B_base, MatrixBase) and \
+            A_base == B_base.inverse():
+            new_exp = A_exp - B_exp
+            new_args[-1] = MatPow(A_base, new_exp).doit(deep=False)
+        else:
+            new_args.append(B)
 
     return newmul(factor, *new_args)
 
diff --git a/sympy/matrices/expressions/matpow.py b/sympy/matrices/expressions/matpow.py
index 14ac5be192..d64746d0fd 100644
--- a/sympy/matrices/expressions/matpow.py
+++ b/sympy/matrices/expressions/matpow.py
@@ -1,10 +1,13 @@
 from __future__ import print_function, division
 
 from sympy.matrices.common import NonSquareMatrixError
-from .matexpr import MatrixExpr, Identity
+from .matexpr import MatrixExpr, Identity, ZeroMatrix
 from sympy.core import S
 from sympy.core.sympify import _sympify
 from sympy.matrices import MatrixBase
+from sympy.matrices.common import NonInvertibleMatrixError
+
+from .permutation import PermutationMatrix
 
 
 class MatPow(MatrixExpr):
@@ -12,13 +15,10 @@ class MatPow(MatrixExpr):
     def __new__(cls, base, exp, evaluate=False, **options):
         base = _sympify(base)
         if not base.is_Matrix:
-            raise TypeError("MatPow base should be a matrix")
-        if not base.is_square:
-            raise NonSquareMatrixError("Power of non-square matrix %s" % base)
-
+            raise TypeError("Function parameter should be a matrix")
         exp = _sympify(exp)
-        obj = super().__new__(cls, base, exp)
 
+        obj = super(MatPow, cls).__new__(cls, base, exp)
         if evaluate:
             obj = obj.doit(deep=False)
 
@@ -41,7 +41,9 @@ def _entry(self, i, j, **kwargs):
         A = self.doit()
         if isinstance(A, MatPow):
             # We still have a MatPow, make an explicit MatMul out of it.
-            if A.exp.is_Integer and A.exp.is_positive:
+            if not A.base.is_square:
+                raise NonSquareMatrixError("Power of non-square matrix %s" % A.base)
+            elif A.exp.is_Integer and A.exp.is_positive:
                 A = MatMul(*[A.base for k in range(A.exp)])
             #elif A.exp.is_Integer and self.exp.is_negative:
             # Note: possible future improvement: in principle we can take
@@ -53,36 +55,42 @@ def _entry(self, i, j, **kwargs):
                 # Leave the expression unevaluated:
                 from sympy.matrices.expressions.matexpr import MatrixElement
                 return MatrixElement(self, i, j)
-        return A[i, j]
+        return A._entry(i, j)
 
     def doit(self, **kwargs):
-        if kwargs.get('deep', True):
-            base, exp = [arg.doit(**kwargs) for arg in self.args]
+        from sympy.matrices.expressions import Inverse
+        deep = kwargs.get('deep', True)
+        if deep:
+            args = [arg.doit(**kwargs) for arg in self.args]
         else:
-            base, exp = self.args
+            args = self.args
 
-        # combine all powers, e.g. (A ** 2) ** 3 -> A ** 6
+        base, exp = args
+        # combine all powers, e.g. (A**2)**3 = A**6
         while isinstance(base, MatPow):
-            exp *= base.args[1]
+            exp = exp*base.args[1]
             base = base.args[0]
 
-        if isinstance(base, MatrixBase):
-            # Delegate
-            return base ** exp
-
-        # Handle simple cases so that _eval_power() in MatrixExpr sub-classes can ignore them
-        if exp == S.One:
+        if exp.is_zero and base.is_square:
+            if isinstance(base, MatrixBase):
+                return base.func(Identity(base.shape[0]))
+            return Identity(base.shape[0])
+        elif isinstance(base, ZeroMatrix) and exp.is_negative:
+            raise NonInvertibleMatrixError("Matrix determinant is 0, not invertible")
+        elif isinstance(base, (Identity, ZeroMatrix)):
             return base
-        if exp == S.Zero:
-            return Identity(base.rows)
-        if exp == S.NegativeOne:
-            from sympy.matrices.expressions import Inverse
+        elif isinstance(base, PermutationMatrix):
+            return PermutationMatrix(base.args[0] ** exp).doit()
+        elif isinstance(base, MatrixBase):
+            if exp is S.One:
+                return base
+            return base**exp
+        # Note: just evaluate cases we know, return unevaluated on others.
+        # E.g., MatrixSymbol('x', n, m) to power 0 is not an error.
+        elif exp is S.NegativeOne and base.is_square:
             return Inverse(base).doit(**kwargs)
-
-        eval_power = getattr(base, '_eval_power', None)
-        if eval_power is not None:
-            return eval_power(exp)
-
+        elif exp is S.One:
+            return base
         return MatPow(base, exp)
 
     def _eval_transpose(self):
@@ -136,6 +144,3 @@ def _eval_derivative_matrix_lines(self, x):
         else:
             raise NotImplementedError("cannot evaluate %s derived by %s" % (self, x))
         return newexpr._eval_derivative_matrix_lines(x)
-
-    def _eval_inverse(self):
-        return MatPow(self.base, -self.exp)
diff --git a/sympy/matrices/expressions/permutation.py b/sympy/matrices/expressions/permutation.py
index e8df0c8d8d..9c51e78e8e 100644
--- a/sympy/matrices/expressions/permutation.py
+++ b/sympy/matrices/expressions/permutation.py
@@ -85,9 +85,6 @@ def _entry(self, i, j, **kwargs):
         perm = self.args[0]
         return KroneckerDelta(perm.apply(i), j)
 
-    def _eval_power(self, exp):
-        return PermutationMatrix(self.args[0] ** exp).doit()
-
     def _eval_inverse(self):
         return PermutationMatrix(self.args[0] ** -1)
 
diff --git a/sympy/polys/polyroots.py b/sympy/polys/polyroots.py
index d3bf7e06e4..f07fb011a8 100644
--- a/sympy/polys/polyroots.py
+++ b/sympy/polys/polyroots.py
@@ -870,7 +870,11 @@ def roots(f, *gens, **flags):
         f = Poly(poly, x, field=True)
     else:
         try:
-            f = Poly(f, *gens, **flags)
+            F = Poly(f, *gens, **flags)
+            if not isinstance(f, Poly) and not F.gen.is_Symbol:
+                raise PolynomialError("generator must be a Symbol")
+            else:
+                f = F
             if f.length == 2 and f.degree() != 1:
                 # check for foo**n factors in the constant
                 n = f.degree()
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index 68b69d1a24..6824c7f410 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -6429,6 +6429,11 @@ def refine_root(f, s, t, eps=None, steps=None, fast=False, check_sqf=False):
     """
     try:
         F = Poly(f)
+        if not isinstance(f, Poly) and not F.gen.is_Symbol:
+            # root of sin(x) + 1 is -1 but when someone
+            # passes an Expr instead of Poly they may not expect
+            # that the generator will be sin(x), not x
+            raise PolynomialError("generator must be a Symbol")
     except GeneratorsNeeded:
         raise PolynomialError(
             "can't refine a root of %s, not a polynomial" % f)
@@ -6458,6 +6463,11 @@ def count_roots(f, inf=None, sup=None):
     """
     try:
         F = Poly(f, greedy=False)
+        if not isinstance(f, Poly) and not F.gen.is_Symbol:
+            # root of sin(x) + 1 is -1 but when someone
+            # passes an Expr instead of Poly they may not expect
+            # that the generator will be sin(x), not x
+            raise PolynomialError("generator must be a Symbol")
     except GeneratorsNeeded:
         raise PolynomialError("can't count roots of %s, not a polynomial" % f)
 
@@ -6477,10 +6487,14 @@ def real_roots(f, multiple=True):
 
     >>> real_roots(2*x**3 - 7*x**2 + 4*x + 4)
     [-1/2, 2, 2]
-
     """
     try:
         F = Poly(f, greedy=False)
+        if not isinstance(f, Poly) and not F.gen.is_Symbol:
+            # root of sin(x) + 1 is -1 but when someone
+            # passes an Expr instead of Poly they may not expect
+            # that the generator will be sin(x), not x
+            raise PolynomialError("generator must be a Symbol")
     except GeneratorsNeeded:
         raise PolynomialError(
             "can't compute real roots of %s, not a polynomial" % f)
@@ -6507,6 +6521,11 @@ def nroots(f, n=15, maxsteps=50, cleanup=True):
     """
     try:
         F = Poly(f, greedy=False)
+        if not isinstance(f, Poly) and not F.gen.is_Symbol:
+            # root of sin(x) + 1 is -1 but when someone
+            # passes an Expr instead of Poly they may not expect
+            # that the generator will be sin(x), not x
+            raise PolynomialError("generator must be a Symbol")
     except GeneratorsNeeded:
         raise PolynomialError(
             "can't compute numerical roots of %s, not a polynomial" % f)
@@ -6533,6 +6552,11 @@ def ground_roots(f, *gens, **args):
 
     try:
         F, opt = poly_from_expr(f, *gens, **args)
+        if not isinstance(f, Poly) and not F.gen.is_Symbol:
+            # root of sin(x) + 1 is -1 but when someone
+            # passes an Expr instead of Poly they may not expect
+            # that the generator will be sin(x), not x
+            raise PolynomialError("generator must be a Symbol")
     except PolificationFailed as exc:
         raise ComputationFailed('ground_roots', 1, exc)
 
@@ -6567,6 +6591,11 @@ def nth_power_roots_poly(f, n, *gens, **args):
 
     try:
         F, opt = poly_from_expr(f, *gens, **args)
+        if not isinstance(f, Poly) and not F.gen.is_Symbol:
+            # root of sin(x) + 1 is -1 but when someone
+            # passes an Expr instead of Poly they may not expect
+            # that the generator will be sin(x), not x
+            raise PolynomialError("generator must be a Symbol")
     except PolificationFailed as exc:
         raise ComputationFailed('nth_power_roots_poly', 1, exc)
 
diff --git a/sympy/polys/rootoftools.py b/sympy/polys/rootoftools.py
index 30980f5290..3a3eaca756 100644
--- a/sympy/polys/rootoftools.py
+++ b/sympy/polys/rootoftools.py
@@ -21,11 +21,13 @@
 from sympy.polys.rootisolation import (
     dup_isolate_complex_roots_sqf,
     dup_isolate_real_roots_sqf)
-from sympy.utilities import lambdify, public, sift
+from sympy.utilities import lambdify, public, sift, numbered_symbols
 
 from mpmath import mpf, mpc, findroot, workprec
 from mpmath.libmp.libmpf import dps_to_prec, prec_to_dps
 
+from itertools import chain
+
 
 __all__ = ['CRootOf']
 
@@ -733,13 +735,22 @@ def _get_roots(cls, method, poly, radicals):
         """Return postprocessed roots of specified kind. """
         if not poly.is_univariate:
             raise PolynomialError("only univariate polynomials are allowed")
-
+        # get rid of gen and it's free symbol
+        d = Dummy()
+        poly = poly.subs(poly.gen, d)
+        x = symbols('x')
+        # see what others are left and select x or a numbered x
+        # that doesn't clash
+        free_names = {str(i) for i in poly.free_symbols}
+        for x in chain((symbols('x'),), numbered_symbols('x')):
+            if x.name not in free_names:
+                poly = poly.xreplace({d: x})
+                break
         coeff, poly = cls._preprocess_roots(poly)
         roots = []
 
         for root in getattr(cls, method)(poly):
             roots.append(coeff*cls._postprocess_root(root, radicals))
-
         return roots
 
     @classmethod
diff --git a/sympy/simplify/cse_opts.py b/sympy/simplify/cse_opts.py
index f02b6e61a7..754394d43f 100644
--- a/sympy/simplify/cse_opts.py
+++ b/sympy/simplify/cse_opts.py
@@ -14,8 +14,7 @@ def sub_pre(e):
     """
     # replacing Add, A, from which -1 can be extracted with -1*-A
     adds = [a for a in e.atoms(Add) if a.could_extract_minus_sign()]
-    reps = dict((a, Mul(S.NegativeOne, -a, evaluate=False)) for a in adds)
-
+    reps = dict((a, Mul._from_args([S.NegativeOne, -a])) for a in adds)
     e = e.xreplace(reps)
 
     # repeat again for persisting Adds but mark these with a leading 1, -1
@@ -26,7 +25,7 @@ def sub_pre(e):
             if a in reps:
                 negs[a] = reps[a]
             elif a.could_extract_minus_sign():
-                negs[a] = Mul(S.One, S.NegativeOne, -a, evaluate=False)
+                negs[a] = Mul._from_args([S.One, S.NegativeOne, -a])
         e = e.xreplace(negs)
     return e
 
@@ -38,7 +37,7 @@ def sub_post(e):
     for node in preorder_traversal(e):
         if isinstance(node, Mul) and \
             node.args[0] is S.One and node.args[1] is S.NegativeOne:
-            replacements.append((node, -Mul(*node.args[2:], evaluate=False)))
+            replacements.append((node, -Mul._from_args(node.args[2:])))
     for node, replacement in replacements:
         e = e.xreplace({node: replacement})
 
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index 6205ad24c3..64c644fca0 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -40,7 +40,7 @@
     powdenest, nsimplify, denom, logcombine, sqrtdenest, fraction,
     separatevars)
 from sympy.simplify.sqrtdenest import sqrt_depth
-from sympy.simplify.fu import TR1
+from sympy.simplify.fu import TR1, TR2i
 from sympy.matrices.common import NonInvertibleMatrixError
 from sympy.matrices import Matrix, zeros
 from sympy.polys import roots, cancel, factor, Poly, degree
@@ -791,8 +791,8 @@ def solve(f, *symbols, **flags):
         >>> from sympy import real_root, S
         >>> eq = root(x, 3) - root(x, 5) + S(1)/7
         >>> solve(eq)  # this gives 2 solutions but misses a 3rd
-        [CRootOf(7*_p**5 - 7*_p**3 + 1, 1)**15,
-        CRootOf(7*_p**5 - 7*_p**3 + 1, 2)**15]
+        [CRootOf(7*x**5 - 7*x**3 + 1, 1)**15,
+        CRootOf(7*x**5 - 7*x**3 + 1, 2)**15]
         >>> sol = solve(eq, check=False)
         >>> [abs(eq.subs(x,i).n(2)) for i in sol]
         [0.48, 0.e-110, 0.e-110, 0.052, 0.052]
@@ -1596,7 +1596,26 @@ def _as_base_q(x):
                     isinstance(_, TrigonometricFunction)])
                 other = funcs - trig
                 if not other and len(funcs.intersection(trig)) > 1:
-                    newf = TR1(f_num).rewrite(tan)
+                    newf = None
+                    if f_num.is_Add and len(f_num.args) == 2:
+                        # check for sin(x)**p = cos(x)**p
+                        _args = f_num.args
+                        t = a, b = [i.atoms(Function).intersection(
+                            trig) for i in _args]
+                        if all(len(i) == 1 for i in t):
+                            a, b = [i.pop() for i in t]
+                            if isinstance(a, cos):
+                                a, b = b, a
+                                _args = _args[::-1]
+                            if isinstance(a, sin) and isinstance(b, cos
+                                    ) and a.args[0] == b.args[0]:
+                                # sin(x) + cos(x) = 0 -> tan(x) + 1 = 0
+                                newf, _d = (TR2i(_args[0]/_args[1]) + 1
+                                    ).as_numer_denom()
+                                if not _d.is_Number:
+                                    newf = None
+                    if newf is None:
+                        newf = TR1(f_num).rewrite(tan)
                     if newf != f_num:
                         # don't check the rewritten form --check
                         # solutions in the un-rewritten form below


Test patch:
diff --git a/sympy/combinatorics/tests/test_perm_groups.py b/sympy/combinatorics/tests/test_perm_groups.py
index 74144014dc..2176ad2723 100644
--- a/sympy/combinatorics/tests/test_perm_groups.py
+++ b/sympy/combinatorics/tests/test_perm_groups.py
@@ -1,5 +1,5 @@
 from sympy.combinatorics.perm_groups import (PermutationGroup,
-    _orbit_transversal, Coset, SymmetricPermutationGroup)
+    _orbit_transversal)
 from sympy.combinatorics.named_groups import SymmetricGroup, CyclicGroup,\
     DihedralGroup, AlternatingGroup, AbelianGroup, RubikGroup
 from sympy.combinatorics.permutations import Permutation
@@ -1116,46 +1116,3 @@ def test_conjugacy_classes():
 
     assert len(expected) == len(computed)
     assert all(e in computed for e in expected)
-
-def test_coset_class():
-    a = Permutation(1, 2)
-    b = Permutation(0, 1)
-    G = PermutationGroup([a, b])
-    #Creating right coset
-    rht_coset = Coset(a, G, dir='+')
-    #Checking whether it is left coset or right coset
-    assert rht_coset.is_right_coset
-    assert not rht_coset.is_left_coset
-    #Creating list representation of coset
-    list_repr = rht_coset.as_list()
-    expected = [Permutation(0, 2), Permutation(0, 2, 1), Permutation(1, 2), Permutation(2), Permutation(2)(0, 1), Permutation(0, 1, 2)]
-    for ele in list_repr:
-        assert ele in expected
-    #Creating left coset
-    left_coset = Coset(a, G, dir='-')
-    #Checking whether it is left coset or right coset
-    assert not left_coset.is_right_coset
-    assert left_coset.is_left_coset
-    #Creating list representation of Coset
-    list_repr = left_coset.as_list()
-    expected = [Permutation(2)(0, 1), Permutation(0, 1, 2), Permutation(1, 2),
-    Permutation(2), Permutation(0, 2), Permutation(0, 2, 1)]
-    for ele in list_repr:
-        assert ele in expected
-
-    G = PermutationGroup(Permutation(1, 2, 3, 4), Permutation(2, 3, 4))
-    H = PermutationGroup(Permutation(1, 2, 3, 4))
-    g = Permutation(1, 3)(2, 4)
-    rht_coset = Coset(g, H, G, dir='+')
-    assert rht_coset.is_right_coset
-    list_repr = rht_coset.as_list()
-    expected = [Permutation(1, 2, 3, 4), Permutation(4), Permutation(1, 3)(2, 4),
-    Permutation(1, 4, 3, 2)]
-    for ele in list_repr:
-        assert ele in expected
-
-def test_symmetricpermutationgroup():
-    a = SymmetricPermutationGroup(5)
-    assert a.degree == 5
-    assert a.order() == 120
-    assert a.identity() == Permutation(4)
diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index bf78acf3c4..5599fcc703 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -4957,15 +4957,3 @@ def test_sympy__integrals__rubi__utility_function__ProductLog():
 def test_sympy__combinatorics__schur_number__SchurNumber():
     from sympy.combinatorics.schur_number import SchurNumber
     assert _test_args(SchurNumber(1))
-
-def test_sympy__combinatorics__perm_groups__SymmetricPermutationGroup():
-    from sympy.combinatorics.perm_groups import SymmetricPermutationGroup
-    assert _test_args(SymmetricPermutationGroup(5))
-
-def test_sympy__combinatorics__perm_groups__Coset():
-    from sympy.combinatorics.permutations import Permutation
-    from sympy.combinatorics.perm_groups import PermutationGroup, Coset
-    a = Permutation(1, 2)
-    b = Permutation(0, 1)
-    G = PermutationGroup([a, b])
-    assert _test_args(Coset(a, G))
diff --git a/sympy/matrices/expressions/tests/test_inverse.py b/sympy/matrices/expressions/tests/test_inverse.py
index 83eb770b61..4fe753870f 100644
--- a/sympy/matrices/expressions/tests/test_inverse.py
+++ b/sympy/matrices/expressions/tests/test_inverse.py
@@ -19,7 +19,6 @@ def test_inverse():
     assert Inverse(A*E).shape == (n, n)
     assert Inverse(E*A).shape == (m, m)
     assert Inverse(C).inverse() == C
-    assert Inverse(Inverse(C)).doit() == C
     assert isinstance(Inverse(Inverse(C)), Inverse)
 
     assert Inverse(*Inverse(E*A).args) == Inverse(E*A)
diff --git a/sympy/matrices/expressions/tests/test_kronecker.py b/sympy/matrices/expressions/tests/test_kronecker.py
index 230df2119b..4a7b142945 100644
--- a/sympy/matrices/expressions/tests/test_kronecker.py
+++ b/sympy/matrices/expressions/tests/test_kronecker.py
@@ -127,8 +127,6 @@ def test_KroneckerProduct_combine_pow():
                              ** 2) == x * KroneckerProduct(X**2, Y**2)
     assert combine_kronecker(
         x * (KroneckerProduct(X, Y)**2) * KroneckerProduct(A, B)) == x * KroneckerProduct(X**2 * A, Y**2 * B)
-    # cannot simplify because of non-square arguments to kronecker product:
-    assert combine_kronecker(KroneckerProduct(A, B.T) ** m) == KroneckerProduct(A, B.T) ** m
 
 
 def test_KroneckerProduct_expand():
diff --git a/sympy/matrices/expressions/tests/test_matexpr.py b/sympy/matrices/expressions/tests/test_matexpr.py
index b75e6ff63d..af0a609e7d 100644
--- a/sympy/matrices/expressions/tests/test_matexpr.py
+++ b/sympy/matrices/expressions/tests/test_matexpr.py
@@ -7,7 +7,7 @@
 from sympy.simplify import simplify
 from sympy.matrices import (Identity, ImmutableMatrix, Inverse, MatAdd, MatMul,
         MatPow, Matrix, MatrixExpr, MatrixSymbol, ShapeError, ZeroMatrix,
-        SparseMatrix, Transpose, Adjoint, NonSquareMatrixError)
+        SparseMatrix, Transpose, Adjoint)
 from sympy.matrices.expressions.matexpr import (MatrixElement,
                                                 GenericZeroMatrix, GenericIdentity, OneMatrix)
 from sympy.testing.pytest import raises, XFAIL
@@ -133,11 +133,9 @@ def test_ZeroMatrix():
     assert Z.conjugate() == Z
 
     assert ZeroMatrix(n, n)**0 == Identity(n)
-    with raises(NonSquareMatrixError):
+    with raises(ShapeError):
         Z**0
-    with raises(NonSquareMatrixError):
-        Z**1
-    with raises(NonSquareMatrixError):
+    with raises(ShapeError):
         Z**2
 
 
@@ -159,11 +157,9 @@ def test_OneMatrix():
     assert U.conjugate() == U
 
     assert OneMatrix(n, n) ** 0 == Identity(n)
-    with raises(NonSquareMatrixError):
+    with raises(ShapeError):
         U ** 0
-    with raises(NonSquareMatrixError):
-        U ** 1
-    with raises(NonSquareMatrixError):
+    with raises(ShapeError):
         U ** 2
     with raises(ShapeError):
         a + U
@@ -294,7 +290,7 @@ def test_MatPow():
     assert (A**2)**3 == A**6
     assert A**S.Half == sqrt(A)
     assert A**Rational(1, 3) == cbrt(A)
-    raises(NonSquareMatrixError, lambda: MatrixSymbol('B', 3, 2)**2)
+    raises(ShapeError, lambda: MatrixSymbol('B', 3, 2)**2)
 
 
 def test_MatrixSymbol():
diff --git a/sympy/matrices/expressions/tests/test_matpow.py b/sympy/matrices/expressions/tests/test_matpow.py
index 4903b3506a..b52d807513 100644
--- a/sympy/matrices/expressions/tests/test_matpow.py
+++ b/sympy/matrices/expressions/tests/test_matpow.py
@@ -1,12 +1,8 @@
-from sympy.simplify.powsimp import powsimp
 from sympy.testing.pytest import raises
-from sympy.core.expr import unchanged
-from sympy.core import symbols, S
-from sympy.matrices import Identity, MatrixSymbol, ImmutableMatrix, ZeroMatrix, OneMatrix
-from sympy.matrices.common import NonSquareMatrixError
+from sympy.core import symbols, pi, S
+from sympy.matrices import Identity, MatrixSymbol, ImmutableMatrix, ZeroMatrix
 from sympy.matrices.expressions import MatPow, MatAdd, MatMul
-from sympy.matrices.expressions.inverse import Inverse
-from sympy.matrices.expressions.matexpr import MatrixElement
+from sympy.matrices.common import ShapeError
 
 n, m, l, k = symbols('n m l k', integer=True)
 A = MatrixSymbol('A', n, m)
@@ -16,22 +12,12 @@
 E = MatrixSymbol('E', m, n)
 
 
-def test_entry_matrix():
-    X = ImmutableMatrix([[1, 2], [3, 4]])
-    assert MatPow(X, 0)[0, 0] == 1
-    assert MatPow(X, 0)[0, 1] == 0
-    assert MatPow(X, 1)[0, 0] == 1
-    assert MatPow(X, 1)[0, 1] == 2
-    assert MatPow(X, 2)[0, 0] == 7
-
-
-def test_entry_symbol():
+def test_entry():
     from sympy.concrete import Sum
+    assert MatPow(A, 1)[0, 0] == A[0, 0]
     assert MatPow(C, 0)[0, 0] == 1
     assert MatPow(C, 0)[0, 1] == 0
-    assert MatPow(C, 1)[0, 0] == C[0, 0]
     assert isinstance(MatPow(C, 2)[0, 0], Sum)
-    assert isinstance(MatPow(C, n)[0, 0], MatrixElement)
 
 
 def test_as_explicit_symbol():
@@ -39,13 +25,16 @@ def test_as_explicit_symbol():
     assert MatPow(X, 0).as_explicit() == ImmutableMatrix(Identity(2))
     assert MatPow(X, 1).as_explicit() == X.as_explicit()
     assert MatPow(X, 2).as_explicit() == (X.as_explicit())**2
-    assert MatPow(X, n).as_explicit() == ImmutableMatrix([
-        [(X ** n)[0, 0], (X ** n)[0, 1]],
-        [(X ** n)[1, 0], (X ** n)[1, 1]],
-    ])
 
 
-def test_as_explicit_matrix():
+def test_as_explicit_nonsquare_symbol():
+    X = MatrixSymbol('X', 2, 3)
+    assert MatPow(X, 1).as_explicit() == X.as_explicit()
+    for r in [0, 2, S.Half, S.Pi]:
+        raises(ShapeError, lambda: MatPow(X, r).as_explicit())
+
+
+def test_as_explicit():
     A = ImmutableMatrix([[1, 2], [3, 4]])
     assert MatPow(A, 0).as_explicit() == ImmutableMatrix(Identity(2))
     assert MatPow(A, 1).as_explicit() == A
@@ -53,19 +42,33 @@ def test_as_explicit_matrix():
     assert MatPow(A, -1).as_explicit() == A.inv()
     assert MatPow(A, -2).as_explicit() == (A.inv())**2
     # less expensive than testing on a 2x2
-    A = ImmutableMatrix([4])
+    A = ImmutableMatrix([4]);
     assert MatPow(A, S.Half).as_explicit() == A**S.Half
 
 
-def test_doit_symbol():
+def test_as_explicit_nonsquare():
+    A = ImmutableMatrix([[1, 2, 3], [4, 5, 6]])
+    assert MatPow(A, 1).as_explicit() == A
+    raises(ShapeError, lambda: MatPow(A, 0).as_explicit())
+    raises(ShapeError, lambda: MatPow(A, 2).as_explicit())
+    raises(ShapeError, lambda: MatPow(A, -1).as_explicit())
+    raises(ValueError, lambda: MatPow(A, pi).as_explicit())
+
+
+def test_doit_nonsquare_MatrixSymbol():
+    assert MatPow(A, 1).doit() == A
+    for r in [0, 2, -1, pi]:
+        assert MatPow(A, r).doit() == MatPow(A, r)
+
+
+def test_doit_square_MatrixSymbol_symsize():
     assert MatPow(C, 0).doit() == Identity(n)
     assert MatPow(C, 1).doit() == C
-    assert MatPow(C, -1).doit() == C.I
-    for r in [2, S.Half, S.Pi, n]:
+    for r in [2, pi]:
         assert MatPow(C, r).doit() == MatPow(C, r)
 
 
-def test_doit_matrix():
+def test_doit_with_MatrixBase():
     X = ImmutableMatrix([[1, 2], [3, 4]])
     assert MatPow(X, 0).doit() == ImmutableMatrix(Identity(2))
     assert MatPow(X, 1).doit() == X
@@ -79,12 +82,13 @@ def test_doit_matrix():
     raises(ValueError, lambda: MatPow(X,-2).doit())
 
 
-def test_nonsquare():
-    A = MatrixSymbol('A', 2, 3)
-    B = ImmutableMatrix([[1, 2, 3], [4, 5, 6]])
-    for r in [-1, 0, 1, 2, S.Half, S.Pi, n]:
-        raises(NonSquareMatrixError, lambda: MatPow(A, r))
-        raises(NonSquareMatrixError, lambda: MatPow(B, r))
+def test_doit_nonsquare():
+    X = ImmutableMatrix([[1, 2, 3], [4, 5, 6]])
+    assert MatPow(X, 1).doit() == X
+    raises(ShapeError, lambda: MatPow(X, 0).doit())
+    raises(ShapeError, lambda: MatPow(X, 2).doit())
+    raises(ShapeError, lambda: MatPow(X, -1).doit())
+    raises(ShapeError, lambda: MatPow(X, pi).doit())
 
 
 def test_doit_equals_pow(): #17179
@@ -126,19 +130,6 @@ def test_zero_power():
     raises(ValueError, lambda:MatPow(z2, -1).doit())
 
 
-def test_OneMatrix_power():
-    o = OneMatrix(3, 3)
-    assert o ** 0 == Identity(3)
-    assert o ** 1 == o
-    assert o * o == o ** 2 == 3 * o
-    assert o * o * o == o ** 3 == 9 * o
-
-    o = OneMatrix(n, n)
-    assert o * o == o ** 2 == n * o
-    # powsimp necessary as n ** (n - 2) * n does not produce n ** (n - 1)
-    assert powsimp(o ** (n - 1) * o) == o ** n == n ** (n - 1) * o
-
-
 def test_transpose_power():
     from sympy.matrices.expressions.transpose import Transpose as TP
 
@@ -154,33 +145,3 @@ def test_transpose_power():
 
     assert ((D*C)**-5).T**-5 == ((D*C)**25).T
     assert (((D*C)**l).T**k).T == (D*C)**(l*k)
-
-
-def test_Inverse():
-    assert Inverse(MatPow(C, 0)).doit() == Identity(n)
-    assert Inverse(MatPow(C, 1)).doit() == Inverse(C)
-    assert Inverse(MatPow(C, 2)).doit() == MatPow(C, -2)
-    assert Inverse(MatPow(C, -1)).doit() == C
-
-    assert MatPow(Inverse(C), 0).doit() == Identity(n)
-    assert MatPow(Inverse(C), 1).doit() == Inverse(C)
-    assert MatPow(Inverse(C), 2).doit() == MatPow(C, -2)
-    assert MatPow(Inverse(C), -1).doit() == C
-
-
-def test_combine_powers():
-    assert (C ** 1) ** 1 == C
-    assert (C ** 2) ** 3 == MatPow(C, 6)
-    assert (C ** -2) ** -3 == MatPow(C, 6)
-    assert (C ** -1) ** -1 == C
-    assert (((C ** 2) ** 3) ** 4) ** 5 == MatPow(C, 120)
-    assert (C ** n) ** n == C ** (n ** 2)
-
-
-def test_unchanged():
-    assert unchanged(MatPow, C, 0)
-    assert unchanged(MatPow, C, 1)
-    assert unchanged(MatPow, Inverse(C), -1)
-    assert unchanged(Inverse, MatPow(C, -1), -1)
-    assert unchanged(MatPow, MatPow(C, -1), -1)
-    assert unchanged(MatPow, MatPow(C, 1), 1)
diff --git a/sympy/polys/tests/test_polyroots.py b/sympy/polys/tests/test_polyroots.py
index 786d534afd..c18d65f48d 100644
--- a/sympy/polys/tests/test_polyroots.py
+++ b/sympy/polys/tests/test_polyroots.py
@@ -11,6 +11,7 @@
     roots_binomial, preprocess_roots, roots)
 
 from sympy.polys.orthopolys import legendre_poly
+from sympy.polys.polyerrors import PolynomialError
 from sympy.polys.polyutils import _nsort
 
 from sympy.utilities.iterables import cartes
@@ -692,3 +693,8 @@ def test_nroots2():
 
 def test_roots_composite():
     assert len(roots(Poly(y**3 + y**2*sqrt(x) + y + x, y, composite=True))) == 3
+
+
+def test_issue_19113():
+    eq = cos(x)**3 - cos(x) + 1
+    raises(PolynomialError, lambda: roots(eq))
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index 174c86c4d8..19a86c46e8 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -3323,17 +3323,30 @@ def test_issue_15669():
         2*2**Rational(4, 5)*x*(-x**2 + sqrt(8*x**2 + (x**2 - 2)**2) + 2)**Rational(3, 5) + 10*x)
     assert factor(expr, deep=True) == x*(x**2 + 2)
 
+
 def test_issue_17988():
     x = Symbol('x')
     p = poly(x - 1)
     M = Matrix([[poly(x + 1), poly(x + 1)]])
     assert p * M == M * p == Matrix([[poly(x**2 - 1), poly(x**2 - 1)]])
 
+
 def test_issue_18205():
     assert cancel((2 + I)*(3 - I)) == 7 + I
     assert cancel((2 + I)*(2 - I)) == 5
 
+
 def test_issue_8695():
     p = (x**2 + 1) * (x - 1)**2 * (x - 2)**3 * (x - 3)**3
     result = (1, [(x**2 + 1, 1), (x - 1, 2), (x**2 - 5*x + 6, 3)])
     assert sqf_list(p) == result
+
+
+def test_issue_19113():
+    eq = sin(x)**3 - sin(x) + 1
+    raises(PolynomialError, lambda: refine_root(eq, 1, 2, 1e-2))
+    raises(PolynomialError, lambda: count_roots(eq, -1, 1))
+    raises(PolynomialError, lambda: real_roots(eq))
+    raises(PolynomialError, lambda: nroots(eq))
+    raises(PolynomialError, lambda: ground_roots(eq))
+    raises(PolynomialError, lambda: nth_power_roots_poly(eq, 2))
diff --git a/sympy/polys/tests/test_rootoftools.py b/sympy/polys/tests/test_rootoftools.py
index 88c84f96b9..5ef5fd96fb 100644
--- a/sympy/polys/tests/test_rootoftools.py
+++ b/sympy/polys/tests/test_rootoftools.py
@@ -574,3 +574,13 @@ def test_issue_15920():
     r = rootof(x**5 - x + 1, 0)
     p = Integral(x, (x, 1, y))
     assert unchanged(Eq, r, p)
+
+
+def test_issue_19113():
+    eq = y**3 - y + 1
+    # generator is a canonical x in RootOf
+    assert str(Poly(eq).real_roots()) == '[CRootOf(x**3 - x + 1, 0)]'
+    assert str(Poly(eq.subs(y, tan(y))).real_roots()
+        ) == '[CRootOf(x**3 - x + 1, 0)]'
+    assert str(Poly(eq.subs(y, tan(x))).real_roots()
+        ) == '[CRootOf(x**3 - x + 1, 0)]'
diff --git a/sympy/simplify/tests/test_cse.py b/sympy/simplify/tests/test_cse.py
index f076d8f533..4e46adafe0 100644
--- a/sympy/simplify/tests/test_cse.py
+++ b/sympy/simplify/tests/test_cse.py
@@ -5,7 +5,7 @@
 from sympy import (
     Add, Mul, Pow, Symbol, exp, sqrt, symbols, sympify, cse,
     Matrix, S, cos, sin, Eq, Function, Tuple, CRootOf,
-    IndexedBase, Idx, Piecewise, O, signsimp
+    IndexedBase, Idx, Piecewise, O
 )
 from sympy.core.function import count_ops
 from sympy.simplify.cse_opts import sub_pre, sub_post
@@ -545,7 +545,3 @@ def test_issue_18203():
 def test_unevaluated_mul():
     eq = Mul(x + y, x + y, evaluate=False)
     assert cse(eq) == ([(x0, x + y)], [x0**2])
-
-def test_issue_18991():
-    A = MatrixSymbol('A', 2, 2)
-    assert signsimp(-A * A - A) == -A * A - A
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index 88bd8d0cb1..b43e690fe9 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -356,7 +356,7 @@ def test_solve_transcendental():
     assert set(solve((a*x + b)*(exp(x) - 3), x)) == set([-b/a, log(3)])
     assert solve(cos(x) - y, x) == [-acos(y) + 2*pi, acos(y)]
     assert solve(2*cos(x) - y, x) == [-acos(y/2) + 2*pi, acos(y/2)]
-    assert solve(Eq(cos(x), sin(x)), x) == [pi*Rational(-3, 4), pi/4]
+    assert solve(Eq(cos(x), sin(x)), x) == [pi/4]
 
     assert set(solve(exp(x) + exp(-x) - y, x)) in [set([
         log(y/2 - sqrt(y**2 - 4)/2),
@@ -510,6 +510,7 @@ def test_solve_for_functions_derivatives():
             f(x).diff(x, 3), 1 - f(x).diff(x,3)], f(x).diff(x, 2), f(x).diff(x, 3))
     assert soln == { f(x).diff(x, 2): 0, f(x).diff(x, 3): 1 }
 
+
 def test_issue_3725():
     f = Function('f')
     F = x**2 + f(x)**2 - 4*x - 1
@@ -734,8 +735,7 @@ def test_issue_4671_4463_4467():
     assert solve(1 - log(a + 4*x**2), x) in (
         [-sqrt(-a + E)/2, sqrt(-a + E)/2],
         [sqrt(-a + E)/2, -sqrt(-a + E)/2],)
-    assert set(solve((
-        a**2 + 1) * (sin(a*x) + cos(a*x)), x)) == set([-pi/(4*a), 3*pi/(4*a)])
+    assert solve((a**2 + 1)*(sin(a*x) + cos(a*x)), x) == [-pi/(4*a)]
     assert solve(3 - (sinh(a*x) + cosh(a*x)), x) == [log(3)/a]
     assert set(solve(3 - (sinh(a*x) + cosh(a*x)**2), x)) == \
         set([log(-2 + sqrt(5))/a, log(-sqrt(2) + 1)/a,
@@ -1160,6 +1160,7 @@ def test_checksol():
     raises(ValueError, lambda: checksol(x, 1))
     raises(ValueError, lambda: checksol([], x, 1))
 
+
 def test__invert():
     assert _invert(x - 2) == (2, x)
     assert _invert(2) == (2, 0)
@@ -1370,6 +1371,7 @@ def test_failing_assumptions():
     'extended_nonpositive': None, 'extended_nonzero': None,
     'extended_positive': None }
 
+
 def test_issue_6056():
     assert solve(tanh(x + 3)*tanh(x - 3) - 1) == []
     assert solve(tanh(x - 1)*tanh(x + 1) + 1) == \
@@ -1692,8 +1694,7 @@ def test_rewrite_trig():
     assert solve(sinh(x) + tanh(x)) == [0, I*pi]
 
     # issue 6157
-    assert solve(2*sin(x) - cos(x), x) == [-2*atan(2 - sqrt(5)),
-                                           -2*atan(2 + sqrt(5))]
+    assert solve(2*sin(x) - cos(x), x) == [atan(S.Half)]
 
 
 @XFAIL
@@ -2057,6 +2058,7 @@ def test_issue_15307():
     eq2 = Eq(-2*x + 8, 2*x - 40)
     assert solve([eq1, eq2]) == {x:12, y:75}
 
+
 def test_issue_15415():
     assert solve(x - 3, x) == [3]
     assert solve([x - 3], x) == {x:3}
@@ -2152,6 +2154,7 @@ def test_issue_17949():
     assert solve(exp(+x-x**2), x) == []
     assert solve(exp(-x-x**2), x) == []
 
+
 def test_issue_10993():
     assert solve(Eq(binomial(x, 2), 3)) == [-2, 3]
     assert solve(Eq(pow(x, 2) + binomial(x, 3), x)) == [-4, 0, 1]
@@ -2160,9 +2163,25 @@ def test_issue_10993():
     assert solve(x-binomial(a, 3) + binomial(y, 2) + sin(a), x) == [-sin(a) + binomial(a, 3) - binomial(y, 2)]
     assert solve((x+1)-binomial(x+1, 3), x) == [-2, -1, 3]
 
+
 def test_issue_11553():
     eq1 = x + y + 1
     eq2 = x + GoldenRatio
     assert solve([eq1, eq2], x, y) == {x: -GoldenRatio, y: -1 + GoldenRatio}
     eq3 = x + 2 + TribonacciConstant
     assert solve([eq1, eq3], x, y) == {x: -2 - TribonacciConstant, y: 1 + TribonacciConstant}
+
+
+def test_issue_19113_19102():
+    t = S(1)/3
+    solve(cos(x)**5-sin(x)**5)
+    assert solve(4*cos(x)**3 - 2*sin(x)**3) == [
+        atan(2**(t)), -atan(2**(t)*(1 - sqrt(3)*I)/2),
+        -atan(2**(t)*(1 + sqrt(3)*I)/2)]
+    h = S.Half
+    assert solve(cos(x)**2 + sin(x)) == [
+        2*atan(-h + sqrt(5)/2 + sqrt(2)*sqrt(1 - sqrt(5))/2),
+        -2*atan(h + sqrt(5)/2 + sqrt(2)*sqrt(1 + sqrt(5))/2),
+        -2*atan(-sqrt(5)/2 + h + sqrt(2)*sqrt(1 - sqrt(5))/2),
+        -2*atan(-sqrt(2)*sqrt(1 + sqrt(5))/2 + h + sqrt(5)/2)]
+    assert solve(3*cos(x) - sin(x)) == [atan(3)]