Repository: sympy/sympy

Issue:
Matrix power is not implemented correctly
Matrix __pow__() operation does not seem to be implemented very rigorously leading to some erroneous cases. This issue has evolved from #17179 and a fix has been suggested by @oscarbenjamin to make matrix power work as such:

1. For `M**n` where n is a non-negative integer the result is always well defined. It can be computed from the Jordan form including for symbolic n (that satisfies the assumptions). The trivial case n=0 always gives the identity.

2. If n is a negative integer then the matrix needs to be invertible but otherwise the result is as in case 1.

3. If n is an integer, not known to be positive or negative then it can evaluate if the matrix is invertible but should remain unevaluated otherwise. Substitution of a negative integer for n should raise. Substitution of a nonnegative integer should evaluate.

4. If n is rational p/q then this is the pth power of the qth root of M (or the qth root of the pth power?). As for integers negative n is only allowed for invertible M. The qth root is not unique and does not always exist e.g. `[[0, 1], [0, 0]]` has no square root so this can raise. Where the roots exist there can be up to q**N possibilities if M is NxN.

4a. If M is positive definite then I think it has a unique positive definite qth root so we use that.

4b. If M has negative or non-real eigenvalues then provided it is not defective we can just take the pth root of each eigenvalue (taking SymPy's definition of the pth root for scalars).

4c. If M is defective then the Jordan form method can still work using pth roots of the eigenvalues.

5. If n is irrational or non-real complex then `M**n = exp(n*log(M))` - `log(M)` exists iff M is invertible and is unique if M is positive definite. If M is not positive definite then different choices to log(M) will not give equivalent results so raise? This can be taken as the definition of correctness although it might still be easiest to compute via the Jordan form method.

6. All of the above assumes that the entries of M are complex-valued. If M has matrices or operators or something as entries then only integer powers can be defined? The result could depend on whether the elements are invertible...

The problem is that some tests as well as probably some code out in the wild is coded to the current particular implementation and will break when we correct the power function.


Discussion:

Patch:
diff --git a/doc/src/tutorial/matrices.rst b/doc/src/tutorial/matrices.rst
index 2c88619adc..cf5dca5f3f 100644
--- a/doc/src/tutorial/matrices.rst
+++ b/doc/src/tutorial/matrices.rst
@@ -151,7 +151,7 @@ raise it to the ``-1`` power.
     >>> N**-1
     Traceback (most recent call last):
     ...
-    ValueError: Matrix det == 0; not invertible.
+    NonInvertibleMatrixError: Matrix det == 0; not invertible.
 
 To take the transpose of a Matrix, use ``T``.
 
diff --git a/sympy/functions/combinatorial/factorials.py b/sympy/functions/combinatorial/factorials.py
index 07e794abc8..8079184824 100644
--- a/sympy/functions/combinatorial/factorials.py
+++ b/sympy/functions/combinatorial/factorials.py
@@ -579,10 +579,6 @@ def eval(cls, x, k):
                                             r*(x - i),
                                             range(1, abs(int(k)) + 1), 1)
 
-        if k.is_integer == False:
-            if x.is_integer and x.is_negative:
-                return S.Zero
-
     def _eval_rewrite_as_gamma(self, x, k, **kwargs):
         from sympy import gamma
         return gamma(x + k) / gamma(x)
diff --git a/sympy/matrices/common.py b/sympy/matrices/common.py
index c5035c9080..881dc8a532 100644
--- a/sympy/matrices/common.py
+++ b/sympy/matrices/common.py
@@ -39,6 +39,11 @@ class NonSquareMatrixError(ShapeError):
     pass
 
 
+class NonInvertibleMatrixError(ValueError, MatrixError):
+    """The matrix in not invertible (division by multidimensional zero error)."""
+    pass
+
+
 class NonPositiveDefiniteMatrixError(ValueError, MatrixError):
     """The matrix is not a positive-definite matrix."""
     pass
@@ -2218,37 +2223,51 @@ def __neg__(self):
         return self._eval_scalar_mul(-1)
 
     @call_highest_priority('__rpow__')
