ASKSAGE: Sage Q&A Forum - Individual question feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Mon, 02 Apr 2018 13:04:15 -0500Bug in general power of a matrixhttps://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/The code to [question 35658](http://ask.sagemath.org/question/35658/ge) gives a wrong answer i.e. for the matrix
A=matrix([[4,1,2],[0,2,-4],[0,1,6]])
Where can I find an improvement?
For $k\in {\bf N}$ the $k$-th power of a matrix $A\in {\bf R}^{n\times n}$ satisfies by definition the recursion $A^{k+1}=A*A^k$ and the initial condition $A^0=I$.
Substitution of $k=0$ and $k=1$ should therefore give the identity matrix $I$ resp. the matrix $A$. The following
[live code](http://sagecell.sagemath.org/?z=eJzLVrBVKEss0lDPVtfk5XK0zU0sKcqs0IiONtEx1DGK1Yk20DHS0TUBMwx1zGJjgaoK8ktsHeOyebmKM_LLNYA8oBiCo1dcmlSskW1roKmJKWioqQkA_MQgTQ==&lang=sage) does (currently) not give the expected answer.
Btw: One of the [M-programs](http://www.wolframalpha.com/input/?i=MatrixPower[{{4,1,2},{0,2,-4},{0,1,6}},k])
gives the correct answer.Mon, 19 Mar 2018 07:45:23 -0500https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/Comment by dan_fulea for <p>The code to <a href="http://ask.sagemath.org/question/35658/ge">question 35658</a> gives a wrong answer i.e. for the matrix</p>
<pre><code>A=matrix([[4,1,2],[0,2,-4],[0,1,6]])
</code></pre>
<p>Where can I find an improvement?</p>
<p>For $k\in {\bf N}$ the $k$-th power of a matrix $A\in {\bf R}^{n\times n}$ satisfies by definition the recursion $A^{k+1}=A*A^k$ and the initial condition $A^0=I$.
Substitution of $k=0$ and $k=1$ should therefore give the identity matrix $I$ resp. the matrix $A$. The following
<a href="http://sagecell.sagemath.org/?z=eJzLVrBVKEss0lDPVtfk5XK0zU0sKcqs0IiONtEx1DGK1Yk20DHS0TUBMwx1zGJjgaoK8ktsHeOyebmKM_LLNYA8oBiCo1dcmlSskW1roKmJKWioqQkA_MQgTQ==&lang=sage">live code</a> does (currently) not give the expected answer.</p>
<p>Btw: One of the <a href="http://www.wolframalpha.com/input/?i=MatrixPower[{{4,1,2},{0,2,-4},{0,1,6}},k]">M-programs</a>
gives the correct answer.</p>
https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41725#post-id-41725I finally got the point of the question. (Too many hidden links in a simple question. Please always give code here, not elsewhere with a to be copied html address...)
Bug:
sage: k = var( 'k' )
sage: A = matrix( QQbar, 3, 3, [[4,1,0],[0,4,0],[0,0,4]] )
sage: A
[4 1 0]
[0 4 0]
[0 0 4]
sage: A.jordan_form()
[4 1|0]
[0 4|0]
[---+-]
[0 0|4]
sage: A^k
[ 4^k 4^(k - 1)*k 0]
[ 0 4^k 0]
[ 0 0 0]
The Jordan blocks are not correctly parsed, the $0$ in the bottom right corner...Thu, 22 Mar 2018 19:01:48 -0500https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41725#post-id-41725Comment by slelievre for <p>The code to <a href="http://ask.sagemath.org/question/35658/ge">question 35658</a> gives a wrong answer i.e. for the matrix</p>
<pre><code>A=matrix([[4,1,2],[0,2,-4],[0,1,6]])
</code></pre>
<p>Where can I find an improvement?</p>
<p>For $k\in {\bf N}$ the $k$-th power of a matrix $A\in {\bf R}^{n\times n}$ satisfies by definition the recursion $A^{k+1}=A*A^k$ and the initial condition $A^0=I$.
Substitution of $k=0$ and $k=1$ should therefore give the identity matrix $I$ resp. the matrix $A$. The following
<a href="http://sagecell.sagemath.org/?z=eJzLVrBVKEss0lDPVtfk5XK0zU0sKcqs0IiONtEx1DGK1Yk20DHS0TUBMwx1zGJjgaoK8ktsHeOyebmKM_LLNYA8oBiCo1dcmlSskW1roKmJKWioqQkA_MQgTQ==&lang=sage">live code</a> does (currently) not give the expected answer.</p>
<p>Btw: One of the <a href="http://www.wolframalpha.com/input/?i=MatrixPower[{{4,1,2},{0,2,-4},{0,1,6}},k]">M-programs</a>
gives the correct answer.</p>
https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41732#post-id-41732Thanks for this simpler example, which put me on the track to pin down the bug.
The problem was with `FJ.set_block(k, k, ...)` in the code of `_matrix_power_symbolic`.
See my answer, with a replacement of `k, k` by the appropriate index where the block needs to go.Fri, 23 Mar 2018 05:09:57 -0500https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41732#post-id-41732Comment by kcrisman for <p>The code to <a href="http://ask.sagemath.org/question/35658/ge">question 35658</a> gives a wrong answer i.e. for the matrix</p>
<pre><code>A=matrix([[4,1,2],[0,2,-4],[0,1,6]])
</code></pre>
<p>Where can I find an improvement?</p>
<p>For $k\in {\bf N}$ the $k$-th power of a matrix $A\in {\bf R}^{n\times n}$ satisfies by definition the recursion $A^{k+1}=A*A^k$ and the initial condition $A^0=I$.
Substitution of $k=0$ and $k=1$ should therefore give the identity matrix $I$ resp. the matrix $A$. The following
<a href="http://sagecell.sagemath.org/?z=eJzLVrBVKEss0lDPVtfk5XK0zU0sKcqs0IiONtEx1DGK1Yk20DHS0TUBMwx1zGJjgaoK8ktsHeOyebmKM_LLNYA8oBiCo1dcmlSskW1roKmJKWioqQkA_MQgTQ==&lang=sage">live code</a> does (currently) not give the expected answer.</p>
<p>Btw: One of the <a href="http://www.wolframalpha.com/input/?i=MatrixPower[{{4,1,2},{0,2,-4},{0,1,6}},k]">M-programs</a>
gives the correct answer.</p>
https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41641#post-id-41641Also, if you read the other ticket carefully, you will see this is already *in* Sage now: see [here for a live example](http://sagecell.sagemath.org/?z=eJzLVrBVKEss0lDPVtfk5XK0zU0sKcqs0IiONtEx1DGK1Yk20DHS0TUBMwx1zGJjQariNIy0srUNNQHvFw8S&lang=sage). Just define as completely normal and take the power.Mon, 19 Mar 2018 12:29:11 -0500https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41641#post-id-41641Comment by kcrisman for <p>The code to <a href="http://ask.sagemath.org/question/35658/ge">question 35658</a> gives a wrong answer i.e. for the matrix</p>
<pre><code>A=matrix([[4,1,2],[0,2,-4],[0,1,6]])
</code></pre>
<p>Where can I find an improvement?</p>
<p>For $k\in {\bf N}$ the $k$-th power of a matrix $A\in {\bf R}^{n\times n}$ satisfies by definition the recursion $A^{k+1}=A*A^k$ and the initial condition $A^0=I$.
Substitution of $k=0$ and $k=1$ should therefore give the identity matrix $I$ resp. the matrix $A$. The following
<a href="http://sagecell.sagemath.org/?z=eJzLVrBVKEss0lDPVtfk5XK0zU0sKcqs0IiONtEx1DGK1Yk20DHS0TUBMwx1zGJjgaoK8ktsHeOyebmKM_LLNYA8oBiCo1dcmlSskW1roKmJKWioqQkA_MQgTQ==&lang=sage">live code</a> does (currently) not give the expected answer.</p>
<p>Btw: One of the <a href="http://www.wolframalpha.com/input/?i=MatrixPower[{{4,1,2},{0,2,-4},{0,1,6}},k]">M-programs</a>
gives the correct answer.</p>
https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41639#post-id-41639I bet this method will work if you do `A=matrix(QQ,[[4,1,2],[0,2,-4],[0,1,6]])`.Mon, 19 Mar 2018 12:25:59 -0500https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41639#post-id-41639Answer by slelievre for <p>The code to <a href="http://ask.sagemath.org/question/35658/ge">question 35658</a> gives a wrong answer i.e. for the matrix</p>
<pre><code>A=matrix([[4,1,2],[0,2,-4],[0,1,6]])
</code></pre>
<p>Where can I find an improvement?</p>
<p>For $k\in {\bf N}$ the $k$-th power of a matrix $A\in {\bf R}^{n\times n}$ satisfies by definition the recursion $A^{k+1}=A*A^k$ and the initial condition $A^0=I$.
Substitution of $k=0$ and $k=1$ should therefore give the identity matrix $I$ resp. the matrix $A$. The following
<a href="http://sagecell.sagemath.org/?z=eJzLVrBVKEss0lDPVtfk5XK0zU0sKcqs0IiONtEx1DGK1Yk20DHS0TUBMwx1zGJjgaoK8ktsHeOyebmKM_LLNYA8oBiCo1dcmlSskW1roKmJKWioqQkA_MQgTQ==&lang=sage">live code</a> does (currently) not give the expected answer.</p>
<p>Btw: One of the <a href="http://www.wolframalpha.com/input/?i=MatrixPower[{{4,1,2},{0,2,-4},{0,1,6}},k]">M-programs</a>
gives the correct answer.</p>
https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?answer=41730#post-id-41730Thank you for pointing out this problem.
The computation of symbolic powers of a matrix was recently
[introduced in Sage](https://trac.sagemath.org/ticket/22523), based on an idea
[proposed on Ask Sage](https://ask.sagemath.org/question/35658/general-power-of-a-matrix/):
- compute the Jordan form, and the transformation matrix
- compute powers of each Jordan block
- form a block-diagonal matrix with these powers
- change back to the old basis using the transformation matrix
The idea is good. Unfortunately there was a glitch in the
implementation: the powers of the Jordan blocks were being
inserted at the wrong place.
To illustrate the bug in an enlightening way, consider the matrix:
sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]
The current implementation in Sage gives:
sage: B = A^n; B
[ 2^n 2^(n - 1)*n 0]
[ 0 3^n 0]
[ 0 0 0]
Block k is inserted starting at position (k, k)
instead of its correct position.
The revised implementation, suggested below, gives:
sage: B = _matrix_power_symbolic(A, n); B
[ 2^n 2^(n - 1)*n 0]
[ 0 2^n 0]
[ 0 0 3^n]
and can deal with the matrix in the question here.
sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4 1 2]
[ 0 2 -4]
[ 0 1 6]
sage: B = _matrix_power_symbolic(A, n); B
[ 4^n 4^(n - 1)*n 2*4^(n - 1)*n]
[ 0 -2*4^(n - 1)*n + 4^n -4*4^(n - 1)*n]
[ 0 4^(n - 1)*n 2*4^(n - 1)*n + 4^n]
Here is a revised implementation with a fix for the glitch,
along with some efficiency improvements.
def _matrix_power_symbolic(A, n):
r"""
Return the symbolic matrix power `A^n`
This function implements the computation of `A^n` for symbolic `n`,
relying on the Jordan normal form of `A`, available for exact rings
as :meth:`jordan_form`. See [Hig2008]_, §1.2, for further details.
INPUT:
- ``A`` -- a square matrix over an exact field
- ``n`` -- the symbolic exponent
OUTPUT:
The matrix `A^n` (with symbolic entries).
EXAMPLES:
Powers of two by two matrix::
sage: n = SR.var('n')
sage: A = matrix(QQ, [[2, -1], [1, 0]])
sage: B = A^n
sage: B
[ n + 1 -n]
[ n -n + 1]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
Powers of a three by three matrix in Jordan form::
sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]
sage: B = A^n; B
[ 2^n 2^(n - 1)*n 0]
[ 0 2^n 0]
[ 0 0 3^n]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
Powers of a three by three matrix not in Jordan form::
sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4 1 2]
[ 0 2 -4]
[ 0 1 6]
sage: B = A^n
sage: B
[ 4^n 4^(n - 1)*n 2*4^(n - 1)*n]
[ 0 -2*4^(n - 1)*n + 4^n -4*4^(n - 1)*n]
[ 0 4^(n - 1)*n 2*4^(n - 1)*n + 4^n]
sage: [B.subs({n: k}) for k in range(4)]
[
[1 0 0] [ 4 1 2] [ 16 8 16] [ 64 48 96]
[0 1 0] [ 0 2 -4] [ 0 0 -32] [ 0 -32 -192]
[0 0 1], [ 0 1 6], [ 0 8 32], [ 0 48 160]
]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
TESTS:
Testing exponentiation in the symbolic ring::
sage: n = var('n')
sage: A = matrix([[pi, e],[0, -2*I]])
sage: (A^n).list()
[pi^n,
-(-2*I)^n/(pi*e^(-1) + 2*I*e^(-1)) + pi^n/(pi*e^(-1) + 2*I*e^(-1)),
0,
(-2*I)^n]
If the base ring is inexact, the Jordan normal form is not available::
sage: A = matrix(RDF, [[2, -1], [1, 0]])
sage: A^n
Traceback (most recent call last):
...
ValueError: Jordan normal form not implemented over inexact rings.
Testing exponentiation in the integer ring::
sage: A = matrix(ZZ, [[1,-1],[-1,1]])
sage: A^(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)]
Check if :trac:`23215` is fixed::
sage: a, b, k = var('a, b, k')
sage: (matrix(2, [a, b, -b, a])^k).list()
[1/2*(a + I*b)^k + 1/2*(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*I*(a - I*b)^k,
1/2*(a + I*b)^k + 1/2*(a - I*b)^k]
"""
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
got_SR = A.base_ring() == SR
# Change to QQbar if possible
try:
A = A.change_ring(QQbar)
except (TypeError, NotImplementedError):
pass
# Get Jordan matrix J and invertible matrix P such that A = P*J*~P
# From that, we will compute M = J^n, and obtain A^n = P*J^n*~P
J, P = A.jordan_form(transformation=True)
# Where each Jordan block starts, and number of blocks
block_start = [0] + J.subdivisions()[0]
num_blocks = len(block_start)
# Prepare matrix M to store `J^n`, computed by Jordan block
M = matrix(SR, J.ncols())
M.subdivide(J.subdivisions())
for k in range(num_blocks):
# Jordan block Jk, its dimension nk, the eigenvalue m
Jk = J.subdivision(k, k)
nk = Jk.ncols()
mk = Jk[0,0]
# First row of block Mk; its entries are of the form
# 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)]
# Form block Mk and insert it in M
Mk = matrix(SR, [[SR.zero()]*i + vk[:nk-i] for i in range(nk)])
M.set_block(block_start[k], block_start[k], Mk)
# Change entries of P and P^-1 into symbolic expressions
if not got_SR:
Pinv = (~P).apply_map(AlgebraicNumber.radical_expression)
P = P.apply_map(AlgebraicNumber.radical_expression)
else:
Pinv = ~P
return P * M * Pinv
I will open a ticket on Sage's issue tracker to get the fix in Sage.
Edit: The issue is now tracked at
[Sage Trac ticket #25082](https://trac.sagemath.org/ticket/25082).Fri, 23 Mar 2018 05:02:24 -0500https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?answer=41730#post-id-41730Comment by slelievre for <p>Thank you for pointing out this problem.</p>
<p>The computation of symbolic powers of a matrix was recently
<a href="https://trac.sagemath.org/ticket/22523">introduced in Sage</a>, based on an idea
<a href="https://ask.sagemath.org/question/35658/general-power-of-a-matrix/">proposed on Ask Sage</a>:</p>
<ul>
<li>compute the Jordan form, and the transformation matrix</li>
<li>compute powers of each Jordan block</li>
<li>form a block-diagonal matrix with these powers</li>
<li>change back to the old basis using the transformation matrix</li>
</ul>
<p>The idea is good. Unfortunately there was a glitch in the
implementation: the powers of the Jordan blocks were being
inserted at the wrong place.</p>
<p>To illustrate the bug in an enlightening way, consider the matrix:</p>
<pre><code>sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]
</code></pre>
<p>The current implementation in Sage gives:</p>
<pre><code>sage: B = A^n; B
[ 2^n 2^(n - 1)*n 0]
[ 0 3^n 0]
[ 0 0 0]
</code></pre>
<p>Block k is inserted starting at position (k, k)
instead of its correct position.</p>
<p>The revised implementation, suggested below, gives:</p>
<pre><code>sage: B = _matrix_power_symbolic(A, n); B
[ 2^n 2^(n - 1)*n 0]
[ 0 2^n 0]
[ 0 0 3^n]
</code></pre>
<p>and can deal with the matrix in the question here.</p>
<pre><code>sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4 1 2]
[ 0 2 -4]
[ 0 1 6]
sage: B = _matrix_power_symbolic(A, n); B
[ 4^n 4^(n - 1)*n 2*4^(n - 1)*n]
[ 0 -2*4^(n - 1)*n + 4^n -4*4^(n - 1)*n]
[ 0 4^(n - 1)*n 2*4^(n - 1)*n + 4^n]
</code></pre>
<p>Here is a revised implementation with a fix for the glitch,
along with some efficiency improvements.</p>
<pre><code>def _matrix_power_symbolic(A, n):
r"""
Return the symbolic matrix power `A^n`
This function implements the computation of `A^n` for symbolic `n`,
relying on the Jordan normal form of `A`, available for exact rings
as :meth:`jordan_form`. See [Hig2008]_, §1.2, for further details.
INPUT:
- ``A`` -- a square matrix over an exact field
- ``n`` -- the symbolic exponent
OUTPUT:
The matrix `A^n` (with symbolic entries).
EXAMPLES:
Powers of two by two matrix::
sage: n = SR.var('n')
sage: A = matrix(QQ, [[2, -1], [1, 0]])
sage: B = A^n
sage: B
[ n + 1 -n]
[ n -n + 1]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
Powers of a three by three matrix in Jordan form::
sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]
sage: B = A^n; B
[ 2^n 2^(n - 1)*n 0]
[ 0 2^n 0]
[ 0 0 3^n]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
Powers of a three by three matrix not in Jordan form::
sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4 1 2]
[ 0 2 -4]
[ 0 1 6]
sage: B = A^n
sage: B
[ 4^n 4^(n - 1)*n 2*4^(n - 1)*n]
[ 0 -2*4^(n - 1)*n + 4^n -4*4^(n - 1)*n]
[ 0 4^(n - 1)*n 2*4^(n - 1)*n + 4^n]
sage: [B.subs({n: k}) for k in range(4)]
[
[1 0 0] [ 4 1 2] [ 16 8 16] [ 64 48 96]
[0 1 0] [ 0 2 -4] [ 0 0 -32] [ 0 -32 -192]
[0 0 1], [ 0 1 6], [ 0 8 32], [ 0 48 160]
]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
TESTS:
Testing exponentiation in the symbolic ring::
sage: n = var('n')
sage: A = matrix([[pi, e],[0, -2*I]])
sage: (A^n).list()
[pi^n,
-(-2*I)^n/(pi*e^(-1) + 2*I*e^(-1)) + pi^n/(pi*e^(-1) + 2*I*e^(-1)),
0,
(-2*I)^n]
If the base ring is inexact, the Jordan normal form is not available::
sage: A = matrix(RDF, [[2, -1], [1, 0]])
sage: A^n
Traceback (most recent call last):
...
ValueError: Jordan normal form not implemented over inexact rings.
Testing exponentiation in the integer ring::
sage: A = matrix(ZZ, [[1,-1],[-1,1]])
sage: A^(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)]
Check if :trac:`23215` is fixed::
sage: a, b, k = var('a, b, k')
sage: (matrix(2, [a, b, -b, a])^k).list()
[1/2*(a + I*b)^k + 1/2*(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*I*(a - I*b)^k,
1/2*(a + I*b)^k + 1/2*(a - I*b)^k]
"""
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
got_SR = A.base_ring() == SR
# Change to QQbar if possible
try:
A = A.change_ring(QQbar)
except (TypeError, NotImplementedError):
pass
# Get Jordan matrix J and invertible matrix P such that A = P*J*~P
# From that, we will compute M = J^n, and obtain A^n = P*J^n*~P
J, P = A.jordan_form(transformation=True)
# Where each Jordan block starts, and number of blocks
block_start = [0] + J.subdivisions()[0]
num_blocks = len(block_start)
# Prepare matrix M to store `J^n`, computed by Jordan block
M = matrix(SR, J.ncols())
M.subdivide(J.subdivisions())
for k in range(num_blocks):
# Jordan block Jk, its dimension nk, the eigenvalue m
Jk = J.subdivision(k, k)
nk = Jk.ncols()
mk = Jk[0,0]
# First row of block Mk; its entries are of the form
# 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)]
# Form block Mk and insert it in M
Mk = matrix(SR, [[SR.zero()]*i + vk[:nk-i] for i in range(nk)])
M.set_block(block_start[k], block_start[k], Mk)
# Change entries of P and P^-1 into symbolic expressions
if not got_SR:
Pinv = (~P).apply_map(AlgebraicNumber.radical_expression)
P = P.apply_map(AlgebraicNumber.radical_expression)
else:
Pinv = ~P
return P * M * Pinv
</code></pre>
<p>I will open a ticket on Sage's issue tracker to get the fix in Sage.</p>
<p>Edit: The issue is now tracked at
<a href="https://trac.sagemath.org/ticket/25082">Sage Trac ticket #25082</a>.</p>
https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41867#post-id-41867Sorry it took me a while. The issue is now tracked at
[Sage Trac ticket #25082](https://trac.sagemath.org/ticket/25082).Mon, 02 Apr 2018 13:04:15 -0500https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41867#post-id-41867Comment by tmonteil for <p>Thank you for pointing out this problem.</p>
<p>The computation of symbolic powers of a matrix was recently
<a href="https://trac.sagemath.org/ticket/22523">introduced in Sage</a>, based on an idea
<a href="https://ask.sagemath.org/question/35658/general-power-of-a-matrix/">proposed on Ask Sage</a>:</p>
<ul>
<li>compute the Jordan form, and the transformation matrix</li>
<li>compute powers of each Jordan block</li>
<li>form a block-diagonal matrix with these powers</li>
<li>change back to the old basis using the transformation matrix</li>
</ul>
<p>The idea is good. Unfortunately there was a glitch in the
implementation: the powers of the Jordan blocks were being
inserted at the wrong place.</p>
<p>To illustrate the bug in an enlightening way, consider the matrix:</p>
<pre><code>sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]
</code></pre>
<p>The current implementation in Sage gives:</p>
<pre><code>sage: B = A^n; B
[ 2^n 2^(n - 1)*n 0]
[ 0 3^n 0]
[ 0 0 0]
</code></pre>
<p>Block k is inserted starting at position (k, k)
instead of its correct position.</p>
<p>The revised implementation, suggested below, gives:</p>
<pre><code>sage: B = _matrix_power_symbolic(A, n); B
[ 2^n 2^(n - 1)*n 0]
[ 0 2^n 0]
[ 0 0 3^n]
</code></pre>
<p>and can deal with the matrix in the question here.</p>
<pre><code>sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4 1 2]
[ 0 2 -4]
[ 0 1 6]
sage: B = _matrix_power_symbolic(A, n); B
[ 4^n 4^(n - 1)*n 2*4^(n - 1)*n]
[ 0 -2*4^(n - 1)*n + 4^n -4*4^(n - 1)*n]
[ 0 4^(n - 1)*n 2*4^(n - 1)*n + 4^n]
</code></pre>
<p>Here is a revised implementation with a fix for the glitch,
along with some efficiency improvements.</p>
<pre><code>def _matrix_power_symbolic(A, n):
r"""
Return the symbolic matrix power `A^n`
This function implements the computation of `A^n` for symbolic `n`,
relying on the Jordan normal form of `A`, available for exact rings
as :meth:`jordan_form`. See [Hig2008]_, §1.2, for further details.
INPUT:
- ``A`` -- a square matrix over an exact field
- ``n`` -- the symbolic exponent
OUTPUT:
The matrix `A^n` (with symbolic entries).
EXAMPLES:
Powers of two by two matrix::
sage: n = SR.var('n')
sage: A = matrix(QQ, [[2, -1], [1, 0]])
sage: B = A^n
sage: B
[ n + 1 -n]
[ n -n + 1]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
Powers of a three by three matrix in Jordan form::
sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]
sage: B = A^n; B
[ 2^n 2^(n - 1)*n 0]
[ 0 2^n 0]
[ 0 0 3^n]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
Powers of a three by three matrix not in Jordan form::
sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4 1 2]
[ 0 2 -4]
[ 0 1 6]
sage: B = A^n
sage: B
[ 4^n 4^(n - 1)*n 2*4^(n - 1)*n]
[ 0 -2*4^(n - 1)*n + 4^n -4*4^(n - 1)*n]
[ 0 4^(n - 1)*n 2*4^(n - 1)*n + 4^n]
sage: [B.subs({n: k}) for k in range(4)]
[
[1 0 0] [ 4 1 2] [ 16 8 16] [ 64 48 96]
[0 1 0] [ 0 2 -4] [ 0 0 -32] [ 0 -32 -192]
[0 0 1], [ 0 1 6], [ 0 8 32], [ 0 48 160]
]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
TESTS:
Testing exponentiation in the symbolic ring::
sage: n = var('n')
sage: A = matrix([[pi, e],[0, -2*I]])
sage: (A^n).list()
[pi^n,
-(-2*I)^n/(pi*e^(-1) + 2*I*e^(-1)) + pi^n/(pi*e^(-1) + 2*I*e^(-1)),
0,
(-2*I)^n]
If the base ring is inexact, the Jordan normal form is not available::
sage: A = matrix(RDF, [[2, -1], [1, 0]])
sage: A^n
Traceback (most recent call last):
...
ValueError: Jordan normal form not implemented over inexact rings.
Testing exponentiation in the integer ring::
sage: A = matrix(ZZ, [[1,-1],[-1,1]])
sage: A^(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)]
Check if :trac:`23215` is fixed::
sage: a, b, k = var('a, b, k')
sage: (matrix(2, [a, b, -b, a])^k).list()
[1/2*(a + I*b)^k + 1/2*(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*I*(a - I*b)^k,
1/2*(a + I*b)^k + 1/2*(a - I*b)^k]
"""
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
got_SR = A.base_ring() == SR
# Change to QQbar if possible
try:
A = A.change_ring(QQbar)
except (TypeError, NotImplementedError):
pass
# Get Jordan matrix J and invertible matrix P such that A = P*J*~P
# From that, we will compute M = J^n, and obtain A^n = P*J^n*~P
J, P = A.jordan_form(transformation=True)
# Where each Jordan block starts, and number of blocks
block_start = [0] + J.subdivisions()[0]
num_blocks = len(block_start)
# Prepare matrix M to store `J^n`, computed by Jordan block
M = matrix(SR, J.ncols())
M.subdivide(J.subdivisions())
for k in range(num_blocks):
# Jordan block Jk, its dimension nk, the eigenvalue m
Jk = J.subdivision(k, k)
nk = Jk.ncols()
mk = Jk[0,0]
# First row of block Mk; its entries are of the form
# 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)]
# Form block Mk and insert it in M
Mk = matrix(SR, [[SR.zero()]*i + vk[:nk-i] for i in range(nk)])
M.set_block(block_start[k], block_start[k], Mk)
# Change entries of P and P^-1 into symbolic expressions
if not got_SR:
Pinv = (~P).apply_map(AlgebraicNumber.radical_expression)
P = P.apply_map(AlgebraicNumber.radical_expression)
else:
Pinv = ~P
return P * M * Pinv
</code></pre>
<p>I will open a ticket on Sage's issue tracker to get the fix in Sage.</p>
<p>Edit: The issue is now tracked at
<a href="https://trac.sagemath.org/ticket/25082">Sage Trac ticket #25082</a>.</p>
https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41854#post-id-41854Which ticket is that ?Sat, 31 Mar 2018 19:37:54 -0500https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41854#post-id-41854Comment by kcrisman for <p>Thank you for pointing out this problem.</p>
<p>The computation of symbolic powers of a matrix was recently
<a href="https://trac.sagemath.org/ticket/22523">introduced in Sage</a>, based on an idea
<a href="https://ask.sagemath.org/question/35658/general-power-of-a-matrix/">proposed on Ask Sage</a>:</p>
<ul>
<li>compute the Jordan form, and the transformation matrix</li>
<li>compute powers of each Jordan block</li>
<li>form a block-diagonal matrix with these powers</li>
<li>change back to the old basis using the transformation matrix</li>
</ul>
<p>The idea is good. Unfortunately there was a glitch in the
implementation: the powers of the Jordan blocks were being
inserted at the wrong place.</p>
<p>To illustrate the bug in an enlightening way, consider the matrix:</p>
<pre><code>sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]
</code></pre>
<p>The current implementation in Sage gives:</p>
<pre><code>sage: B = A^n; B
[ 2^n 2^(n - 1)*n 0]
[ 0 3^n 0]
[ 0 0 0]
</code></pre>
<p>Block k is inserted starting at position (k, k)
instead of its correct position.</p>
<p>The revised implementation, suggested below, gives:</p>
<pre><code>sage: B = _matrix_power_symbolic(A, n); B
[ 2^n 2^(n - 1)*n 0]
[ 0 2^n 0]
[ 0 0 3^n]
</code></pre>
<p>and can deal with the matrix in the question here.</p>
<pre><code>sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4 1 2]
[ 0 2 -4]
[ 0 1 6]
sage: B = _matrix_power_symbolic(A, n); B
[ 4^n 4^(n - 1)*n 2*4^(n - 1)*n]
[ 0 -2*4^(n - 1)*n + 4^n -4*4^(n - 1)*n]
[ 0 4^(n - 1)*n 2*4^(n - 1)*n + 4^n]
</code></pre>
<p>Here is a revised implementation with a fix for the glitch,
along with some efficiency improvements.</p>
<pre><code>def _matrix_power_symbolic(A, n):
r"""
Return the symbolic matrix power `A^n`
This function implements the computation of `A^n` for symbolic `n`,
relying on the Jordan normal form of `A`, available for exact rings
as :meth:`jordan_form`. See [Hig2008]_, §1.2, for further details.
INPUT:
- ``A`` -- a square matrix over an exact field
- ``n`` -- the symbolic exponent
OUTPUT:
The matrix `A^n` (with symbolic entries).
EXAMPLES:
Powers of two by two matrix::
sage: n = SR.var('n')
sage: A = matrix(QQ, [[2, -1], [1, 0]])
sage: B = A^n
sage: B
[ n + 1 -n]
[ n -n + 1]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
Powers of a three by three matrix in Jordan form::
sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]
sage: B = A^n; B
[ 2^n 2^(n - 1)*n 0]
[ 0 2^n 0]
[ 0 0 3^n]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
Powers of a three by three matrix not in Jordan form::
sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4 1 2]
[ 0 2 -4]
[ 0 1 6]
sage: B = A^n
sage: B
[ 4^n 4^(n - 1)*n 2*4^(n - 1)*n]
[ 0 -2*4^(n - 1)*n + 4^n -4*4^(n - 1)*n]
[ 0 4^(n - 1)*n 2*4^(n - 1)*n + 4^n]
sage: [B.subs({n: k}) for k in range(4)]
[
[1 0 0] [ 4 1 2] [ 16 8 16] [ 64 48 96]
[0 1 0] [ 0 2 -4] [ 0 0 -32] [ 0 -32 -192]
[0 0 1], [ 0 1 6], [ 0 8 32], [ 0 48 160]
]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
TESTS:
Testing exponentiation in the symbolic ring::
sage: n = var('n')
sage: A = matrix([[pi, e],[0, -2*I]])
sage: (A^n).list()
[pi^n,
-(-2*I)^n/(pi*e^(-1) + 2*I*e^(-1)) + pi^n/(pi*e^(-1) + 2*I*e^(-1)),
0,
(-2*I)^n]
If the base ring is inexact, the Jordan normal form is not available::
sage: A = matrix(RDF, [[2, -1], [1, 0]])
sage: A^n
Traceback (most recent call last):
...
ValueError: Jordan normal form not implemented over inexact rings.
Testing exponentiation in the integer ring::
sage: A = matrix(ZZ, [[1,-1],[-1,1]])
sage: A^(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)]
Check if :trac:`23215` is fixed::
sage: a, b, k = var('a, b, k')
sage: (matrix(2, [a, b, -b, a])^k).list()
[1/2*(a + I*b)^k + 1/2*(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*I*(a - I*b)^k,
1/2*(a + I*b)^k + 1/2*(a - I*b)^k]
"""
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
got_SR = A.base_ring() == SR
# Change to QQbar if possible
try:
A = A.change_ring(QQbar)
except (TypeError, NotImplementedError):
pass
# Get Jordan matrix J and invertible matrix P such that A = P*J*~P
# From that, we will compute M = J^n, and obtain A^n = P*J^n*~P
J, P = A.jordan_form(transformation=True)
# Where each Jordan block starts, and number of blocks
block_start = [0] + J.subdivisions()[0]
num_blocks = len(block_start)
# Prepare matrix M to store `J^n`, computed by Jordan block
M = matrix(SR, J.ncols())
M.subdivide(J.subdivisions())
for k in range(num_blocks):
# Jordan block Jk, its dimension nk, the eigenvalue m
Jk = J.subdivision(k, k)
nk = Jk.ncols()
mk = Jk[0,0]
# First row of block Mk; its entries are of the form
# 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)]
# Form block Mk and insert it in M
Mk = matrix(SR, [[SR.zero()]*i + vk[:nk-i] for i in range(nk)])
M.set_block(block_start[k], block_start[k], Mk)
# Change entries of P and P^-1 into symbolic expressions
if not got_SR:
Pinv = (~P).apply_map(AlgebraicNumber.radical_expression)
P = P.apply_map(AlgebraicNumber.radical_expression)
else:
Pinv = ~P
return P * M * Pinv
</code></pre>
<p>I will open a ticket on Sage's issue tracker to get the fix in Sage.</p>
<p>Edit: The issue is now tracked at
<a href="https://trac.sagemath.org/ticket/25082">Sage Trac ticket #25082</a>.</p>
https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41799#post-id-41799Yes, great detective work. Sorry for confusing the issue earlier since I thought it was just that an error was raised.Tue, 27 Mar 2018 20:56:44 -0500https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41799#post-id-41799Comment by mforets for <p>Thank you for pointing out this problem.</p>
<p>The computation of symbolic powers of a matrix was recently
<a href="https://trac.sagemath.org/ticket/22523">introduced in Sage</a>, based on an idea
<a href="https://ask.sagemath.org/question/35658/general-power-of-a-matrix/">proposed on Ask Sage</a>:</p>
<ul>
<li>compute the Jordan form, and the transformation matrix</li>
<li>compute powers of each Jordan block</li>
<li>form a block-diagonal matrix with these powers</li>
<li>change back to the old basis using the transformation matrix</li>
</ul>
<p>The idea is good. Unfortunately there was a glitch in the
implementation: the powers of the Jordan blocks were being
inserted at the wrong place.</p>
<p>To illustrate the bug in an enlightening way, consider the matrix:</p>
<pre><code>sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]
</code></pre>
<p>The current implementation in Sage gives:</p>
<pre><code>sage: B = A^n; B
[ 2^n 2^(n - 1)*n 0]
[ 0 3^n 0]
[ 0 0 0]
</code></pre>
<p>Block k is inserted starting at position (k, k)
instead of its correct position.</p>
<p>The revised implementation, suggested below, gives:</p>
<pre><code>sage: B = _matrix_power_symbolic(A, n); B
[ 2^n 2^(n - 1)*n 0]
[ 0 2^n 0]
[ 0 0 3^n]
</code></pre>
<p>and can deal with the matrix in the question here.</p>
<pre><code>sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4 1 2]
[ 0 2 -4]
[ 0 1 6]
sage: B = _matrix_power_symbolic(A, n); B
[ 4^n 4^(n - 1)*n 2*4^(n - 1)*n]
[ 0 -2*4^(n - 1)*n + 4^n -4*4^(n - 1)*n]
[ 0 4^(n - 1)*n 2*4^(n - 1)*n + 4^n]
</code></pre>
<p>Here is a revised implementation with a fix for the glitch,
along with some efficiency improvements.</p>
<pre><code>def _matrix_power_symbolic(A, n):
r"""
Return the symbolic matrix power `A^n`
This function implements the computation of `A^n` for symbolic `n`,
relying on the Jordan normal form of `A`, available for exact rings
as :meth:`jordan_form`. See [Hig2008]_, §1.2, for further details.
INPUT:
- ``A`` -- a square matrix over an exact field
- ``n`` -- the symbolic exponent
OUTPUT:
The matrix `A^n` (with symbolic entries).
EXAMPLES:
Powers of two by two matrix::
sage: n = SR.var('n')
sage: A = matrix(QQ, [[2, -1], [1, 0]])
sage: B = A^n
sage: B
[ n + 1 -n]
[ n -n + 1]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
Powers of a three by three matrix in Jordan form::
sage: n = SR.var('n')
sage: A = matrix(QQ, 3, [[2, 1, 0], [0, 2, 0], [0, 0, 3]])
sage: A
[2 1 0]
[0 2 0]
[0 0 3]
sage: B = A^n; B
[ 2^n 2^(n - 1)*n 0]
[ 0 2^n 0]
[ 0 0 3^n]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
Powers of a three by three matrix not in Jordan form::
sage: A = matrix([[4, 1, 2], [0, 2, -4], [0, 1, 6]])
sage: A
[ 4 1 2]
[ 0 2 -4]
[ 0 1 6]
sage: B = A^n
sage: B
[ 4^n 4^(n - 1)*n 2*4^(n - 1)*n]
[ 0 -2*4^(n - 1)*n + 4^n -4*4^(n - 1)*n]
[ 0 4^(n - 1)*n 2*4^(n - 1)*n + 4^n]
sage: [B.subs({n: k}) for k in range(4)]
[
[1 0 0] [ 4 1 2] [ 16 8 16] [ 64 48 96]
[0 1 0] [ 0 2 -4] [ 0 0 -32] [ 0 -32 -192]
[0 0 1], [ 0 1 6], [ 0 8 32], [ 0 48 160]
]
sage: all(A^k == B.subs({n: k}) for k in range(8))
True
TESTS:
Testing exponentiation in the symbolic ring::
sage: n = var('n')
sage: A = matrix([[pi, e],[0, -2*I]])
sage: (A^n).list()
[pi^n,
-(-2*I)^n/(pi*e^(-1) + 2*I*e^(-1)) + pi^n/(pi*e^(-1) + 2*I*e^(-1)),
0,
(-2*I)^n]
If the base ring is inexact, the Jordan normal form is not available::
sage: A = matrix(RDF, [[2, -1], [1, 0]])
sage: A^n
Traceback (most recent call last):
...
ValueError: Jordan normal form not implemented over inexact rings.
Testing exponentiation in the integer ring::
sage: A = matrix(ZZ, [[1,-1],[-1,1]])
sage: A^(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)]
Check if :trac:`23215` is fixed::
sage: a, b, k = var('a, b, k')
sage: (matrix(2, [a, b, -b, a])^k).list()
[1/2*(a + I*b)^k + 1/2*(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*I*(a - I*b)^k,
1/2*(a + I*b)^k + 1/2*(a - I*b)^k]
"""
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
got_SR = A.base_ring() == SR
# Change to QQbar if possible
try:
A = A.change_ring(QQbar)
except (TypeError, NotImplementedError):
pass
# Get Jordan matrix J and invertible matrix P such that A = P*J*~P
# From that, we will compute M = J^n, and obtain A^n = P*J^n*~P
J, P = A.jordan_form(transformation=True)
# Where each Jordan block starts, and number of blocks
block_start = [0] + J.subdivisions()[0]
num_blocks = len(block_start)
# Prepare matrix M to store `J^n`, computed by Jordan block
M = matrix(SR, J.ncols())
M.subdivide(J.subdivisions())
for k in range(num_blocks):
# Jordan block Jk, its dimension nk, the eigenvalue m
Jk = J.subdivision(k, k)
nk = Jk.ncols()
mk = Jk[0,0]
# First row of block Mk; its entries are of the form
# 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)]
# Form block Mk and insert it in M
Mk = matrix(SR, [[SR.zero()]*i + vk[:nk-i] for i in range(nk)])
M.set_block(block_start[k], block_start[k], Mk)
# Change entries of P and P^-1 into symbolic expressions
if not got_SR:
Pinv = (~P).apply_map(AlgebraicNumber.radical_expression)
P = P.apply_map(AlgebraicNumber.radical_expression)
else:
Pinv = ~P
return P * M * Pinv
</code></pre>
<p>I will open a ticket on Sage's issue tracker to get the fix in Sage.</p>
<p>Edit: The issue is now tracked at
<a href="https://trac.sagemath.org/ticket/25082">Sage Trac ticket #25082</a>.</p>
https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41795#post-id-41795Thanks for finding and fixing the bug. Feel free to CC me when you create the ticket.Tue, 27 Mar 2018 13:48:55 -0500https://ask.sagemath.org/question/41622/bug-in-general-power-of-a-matrix/?comment=41795#post-id-41795