sagemathgh-36845: Give more precise answer on symbolic power of a matrix
    
- Changed _matrix_power_symbolic function to consider mk=0 case
- Changed the _matrix_power_symbolic function to handle all the cases
- Created tests covering the changes

This PR improves the answer given by _matrix_power_symbolic() function,
which involves finding the symbolic power of a matrix, and also gives
more precise answer on evalution. eg, before, when we were evaluating
the zeroth power of a matrix sometimes it was giving wrong answer like
below,
```python
        sage: n = var("n")
        sage: A = matrix(ZZ, [[1,-1], [-1,1]])
        sage: An = A^(n)
        [ 1/2*2^n -1/2*2^n]
        [-1/2*2^n  1/2*2^n]
        sage: An({n:0})
        [ 1/2 -1/2]
        [-1/2  1/2]
```
Here, for n=0, it is not giving an identity matrix which is wrong, but
after the change it will give answer as below,
```python
        sage: n = var("n")
        sage: A = matrix(ZZ, [[1,-1], [-1,1]])
        sage: An = A^(n)
        [ 1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1) -1/2*2^(2*n
+ 1) + 1/2*kronecker_delta(0, 2*n + 1)]
        [-1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)  1/2*2^(2*n
+ 1) + 1/2*kronecker_delta(0, 2*n + 1)]
        sage: An({n:0})
        [1 0]
        [0 1]
```

Also, This patch fixes sagemath#36838 . The function, _matrix_power_symbolic()
was giving symbolic power
of a nilpotent matrix as a null matrix but the correct answer should be
represented
in the form of kronecker_delta() function. In the case of mk=0,
simplifying only binomial
term should work because afterall the form 0^(n-i) should be simplifed
to kronecker_delta() function and no further simplification is needed.
And in every other case it is simplifying whole term so it will handle
all other cases.
<!-- ^^^^^
Please provide a concise, informative and self-explanatory title.
Don't put issue numbers in there, do this in the PR body below.
For example, instead of "Fixes sagemath#1234" use "Introduce new method to
calculate 1+1"
-->
<!-- Describe your changes here in detail -->

<!-- Why is this change required? What problem does it solve? -->
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes sagemath#12345". -->
<!-- If your change requires a documentation PR, please link it
appropriately. -->

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->
<!-- If your change requires a documentation PR, please link it
appropriately -->
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
<!-- Feel free to remove irrelevant items. -->

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on
- sagemath#12345: short description why this is a dependency
- sagemath#34567: ...
-->

<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
None
    
URL: sagemath#36845
Reported by: Ruchit Jagodara
Reviewer(s): Dima Pasechnik, phul-ste, Ruchit Jagodara

Author: Release Manager

src/sage/matrix/matrix2.pyx

@@ -216,10 +216,10 @@ cdef class Matrix(Matrix1):
         from sage.matrix.constructor import matrix
         if self.is_sparse():
             return matrix({ij: self[ij].subs(*args, **kwds) for ij in self.nonzero_positions()},
-                    nrows=self._nrows, ncols=self._ncols, sparse=True)
+                          nrows=self._nrows, ncols=self._ncols, sparse=True)
         else:
             return matrix([a.subs(*args, **kwds) for a in self.list()],
-                        nrows=self._nrows, ncols=self._ncols, sparse=False)
+                          nrows=self._nrows, ncols=self._ncols, sparse=False)
 
     def solve_left(self, B, check=True):
         """
@@ -3183,14 +3183,12 @@ cdef class Matrix(Matrix1):
         """
 
         # 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
-        #
         cdef Matrix M  = <Matrix> self
         n  = M._ncols
         R  = M._base_ring
@@ -3199,7 +3197,6 @@ cdef class Matrix(Matrix1):
         # 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.
-        #
         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
@@ -18607,8 +18598,8 @@ def _matrix_power_symbolic(A, n):
 
         sage: A = matrix(ZZ, [[1,-1], [-1,1]])
         sage: A^(2*n+1)                                                                 # needs sage.symbolic
-        [ 1/2*2^(2*n + 1) -1/2*2^(2*n + 1)]
-        [-1/2*2^(2*n + 1)  1/2*2^(2*n + 1)]
+        [ 1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1) -1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)]
+        [-1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)  1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)]
 
     Check if :trac:`23215` is fixed::
 
@@ -18618,12 +18609,25 @@ def _matrix_power_symbolic(A, n):
          -1/2*I*(a + I*b)^k + 1/2*I*(a - I*b)^k,
          1/2*I*(a + I*b)^k - 1/2*I*(a - I*b)^k,
          1/2*(a + I*b)^k + 1/2*(a - I*b)^k]
+
+    Check if :trac:`36838` is fixed:Checking symbolic power of
+    nilpotent matrix::
+
+        sage: A = matrix([[0,1],[0,0]]); A
+        [0 1]
+        [0 0]
+        sage: n = var('n'); n
+        n
+        sage: B = A^n; B
+        [  kronecker_delta(0, n) n*kronecker_delta(1, n)]
+        [                      0   kronecker_delta(0, n)]
     """
     from sage.rings.qqbar import AlgebraicNumber
     from sage.matrix.constructor import matrix
     from sage.functions.other import binomial
     from sage.symbolic.ring import SR
     from sage.rings.qqbar import QQbar
+    from sage.functions.generalized import kronecker_delta
 
     got_SR = A.base_ring() == SR
 
@@ -18656,8 +18660,17 @@ def _matrix_power_symbolic(A, n):
         # D^i(f) / i! with f = x^n and D = differentiation wrt x
         if hasattr(mk, 'radical_expression'):
             mk = mk.radical_expression()
-        vk = [(binomial(n, i) * mk**(n-i)).simplify_full()
-              for i in range(nk)]
+
+
+        # When the variable "mk" is equal to zero, it is advisable to employ the Kronecker delta function
+        # instead of utilizing the numerical value zero. This choice is made to encompass scenarios where
+        # the power of zero is also equal to zero.
+        if mk:
+            vk = [(binomial(n, i) * mk._pow_(n-i)).simplify_full()
+                  for i in range(nk)]
+        else:
+            vk = [(binomial(n, i).simplify_full() * kronecker_delta(n,i))
+                  for i in range(nk)]
 
         # Form block Mk and insert it in M
         Mk = matrix(SR, [[SR.zero()]*i + vk[:nk-i] for i in range(nk)])