Skip to content

Commit

Permalink
gh-36845: Give more precise answer on symbolic power of a matrix
Browse files Browse the repository at this point in the history
    
- 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 #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 #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 #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
- #12345: short description why this is a dependency
- #34567: ...
-->

<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
None
    
URL: #36845
Reported by: Ruchit Jagodara
Reviewer(s): Dima Pasechnik, phul-ste, Ruchit Jagodara
  • Loading branch information
Release Manager committed Dec 25, 2023
2 parents 003c1d9 + f95793b commit 58e899c
Showing 1 changed file with 30 additions and 17 deletions.
47 changes: 30 additions & 17 deletions src/sage/matrix/matrix2.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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
Expand All @@ -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()

Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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::

Expand All @@ -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

Expand Down Expand Up @@ -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)])
Expand Down

0 comments on commit 58e899c

Please sign in to comment.