-    def __pow__(self, num):
+    def __pow__(self, exp):
         if self.rows != self.cols:
             raise NonSquareMatrixError()
         a = self
         jordan_pow = getattr(a, '_matrix_pow_by_jordan_blocks', None)
-        num = sympify(num)
-        if num.is_Number and num % 1 == 0:
+        exp = sympify(exp)
+
+        if exp.is_zero:
+            return a._new(a.rows, a.cols, lambda i, j: int(i == j))
+        if exp == 1:
+            return a
+
+        diagonal = getattr(a, 'is_diagonal', None)
+        if diagonal is not None and diagonal():
+            return a._new(a.rows, a.cols, lambda i, j: a[i,j]**exp if i == j else 0)
+
+        if exp.is_Number and exp % 1 == 0:
             if a.rows == 1:
-                return a._new([[a[0]**num]])
-            if num == 0:
-                return self._new(self.rows, self.cols, lambda i, j: int(i == j))
-            if num < 0:
-                num = -num
+                return a._new([[a[0]**exp]])
+            if exp < 0:
+                exp = -exp
                 a = a.inv()
             # When certain conditions are met,
             # Jordan block algorithm is faster than
             # computation by recursion.
-            elif a.rows == 2 and num > 100000 and jordan_pow is not None:
+            elif a.rows == 2 and exp > 100000 and jordan_pow is not None:
                 try:
-                    return jordan_pow(num)
+                    return jordan_pow(exp)
                 except MatrixError:
                     pass
-            return a._eval_pow_by_recursion(num)
-        elif not num.is_Number and num.is_negative is None and a.det() == 0:
-            from sympy.matrices.expressions import MatPow
-            return MatPow(a, num)
-        elif isinstance(num, (Expr, float)):
-            return jordan_pow(num)
-        else:
-            raise TypeError(
-                "Only SymPy expressions or integers are supported as exponent for matrices")
+            return a._eval_pow_by_recursion(exp)
+
+        if jordan_pow:
+            try:
+                return jordan_pow(exp)
+            except NonInvertibleMatrixError:
+                # Raised by jordan_pow on zero determinant matrix unless exp is
+                # definitely known to be a non-negative integer.
+                # Here we raise if n is definitely not a non-negative integer
+                # but otherwise we can leave this as an unevaluated MatPow.
+                if exp.is_integer is False or exp.is_nonnegative is False:
+                    raise
+
+        from sympy.matrices.expressions import MatPow
+        return MatPow(a, exp)
 
     @call_highest_priority('__add__')
     def __radd__(self, other):
diff --git a/sympy/matrices/expressions/matpow.py b/sympy/matrices/expressions/matpow.py
index d079314946..d34ddb13bf 100644
--- a/sympy/matrices/expressions/matpow.py
+++ b/sympy/matrices/expressions/matpow.py
@@ -71,7 +71,7 @@ def doit(self, **kwargs):
             raise ValueError("Matrix determinant is 0, not invertible.")
         elif isinstance(base, (Identity, ZeroMatrix)):
             return base
-        elif isinstance(base, MatrixBase) and exp.is_number:
+        elif isinstance(base, MatrixBase):
             if exp is S.One:
                 return base
             return base**exp
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index f7cb139ef0..62a0a54877 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -20,6 +20,7 @@
 from sympy.core.sympify import sympify
 from sympy.functions import exp, factorial
 from sympy.functions.elementary.miscellaneous import Max, Min, sqrt
+from sympy.functions.special.tensor_functions import KroneckerDelta
 from sympy.polys import PurePoly, cancel, roots
 from sympy.printing import sstr
 from sympy.simplify import nsimplify
@@ -29,8 +30,8 @@
 from sympy.utilities.misc import filldedent
 
 from .common import (
-    MatrixCommon, MatrixError, NonSquareMatrixError, ShapeError,
-    NonPositiveDefiniteMatrixError)
+    MatrixCommon, MatrixError, NonSquareMatrixError, NonInvertibleMatrixError,
+    ShapeError, NonPositiveDefiniteMatrixError)
 
 
 def _iszero(x):
