# Revision history [back]

Thank you for pointing out this problem.

The computation of symbolic powers of a matrix was recently introduced in Sage, based on an idea proposed on Ask Sage:

• 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
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[:-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:
else:
Pinv = ~P

return P * M * Pinv


I will open a ticket on Sage's issue tracker to get the fix in Sage.

Thank you for pointing out this problem.

The computation of symbolic powers of a matrix was recently introduced in Sage, Sage, based on an idea idea proposed on Ask Sage:Sage:

• 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
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[:-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: