Give more precise answer by using kroncker_delta function

Instead of returning 0^(n-i), it would be more precise if
we reutrn value of  kroncker_delta function, it will be helpful when
we try to evaluate the final matrix using the value of n.

src/sage/matrix/matrix2.pyx

@@ -3182,24 +3182,21 @@ cdef class Matrix(Matrix1):
 
         """
 
-        # Validate assertions
-        #
+        # Validate assertions        
         if not self.is_square():
             raise ValueError("self must be a square matrix")
 
         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 
-        # Extract parameters
-        #
+        # Extract parameters        
         cdef Matrix M  = <Matrix> self
         n  = M._ncols
         R  = M._base_ring
         S  = PolynomialRing(R, var)
 
         # Corner cases
         # N.B.  We already tested for M to be square, hence we do not need to
-        # test for 0 x n or m x 0 matrices.
-        #
+        # test for 0 x n or m x 0 matrices.      
         if n == 0:
             return S.one()
 
@@ -3215,7 +3212,6 @@ cdef class Matrix(Matrix1):
         #
         # N.B.  The documentation is still 1-based, although the code, after
         # having been ported from Magma to Sage, is 0-based.
-        #
         from sage.matrix.constructor import matrix
 
         F = [R.zero()] * n
@@ -3226,30 +3222,25 @@ cdef class Matrix(Matrix1):
         for t in range(1,n):
 
             # Set a(1, t) to be M(<=t, t)
-            #
             for i in range(t+1):
                 a.set_unsafe(0, i, M.get_unsafe(i, t))
 
             # Set A[1, t] to be the (t)th entry in a[1, t]
-            #
             A[0] = M.get_unsafe(t, t)
 
             for p in range(1, t):
 
                 # Set a(p, t) to the product of M[<=t, <=t] * a(p-1, t)
-                #
                 for i in range(t+1):
                     s = R.zero()
                     for j in range(t+1):
                         s = s + M.get_unsafe(i, j) * a.get_unsafe(p-1, j)
                     a.set_unsafe(p, i, s)
 
                 # Set A[p, t] to be the (t)th entry in a[p, t]
-                #
                 A[p] = a.get_unsafe(p, t)
 
             # Set A[t, t] to be M[t, <=t] * a(p-1, t)
-            #
             s = R.zero()
             for j in range(t+1):
                 s = s + M.get_unsafe(t, j) * a.get_unsafe(t-1, j)
@@ -3262,7 +3253,7 @@ cdef class Matrix(Matrix1):
                 F[p] = s - A[p]
 
         X = S.gen(0)
-        f = X ** n + S(list(reversed(F)))
+        f = X**n + S(list(reversed(F)))
 
         return f
 
@@ -3634,7 +3625,7 @@ cdef class Matrix(Matrix1):
         # using the rows of a matrix.)  Also see Cohen's first GTM,
         # Algorithm 2.2.9.
 
-        cdef Py_ssize_t i, m, n,
+        cdef Py_ssize_t i, m, n
         n = self._nrows
 
         cdef Matrix c
@@ -18636,7 +18627,7 @@ def _matrix_power_symbolic(A, n):
     from sage.functions.other import binomial
     from sage.symbolic.ring import SR
     from sage.rings.qqbar import QQbar
-    from sage.symbolic.expression import Expression
+    from sage.functions.generalized import kronecker_delta
 
     got_SR = A.base_ring() == SR
 
@@ -18675,7 +18666,7 @@ def _matrix_power_symbolic(A, n):
             vk = [(binomial(n, i) * mk._pow_(n-i)).simplify_full()
                   for i in range(nk)]
         else:
-            vk = [(binomial(n, i).simplify_full() * mk._pow_(n-i))
+            vk = [(binomial(n, i).simplify_full() * kronecker_delta(n,i))
                   for i in range(nk)]
 
         # Form block Mk and insert it in M