@@ -2362,17 +2363,24 @@ def _matrix_pow_by_jordan_blocks(self, num):
 
         def jordan_cell_power(jc, n):
             N = jc.shape[0]
-            l = jc[0, 0]
-            if l == 0 and (n < N - 1) != False:
-                raise ValueError("Matrix det == 0; not invertible")
-            elif l == 0 and N > 1 and n % 1 != 0:
-                raise ValueError("Non-integer power cannot be evaluated")
-            for i in range(N):
-                for j in range(N-i):
+            l = jc[0,0]
+            if l.is_zero:
+                if N == 1 and n.is_nonnegative:
+                    jc[0,0] = l**n
+                elif not (n.is_integer and n.is_nonnegative):
+                    raise NonInvertibleMatrixError("Non-invertible matrix can only be raised to a nonnegative integer")
+                else:
+                    for i in range(N):
+                        jc[0,i] = KroneckerDelta(i, n)
+            else:
+                for i in range(N):
                     bn = binomial(n, i)
                     if isinstance(bn, binomial):
                         bn = bn._eval_expand_func()
-                    jc[j, i+j] = l**(n-i)*bn
+                    jc[0,i] = l**(n-i)*bn
+            for i in range(N):
+                for j in range(1, N-i):
+                    jc[j,i+j] = jc [j-1,i+j-1]
 
         P, J = self.jordan_form()
         jordan_cells = J.get_diag_blocks()
@@ -3364,7 +3372,7 @@ def inv_mod(self, m):
         try:
             det_inv = mod_inverse(det_K, m)
         except ValueError:
