From afa6a7ca353f16473528bcada997138423755844 Mon Sep 17 00:00:00 2001 From: RuchitJagodara Date: Sun, 10 Dec 2023 01:44:04 +0530 Subject: [PATCH 1/9] Changed _matrix_power_symbolic function to condier mk=0 case This fixes #36838. Here, I have added a condition to check whether mk=0 or not. Because whenever mk=0 we should give mk^(n-i) (i.e. 0^(n-i)) instead of only 0 considering (n-i) can be equal to zero and in this case 0^(n-i) will be more accurate than only 0. --- src/sage/matrix/matrix2.pyx | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 4acd46db4de..020a1660028 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -18656,8 +18656,14 @@ 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)] + + # Return mk^(n-i) instead of mk**(n-i) if mk=0 + if mk: + vk = [(binomial(n, i) * mk**(n-i)).simplify_full() + for i in range(nk)] + else: + vk = [(binomial(n, i)).simplify_full() * (mk^(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)]) From dc79089d51ac56d2651960660b89eb12e389fd4c Mon Sep 17 00:00:00 2001 From: RuchitJagodara Date: Sun, 10 Dec 2023 02:46:02 +0530 Subject: [PATCH 2/9] Changed the _matrix_power_symbolic function to handle all the cases This change handles all the errors which were occuring in previous commit and this commit handles the case of mx=0 very effectively. --- src/sage/matrix/matrix2.pyx | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 020a1660028..ec120870785 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -18624,6 +18624,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 got_SR = A.base_ring() == SR @@ -18659,11 +18660,11 @@ def _matrix_power_symbolic(A, n): # Return mk^(n-i) instead of mk**(n-i) if mk=0 if mk: - vk = [(binomial(n, i) * mk**(n-i)).simplify_full() + vk = [(binomial(n, i) * mk._pow_(n-i)).simplify_full() for i in range(nk)] else: - vk = [(binomial(n, i)).simplify_full() * (mk^(n-i)) - for i in range(nk)] + vk = [(binomial(n, i).simplify_full() * mk._pow_(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)]) From 7139bb86a7d0c574a92a1ffc176adae9d57aa2a8 Mon Sep 17 00:00:00 2001 From: RuchitJagodara Date: Sun, 10 Dec 2023 02:51:15 +0530 Subject: [PATCH 3/9] Created tests covering the changes Created tests covering the changes and checking whether the trac:`36838` is fixed or not. --- src/sage/matrix/matrix2.pyx | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index ec120870785..eade57082af 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -18618,6 +18618,18 @@ 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: An = A^n; An + [ 0^n 0^(n - 1)*n] + [ 0 0^n] """ from sage.rings.qqbar import AlgebraicNumber from sage.matrix.constructor import matrix From 42454c0fd387bab302907e68b7ed695b9563af17 Mon Sep 17 00:00:00 2001 From: RuchitJagodara Date: Sun, 10 Dec 2023 09:48:29 +0530 Subject: [PATCH 4/9] 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 | 23 +++++++---------------- 1 file changed, 7 insertions(+), 16 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index eade57082af..0f69ef06f02 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -3182,15 +3182,13 @@ 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 = self n = M._ncols R = M._base_ring @@ -3198,8 +3196,7 @@ 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. - # + # 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,18 +3222,15 @@ 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): @@ -3245,11 +3238,9 @@ cdef class Matrix(Matrix1): 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 From 96442fb9494cab668b7b0f227a204a758b2ccc8f Mon Sep 17 00:00:00 2001 From: RuchitJagodara Date: Sun, 10 Dec 2023 09:55:38 +0530 Subject: [PATCH 5/9] Rewriting tests for the PR Changed the final answer given by test of this PR. --- src/sage/matrix/matrix2.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 0f69ef06f02..cb08e47dc51 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -18619,8 +18619,8 @@ def _matrix_power_symbolic(A, n): sage: n = var('n'); n n sage: An = A^n; An - [ 0^n 0^(n - 1)*n] - [ 0 0^n] + [ 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 d40d348e5ed2198f1a79795d9065699e1d4ec7d1 Mon Sep 17 00:00:00 2001 From: RuchitJagodara Date: Thu, 14 Dec 2023 00:42:42 +0530 Subject: [PATCH 6/9] Correct answers of the doctest for _matrix_power_symbolic function Changed the answer of doctest according to new changes. --- src/sage/matrix/matrix2.pyx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index cb08e47dc51..789eed0ef8c 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -18598,8 +18598,9 @@ 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:: From 935a42095647e3e51aaaa07e88b2f4114fa047d8 Mon Sep 17 00:00:00 2001 From: RuchitJagodara Date: Thu, 14 Dec 2023 02:05:03 +0530 Subject: [PATCH 7/9] Corrected the doctest Corrected the doctest and improved some code styles, which were not correct according to guidelines for python coding for sage. --- src/sage/matrix/matrix2.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 789eed0ef8c..970a89893e0 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/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): """ From 7955a51b06310e09a55f6f9f891b52b1a45fb1b3 Mon Sep 17 00:00:00 2001 From: RuchitJagodara Date: Thu, 14 Dec 2023 13:50:37 +0530 Subject: [PATCH 8/9] Modified the doctest and changed the comment Updated the old doctests which were failing before and updated the comment to make it more readable. --- src/sage/matrix/matrix2.pyx | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 970a89893e0..488098f25e8 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -18601,7 +18601,6 @@ def _matrix_power_symbolic(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)] - Check if :trac:`23215` is fixed:: sage: a, b, k = var('a, b, k') # needs sage.symbolic @@ -18619,7 +18618,7 @@ def _matrix_power_symbolic(A, n): [0 0] sage: n = var('n'); n n - sage: An = A^n; An + sage: B = A^n; B [ kronecker_delta(0, n) n*kronecker_delta(1, n)] [ 0 kronecker_delta(0, n)] """ @@ -18662,7 +18661,10 @@ def _matrix_power_symbolic(A, n): if hasattr(mk, 'radical_expression'): mk = mk.radical_expression() - # Return mk^(n-i) instead of mk**(n-i) if mk=0 + + # 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)] From ed4e71208b7822900e20dfcfbc0581292ded0c5e Mon Sep 17 00:00:00 2001 From: Ruchit Jagodara Date: Sun, 17 Dec 2023 23:00:01 +0530 Subject: [PATCH 9/9] Corrected lint errors --- src/sage/matrix/matrix2.pyx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 488098f25e8..0c2aa9eda4b 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -3182,13 +3182,13 @@ 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 = self n = M._ncols R = M._base_ring @@ -3196,7 +3196,7 @@ 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. + # test for 0 x n or m x 0 matrices. if n == 0: return S.one() @@ -18662,8 +18662,8 @@ def _matrix_power_symbolic(A, n): mk = mk.radical_expression() - # 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 + # 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()