Repository: sympy/sympy

Issue:
Expression blowup in matrix operations
Coming from #16823 I'm opening a new issue because this topic doesn't exactly deal with the previous. I'm going to send up a first quick ROUGH DRAFT example of the idea. The simplification is applied only to the power function currently but if it works the idea is apply it to all the matrix operations that blow up. The function is a copy of certain sections of the simplify function which seem effective at taming the blowup. I don't like the copy-and-paste reuse, maybe more flags could be added to the simplify function for more granular control of what simplifications get tried and then that is used instead of a dedicated matrix member.

Other issues: I tried the function in `rref()` as suggested but it actually came out worse than doing it at the end. The exact positioning of the simplification is still to be decided but within the rref seems to be too deep to be able to really cut down the expression, also it could be expensive there. Right now judging by A FEW test cases this function is 2-5x faster than a full simplify and in one case even significantly faster that `cancel()`, but that will depend on the content of the matrix and could vary widely.

Discussion:

Patch:
diff --git a/doc/src/tutorial/matrices.rst b/doc/src/tutorial/matrices.rst
index c09c02e6ed..4ad292efbb 100644
--- a/doc/src/tutorial/matrices.rst
+++ b/doc/src/tutorial/matrices.rst
@@ -412,6 +412,81 @@ expensive to calculate.
 Possible Issues
 ===============
 
+Controlling matrix expression blowup
+------------------------------------
+
+Normally a matrix ``*`` multiplication or ``**`` power will not simplify
+intermediate terms as they are calcualted to give the best general performance.
+Unfortunately this can lead to an explosion in the size and complexity of matrix
+element symbolic expressions (including complex number expressions), and even
+impair the ability of the calculation to finish, though in the majority of small
+matrices it is much faster. This can be a problem however for larger matrices
+and can even cause multiplication to never end or may be impossible to simplify
+afterwards.
+
+For this reason two separate functions - ``multiply`` and ``pow`` for
+multiplication and exponentiation have been added with a flag which allows you
+to specify that an optimized intermediate simplification step is to be performed
+at each step of calculation. The matrix ``exp`` function has also been modified
+to take this flag. By default these functions have the flag ``dotprodsimp`` set
+to None which means that intermediate simplification is not done, set this to
+``True`` to use the simplification. This will leave the matrix in a relatively
+simplified state after the operation and permit calculations on some large
+matrices which were not possible or extremely slow previously.
+
+This flag has also been added to the majority of matrix functions to enable
+this intermedate simplification in order to speed up calculation in many
+instances and return a simplified result. NOTE: It does not work in all cases,
+but so far testing has shown that it helps with the majority.
+
+    >>> x = Symbol('x')
+    >>> M = Matrix([[1+x, 1-x], [1-x, 1+x]])
+    >>> M*M
+    ⎡       2          2                     ⎤
+    ⎢(1 - x)  + (x + 1)    2⋅(1 - x)⋅(x + 1) ⎥
+    ⎢                                        ⎥
+    ⎢                            2          2⎥
+    ⎣ 2⋅(1 - x)⋅(x + 1)   (1 - x)  + (x + 1) ⎦
+    >>> M.multiply(M, dotprodsimp=True)
+    ⎡   2             2⎤
+    ⎢2⋅x  + 2  2 - 2⋅x ⎥
+    ⎢                  ⎥
+    ⎢       2     2    ⎥
+    ⎣2 - 2⋅x   2⋅x  + 2⎦
+    >>> M**2
+    ⎡       2          2                     ⎤
+    ⎢(1 - x)  + (x + 1)    2⋅(1 - x)⋅(x + 1) ⎥
+    ⎢                                        ⎥
+    ⎢                            2          2⎥
+    ⎣ 2⋅(1 - x)⋅(x + 1)   (1 - x)  + (x + 1) ⎦
+    >>> M.pow(8, dotprodsimp=True)
+    ⎡     8                   8⎤
+    ⎢128⋅x  + 128  128 - 128⋅x ⎥
+    ⎢                          ⎥
+    ⎢           8       8      ⎥
+    ⎣128 - 128⋅x   128⋅x  + 128⎦
+    >>> M.exp()
+    ⎡                               ⎛  1 - x      ⎞  2               ⎤
+    ⎢           2    2⋅x    (1 - x)⋅⎜───────── + 1⎟⋅ℯ             2⋅x⎥
+    ⎢  (1 - x)⋅ℯ    ℯ               ⎝2⋅(x - 1)    ⎠      (1 - x)⋅ℯ   ⎥
+    ⎢- ────────── + ────  - ────────────────────────── + ────────────⎥
+    ⎢  2⋅(x - 1)     2                x - 1               2⋅(x - 1)  ⎥
+    ⎢                                                                ⎥
+    ⎢       2⋅x    2                     2⋅x                         ⎥
+    ⎢      ℯ      ℯ             (1 - x)⋅ℯ      ⎛  1 - x      ⎞  2    ⎥
+    ⎢    - ──── + ──          - ──────────── + ⎜───────── + 1⎟⋅ℯ     ⎥
+    ⎣       2     2              2⋅(x - 1)     ⎝2⋅(x - 1)    ⎠       ⎦
+    >>> M.exp(dotprodsimp=True)
+    ⎡  2⋅x    2      2⋅x    2⎤
+    ⎢ ℯ      ℯ      ℯ      ℯ ⎥
+    ⎢ ──── + ──   - ──── + ──⎥
+    ⎢  2     2       2     2 ⎥
+    ⎢                        ⎥
+    ⎢   2⋅x    2    2⋅x    2 ⎥
+    ⎢  ℯ      ℯ    ℯ      ℯ  ⎥
+    ⎢- ──── + ──   ──── + ── ⎥
+    ⎣   2     2     2     2  ⎦
+
 Zero Testing
 ------------
 
diff --git a/sympy/matrices/common.py b/sympy/matrices/common.py
index faf27e2cb6..996a519378 100644
--- a/sympy/matrices/common.py
+++ b/sympy/matrices/common.py
@@ -10,15 +10,18 @@
 from inspect import isfunction
 
 from sympy.assumptions.refine import refine
+from sympy.core import SympifyError, Add
 from sympy.core.basic import Atom
 from sympy.core.compatibility import (
     Iterable, as_int, is_sequence, range, reduce)
 from sympy.core.decorators import call_highest_priority
+from sympy.core.expr import Expr
 from sympy.core.singleton import S
 from sympy.core.symbol import Symbol
 from sympy.core.sympify import sympify
 from sympy.functions import Abs
-from sympy.simplify import simplify as _simplify
+from sympy.polys import cancel, together
+from sympy.simplify import simplify as _simplify, dotprodsimp as _dotprodsimp
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 from sympy.utilities.iterables import flatten
 from sympy.utilities.misc import filldedent
@@ -2139,18 +2142,14 @@ def _eval_add(self, other):
 
     def _eval_matrix_mul(self, other):
         def entry(i, j):
+            vec = [self[i,k]*other[k,j] for k in range(self.cols)]
             try:
-                return sum(self[i,k]*other[k,j] for k in range(self.cols))
-            except TypeError:
-                # Block matrices don't work with `sum` or `Add` (ISSUE #11599)
+                return Add(*vec)
+            except (TypeError, SympifyError):
+                # Some matrices don't work with `sum` or `Add`
                 # They don't work with `sum` because `sum` tries to add `0`
-                # initially, and for a matrix, that is a mix of a scalar and
-                # a matrix, which raises a TypeError. Fall back to a
-                # block-matrix-safe way to multiply if the `sum` fails.
-                ret = self[i, 0]*other[0, j]
-                for k in range(1, self.cols):
-                    ret += self[i, k]*other[k, j]
-                return ret
+                # Fall back to a safe way to multiply if the `Add` fails.
+                return reduce(lambda a, b: a + b, vec)
 
         return self._new(self.rows, other.cols, entry)
 
@@ -2165,10 +2164,43 @@ def entry(i, j):
     def _eval_pow_by_recursion(self, num):
         if num == 1:
             return self
+
         if num % 2 == 1:
-            return self * self._eval_pow_by_recursion(num - 1)
-        ret = self._eval_pow_by_recursion(num // 2)
-        return ret * ret
+            a, b = self, self._eval_pow_by_recursion(num - 1)
+        else:
+            a = b = self._eval_pow_by_recursion(num // 2)
+
+        return a.multiply(b)
+
+    def _eval_pow_by_recursion_mulsimp(self, num, dotprodsimp=None, prevsimp=None):
+        if dotprodsimp and prevsimp is None:
+            prevsimp = [True]*len(self)
+
+        if num == 1:
+            return self
+
+        if num % 2 == 1:
+            a, b = self, self._eval_pow_by_recursion_mulsimp(num - 1,
+                    dotprodsimp=dotprodsimp, prevsimp=prevsimp)
+        else:
+            a = b = self._eval_pow_by_recursion_mulsimp(num // 2,
+                    dotprodsimp=dotprodsimp, prevsimp=prevsimp)
+
+        m = a.multiply(b, dotprodsimp=False)
+
+        if not dotprodsimp:
+            return m
+
+        lenm  = len(m)
+        elems = [None]*lenm
+
+        for i in range(lenm):
+            if prevsimp[i]:
+                elems[i], prevsimp[i] = _dotprodsimp(m[i], withsimp=True)
+            else:
+                elems[i] = m[i]
+
+        return m._new(m.rows, m.cols, elems)
 
     def _eval_scalar_mul(self, other):
         return self._new(self.rows, self.cols, lambda i, j: self[i,j]*other)
@@ -2252,6 +2284,21 @@ def __mul__(self, other):
 
         matrix_multiply_elementwise
         """
+
+        return self.multiply(other)
+
+    def multiply(self, other, dotprodsimp=None):
+        """Same as __mul__() but with optional simplification.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+        """
+
         other = _matrixify(other)
         # matrix-like objects can have shapes.  This is
         # our first sanity check.
@@ -2262,7 +2309,11 @@ def __mul__(self, other):
 
         # honest sympy matrices defer to their class's routine
         if getattr(other, 'is_Matrix', False):
-            return self._eval_matrix_mul(other)
+            m = self._eval_matrix_mul(other)
+            if dotprodsimp:
+                return m.applyfunc(_dotprodsimp)
+            return m
+
         # Matrix-like objects can be passed to CommonMatrix routines directly.
         if getattr(other, 'is_MatrixLike', False):
             return MatrixArithmetic._eval_matrix_mul(self, other)
@@ -2276,11 +2327,59 @@ def __mul__(self, other):
 
         return NotImplemented
 
+    def multiply_elementwise(self, other):
+        """Return the Hadamard product (elementwise product) of A and B
+
+        Examples
+        ========
+
+        >>> from sympy.matrices import Matrix
+        >>> A = Matrix([[0, 1, 2], [3, 4, 5]])
+        >>> B = Matrix([[1, 10, 100], [100, 10, 1]])
+        >>> A.multiply_elementwise(B)
+        Matrix([
+        [  0, 10, 200],
+        [300, 40,   5]])
+
+        See Also
+        ========
+
+        sympy.matrices.matrices.MatrixBase.cross
+        sympy.matrices.matrices.MatrixBase.dot
+        multiply
+        """
+        if self.shape != other.shape:
+            raise ShapeError("Matrix shapes must agree {} != {}".format(self.shape, other.shape))
+
+        return self._eval_matrix_mul_elementwise(other)
+
     def __neg__(self):
         return self._eval_scalar_mul(-1)
 
     @call_highest_priority('__rpow__')
     def __pow__(self, exp):
+        """Return self**exp a scalar or symbol."""
+
+        return self.pow(exp)
+
+    def pow(self, exp, dotprodsimp=None, jordan=None):
+        """Return self**exp a scalar or symbol.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
+        jordan : bool, optional
+            If left as None then Jordan form exponentiation will be used under
+            certain conditions, True specifies that jordan_pow should always
+            be used if possible and False means it should not be used unless
+            it is the only way to calculate the power.
+        """
+
         if self.rows != self.cols:
             raise NonSquareMatrixError()
         a = self
@@ -2305,16 +2404,22 @@ def __pow__(self, exp):
             # When certain conditions are met,
             # Jordan block algorithm is faster than
             # computation by recursion.
-            elif a.rows == 2 and exp > 100000 and jordan_pow is not None:
+            elif jordan_pow is not None and (jordan or \
+                    (jordan is not False and a.rows == 2 and exp > 100000)):
                 try:
-                    return jordan_pow(exp)
+                    return jordan_pow(exp, dotprodsimp=dotprodsimp)
                 except MatrixError:
-                    pass
-            return a._eval_pow_by_recursion(exp)
+                    if jordan:
+                        raise
+
+            if dotprodsimp is not None:
+                return a._eval_pow_by_recursion_mulsimp(exp, dotprodsimp=dotprodsimp)
+            else:
+                return a._eval_pow_by_recursion(exp)
 
         if jordan_pow:
             try:
-                return jordan_pow(exp)
+                return jordan_pow(exp, dotprodsimp=dotprodsimp)
             except NonInvertibleMatrixError:
                 # Raised by jordan_pow on zero determinant matrix unless exp is
                 # definitely known to be a non-negative integer.
@@ -2375,32 +2480,6 @@ def __sub__(self, a):
     def __truediv__(self, other):
         return self.__div__(other)
 
-    def multiply_elementwise(self, other):
-        """Return the Hadamard product (elementwise product) of A and B
-
-        Examples
-        ========
-
-        >>> from sympy.matrices import Matrix
-        >>> A = Matrix([[0, 1, 2], [3, 4, 5]])
-        >>> B = Matrix([[1, 10, 100], [100, 10, 1]])
-        >>> A.multiply_elementwise(B)
-        Matrix([
-        [  0, 10, 200],
-        [300, 40,   5]])
-
-        See Also
-        ========
-
-        sympy.matrices.matrices.MatrixBase.cross
-        sympy.matrices.matrices.MatrixBase.dot
-        multiply
-        """
-        if self.shape != other.shape:
-            raise ShapeError("Matrix shapes must agree {} != {}".format(self.shape, other.shape))
-
-        return self._eval_matrix_mul_elementwise(other)
-
 
 class MatrixCommon(MatrixArithmetic, MatrixOperations, MatrixProperties,
                   MatrixSpecial, MatrixShaping):
diff --git a/sympy/matrices/dense.py b/sympy/matrices/dense.py
index 6b42b7ab8e..018d60b3a4 100644
--- a/sympy/matrices/dense.py
+++ b/sympy/matrices/dense.py
@@ -2,11 +2,11 @@
 
 import random
 
-from sympy.core import SympifyError
+from sympy.core import SympifyError, Add
 from sympy.core.basic import Basic
 from sympy.core.compatibility import is_sequence, range, reduce
 from sympy.core.expr import Expr
-from sympy.core.function import expand_mul
+from sympy.core.function import count_ops, expand_mul
 from sympy.core.singleton import S
 from sympy.core.symbol import Symbol
 from sympy.core.sympify import sympify
@@ -172,39 +172,30 @@ def _eval_extract(self, rowsList, colsList):
                          list(mat[i] for i in indices), copy=False)
 
     def _eval_matrix_mul(self, other):
-        from sympy import Add
-        # cache attributes for faster access
-        self_cols = self.cols
-        other_rows, other_cols = other.rows, other.cols
-        other_len = other_rows * other_cols
-        new_mat_rows = self.rows
-        new_mat_cols = other.cols
-
-        # preallocate the array
-        new_mat = [self.zero]*new_mat_rows*new_mat_cols
+        other_len = other.rows*other.cols
+        new_len = self.rows*other.cols
+        new_mat = [self.zero]*new_len
 
         # if we multiply an n x 0 with a 0 x m, the
         # expected behavior is to produce an n x m matrix of zeros
         if self.cols != 0 and other.rows != 0:
-            # cache self._mat and other._mat for performance
+            self_cols = self.cols
             mat = self._mat
             other_mat = other._mat
-            for i in range(len(new_mat)):
-                row, col = i // new_mat_cols, i % new_mat_cols
+            for i in range(new_len):
+                row, col = i // other.cols, i % other.cols
                 row_indices = range(self_cols*row, self_cols*(row+1))
-                col_indices = range(col, other_len, other_cols)
-                vec = (mat[a]*other_mat[b] for a,b in zip(row_indices, col_indices))
+                col_indices = range(col, other_len, other.cols)
+                vec = [mat[a]*other_mat[b] for a, b in zip(row_indices, col_indices)]
                 try:
                     new_mat[i] = Add(*vec)
                 except (TypeError, SympifyError):
-                    # Block matrices don't work with `sum` or `Add` (ISSUE #11599)
+                    # Some matrices don't work with `sum` or `Add`
                     # They don't work with `sum` because `sum` tries to add `0`
-                    # initially, and for a matrix, that is a mix of a scalar and
-                    # a matrix, which raises a TypeError. Fall back to a
-                    # block-matrix-safe way to multiply if the `sum` fails.
-                    vec = (mat[a]*other_mat[b] for a,b in zip(row_indices, col_indices))
-                    new_mat[i] = reduce(lambda a,b: a + b, vec)
-        return classof(self, other)._new(new_mat_rows, new_mat_cols, new_mat, copy=False)
+                    # Fall back to a safe way to multiply if the `Add` fails.
+                    new_mat[i] = reduce(lambda a, b: a + b, vec)
+
+        return classof(self, other)._new(self.rows, other.cols, new_mat, copy=False)
 
     def _eval_matrix_mul_elementwise(self, other):
         mat = [a*b for a,b in zip(self._mat, other._mat)]
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index 270f5fa464..6f3bd2d2a9 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -24,7 +24,7 @@
 from sympy.polys import PurePoly, cancel, roots
 from sympy.printing import sstr
 from sympy.simplify import nsimplify
-from sympy.simplify import simplify as _simplify
+from sympy.simplify import simplify as _simplify, dotprodsimp as _dotprodsimp
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 from sympy.utilities.iterables import flatten, numbered_symbols
 from sympy.utilities.misc import filldedent
@@ -80,9 +80,18 @@ class MatrixDeterminant(MatrixCommon):
     """Provides basic matrix determinant operations.
     Should not be instantiated directly."""
 
-    def _eval_berkowitz_toeplitz_matrix(self):
+    def _eval_berkowitz_toeplitz_matrix(self, dotprodsimp=None):
         """Return (A,T) where T the Toeplitz matrix used in the Berkowitz algorithm
-        corresponding to ``self`` and A is the first principal submatrix."""
+        corresponding to ``self`` and A is the first principal submatrix.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+        """
 
         # the 0 x 0 case is trivial
         if self.rows == 0 and self.cols == 0:
@@ -112,8 +121,8 @@ def _eval_berkowitz_toeplitz_matrix(self):
         # compute -R * A**n * C.
         diags = [C]
         for i in range(self.rows - 2):
-            diags.append(A * diags[i])
-        diags = [(-R*d)[0, 0] for d in diags]
+            diags.append(A.multiply(diags[i], dotprodsimp=dotprodsimp))
+        diags = [(-R).multiply(d, dotprodsimp=dotprodsimp)[0, 0] for d in diags]
         diags = [self.one, -a] + diags
 
         def entry(i,j):
@@ -124,7 +133,7 @@ def entry(i,j):
         toeplitz = self._new(self.cols + 1, self.rows, entry)
         return (A, toeplitz)
 
-    def _eval_berkowitz_vector(self):
+    def _eval_berkowitz_vector(self, dotprodsimp=None):
         """ Run the Berkowitz algorithm and return a vector whose entries
             are the coefficients of the characteristic polynomial of ``self``.
 
@@ -148,6 +157,14 @@ def _eval_berkowitz_vector(self):
             As a consequence, all polynomials generated by Berkowitz
             algorithm are monic.
 
+            Parameters
+            ==========
+
+            dotprodsimp : bool, optional
+                Specifies whether intermediate term algebraic simplification is
+                used during matrix multiplications to control expression blowup
+                and thus speed up calculation.
+
            For more information on the implemented algorithm refer to:
 
            [1] S.J. Berkowitz, On computing the determinant in small
@@ -165,10 +182,12 @@ def _eval_berkowitz_vector(self):
         elif self.rows == 1 and self.cols == 1:
             return self._new(2, 1, [self.one, -self[0,0]])
 
-        submat, toeplitz = self._eval_berkowitz_toeplitz_matrix()
-        return toeplitz * submat._eval_berkowitz_vector()
+        submat, toeplitz = self._eval_berkowitz_toeplitz_matrix(dotprodsimp=dotprodsimp)
+        return toeplitz.multiply(submat._eval_berkowitz_vector(dotprodsimp=dotprodsimp),
+                dotprodsimp=dotprodsimp)
 
-    def _eval_det_bareiss(self, iszerofunc=_is_zero_after_expand_mul):
+    def _eval_det_bareiss(self, iszerofunc=_is_zero_after_expand_mul,
+            dotprodsimp=None):
         """Compute matrix determinant using Bareiss' fraction-free
         algorithm which is an extension of the well known Gaussian
         elimination method. This approach is best suited for dense
@@ -176,6 +195,14 @@ def _eval_det_bareiss(self, iszerofunc=_is_zero_after_expand_mul):
         minimal number of fractions. It means that less term
         rewriting is needed on resulting formulae.
 
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         TODO: Implement algorithm for sparse matrices (SFF),
         http://www.eecis.udel.edu/~saunders/papers/sffge/it5.ps.
         """
@@ -208,20 +235,32 @@ def bareiss(mat, cumm=1):
 
             def entry(i, j):
                 ret = (pivot_val*tmp_mat[i, j + 1] - mat[pivot_pos, j + 1]*tmp_mat[i, 0]) / cumm
-                if not ret.is_Atom:
+                if dotprodsimp:
+                    return _dotprodsimp(ret)
+                elif not ret.is_Atom:
                     return cancel(ret)
                 return ret
 
             return sign*bareiss(self._new(mat.rows - 1, mat.cols - 1, entry), pivot_val)
 
-        return cancel(bareiss(self))
+        return bareiss(self)
+
+    def _eval_det_berkowitz(self, dotprodsimp=None):
+        """ Use the Berkowitz algorithm to compute the determinant.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+        """
 
-    def _eval_det_berkowitz(self):
-        """ Use the Berkowitz algorithm to compute the determinant."""
-        berk_vector = self._eval_berkowitz_vector()
+        berk_vector = self._eval_berkowitz_vector(dotprodsimp=dotprodsimp)
         return (-1)**(len(berk_vector) - 1) * berk_vector[-1]
 
-    def _eval_det_lu(self, iszerofunc=_iszero, simpfunc=None):
+    def _eval_det_lu(self, iszerofunc=_iszero, simpfunc=None, dotprodsimp=None):
         """ Computes the determinant of a matrix from its LU decomposition.
         This function uses the LU decomposition computed by
         LUDecomposition_Simple().
@@ -234,7 +273,16 @@ def _eval_det_lu(self, iszerofunc=_iszero, simpfunc=None):
         The default is simpfunc=None, which indicate that the pivot search
         algorithm should not attempt to simplify any candidate pivots.
         If simpfunc fails to simplify its input, then it must return its input
-        instead of a copy."""
+        instead of a copy.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+        """
 
         if self.rows == 0:
             return self.one
@@ -242,7 +290,8 @@ def _eval_det_lu(self, iszerofunc=_iszero, simpfunc=None):
             # suggests that the determinant of a 0 x 0 matrix is one, by
             # convention.
 
-        lu, row_swaps = self.LUdecomposition_Simple(iszerofunc=iszerofunc, simpfunc=None)
+        lu, row_swaps = self.LUdecomposition_Simple(iszerofunc=iszerofunc, simpfunc=None,
+                dotprodsimp=dotprodsimp)
         # P*A = L*U => det(A) = det(L)*det(U)/det(P) = det(P)*det(U).
         # Lower triangular factor L encoded in lu has unit diagonal => det(L) = 1.
         # P is a permutation matrix => det(P) in {-1, 1} => 1/det(P) = det(P).
@@ -288,7 +337,7 @@ def adjugate(self, method="berkowitz"):
         """
         return self.cofactor_matrix(method).transpose()
 
-    def charpoly(self, x='lambda', simplify=_simplify):
+    def charpoly(self, x='lambda', simplify=_simplify, dotprodsimp=None):
         """Computes characteristic polynomial det(x*I - self) where I is
         the identity matrix.
 
@@ -326,6 +375,14 @@ def charpoly(self, x='lambda', simplify=_simplify):
         >>> A.charpoly(x).gen == x
         False
 
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            to control expression blowup during matrix multiplication. If this
+            is true then the simplify function is not used.
+
         Notes
         =====
 
@@ -343,13 +400,24 @@ def charpoly(self, x='lambda', simplify=_simplify):
         if not self.is_square:
             raise NonSquareMatrixError()
 
-        berk_vector = self._eval_berkowitz_vector()
+        if dotprodsimp:
+            simplify = lambda e: e
+
+        berk_vector = self._eval_berkowitz_vector(dotprodsimp=dotprodsimp)
         x = _uniquely_named_symbol(x, berk_vector)
         return PurePoly([simplify(a) for a in berk_vector], x)
 
-    def cofactor(self, i, j, method="berkowitz"):
+    def cofactor(self, i, j, method="berkowitz", dotprodsimp=None):
         """Calculate the cofactor of an element.
 
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            to control expression blowup during matrix multiplication. If this
+            is true then the simplify function is not used.
+
         See Also
         ========
 
@@ -361,11 +429,19 @@ def cofactor(self, i, j, method="berkowitz"):
         if not self.is_square or self.rows < 1:
             raise NonSquareMatrixError()
 
-        return (-1)**((i + j) % 2) * self.minor(i, j, method)
+        return (-1)**((i + j) % 2) * self.minor(i, j, method, dotprodsimp=dotprodsimp)
 
-    def cofactor_matrix(self, method="berkowitz"):
+    def cofactor_matrix(self, method="berkowitz", dotprodsimp=None):
         """Return a matrix containing the cofactor of each element.
 
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            to control expression blowup during matrix multiplication. If this
+            is true then the simplify function is not used.
+
         See Also
         ========
 
@@ -379,9 +455,9 @@ def cofactor_matrix(self, method="berkowitz"):
             raise NonSquareMatrixError()
 
         return self._new(self.rows, self.cols,
-                         lambda i, j: self.cofactor(i, j, method))
+                         lambda i, j: self.cofactor(i, j, method, dotprodsimp=dotprodsimp))
 
-    def det(self, method="bareiss", iszerofunc=None):
+    def det(self, method="bareiss", iszerofunc=None, dotprodsimp=None):
         """Computes the determinant of a matrix.
 
         Parameters
@@ -418,6 +494,11 @@ def det(self, method="bareiss", iszerofunc=None):
             and returns ``True`` if it is tested as zero and ``False`` if it
             tested as non-zero, and also ``None`` if it is undecidable.
 
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         Returns
         =======
 
@@ -463,27 +544,37 @@ def det(self, method="bareiss", iszerofunc=None):
         elif n == 1:
             return self[0,0]
         elif n == 2:
-            return self[0, 0] * self[1, 1] - self[0, 1] * self[1, 0]
+            m = self[0, 0] * self[1, 1] - self[0, 1] * self[1, 0]
+            return _dotprodsimp(m) if dotprodsimp else m
         elif n == 3:
-            return  (self[0, 0] * self[1, 1] * self[2, 2]
-                   + self[0, 1] * self[1, 2] * self[2, 0]
-                   + self[0, 2] * self[1, 0] * self[2, 1]
-                   - self[0, 2] * self[1, 1] * self[2, 0]
-                   - self[0, 0] * self[1, 2] * self[2, 1]
-                   - self[0, 1] * self[1, 0] * self[2, 2])
+            m =  (self[0, 0] * self[1, 1] * self[2, 2]
+                + self[0, 1] * self[1, 2] * self[2, 0]
+                + self[0, 2] * self[1, 0] * self[2, 1]
+                - self[0, 2] * self[1, 1] * self[2, 0]
+                - self[0, 0] * self[1, 2] * self[2, 1]
+                - self[0, 1] * self[1, 0] * self[2, 2])
+            return _dotprodsimp(m) if dotprodsimp else m
 
         if method == "bareiss":
-            return self._eval_det_bareiss(iszerofunc=iszerofunc)
+            return self._eval_det_bareiss(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
         elif method == "berkowitz":
-            return self._eval_det_berkowitz()
+            return self._eval_det_berkowitz(dotprodsimp=dotprodsimp)
         elif method == "lu":
-            return self._eval_det_lu(iszerofunc=iszerofunc)
+            return self._eval_det_lu(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
 
-    def minor(self, i, j, method="berkowitz"):
+    def minor(self, i, j, method="berkowitz", dotprodsimp=None):
         """Return the (i,j) minor of ``self``.  That is,
         return the determinant of the matrix obtained by deleting
         the `i`th row and `j`th column from ``self``.
 
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         See Also
         ========
 
@@ -495,7 +586,7 @@ def minor(self, i, j, method="berkowitz"):
         if not self.is_square or self.rows < 1:
             raise NonSquareMatrixError()
 
-        return self.minor_submatrix(i, j).det(method=method)
+        return self.minor_submatrix(i, j).det(method=method, dotprodsimp=dotprodsimp)
 
     def minor_submatrix(self, i, j):
         """Return the submatrix obtained by removing the `i`th row
@@ -572,13 +663,14 @@ def entry(i, j):
             return self[i, j]
         return self._new(self.rows, self.cols, entry)
 
-    def _eval_echelon_form(self, iszerofunc, simpfunc):
+    def _eval_echelon_form(self, iszerofunc, simpfunc, dotprodsimp=None):
         """Returns (mat, swaps) where ``mat`` is a row-equivalent matrix
         in echelon form and ``swaps`` is a list of row-swaps performed."""
         reduced, pivot_cols, swaps = self._row_reduce(iszerofunc, simpfunc,
                                                       normalize_last=True,
                                                       normalize=False,
-                                                      zero_above=False)
+                                                      zero_above=False,
+                                                      dotprodsimp=dotprodsimp)
         return reduced, pivot_cols, swaps
 
     def _eval_is_echelon(self, iszerofunc):
@@ -589,10 +681,11 @@ def _eval_is_echelon(self, iszerofunc):
             return zeros_below and self[:, 1:]._eval_is_echelon(iszerofunc)
         return zeros_below and self[1:, 1:]._eval_is_echelon(iszerofunc)
 
-    def _eval_rref(self, iszerofunc, simpfunc, normalize_last=True):
+    def _eval_rref(self, iszerofunc, simpfunc, normalize_last=True, dotprodsimp=None):
         reduced, pivot_cols, swaps = self._row_reduce(iszerofunc, simpfunc,
                                                       normalize_last, normalize=True,
-                                                      zero_above=True)
+                                                      zero_above=True,
+                                                      dotprodsimp=dotprodsimp)
         return reduced, pivot_cols
 
     def _normalize_op_args(self, op, col, k, col1, col2, error_str="col"):
@@ -667,7 +760,7 @@ def complexity(i):
         return (self.permute(perm, orientation='cols'), perm)
 
     def _row_reduce(self, iszerofunc, simpfunc, normalize_last=True,
-                    normalize=True, zero_above=True):
+                    normalize=True, zero_above=True, dotprodsimp=None):
         """Row reduce ``self`` and return a tuple (rref_matrix,
         pivot_cols, swaps) where pivot_cols are the pivot columns
         and swaps are any row swaps that were used in the process
@@ -677,17 +770,27 @@ def _row_reduce(self, iszerofunc, simpfunc, normalize_last=True,
         ==========
 
         iszerofunc : determines if an entry can be used as a pivot
+
         simpfunc : used to simplify elements and test if they are
             zero if ``iszerofunc`` returns `None`
+
         normalize_last : indicates where all row reduction should
             happen in a fraction-free manner and then the rows are
             normalized (so that the pivots are 1), or whether
             rows should be normalized along the way (like the naive
             row reduction algorithm)
+
         normalize : whether pivot rows should be normalized so that
             the pivot value is 1
+
         zero_above : whether entries above the pivot should be zeroed.
             If ``zero_above=False``, an echelon matrix will be returned.
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         """
         rows, cols = self.rows, self.cols
         mat = list(self)
@@ -702,7 +805,8 @@ def cross_cancel(a, i, b, j):
             """Does the row op row[i] = a*row[i] - b*row[j]"""
             q = (j - i)*cols
             for p in range(i*cols, (i + 1)*cols):
-                mat[p] = a*mat[p] - b*mat[p + q]
+                m = a*mat[p] - b*mat[p + q]
+                mat[p] = _dotprodsimp(m) if dotprodsimp else m
 
         piv_row, piv_col = 0, 0
         pivot_cols = []
@@ -734,7 +838,8 @@ def cross_cancel(a, i, b, j):
                 i, j = piv_row, piv_col
                 mat[i*cols + j] = self.one
                 for p in range(i*cols + j + 1, (i + 1)*cols):
-                    mat[p] = mat[p] / pivot_val
+                    m = mat[p] / pivot_val
+                    mat[p] = _dotprodsimp(m) if dotprodsimp else m
                 # after normalizing, the pivot value is 1
                 pivot_val = self.one
 
@@ -760,19 +865,31 @@ def cross_cancel(a, i, b, j):
                 pivot_val = mat[piv_i*cols + piv_j]
                 mat[piv_i*cols + piv_j] = self.one
                 for p in range(piv_i*cols + piv_j + 1, (piv_i + 1)*cols):
-                    mat[p] = mat[p] / pivot_val
+                    m = mat[p] / pivot_val
+                    mat[p] = _dotprodsimp(m) if dotprodsimp else m
 
         return self._new(self.rows, self.cols, mat), tuple(pivot_cols), tuple(swaps)
 
-    def echelon_form(self, iszerofunc=_iszero, simplify=False, with_pivots=False):
+    def echelon_form(self, iszerofunc=_iszero, simplify=False, with_pivots=False,
+            dotprodsimp=None):
         """Returns a matrix row-equivalent to ``self`` that is
         in echelon form.  Note that echelon form of a matrix
         is *not* unique, however, properties like the row
-        space and the null space are preserved."""
+        space and the null space are preserved.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+        """
         simpfunc = simplify if isinstance(
             simplify, FunctionType) else _simplify
 
-        mat, pivots, swaps = self._eval_echelon_form(iszerofunc, simpfunc)
+        mat, pivots, swaps = self._eval_echelon_form(iszerofunc, simpfunc,
+                dotprodsimp=dotprodsimp)
         if with_pivots:
             return mat, pivots
         return mat
@@ -884,7 +1001,8 @@ def rank(self, iszerofunc=_iszero, simplify=False):
         echelon_form, pivots, swaps = mat._eval_echelon_form(iszerofunc=iszerofunc, simpfunc=simpfunc)
         return len(pivots)
 
-    def rref(self, iszerofunc=_iszero, simplify=False, pivots=True, normalize_last=True):
+    def rref(self, iszerofunc=_iszero, simplify=False, pivots=True,
+            normalize_last=True, dotprodsimp=None):
         """Return reduced row-echelon form of matrix and indices of pivot vars.
 
         Parameters
@@ -893,13 +1011,16 @@ def rref(self, iszerofunc=_iszero, simplify=False, pivots=True, normalize_last=T
         iszerofunc : Function
             A function used for detecting whether an element can
             act as a pivot.  ``lambda x: x.is_zero`` is used by default.
+
         simplify : Function
             A function used to simplify elements when looking for a pivot.
             By default SymPy's ``simplify`` is used.
+
         pivots : True or False
             If ``True``, a tuple containing the row-reduced matrix and a tuple
             of pivot columns is returned.  If ``False`` just the row-reduced
             matrix is returned.
+
         normalize_last : True or False
             If ``True``, no pivots are normalized to `1` until after all
             entries above and below each pivot are zeroed.  This means the row
@@ -908,6 +1029,11 @@ def rref(self, iszerofunc=_iszero, simplify=False, pivots=True, normalize_last=T
             each pivot is normalized to be `1` before row operations are
             used to zero above and below the pivot.
 
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         Notes
         =====
 
@@ -935,12 +1061,12 @@ def rref(self, iszerofunc=_iszero, simplify=False, pivots=True, normalize_last=T
         >>> rref_pivots
         (0, 1)
         """
-        simpfunc = simplify if isinstance(
-            simplify, FunctionType) else _simplify
+        simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
 
         ret, pivot_cols = self._eval_rref(iszerofunc=iszerofunc,
                                           simpfunc=simpfunc,
-                                          normalize_last=normalize_last)
+                                          normalize_last=normalize_last,
+                                          dotprodsimp=dotprodsimp)
         if pivots:
             ret = (ret, pivot_cols)
         return ret
@@ -950,9 +1076,17 @@ class MatrixSubspaces(MatrixReductions):
     """Provides methods relating to the fundamental subspaces
     of a matrix.  Should not be instantiated directly."""
 
-    def columnspace(self, simplify=False):
+    def columnspace(self, simplify=False, dotprodsimp=None):
         """Returns a list of vectors (Matrix objects) that span columnspace of ``self``
 
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         Examples
         ========
 
@@ -978,13 +1112,22 @@ def columnspace(self, simplify=False):
         nullspace
         rowspace
         """
-        reduced, pivots = self.echelon_form(simplify=simplify, with_pivots=True)
+        reduced, pivots = self.echelon_form(simplify=simplify, with_pivots=True,
+                dotprodsimp=dotprodsimp)
 
         return [self.col(i) for i in pivots]
 
-    def nullspace(self, simplify=False, iszerofunc=_iszero):
+    def nullspace(self, simplify=False, iszerofunc=_iszero, dotprodsimp=None):
         """Returns list of vectors (Matrix objects) that span nullspace of ``self``
 
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         Examples
         ========
 
@@ -1008,7 +1151,8 @@ def nullspace(self, simplify=False, iszerofunc=_iszero):
         rowspace
         """
 
-        reduced, pivots = self.rref(iszerofunc=iszerofunc, simplify=simplify)
+        reduced, pivots = self.rref(iszerofunc=iszerofunc, simplify=simplify,
+                dotprodsimp=dotprodsimp)
 
         free_vars = [i for i in range(self.cols) if i not in pivots]
 
@@ -1024,10 +1168,20 @@ def nullspace(self, simplify=False, iszerofunc=_iszero):
 
         return [self._new(self.cols, 1, b) for b in basis]
 
-    def rowspace(self, simplify=False):
-        """Returns a list of vectors that span the row space of ``self``."""
+    def rowspace(self, simplify=False, dotprodsimp=None):
+        """Returns a list of vectors that span the row space of ``self``.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+        """
 
-        reduced, pivots = self.echelon_form(simplify=simplify, with_pivots=True)
+        reduced, pivots = self.echelon_form(simplify=simplify, with_pivots=True,
+                dotprodsimp=dotprodsimp)
 
         return [reduced.row(i) for i in range(len(pivots))]
 
@@ -1110,7 +1264,8 @@ class MatrixEigen(MatrixSubspaces):
     """Provides basic matrix eigenvalue/vector operations.
     Should not be instantiated directly."""
 
-    def diagonalize(self, reals_only=False, sort=False, normalize=False):
+    def diagonalize(self, reals_only=False, sort=False, normalize=False,
+            dotprodsimp=None):
         """
         Return (P, D), where D is diagonal and
 
@@ -1123,9 +1278,16 @@ def diagonalize(self, reals_only=False, sort=False, normalize=False):
 
         reals_only : bool. Whether to throw an error if complex numbers are need
                      to diagonalize. (Default: False)
+
         sort : bool. Sort the eigenvalues along the diagonal. (Default: False)
+
         normalize : bool. If True, normalize the columns of P. (Default: False)
 
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         Examples
         ========
 
@@ -1163,10 +1325,11 @@ def diagonalize(self, reals_only=False, sort=False, normalize=False):
         if not self.is_square:
             raise NonSquareMatrixError()
 
-        if not self.is_diagonalizable(reals_only=reals_only):
-            raise MatrixError("Matrix is not diagonalizable")
+        is_diagonalizable, eigenvecs = self.is_diagonalizable_with_eigen(
+                    reals_only=reals_only, dotprodsimp=dotprodsimp)
 
-        eigenvecs = self.eigenvects(simplify=True)
+        if not is_diagonalizable:
+            raise MatrixError("Matrix is not diagonalizable")
 
         if sort:
             eigenvecs = sorted(eigenvecs, key=default_sort_key)
@@ -1181,7 +1344,7 @@ def diagonalize(self, reals_only=False, sort=False, normalize=False):
 
         return self.hstack(*p_cols), self.diag(*diag)
 
-    def eigenvals(self, error_when_incomplete=True, **flags):
+    def eigenvals(self, error_when_incomplete=True, dotprodsimp=None, **flags):
         r"""Return eigenvalues using the Berkowitz agorithm to compute
         the characteristic polynomial.
 
@@ -1193,6 +1356,11 @@ def eigenvals(self, error_when_incomplete=True, **flags):
             eigenvalues are computed. This is caused by ``roots`` not returning
             a full list of eigenvalues.
 
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         simplify : bool or function, optional
             If it is set to ``True``, it attempts to return the most
             simplified form of expressions returned by applying default
@@ -1273,9 +1441,10 @@ def eigenvals(self, error_when_incomplete=True, **flags):
         else:
             flags.pop('simplify', None)  # pop unsupported flag
             if isinstance(simplify, FunctionType):
-                eigs = roots(mat.charpoly(x=Dummy('x'), simplify=simplify), **flags)
+                eigs = roots(mat.charpoly(x=Dummy('x'), simplify=simplify,
+                        dotprodsimp=dotprodsimp), **flags)
             else:
-                eigs = roots(mat.charpoly(x=Dummy('x')), **flags)
+                eigs = roots(mat.charpoly(x=Dummy('x'), dotprodsimp=dotprodsimp), **flags)
 
         # make sure the algebraic multiplicity sums to the
         # size of the matrix
@@ -1289,6 +1458,7 @@ def eigenvals(self, error_when_incomplete=True, **flags):
             return eigs
         if not isinstance(simplify, FunctionType):
             simplify = _simplify
+
         # With 'multiple' flag set true, simplify() will be mapped for the list
         # Otherwise, simplify() will be mapped for the keys of the dictionary
         if not multiple:
@@ -1296,7 +1466,8 @@ def eigenvals(self, error_when_incomplete=True, **flags):
         else:
             return [simplify(value) for value in eigs]
 
-    def eigenvects(self, error_when_incomplete=True, iszerofunc=_iszero, **flags):
+    def eigenvects(self, error_when_incomplete=True, iszerofunc=_iszero,
+            dotprodsimp=None, **flags):
         """Return list of triples (eigenval, multiplicity, eigenspace).
 
         Parameters
@@ -1317,6 +1488,11 @@ def eigenvects(self, error_when_incomplete=True, iszerofunc=_iszero, **flags):
             and returns ``True`` if it is tested as zero and ``False`` if it
             is tested as non-zero, and ``None`` if it is undecidable.
 
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         simplify : bool or function, optional
             If ``True``, ``as_content_primitive()`` will be used to tidy up
             normalization artifacts.
@@ -1372,12 +1548,12 @@ def eigenvects(self, error_when_incomplete=True, iszerofunc=_iszero, **flags):
         def eigenspace(eigenval):
             """Get a basis for the eigenspace for a particular eigenvalue"""
             m = mat - self.eye(mat.rows) * eigenval
-            ret = m.nullspace(iszerofunc=iszerofunc)
+            ret = m.nullspace(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
             # the nullspace for a real eigenvalue should be
             # non-trivial.  If we didn't find an eigenvector, try once
             # more a little harder
             if len(ret) == 0 and simplify:
-                ret = m.nullspace(iszerofunc=iszerofunc, simplify=True)
+                ret = m.nullspace(iszerofunc=iszerofunc, simplify=True, dotprodsimp=dotprodsimp)
             if len(ret) == 0:
                 raise NotImplementedError(
                         "Can't evaluate eigenvector for eigenvalue %s" % eigenval)
@@ -1385,6 +1561,7 @@ def eigenspace(eigenval):
 
         eigenvals = mat.eigenvals(rational=False,
                                   error_when_incomplete=error_when_incomplete,
+                                  dotprodsimp=dotprodsimp,
                                   **flags)
         ret = [(val, mult, eigenspace(val)) for val, mult in
                     sorted(eigenvals.items(), key=default_sort_key)]
@@ -1400,7 +1577,31 @@ def denom_clean(l):
             ret = [(val.evalf(chop=chop), mult, [v.evalf(chop=chop) for v in es]) for val, mult, es in ret]
         return ret
 
-    def is_diagonalizable(self, reals_only=False, **kwargs):
+    def is_diagonalizable_with_eigen(self, reals_only=False, dotprodsimp=None):
+        """See is_diagonalizable. This function returns the bool along with the
+        eigenvectors to avoid calculating them again in functions like
+        ``diagonalize``."""
+
+        if not self.is_square:
+            return False, []
+
+        eigenvecs = self.eigenvects(simplify=True, dotprodsimp=dotprodsimp)
+
+        if all(e.is_real for e in self) and self.is_symmetric():
+            # every real symmetric matrix is real diagonalizable
+            return True, eigenvecs
+
+        for val, mult, basis in eigenvecs:
+            # if we have a complex eigenvalue
+            if reals_only and not val.is_real:
+                return False, eigenvecs
+            # if the geometric multiplicity doesn't equal the algebraic
+            if mult != len(basis):
+                return False, eigenvecs
+
+        return True, eigenvecs
+
+    def is_diagonalizable(self, reals_only=False, dotprodsimp=None, **kwargs):
         """Returns true if a matrix is diagonalizable.
 
         Parameters
@@ -1409,6 +1610,11 @@ def is_diagonalizable(self, reals_only=False, **kwargs):
         reals_only : bool. If reals_only=True, determine whether the matrix can be
                      diagonalized without complex numbers. (Default: False)
 
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         kwargs
         ======
 
@@ -1450,6 +1656,7 @@ def is_diagonalizable(self, reals_only=False, **kwargs):
         is_diagonal
         diagonalize
         """
+
         if 'clear_cache' in kwargs:
             SymPyDeprecationWarning(
                 feature='clear_cache',
@@ -1463,24 +1670,8 @@ def is_diagonalizable(self, reals_only=False, **kwargs):
                 issue=15887
             ).warn()
 
-        if not self.is_square:
-            return False
-
-        if all(e.is_real for e in self) and self.is_symmetric():
-            # every real symmetric matrix is real diagonalizable
-            return True
-
-        eigenvecs = self.eigenvects(simplify=True)
-
-        ret = True
-        for val, mult, basis in eigenvecs:
-            # if we have a complex eigenvalue
-            if reals_only and not val.is_real:
-                ret = False
-            # if the geometric multiplicity doesn't equal the algebraic
-            if mult != len(basis):
-                ret = False
-        return ret
+        return self.is_diagonalizable_with_eigen(reals_only=reals_only,
+                dotprodsimp=dotprodsimp)[0]
 
     def _eval_is_positive_definite(self, method="eigen"):
         """Algorithm dump for computing positive-definiteness of a
@@ -2366,7 +2557,7 @@ def _diagonal_solve(self, rhs):
         return self._new(
             rhs.rows, rhs.cols, lambda i, j: rhs[i, j] / self[i, i])
 
-    def _matrix_pow_by_jordan_blocks(self, num):
+    def _matrix_pow_by_jordan_blocks(self, num, dotprodsimp=None):
         from sympy.matrices import diag, MutableMatrix
         from sympy import binomial
 
@@ -2397,7 +2588,8 @@ def jordan_cell_power(jc, n):
         jordan_cells = [MutableMatrix(j) for j in jordan_cells]
         for j in jordan_cells:
             jordan_cell_power(j, num)
-        return self._new(P*diag(*jordan_cells)*P.inv())
+        return self._new(P.multiply(diag(*jordan_cells), dotprodsimp=dotprodsimp).multiply(
+                P.inv(), dotprodsimp=dotprodsimp))
 
     def __repr__(self):
         return sstr(self)
@@ -3186,7 +3378,7 @@ def _eval_matrix_exp_jblock(self):
         from .sparsetools import banded
         return self.__class__(banded(size, bands))
 
-    def exp(self):
+    def exp(self, dotprodsimp=None):
         """Return the exponential of a square matrix
 
         Examples
@@ -3200,6 +3392,14 @@ def exp(self):
         Matrix([
         [    exp(I*t)/2 + exp(-I*t)/2, -I*exp(I*t)/2 + I*exp(-I*t)/2],
         [I*exp(I*t)/2 - I*exp(-I*t)/2,      exp(I*t)/2 + exp(-I*t)/2]])
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
         """
         if not self.is_square:
             raise NonSquareMatrixError(
@@ -3216,7 +3416,7 @@ def exp(self):
         from sympy import re
         eJ = diag(*blocks)
         # n = self.rows
-        ret = P * eJ * P.inv()
+        ret = P.multiply(eJ, dotprodsimp=dotprodsimp).multiply(P.inv(), dotprodsimp=dotprodsimp)
         if all(value.is_real for value in self.values()):
             return type(self)(re(ret))
         else:
@@ -3671,12 +3871,20 @@ def inv(self, method=None, **kwargs):
             kwargs['method'] = method
         return self._eval_inverse(**kwargs)
 
-    def is_nilpotent(self):
+    def is_nilpotent(self, dotprodsimp=None):
         """Checks if a matrix is nilpotent.
 
         A matrix B is nilpotent if for some integer k, B**k is
         a zero matrix.
 
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         Examples
         ========
 
@@ -3695,7 +3903,7 @@ def is_nilpotent(self):
             raise NonSquareMatrixError(
                 "Nilpotency is valid only for square matrices")
         x = _uniquely_named_symbol('x', self)
-        p = self.charpoly(x)
+        p = self.charpoly(x, dotprodsimp=dotprodsimp)
         if p.args[0] == x ** self.rows:
             return True
         return False
@@ -3958,7 +4166,8 @@ def entry_U(i, j):
     def LUdecomposition_Simple(self,
                                iszerofunc=_iszero,
                                simpfunc=None,
-                               rankcheck=False):
+                               rankcheck=False,
+                               dotprodsimp=None):
         """Compute an lu decomposition of m x n matrix A, where P*A = L*U
 
         * L is m x m lower triangular with unit diagonal
@@ -4030,6 +4239,14 @@ def LUdecomposition_Simple(self,
         Future versions of ``LUdecomposition_simple()`` may use
         ``_find_reasonable_pivot()``.
 
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
         See Also
         ========
 
@@ -4134,7 +4351,8 @@ def LUdecomposition_Simple(self,
                 # in sympy/matrices/tests/test_sparse.py.
                 # c = pivot_row + 1 if pivot_row == pivot_col else pivot_col
                 for c in range(start_col, lu.cols):
-                    lu[row, c] = lu[row, c] - lu[row, pivot_row]*lu[pivot_row, c]
+                    e = lu[row, c] - lu[row, pivot_row]*lu[pivot_row, c]
+                    lu[row, c] = _dotprodsimp(e) if dotprodsimp else e
 
             if pivot_row != pivot_col:
                 # matrix rank < min(num rows, num cols),
@@ -4266,18 +4484,6 @@ def LUsolve(self, rhs, iszerofunc=_iszero):
             b.row_op(i, lambda x, _: x / scale)
         return rhs.__class__(b)
 
-    def multiply(self, b):
-        """Returns ``self*b``
-
-        See Also
-        ========
-
-        dot
-        cross
-        multiply_elementwise
-        """
-        return self * b
-
     def normalized(self, iszerofunc=_iszero):
         """Return the normalized version of ``self``.
 
diff --git a/sympy/matrices/sparse.py b/sympy/matrices/sparse.py
index f02863d13a..19177a38a0 100644
--- a/sympy/matrices/sparse.py
+++ b/sympy/matrices/sparse.py
@@ -3,7 +3,9 @@
 import copy
 from collections import defaultdict
 
-from sympy.core.compatibility import Callable, as_int, is_sequence, range
+from sympy.core import SympifyError, Add
+from sympy.core.compatibility import Callable, as_int, is_sequence, range, \
+    reduce
 from sympy.core.containers import Dict
 from sympy.core.expr import Expr
 from sympy.core.singleton import S
@@ -461,8 +463,9 @@ def _eval_is_symmetric(self, simpfunc):
 
     def _eval_matrix_mul(self, other):
         """Fast multiplication exploiting the sparsity of the matrix."""
+
         if not isinstance(other, SparseMatrix):
-            return self*self._new(other)
+            other = self._new(other)
 
         # if we made it here, we're both sparse matrices
         # create quick lookups for rows and cols
@@ -480,8 +483,15 @@ def _eval_matrix_mul(self, other):
                 # these are the only things that need to be multiplied.
                 indices = set(col_lookup[col].keys()) & set(row_lookup[row].keys())
                 if indices:
-                    val = sum(row_lookup[row][k]*col_lookup[col][k] for k in indices)
-                    smat[row, col] = val
+                    vec = [row_lookup[row][k]*col_lookup[col][k] for k in indices]
+                    try:
+                        smat[row, col] = Add(*vec)
+                    except (TypeError, SympifyError):
+                        # Some matrices don't work with `sum` or `Add`
+                        # They don't work with `sum` because `sum` tries to add `0`
+                        # Fall back to a safe way to multiply if the `Add` fails.
+                        smat[row, col] = reduce(lambda a, b: a + b, vec)
+
         return self._new(self.rows, other.cols, smat)
 
     def _eval_row_insert(self, irow, other):
diff --git a/sympy/simplify/__init__.py b/sympy/simplify/__init__.py
index 15da10c9aa..7e5939cc81 100644
--- a/sympy/simplify/__init__.py
+++ b/sympy/simplify/__init__.py
@@ -6,7 +6,7 @@
 """
 from .simplify import (simplify, hypersimp, hypersimilar,
     logcombine, separatevars, posify, besselsimp, kroneckersimp,
-    signsimp, bottom_up, nsimplify)
+    signsimp, bottom_up, nsimplify, dotprodsimp)
 
 from .fu import FU, fu
 
diff --git a/sympy/simplify/simplify.py b/sympy/simplify/simplify.py
index 347e5f6332..ac551a89db 100644
--- a/sympy/simplify/simplify.py
+++ b/sympy/simplify/simplify.py
@@ -7,7 +7,8 @@
                         expand_power_exp, Eq)
 from sympy.core.compatibility import iterable, ordered, range, as_int
 from sympy.core.evaluate import global_evaluate
-from sympy.core.function import expand_log, count_ops, _mexpand, _coeff_isneg, nfloat
+from sympy.core.function import expand_log, count_ops, _mexpand, _coeff_isneg, \
+    nfloat, expand_mul, expand_multinomial
 from sympy.core.numbers import Float, I, pi, Rational, Integer
 from sympy.core.rules import Transform
 from sympy.core.sympify import _sympify
@@ -1960,3 +1961,157 @@ def compare(s, alt_s):
     else:
         simp = simp.doit(inv_expand=False)
     return simp
+
+
+def dotprodsimp(expr, withsimp=False):
+    """Simplification for a sum of products targeted at the kind of blowup that
+    occurs during summation of products. Intended to reduce expression blowup
+    during matrix multiplication or other similar operations.
+
+    Parameters
+    ==========
+
+    withsimp : bool, optional
+        Specifies whether a flag should be returned along with the expression
+        to indicate roughly whether simplification was successful. It is used
+        in ``MatrixArithmetic._eval_pow_by_recursion`` to avoid attempting to
+        simplify an expression repetitively which does not simplify.
+    """
+
+    from sympy.core.operations import LatticeOp
+
+    def _count_ops_alg(expr):
+        """Optimized count algebraic operations with no recursion into
+        non-algebraic args that ``core.function.count_ops`` does.
+        """
+
+        ops  = 0
+        args = [expr]
+
+        while args:
+            a = args.pop()
+
+            if a.is_Rational:
+                if a is not S.One: # -1/3 = NEG + DIV
+                    ops += bool (a.p < 0) + bool (a.q != 1)
+
+            elif a.is_Mul:
+                if _coeff_isneg(a):
+                    ops += 1
+                    if a.args[0] is S.NegativeOne:
+                        a = a.as_two_terms()[1]
+                    else:
+                        a = -a
+
+                n, d = fraction(a)
+
+                if n.is_Integer:
+                    ops += 1 + bool (n < 0)
+                    args.append(d) # won't be -Mul but could be Add
+
+                elif d is not S.One:
+                    if not d.is_Integer:
+                        args.append(d)
+                    ops += 1
+                    args.append(n) # could be -Mul
+
+                else:
+                    ops += len(a.args) - 1
+                    args.extend(a.args)
+
+            elif a.is_Add:
+                laargs = len(a.args)
+                negs   = 0
+
+                for ai in a.args:
+                    if _coeff_isneg(ai):
+                        negs += 1
+                        ai    = -ai
+                    args.append(ai)
+
+                ops += laargs - (negs != laargs) # -x - y = NEG + SUB
+
+            elif a.is_Pow:
+                ops += 1
+                args.append(a.base)
+
+            elif isinstance(a, LatticeOp):
+                ops += len(a.args) - 1
+                args.extend(a.args)
+
+        return ops
+
+    def _nonalg_subs_dummies(expr, dummies):
+        """Substitute dummy variables for non-algebraic expressions to avoid
+        evaluation of non-algebraic terms that ``polys.polytools.cancel`` does.
+        """
+
+        if not expr.args:
+            return expr
+
+        if expr.is_Add or expr.is_Mul or expr.is_Pow:
+            args = None
+
+            for i, a in enumerate(expr.args):
+                c = _nonalg_subs_dummies(a, dummies)
+
+                if c is a:
+                    continue
+
+                if args is None:
+                    args = list(expr.args)
+
+                args[i] = c
+
+            if args is None:
+                return expr
+
+            return expr.func(*args)
+
+        return dummies.setdefault(expr, Dummy())
+
+    simplified = False # doesn't really mean simplified, rather "can simplify again"
+
+    if isinstance(expr, Expr) and not expr.is_Relational:
+        # XXX: Should really count ops of initial expression and compare to
+        # expand_mul to see if simplification was done but that can be good bit
+        # slower for large initial expressions and this seems to work.
+        expr     = expand_mul(expr) # this and the following expand should not be combined
+        exprops  = _count_ops_alg(expr)
+        expr2    = expand_multinomial(expr)
+        expr2ops = _count_ops_alg(expr2)
+
+        if expr2ops < exprops:
+            expr       = expr2
+            exprops    = expr2ops
+            simplified = True
+
+        if exprops >= 6: # empirically tested cutoff for expensive simplification
+            dummies = {}
+            expr2   = _nonalg_subs_dummies(expr, dummies)
+
+            if expr2 is expr or _count_ops_alg(expr2) >= 6: # check again after substitution
+                expr3 = cancel(expr2)
+
+                if expr3 != expr2:
+                    expr       = expr3.subs([(d, e) for e, d in dummies.items()])
+                    simplified = True
+
+        # very special case: x/(x-1) - 1/(x-1) -> 1
+        elif (exprops == 5 and expr.is_Add and expr.args [0].is_Mul and
+                expr.args [1].is_Mul and expr.args [0].args [-1].is_Pow and
+                expr.args [1].args [-1].is_Pow and
+                expr.args [0].args [-1].exp is S.NegativeOne and
+                expr.args [1].args [-1].exp is S.NegativeOne):
+
+            expr2    = together (expr)
+            expr2ops = _count_ops_alg(expr2)
+
+            if expr2ops < exprops:
+                expr       = expr2
+                simplified = True
+
+        else:
+            simplified = True
+
+    return (expr, simplified) if withsimp else expr


Test patch:
diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py
index 3395999fe7..895aa6ccfa 100644
--- a/sympy/matrices/tests/test_matrices.py
+++ b/sympy/matrices/tests/test_matrices.py
@@ -285,29 +285,25 @@ def test_power():
     assert A**S.Half == Matrix ([[S.Half, S.Half], [3*S.Half, 3*S.Half]])
 
 
-@XFAIL
 def test_issue_17247_expression_blowup_1():
     M = Matrix([[1+x, 1-x], [1-x, 1+x]])
-    assert exp(M).expand() == Matrix([
+    assert M.exp(dotprodsimp=True).expand() == Matrix([
         [ (exp(2*x) + exp(2))/2, (-exp(2*x) + exp(2))/2],
         [(-exp(2*x) + exp(2))/2,  (exp(2*x) + exp(2))/2]])
 
-@XFAIL
 def test_issue_17247_expression_blowup_2():
     M = Matrix([[1+x, 1-x], [1-x, 1+x]])
     P, J = M.jordan_form ()
-    assert P*J*P.inv() == M
+    assert P.multiply(J, dotprodsimp=True).multiply(P.inv(), dotprodsimp=True)
 
-@XFAIL
 def test_issue_17247_expression_blowup_3():
     M = Matrix([[1+x, 1-x], [1-x, 1+x]])
-    assert M**100 == Matrix([
+    assert M.pow(100, dotprodsimp=True) == Matrix([
         [633825300114114700748351602688*x**100 + 633825300114114700748351602688, 633825300114114700748351602688 - 633825300114114700748351602688*x**100],
         [633825300114114700748351602688 - 633825300114114700748351602688*x**100, 633825300114114700748351602688*x**100 + 633825300114114700748351602688]])
 
-# This test commented out because it takes extremely long on current master,
-# it is here for testing when eventually matrix multiplication gets optimized.
-# def test_issue_17247_expression_blowup_4():
+def test_issue_17247_expression_blowup_4():
+# This matrix takes extremely long on current master even with intermediate simplification so an abbreviated version is used. It is left here for test in case of future optimizations.
 #     M = Matrix(S('''[
 #         [             -3/4,       45/32 - 37*I/16,         1/4 + I/2,      -129/64 - 9*I/64,      1/4 - 5*I/16,      65/128 + 87*I/64,         -9/32 - I/16,      183/256 - 97*I/128,       3/64 + 13*I/64,         -23/32 - 59*I/256,      15/128 - 3*I/32,        19/256 + 551*I/1024],
 #         [-149/64 + 49*I/32, -177/128 - 1369*I/128,  125/64 + 87*I/64, -2063/256 + 541*I/128,  85/256 - 33*I/16,  805/128 + 2415*I/512, -219/128 + 115*I/256, 6301/4096 - 6609*I/1024,  119/128 + 143*I/128, -10879/2048 + 4343*I/4096,  129/256 - 549*I/512, 42533/16384 + 29103*I/8192],
@@ -321,7 +317,7 @@ def test_issue_17247_expression_blowup_3():
 #         [      1/4 + 5*I/2,       -23/8 - 57*I/16,      1/4 + 13*I/4,    -825/64 - 147*I/32,          21/8 + I,    -537/64 + 143*I/16,       -5/8 - 39*I/16,     2473/256 + 137*I/64,    -149/64 + 49*I/32,     -177/128 - 1369*I/128,     125/64 + 87*I/64,      -2063/256 + 541*I/128],
 #         [               -4,               9 - 5*I,              -4*I,            27/2 + 6*I,                -2,         17/4 - 13*I/2,                1 + I,           -19/4 + 5*I/4,              1/2 - I,             9/4 + 55*I/16,                 -3/4,            45/32 - 37*I/16],
 #         [             -2*I,        119/8 + 29*I/4,       1/4 + 5*I/2,       -23/8 - 57*I/16,      1/4 + 13*I/4,    -825/64 - 147*I/32,             21/8 + I,      -537/64 + 143*I/16,       -5/8 - 39*I/16,       2473/256 + 137*I/64,    -149/64 + 49*I/32,      -177/128 - 1369*I/128]]'''))
-#     assert (M**10).expand() == Matrix([
+#     assert M.pow(10, dotprodsimp=True).expand() == Matrix([
 #         [    7*(-221393644768594642173548179825793834595 - 1861633166167425978847110897013541127952*I)/9671406556917033397649408,      15*(31670992489131684885307005100073928751695 + 10329090958303458811115024718207404523808*I)/77371252455336267181195264,   7*(-3710978679372178839237291049477017392703 + 1377706064483132637295566581525806894169*I)/19342813113834066795298816,            (9727707023582419994616144751727760051598 - 59261571067013123836477348473611225724433*I)/9671406556917033397649408,      (31896723509506857062605551443641668183707 + 54643444538699269118869436271152084599580*I)/38685626227668133590597632,       (-2024044860947539028275487595741003997397402 + 130959428791783397562960461903698670485863*I)/309485009821345068724781056,     3*(26190251453797590396533756519358368860907 - 27221191754180839338002754608545400941638*I)/77371252455336267181195264,      (1154643595139959842768960128434994698330461 + 3385496216250226964322872072260446072295634*I)/618970019642690137449562112,     3*(-31849347263064464698310044805285774295286 - 11877437776464148281991240541742691164309*I)/77371252455336267181195264,     (4661330392283532534549306589669150228040221 - 4171259766019818631067810706563064103956871*I)/1237940039285380274899124224,        (9598353794289061833850770474812760144506 + 358027153990999990968244906482319780943983*I)/309485009821345068724781056,     (-9755135335127734571547571921702373498554177 - 4837981372692695195747379349593041939686540*I)/2475880078570760549798248448],
 #         [(-379516731607474268954110071392894274962069 - 422272153179747548473724096872271700878296*I)/77371252455336267181195264, (41324748029613152354787280677832014263339501 - 12715121258662668420833935373453570749288074*I)/1237940039285380274899124224, (-339216903907423793947110742819264306542397 + 494174755147303922029979279454787373566517*I)/77371252455336267181195264, (-18121350839962855576667529908850640619878381 - 37413012454129786092962531597292531089199003*I)/1237940039285380274899124224, (2489661087330511608618880408199633556675926 + 1137821536550153872137379935240732287260863*I)/309485009821345068724781056, (-136644109701594123227587016790354220062972119 + 110130123468183660555391413889600443583585272*I)/4951760157141521099596496896, (1488043981274920070468141664150073426459593 - 9691968079933445130866371609614474474327650*I)/1237940039285380274899124224,  27*(4636797403026872518131756991410164760195942 + 3369103221138229204457272860484005850416533*I)/4951760157141521099596496896, (-8534279107365915284081669381642269800472363 + 2241118846262661434336333368511372725482742*I)/1237940039285380274899124224,  (60923350128174260992536531692058086830950875 - 263673488093551053385865699805250505661590126*I)/9903520314283042199192993792, (18520943561240714459282253753348921824172569 + 24846649186468656345966986622110971925703604*I)/4951760157141521099596496896,  (-232781130692604829085973604213529649638644431 + 35981505277760667933017117949103953338570617*I)/9903520314283042199192993792],
 #         [      (8742968295129404279528270438201520488950 + 3061473358639249112126847237482570858327*I)/4835703278458516698824704,      (-245657313712011778432792959787098074935273 + 253113767861878869678042729088355086740856*I)/38685626227668133590597632,      (1947031161734702327107371192008011621193 - 19462330079296259148177542369999791122762*I)/9671406556917033397649408,        (552856485625209001527688949522750288619217 + 392928441196156725372494335248099016686580*I)/77371252455336267181195264,      (-44542866621905323121630214897126343414629 + 3265340021421335059323962377647649632959*I)/19342813113834066795298816,          (136272594005759723105646069956434264218730 - 330975364731707309489523680957584684763587*I)/38685626227668133590597632,       (27392593965554149283318732469825168894401 + 75157071243800133880129376047131061115278*I)/38685626227668133590597632,      7*(-357821652913266734749960136017214096276154 - 45509144466378076475315751988405961498243*I)/309485009821345068724781056,       (104485001373574280824835174390219397141149 - 99041000529599568255829489765415726168162*I)/77371252455336267181195264,      (1198066993119982409323525798509037696321291 + 4249784165667887866939369628840569844519936*I)/618970019642690137449562112,       (-114985392587849953209115599084503853611014 - 52510376847189529234864487459476242883449*I)/77371252455336267181195264,      (6094620517051332877965959223269600650951573 - 4683469779240530439185019982269137976201163*I)/1237940039285380274899124224],
@@ -334,6 +330,141 @@ def test_issue_17247_expression_blowup_3():
 #         [    (27127513117071487872628354831658811211795 - 37765296987901990355760582016892124833857*I)/4835703278458516698824704,     (1741779916057680444272938534338833170625435 + 3083041729779495966997526404685535449810378*I)/77371252455336267181195264, 3*(-60642236251815783728374561836962709533401 - 24630301165439580049891518846174101510744*I)/19342813113834066795298816,      3*(445885207364591681637745678755008757483408 - 350948497734812895032502179455610024541643*I)/38685626227668133590597632,    (-47373295621391195484367368282471381775684 + 219122969294089357477027867028071400054973*I)/19342813113834066795298816,       (-2801565819673198722993348253876353741520438 - 2250142129822658548391697042460298703335701*I)/77371252455336267181195264,      (801448252275607253266997552356128790317119 - 50890367688077858227059515894356594900558*I)/77371252455336267181195264,    (-5082187758525931944557763799137987573501207 + 11610432359082071866576699236013484487676124*I)/309485009821345068724781056,     (-328925127096560623794883760398247685166830 - 643447969697471610060622160899409680422019*I)/77371252455336267181195264,    15*(2954944669454003684028194956846659916299765 + 33434406416888505837444969347824812608566*I)/1237940039285380274899124224,      (-415749104352001509942256567958449835766827 + 479330966144175743357171151440020955412219*I)/77371252455336267181195264,  3*(-4639987285852134369449873547637372282914255 - 11994411888966030153196659207284951579243273*I)/1237940039285380274899124224],
 #         [       (-478846096206269117345024348666145495601 + 1249092488629201351470551186322814883283*I)/302231454903657293676544,         (-17749319421930878799354766626365926894989 - 18264580106418628161818752318217357231971*I)/1208925819614629174706176,         (2801110795431528876849623279389579072819 + 363258850073786330770713557775566973248*I)/604462909807314587353088,          (-59053496693129013745775512127095650616252 + 78143588734197260279248498898321500167517*I)/4835703278458516698824704,         (-283186724922498212468162690097101115349 - 6443437753863179883794497936345437398276*I)/1208925819614629174706176,            (188799118826748909206887165661384998787543 + 84274736720556630026311383931055307398820*I)/9671406556917033397649408,         (-5482217151670072904078758141270295025989 + 1818284338672191024475557065444481298568*I)/1208925819614629174706176,          (56564463395350195513805521309731217952281 - 360208541416798112109946262159695452898431*I)/19342813113834066795298816,        11*(1259539805728870739006416869463689438068 + 1409136581547898074455004171305324917387*I)/4835703278458516698824704,       5*(-123701190701414554945251071190688818343325 + 30997157322590424677294553832111902279712*I)/38685626227668133590597632,          (16130917381301373033736295883982414239781 - 32752041297570919727145380131926943374516*I)/9671406556917033397649408,          (650301385108223834347093740500375498354925 + 899526407681131828596801223402866051809258*I)/77371252455336267181195264],
 #         [      (9011388245256140876590294262420614839483 + 8167917972423946282513000869327525382672*I)/1208925819614629174706176,       (-426393174084720190126376382194036323028924 + 180692224825757525982858693158209545430621*I)/9671406556917033397649408,     (24588556702197802674765733448108154175535 - 45091766022876486566421953254051868331066*I)/4835703278458516698824704,      (1872113939365285277373877183750416985089691 + 3030392393733212574744122057679633775773130*I)/77371252455336267181195264,    (-222173405538046189185754954524429864167549 - 75193157893478637039381059488387511299116*I)/19342813113834066795298816,        (2670821320766222522963689317316937579844558 - 2645837121493554383087981511645435472169191*I)/77371252455336267181195264,     5*(-2100110309556476773796963197283876204940 + 41957457246479840487980315496957337371937*I)/19342813113834066795298816,     (-5733743755499084165382383818991531258980593 - 3328949988392698205198574824396695027195732*I)/154742504910672534362390528,      (707827994365259025461378911159398206329247 - 265730616623227695108042528694302299777294*I)/77371252455336267181195264,    (-1442501604682933002895864804409322823788319 + 11504137805563265043376405214378288793343879*I)/309485009821345068724781056,         (-56130472299445561499538726459719629522285 - 61117552419727805035810982426639329818864*I)/9671406556917033397649408,    (39053692321126079849054272431599539429908717 - 10209127700342570953247177602860848130710666*I)/1237940039285380274899124224]])
+    M = Matrix(S('''[
+        [             -3/4,       45/32 - 37*I/16,         1/4 + I/2,      -129/64 - 9*I/64,      1/4 - 5*I/16,      65/128 + 87*I/64],
+        [-149/64 + 49*I/32, -177/128 - 1369*I/128,  125/64 + 87*I/64, -2063/256 + 541*I/128,  85/256 - 33*I/16,  805/128 + 2415*I/512],
+        [          1/2 - I,         9/4 + 55*I/16,              -3/4,       45/32 - 37*I/16,         1/4 + I/2,      -129/64 - 9*I/64],
+        [   -5/8 - 39*I/16,   2473/256 + 137*I/64, -149/64 + 49*I/32, -177/128 - 1369*I/128,  125/64 + 87*I/64, -2063/256 + 541*I/128],
+        [            1 + I,         -19/4 + 5*I/4,           1/2 - I,         9/4 + 55*I/16,              -3/4,       45/32 - 37*I/16],
+        [         21/8 + I,    -537/64 + 143*I/16,    -5/8 - 39*I/16,   2473/256 + 137*I/64, -149/64 + 49*I/32, -177/128 - 1369*I/128]]'''))
+    assert M.pow(10, dotprodsimp=True) == Matrix(S('''[
+        [     7369525394972778926719607798014571861/604462909807314587353088 - 229284202061790301477392339912557559*I/151115727451828646838272,   -19704281515163975949388435612632058035/1208925819614629174706176 + 14319858347987648723768698170712102887*I/302231454903657293676544,      -3623281909451783042932142262164941211/604462909807314587353088 - 6039240602494288615094338643452320495*I/604462909807314587353088,    109260497799140408739847239685705357695/2417851639229258349412352 - 7427566006564572463236368211555511431*I/2417851639229258349412352, -16095803767674394244695716092817006641/2417851639229258349412352 + 10336681897356760057393429626719177583*I/1208925819614629174706176,    -42207883340488041844332828574359769743/2417851639229258349412352 - 182332262671671273188016400290188468499*I/4835703278458516698824704],
+        [50566491050825573392726324995779608259/1208925819614629174706176 - 90047007594468146222002432884052362145*I/2417851639229258349412352,  74273703462900000967697427843983822011/1208925819614629174706176 + 265947522682943571171988741842776095421*I/1208925819614629174706176, -116900341394390200556829767923360888429/2417851639229258349412352 - 53153263356679268823910621474478756845*I/2417851639229258349412352, 195407378023867871243426523048612490249/1208925819614629174706176 - 1242417915995360200584837585002906728929*I/9671406556917033397649408,   -863597594389821970177319682495878193/302231454903657293676544 + 476936100741548328800725360758734300481*I/9671406556917033397649408, -3154451590535653853562472176601754835575/19342813113834066795298816 - 232909875490506237386836489998407329215*I/2417851639229258349412352],
+        [   -1715444997702484578716037230949868543/302231454903657293676544 + 5009695651321306866158517287924120777*I/302231454903657293676544,     -30551582497996879620371947949342101301/604462909807314587353088 - 7632518367986526187139161303331519629*I/151115727451828646838272,           312680739924495153190604170938220575/18889465931478580854784 - 108664334509328818765959789219208459*I/75557863725914323419136,    -14693696966703036206178521686918865509/604462909807314587353088 + 72345386220900843930147151999899692401*I/1208925819614629174706176,  -8218872496728882299722894680635296519/1208925819614629174706176 - 16776782833358893712645864791807664983*I/1208925819614629174706176,      143237839169380078671242929143670635137/2417851639229258349412352 + 2883817094806115974748882735218469447*I/2417851639229258349412352],
+        [   3087979417831061365023111800749855987/151115727451828646838272 + 34441942370802869368851419102423997089*I/604462909807314587353088, -148309181940158040917731426845476175667/604462909807314587353088 - 263987151804109387844966835369350904919*I/9671406556917033397649408,   50259518594816377378747711930008883165/1208925819614629174706176 - 95713974916869240305450001443767979653*I/2417851639229258349412352,  153466447023875527996457943521467271119/2417851639229258349412352 + 517285524891117105834922278517084871349*I/2417851639229258349412352,  -29184653615412989036678939366291205575/604462909807314587353088 - 27551322282526322041080173287022121083*I/1208925819614629174706176,   196404220110085511863671393922447671649/1208925819614629174706176 - 1204712019400186021982272049902206202145*I/9671406556917033397649408],
+        [     -2632581805949645784625606590600098779/151115727451828646838272 - 589957435912868015140272627522612771*I/37778931862957161709568,     26727850893953715274702844733506310247/302231454903657293676544 - 10825791956782128799168209600694020481*I/302231454903657293676544,      -1036348763702366164044671908440791295/151115727451828646838272 + 3188624571414467767868303105288107375*I/151115727451828646838272,     -36814959939970644875593411585393242449/604462909807314587353088 - 18457555789119782404850043842902832647*I/302231454903657293676544,      12454491297984637815063964572803058647/604462909807314587353088 - 340489532842249733975074349495329171*I/302231454903657293676544,      -19547211751145597258386735573258916681/604462909807314587353088 + 87299583775782199663414539883938008933*I/1208925819614629174706176],
+        [  -40281994229560039213253423262678393183/604462909807314587353088 - 2939986850065527327299273003299736641*I/604462909807314587353088, 331940684638052085845743020267462794181/2417851639229258349412352 - 284574901963624403933361315517248458969*I/1208925819614629174706176,      6453843623051745485064693628073010961/302231454903657293676544 + 36062454107479732681350914931391590957*I/604462909807314587353088,  -147665869053634695632880753646441962067/604462909807314587353088 - 305987938660447291246597544085345123927*I/9671406556917033397649408,  107821369195275772166593879711259469423/2417851639229258349412352 - 11645185518211204108659001435013326687*I/302231454903657293676544,     64121228424717666402009446088588091619/1208925819614629174706176 + 265557133337095047883844369272389762133*I/1208925819614629174706176]]'''))
+
+def test_issue_17247_expression_blowup_5():
+    M = Matrix(6, 6, lambda i, j: 1 + (-1)**(i+j)*I)
+    assert M.charpoly('x', dotprodsimp=True) == PurePoly(x**6 + (-6 - 6*I)*x**5 + 36*I*x**4, x, domain='EX')
+
+def test_issue_17247_expression_blowup_6():
+    M = Matrix(8, 8, [x+i for i in range (64)])
+    assert M.det('bareiss', dotprodsimp=True) == 0
+
+def test_issue_17247_expression_blowup_7():
+    M = Matrix(6, 6, lambda i, j: 1 + (-1)**(i+j)*I)
+    assert M.det('berkowitz', dotprodsimp=True) == 0
+
+def test_issue_17247_expression_blowup_8():
+    M = Matrix(8, 8, [x+i for i in range (64)])
+    assert M.det('lu', dotprodsimp=True) == 0
+
+def test_issue_17247_expression_blowup_9():
+    M = Matrix(8, 8, [x+i for i in range (64)])
+    assert M.rref(dotprodsimp=True) == (Matrix([
+        [1, 0, -1, -2, -3, -4, -5, -6],
+        [0, 1,  2,  3,  4,  5,  6,  7],
+        [0, 0,  0,  0,  0,  0,  0,  0],
+        [0, 0,  0,  0,  0,  0,  0,  0],
+        [0, 0,  0,  0,  0,  0,  0,  0],
+        [0, 0,  0,  0,  0,  0,  0,  0],
+        [0, 0,  0,  0,  0,  0,  0,  0],
+        [0, 0,  0,  0,  0,  0,  0,  0]]), (0, 1))
+
+def test_issue_17247_expression_blowup_10():
+    M = Matrix(6, 6, lambda i, j: 1 + (-1)**(i+j)*I)
+    assert M.cofactor(0, 0, dotprodsimp=True) == 0
+
+def test_issue_17247_expression_blowup_11():
+    M = Matrix(6, 6, lambda i, j: 1 + (-1)**(i+j)*I)
+    assert M.cofactor_matrix(dotprodsimp=True) == Matrix(6, 6, [0]*36)
+
+def test_issue_17247_expression_blowup_12():
+    M = Matrix(6, 6, lambda i, j: 1 + (-1)**(i+j)*I)
+    assert M.eigenvals(dotprodsimp=True) == {6: 1, 6*I: 1, 0: 4}
+
+def test_issue_17247_expression_blowup_13():
+    M = Matrix([
+        [    0, 1 - x, x + 1, 1 - x],
+        [1 - x, x + 1,     0, x + 1],
+        [    0, 1 - x, x + 1, 1 - x],
+        [    0,     0,     1 - x, 0]])
+    ev = M.eigenvects(dotprodsimp=True)
+    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]])
+
+def test_issue_17247_expression_blowup_14():
+    M = Matrix(8, 8, ([1+x, 1-x]*4 + [1-x, 1+x]*4)*4)
+    assert M.echelon_form(dotprodsimp=True) == Matrix([
+        [x + 1, 1 - x, x + 1, 1 - x, x + 1, 1 - x, x + 1, 1 - x],
+        [    0,   4*x,     0,   4*x,     0,   4*x,     0,   4*x],
+        [    0,     0,     0,     0,     0,     0,     0,     0],
+        [    0,     0,     0,     0,     0,     0,     0,     0],
+        [    0,     0,     0,     0,     0,     0,     0,     0],
+        [    0,     0,     0,     0,     0,     0,     0,     0],
+        [    0,     0,     0,     0,     0,     0,     0,     0],
+        [    0,     0,     0,     0,     0,     0,     0,     0]])
+
+def test_issue_17247_expression_blowup_15():
+    M = Matrix(8, 8, ([1+x, 1-x]*4 + [1-x, 1+x]*4)*4)
+    assert M.rowspace(dotprodsimp=True) == [Matrix([[x + 1, 1 - x, x + 1, 1 - x, x + 1, 1 - x, x + 1, 1 - x]]), Matrix([[0, 4*x, 0, 4*x, 0, 4*x, 0, 4*x]])]
+
+def test_issue_17247_expression_blowup_16():
+    M = Matrix(8, 8, ([1+x, 1-x]*4 + [1-x, 1+x]*4)*4)
+    assert M.columnspace(dotprodsimp=True) == [Matrix([[x + 1],[1 - x],[x + 1],[1 - x],[x + 1],[1 - x],[x + 1],[1 - x]]), Matrix([[1 - x],[x + 1],[1 - x],[x + 1],[1 - x],[x + 1],[1 - x],[x + 1]])]
+
+def test_issue_17247_expression_blowup_17():
+    M = Matrix(8, 8, [x+i for i in range (64)])
+    assert M.nullspace(dotprodsimp=True) == [
+        Matrix([[1],[-2],[1],[0],[0],[0],[0],[0]]),
+        Matrix([[2],[-3],[0],[1],[0],[0],[0],[0]]),
+        Matrix([[3],[-4],[0],[0],[1],[0],[0],[0]]),
+        Matrix([[4],[-5],[0],[0],[0],[1],[0],[0]]),
+        Matrix([[5],[-6],[0],[0],[0],[0],[1],[0]]),
+        Matrix([[6],[-7],[0],[0],[0],[0],[0],[1]])]
+
+def test_issue_17247_expression_blowup_18():
+    M = Matrix(6, 6, ([1+x, 1-x]*3 + [1-x, 1+x]*3)*3)
+    assert not M.is_nilpotent(dotprodsimp=True)
+
+def test_issue_17247_expression_blowup_19():
+    M = Matrix(S('''[
+        [             -3/4,                     0,         1/4 + I/2,                     0],
+        [                0, -177/128 - 1369*I/128,                 0, -2063/256 + 541*I/128],
+        [          1/2 - I,                     0,                 0,                     0],
+        [                0,                     0,                 0, -177/128 - 1369*I/128]]'''))
+    assert not M.is_diagonalizable(dotprodsimp=True)
+
+def test_issue_17247_expression_blowup_20():
+    M = Matrix([
+    [x + 1,  1 - x,      0,      0],
+    [1 - x,  x + 1,      0,  x + 1],
+    [    0,  1 - x,  x + 1,      0],
+    [    0,      0,      0,  x + 1]])
+    assert M.diagonalize(dotprodsimp=True) == (Matrix([
+        [1,  1, 0, (x + 1)/(x - 1)],
+        [1, -1, 0,               0],
+        [1,  1, 1,               0],
+        [0,  0, 0,               1]]),
+        Matrix([
+        [2,   0,     0,     0],
+        [0, 2*x,     0,     0],
+        [0,   0, x + 1,     0],
+        [0,   0,     0, x + 1]]))
 
 
 def test_creation():