-            raise ValueError('Matrix is not invertible (mod %d)' % m)
+            raise NonInvertibleMatrixError('Matrix is not invertible (mod %d)' % m)
 
         K_adj = self.adjugate()
         K_inv = self.__class__(N, N,
@@ -3392,7 +3400,7 @@ def inverse_ADJ(self, iszerofunc=_iszero):
             ok = self.rref(simplify=True)[0]
             zero = any(iszerofunc(ok[j, j]) for j in range(ok.rows))
         if zero:
-            raise ValueError("Matrix det == 0; not invertible.")
+            raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
 
         return self.adjugate() / d
 
@@ -3413,7 +3421,7 @@ def inverse_GE(self, iszerofunc=_iszero):
         big = Matrix.hstack(self.as_mutable(), Matrix.eye(self.rows))
         red = big.rref(iszerofunc=iszerofunc, simplify=True)[0]
         if any(iszerofunc(red[j, j]) for j in range(red.rows)):
-            raise ValueError("Matrix det == 0; not invertible.")
+            raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
 
         return self._new(red[:, big.rows:])
 
@@ -3432,7 +3440,7 @@ def inverse_LU(self, iszerofunc=_iszero):
 
         ok = self.rref(simplify=True)[0]
         if any(iszerofunc(ok[j, j]) for j in range(ok.rows)):
-            raise ValueError("Matrix det == 0; not invertible.")
+            raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
 
         return self.LUsolve(self.eye(self.rows), iszerofunc=_iszero)
 
@@ -4876,7 +4884,7 @@ def solve(self, rhs, method='GJ'):
             try:
                 soln, param = self.gauss_jordan_solve(rhs)
                 if param:
-                    raise ValueError("Matrix det == 0; not invertible. "
+                    raise NonInvertibleMatrixError("Matrix det == 0; not invertible. "
                     "Try ``self.gauss_jordan_solve(rhs)`` to obtain a parametric solution.")
             except ValueError:
                 # raise same error as in inv:


Test patch:
diff --git a/sympy/functions/combinatorial/tests/test_comb_factorials.py b/sympy/functions/combinatorial/tests/test_comb_factorials.py
index 49e5d84b3a..859474ee86 100644
--- a/sympy/functions/combinatorial/tests/test_comb_factorials.py
+++ b/sympy/functions/combinatorial/tests/test_comb_factorials.py
@@ -1,7 +1,7 @@
 from sympy import (S, Symbol, symbols, factorial, factorial2, Float, binomial,
                    rf, ff, gamma, polygamma, EulerGamma, O, pi, nan,
                    oo, zoo, simplify, expand_func, Product, Mul, Piecewise,
-                   Mod, Eq, sqrt, Poly, Dummy, I)
+                   Mod, Eq, sqrt, Poly, Dummy)
 from sympy.core.expr import unchanged
 from sympy.core.function import ArgumentIndexError
 from sympy.functions.combinatorial.factorials import subfactorial
@@ -30,14 +30,6 @@ def test_rf_eval_apply():
     assert rf(oo, -6) == oo
     assert rf(-oo, -7) == oo
 
-    assert rf(-1, pi) == 0
-    assert rf(-5, 1 + I) == 0
-
-    assert unchanged(rf, -3, k)
-    assert unchanged(rf, x, Symbol('k', integer=False))
-    assert rf(-3, Symbol('k', integer=False)) == 0
-    assert rf(Symbol('x', negative=True, integer=True), Symbol('k', integer=False)) == 0
-
     assert rf(x, 0) == 1
     assert rf(x, 1) == x
     assert rf(x, 2) == x*(x + 1)
diff --git a/sympy/matrices/expressions/tests/test_matpow.py b/sympy/matrices/expressions/tests/test_matpow.py
index 0fecc79082..906ff05c79 100644
--- a/sympy/matrices/expressions/tests/test_matpow.py
+++ b/sympy/matrices/expressions/tests/test_matpow.py
@@ -77,6 +77,9 @@ def test_doit_with_MatrixBase():
     assert MatPow(X, -2).doit() == (X.inv())**2
     # less expensive than testing on a 2x2
     assert MatPow(ImmutableMatrix([4]), S.Half).doit() == ImmutableMatrix([2])
+    X = ImmutableMatrix([[0, 2], [0, 4]]) # det() == 0
+    raises(ValueError, lambda: MatPow(X,-1).doit())
+    raises(ValueError, lambda: MatPow(X,-2).doit())
 
 
 def test_doit_nonsquare():
@@ -88,6 +91,11 @@ def test_doit_nonsquare():
     raises(ShapeError, lambda: MatPow(X, pi).doit())
 
 
+def test_doit_equals_pow(): #17179
+    X = ImmutableMatrix ([[1,0],[0,1]])
+    assert MatPow(X, n).doit() == X**n == X
+
+
 def test_doit_nested_MatrixExpr():
     X = ImmutableMatrix([[1, 2], [3, 4]])
     Y = ImmutableMatrix([[2, 3], [4, 5]])
diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py
index bb0f04ad79..c7e7b8555d 100644
--- a/sympy/matrices/tests/test_matrices.py
+++ b/sympy/matrices/tests/test_matrices.py
@@ -14,6 +14,7 @@
     rot_axis3, wronskian, zeros, MutableDenseMatrix, ImmutableDenseMatrix, MatrixSymbol)
 from sympy.core.compatibility import long, iterable, range, Hashable
 from sympy.core import Tuple, Wild
+from sympy.functions.special.tensor_functions import KroneckerDelta
 from sympy.utilities.iterables import flatten, capture
 from sympy.utilities.pytest import raises, XFAIL, slow, skip, warns_deprecated_sympy
 from sympy.solvers import solve
@@ -222,17 +223,17 @@ def test_power():
         [0, 0, b**n]])
 
     A = Matrix([[1, 0], [1, 7]])
-    assert A._matrix_pow_by_jordan_blocks(3) == A._eval_pow_by_recursion(3)
+    assert A._matrix_pow_by_jordan_blocks(S(3)) == A._eval_pow_by_recursion(3)
     A = Matrix([[2]])
-    assert A**10 == Matrix([[2**10]]) == A._matrix_pow_by_jordan_blocks(10) == \
+    assert A**10 == Matrix([[2**10]]) == A._matrix_pow_by_jordan_blocks(S(10)) == \
         A._eval_pow_by_recursion(10)
 
     # testing a matrix that cannot be jordan blocked issue 11766
     m = Matrix([[3, 0, 0, 0, -3], [0, -3, -3, 0, 3], [0, 3, 0, 3, 0], [0, 0, 3, 0, 3], [3, 0, 0, 3, 0]])
-    raises(MatrixError, lambda: m._matrix_pow_by_jordan_blocks(10))
+    raises(MatrixError, lambda: m._matrix_pow_by_jordan_blocks(S(10)))
 
     # test issue 11964
-    raises(ValueError, lambda: Matrix([[1, 1], [3, 3]])._matrix_pow_by_jordan_blocks(-10))
+    raises(MatrixError, lambda: Matrix([[1, 1], [3, 3]])._matrix_pow_by_jordan_blocks(S(-10)))
     A = Matrix([[0, 1, 0], [0, 0, 1], [0, 0, 0]])  # Nilpotent jordan block size 3
     assert A**10.0 == Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
     raises(ValueError, lambda: A**2.1)
@@ -240,13 +241,18 @@ def test_power():
     A = Matrix([[8, 1], [3, 2]])
     assert A**10.0 == Matrix([[1760744107, 272388050], [817164150, 126415807]])
     A = Matrix([[0, 0, 1], [0, 0, 1], [0, 0, 1]])  # Nilpotent jordan block size 1
-    assert A**10.2 == Matrix([[0, 0, 1], [0, 0, 1], [0, 0, 1]])
+    assert A**10.0 == Matrix([[0, 0, 1], [0, 0, 1], [0, 0, 1]])
     A = Matrix([[0, 1, 0], [0, 0, 1], [0, 0, 1]])  # Nilpotent jordan block size 2
     assert A**10.0 == Matrix([[0, 0, 1], [0, 0, 1], [0, 0, 1]])
     n = Symbol('n', integer=True)
     assert isinstance(A**n, MatPow)
-    n = Symbol('n', integer=True, nonnegative=True)
+    n = Symbol('n', integer=True, negative=True)
     raises(ValueError, lambda: A**n)
+    n = Symbol('n', integer=True, nonnegative=True)
+    assert A**n == Matrix([
+        [KroneckerDelta(0, n), KroneckerDelta(1, n), -KroneckerDelta(0, n) - KroneckerDelta(1, n) + 1],
+        [                   0, KroneckerDelta(0, n),                         1 - KroneckerDelta(0, n)],
+        [                   0,                    0,                                                1]])
     assert A**(n + 2) == Matrix([[0, 0, 1], [0, 0, 1], [0, 0, 1]])
     raises(ValueError, lambda: A**(S(3)/2))
     A = Matrix([[0, 0, 1], [3, 0, 1], [4, 3, 1]])
@@ -259,6 +265,25 @@ def test_power():
     raises(ValueError, lambda: An.subs(n, -2).doit())
     assert An * An == A**(2*n)
 
+    # concretizing behavior for non-integer and complex powers
+    A = Matrix([[0,0,0],[0,0,0],[0,0,0]])
+    n = Symbol('n', integer=True, positive=True)
+    assert A**n == A
+    n = Symbol('n', integer=True, nonnegative=True)
+    assert A**n == diag(0**n, 0**n, 0**n)
+    assert (A**n).subs(n, 0) == eye(3)
+    assert (A**n).subs(n, 1) == zeros(3)
+    A = Matrix ([[2,0,0],[0,2,0],[0,0,2]])
+    assert A**2.1 == diag (2**2.1, 2**2.1, 2**2.1)
+    assert A**I == diag (2**I, 2**I, 2**I)
+    A = Matrix([[0, 1, 0], [0, 0, 1], [0, 0, 1]])
+    raises(ValueError, lambda: A**2.1)
+    raises(ValueError, lambda: A**I)
+    A = Matrix([[S.Half, S.Half], [S.Half, S.Half]])
+    assert A**S.Half == A
+    A = Matrix([[1, 1],[3, 3]])
+    assert A**S.Half == Matrix ([[S.Half, S.Half], [3*S.Half, 3*S.Half]])
+
 
 def test_creation():
     raises(ValueError, lambda: Matrix(5, 5, range(20)))