Matrix-exponentiation operator

Most programming languages have a built-in implementation of exponentiation for integers and reals only.

Demonstrate how to implement matrix exponentiation as an operator.

Ada

This is a generic solution for any natural power exponent. It will work with any type that has +,*, additive and multiplicative 0s. The implementation factors out powers A^2ⁿ: <lang ada> with Ada.Text_IO; use Ada.Text_IO;

procedure Test_Matrix is

  generic
     type Element is private;
     Zero : Element;
     One  : Element;
     with function "+" (A, B : Element) return Element is <>;
     with function "*" (A, B : Element) return Element is <>;
     with function Image (X : Element) return String is <>;
  package Matrices is
     type Matrix is array (Integer range <>, Integer range <>) of Element;
     function "*" (A, B : Matrix) return Matrix;
     function "**" (A : Matrix; Power : Natural) return Matrix;
     procedure Put (A : Matrix);
  end Matrices;

  package body Matrices is
     function "*" (A, B : Matrix) return Matrix is
        R   : Matrix (A'Range (1), B'Range (2));
        Sum : Element := Zero;
     begin
        for I in R'Range (1) loop
           for J in R'Range (2) loop
              Sum := Zero;
              for K in A'Range (2) loop
                 Sum := Sum + A (I, K) * B (K, J);
              end loop;
              R (I, J) := Sum;
           end loop;
        end loop;
        return R;
     end "*";

     function "**" (A : Matrix; Power : Natural) return Matrix is
     begin
        if Power = 1 then
           return A;
        end if;
        declare
           R : Matrix (A'Range (1), A'Range (2)) := (others => (others => Zero));
           P : Matrix  := A;
           E : Natural := Power;
        begin
           for I in P'Range (1) loop -- R is identity matrix
              R (I, I) := One;
           end loop;
           if E = 0 then
              return R;
           end if;
           loop
              if E mod 2 /= 0 then
                 R := R * P;
              end if;
              E := E / 2;
              exit when E = 0;
              P := P * P;
           end loop;
           return R;
        end;
     end "**";
     
     procedure Put (A : Matrix) is
     begin
        for I in A'Range (1) loop
           for J in A'Range (1) loop
              Put (Image (A (I, J)));
           end loop;
           New_Line;
        end loop;
     end Put;
  end Matrices;
  
  package Integer_Matrices is new Matrices (Integer, 0, 1, Image => Integer'Image);
  use Integer_Matrices;
  
  M : Matrix (1..2, 1..2) := ((3,2),(2,1));

begin

  Put_Line ("M =");       Put (M);
  Put_Line ("M**0 =");    Put (M**0);
  Put_Line ("M**1 =");    Put (M**1);
  Put_Line ("M**2 =");    Put (M**2);
  Put_Line ("M*M =");     Put (M*M);
  Put_Line ("M**3 =");    Put (M**3);
  Put_Line ("M*M*M =");   Put (M*M*M);
  Put_Line ("M**4 =");    Put (M**3);
  Put_Line ("M*M*M*M ="); Put (M*M*M*M);
  Put_Line ("M**10 =");   Put (M**10);
  Put_Line ("M*M*M*M*M*M*M*M*M*M ="); Put (M*M*M*M*M*M*M*M*M*M);

end Test_Matrix; </lang> Sample output:

M =
 3 2
 2 1
M**0 =
 1 0
 0 1
M**1 =
 3 2
 2 1
M**2 =
 13 8
 8 5
M*M =
 13 8
 8 5
M**3 =
 55 34
 34 21
M*M*M =
 55 34
 34 21
M**4 =
 55 34
 34 21
M*M*M*M =
 233 144
 144 89
M**10 =
 1346269 832040
 832040 514229
M*M*M*M*M*M*M*M*M*M =
 1346269 832040
 832040 514229

The following program implements exponentiation of a square Hermitian complex matrix by any complex power. The limitation to be Hermitian is not essential and comes for the limitation of the standard Ada linear algebra library. <lang ada> with Ada.Text_IO; use Ada.Text_IO; with Ada.Complex_Text_IO; use Ada.Complex_Text_IO; with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types; with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays; with Ada.Numerics.Complex_Arrays; use Ada.Numerics.Complex_Arrays; with Ada.Numerics.Complex_Elementary_Functions; use Ada.Numerics.Complex_Elementary_Functions;

procedure Test_Matrix is

  function "**" (A : Complex_Matrix; Power : Complex) return Complex_Matrix is
     L  : Real_Vector (A'Range (1));
     X  : Complex_Matrix (A'Range (1), A'Range (2));
     R  : Complex_Matrix (A'Range (1), A'Range (2));
     RL : Complex_Vector (A'Range (1));
  begin
     Eigensystem (A, L, X);
     for I in L'Range loop
        RL (I) := (L (I), 0.0) ** Power;
     end loop;
     for I in R'Range (1) loop
        for J in R'Range (2) loop
           declare
              Sum : Complex := (0.0, 0.0);
           begin
              for K in RL'Range (1) loop
                 Sum := Sum + X (K, I) * RL (K) * X (K, J);
              end loop;
              R (I, J) := Sum;
           end;
        end loop;
     end loop;
     return R;
  end "**";
  procedure Put (A : Complex_Matrix) is
  begin
     for I in A'Range (1) loop
       for J in A'Range (1) loop
          Put (A (I, J));
       end loop;
       New_Line;
     end loop;
  end Put;
  M : Complex_Matrix (1..2, 1..2) := (((3.0,0.0),(2.0,1.0)),((2.0,-1.0),(1.0,0.0)));

begin

  Put_Line ("M =");      Put (M);
  Put_Line ("M**0 =");   Put (M**(0.0,0.0));
  Put_Line ("M**1 =");   Put (M**(1.0,0.0));
  Put_Line ("M**0.5 ="); Put (M**(0.5,0.0));

end Test_Matrix; </lang> This solution is not tested, because the available version of GNAT GPL Ada compiler (20070405-41) does not provide an implementation of the standard library.

ALGOL 68

main:(

  INT default upb=3;
  MODE VEC = [default upb]COMPL;
  MODE MAT = [default upb,default upb]COMPL;
 
  OP * = (VEC a,b)COMPL: (
      COMPL result:=0;
      FOR i FROM LWB a TO UPB a DO result+:= a[i]*b[i] OD;
      result
    );
 
  OP * = (VEC a, MAT b)VEC: ( # overload vec times matrix #
      [2 LWB b:2 UPB b]COMPL result;
      FOR j FROM 2 LWB b TO 2 UPB b DO result[j]:=a*b[,j] OD;
      result
    );
 
  OP * = (MAT a, b)MAT: ( # overload matrix times matrix #
      [LWB a:UPB a, 2 LWB b:2 UPB b]COMPL result;
      FOR k FROM LWB result TO UPB result DO result[k,]:=a[k,]*b OD;
      result
    );

  OP IDENTITY = (INT upb)MAT:(
    [upb,upb] COMPL out;
    FOR i TO upb DO 
      FOR j TO upb DO
        out[i,j]:= ( i=j |1|0)
      OD
    OD;
    out
  );

   # This is the task part. #
   OP ** = (MAT base, INT exponent)MAT: (
     BITS binary exponent:=BIN exponent ;
     MAT out := IF bits width ELEM binary exponent THEN base ELSE IDENTITY UPB base FI;
     MAT sq:=base;
 
     WHILE 
       binary exponent := binary exponent SHR 1;
       binary exponent /= BIN 0 
     DO
       sq := sq * sq; 
       IF bits width ELEM binary exponent THEN out := out * sq FI
     OD;
     out
   );

  PROC compl matrix printf= (FORMAT compl fmt, MAT m)VOID:(
    FORMAT vec fmt = $"("n(2 UPB m-1)(f(compl fmt)",")f(compl fmt)")"$;
    FORMAT matrix fmt = $x"("n(UPB m-1)(f(vec fmt)","lxx)f(vec fmt)");"$;
    # finally print the result #
    printf((matrix fmt,m))
  );
 
  FORMAT compl fmt = $-z.z,+z.z"i"$; # width of 4, with no leading '+' sign, 1 decimals #
  MAT matrix=((sqrt(0.5)I0         , sqrt(0.5)I0        , 0I0),
              (        0I-sqrt(0.5),         0Isqrt(0.5), 0I0),
              (        0I0         ,         0I0        , 0I1));

  printf(($" matrix ** "g(0)":"l$,24));
  compl matrix printf(compl fmt, matrix**24); print(newline)
)

Output:

matrix ** 24:
(( 1.0+0.0i, 0.0+0.0i, 0.0+0.0i),
 ( 0.0+0.0i, 1.0+0.0i, 0.0+0.0i),
 ( 0.0+0.0i, 0.0+0.0i, 1.0+0.0i));

C++

This is an implementation in C++. <lang cpp>#include <complex>

include <cmath>
include <iostream>

using namespace std;

template<int MSize = 3, class T = complex<double> > class SqMx {

 typedef T Ax[MSize][MSize];
 typedef SqMx<MSize, T> Mx;

private:

 Ax a;
 SqMx() { }

public:

 SqMx(const Ax &_a) { // constructor with pre-defined array
   for (int r = 0; r < MSize; r++)
     for (int c = 0; c < MSize; c++)
       a[r][c] = _a[r][c];
 }

 static Mx identity() {
   Mx m;
   for (int r = 0; r < MSize; r++)
     for (int c = 0; c < MSize; c++)
       m.a[r][c] = (r == c ? 1 : 0);
   return m;
 }

 friend ostream &operator<<(ostream& os, const Mx &p)
 { // ugly print
   for (int i = 0; i < MSize; i++) {
     for (int j = 0; j < MSize; j++)
       os << p.a[i][j] << ",";
     os << endl;
   }
   return os;
 }

 Mx operator*(const Mx &b) {
   Mx d;
   for (int r = 0; r < MSize; r++)
     for (int c = 0; c < MSize; c++) {
       d.a[r][c] = 0;
       for (int k = 0; k < MSize; k++)
         d.a[r][c] += a[r][k] * b.a[k][c];
     }
   return d;
 }</lang>

This is the task part. <lang cpp>

 // C++ does not have a ** operator, instead, ^ (bitwise Xor) is used.
 Mx operator^(int n) {
   if (n < 0)
     throw "Negative exponent not implemented";

   Mx d = identity();
   for (Mx sq = *this; n > 0; sq = sq * sq, n /= 2)
     if (n % 2 != 0)
       d = d * sq;
   return d;
 }

};

typedef SqMx<> M3; typedef complex<double> creal;

int main() {

 double q = sqrt(0.5);
 creal array[3][3] =
   {{creal(q,  0), creal(q, 0), creal(0, 0)},
    {creal(0, -q), creal(0, q), creal(0, 0)},
    {creal(0,  0), creal(0, 0), creal(0, 1)}};
 M3 m(array);

 cout << "m ^ 23=" << endl
      << (m ^ 23) << endl;

 return 0;

}</lang> Output:

m ^ 23=
(0.707107,0),(0,0.707107),(0,0),
(0.707107,0),(0,-0.707107),(0,0),
(0,0),(0,0),(0,-1),

An alternative way would be to implement operator*= and conversion from number (giving multiples of the identity matrix) for the matrix and use the generic code from Exponentiation operator#C++ with support for negative exponents removed (or alternatively, implement matrix inversion as well, implement /= in terms of it, and use the generic code unchanged). Note that the algorithm used there is much faster as well.

D

This is a implementation by D. <lang d>module mxpow ; import std.stdio ; import std.string ; import std.math ;

struct SqMx(int MSize = 3, T = creal) {

 alias T[MSize][MSize] Ax ;
 alias SqMx!(MSize, T) Mx ;
 static string fmt = "%8.3f" ;
 private Ax a ;
 static Mx opCall(Ax a){
   Mx m ;
   m.a[] = a[] ;
   return m ;
 }
 static Mx Identity() {
   Mx m ;
   for(int r = 0; r < MSize ; r++)
     for(int c = 0 ; c < MSize ; c++)
       m.a[r][c] = cast(T) (r == c ? 1 : 0) ;
   return m ;
 }
 string toString() { // pretty print
   string[MSize] s, t ;
   foreach(i, r; a) {
     foreach (j , e ; r)   
       s[j] = format(fmt, e) ;
     t[i] = join(s, ",") ;
   }
   return "<" ~ join(t,"\n ") ~ ">" ;
 }     
 Mx opMul(Mx b) {
   Mx d ;
   for(int r = 0 ; r < MSize ; r++)
     for(int c = 0 ; c < MSize ; c++) {
       d.a[r][c] = cast(T) 0 ;
       for(int k = 0 ; k < MSize ; k++)
         d.a[r][c] += a[r][k]*b.a[k][c] ;
     }
   return d ;
 }

</lang> This is the task part. <lang d>

 // D does not have a ** operator, instead, ^ (bitwise Xor) is used.
 Mx opXor(int n){
   Mx d , sq ;

   if(n < 0)
     throw new Exception("Negative exponent not implemented") ;

   sq.a[] = this.a[] ; 
   d = Mx.Identity ;
   for( ;n > 0 ; sq = sq * sq, n >>= 1)
     if (n & 1)
       d = d * sq ;
   return d ;
 } 
 alias opXor pow ;

}

alias SqMx!() M3 ;

void main() {

 real q = sqrt(0.5) ;
 M3 m = M3(cast(M3.Ax)
           [   q + 0*1.0Li,   q + 0*1.0Li, 0.0L + 0.0Li,
            0.0L - q*1.0Li,0.0L + q*1.0Li, 0.0L + 0.0Li,
            0.0L +   0.0Li,0.0L +   0.0Li, 0.0L + 1.0Li]) ;
 M3.fmt = "%5.2f" ;
 writefln("m ^ 23 =\n", m.pow(23)) ;  
 writefln("m ^ 24 =\n", m ^ 24) ;

}</lang> Output:

m ^ 23 =
< 0.71+ 0.00i, 0.00+ 0.71i, 0.00+ 0.00i
  0.71+ 0.00i, 0.00+-0.71i, 0.00+ 0.00i
  0.00+ 0.00i, 0.00+ 0.00i, 0.00+-1.00i>
m ^ 24 =
< 1.00+ 0.00i, 0.00+ 0.00i, 0.00+ 0.00i
  0.00+ 0.00i, 1.00+ 0.00i, 0.00+ 0.00i
  0.00+ 0.00i, 0.00+ 0.00i, 1.00+ 0.00i>

NOTE: In D, the commutativity of binary operator to be overloading is preset. For instance, arithmetic + * , bitwise & ^ | operators are commutative, - / % >> << >>> is non-commutative.
The exponential operator ^ chose in previous code happened to be commutative, which allow expression like 24 ^ m to be legal. If such expression is not allowed, either a non-commutative operator should be chose, or implement a corresponding opXXX_r overloading that may throw static assert/error.
Details see Operator Loading in D

Haskell

Instead of writing it directly, we can re-use the overloaded exponentiation operator if we declare matrices as an instance of Num, using matrix multiplication (and addition). For simplicity, we use the inefficient representation as list of lists. Note that we don't check the dimensions (there are several ways to do that on the type-level, for example with phantom types).

<lang haskell>import Data.List

xs <+> ys = zipWith (+) xs ys xs <*> ys = sum $ zipWith (*) xs ys

newtype Mat a = Mat {unMat :: a} deriving Eq

instance Show a => Show (Mat a) where

 show xm = "Mat " ++ show (unMat xm)

instance Num a => Num (Mat a) where

 negate xm = Mat $ map (map negate) $ unMat xm
 xm + ym   = Mat $ zipWith (<+>) (unMat xm) (unMat ym)
 xm * ym   = Mat [[xs <*> ys | ys <- transpose $ unMat ym] | xs <- unMat xm]
 fromInteger n = Mat fromInteger n

</lang>

Output:

*Main> Mat [[1,2],[0,1]]^4
Mat [[1,8],[0,1]]

This will work for matrices over any numeric type, including complex numbers. The implementation of (^) uses the fast binary algorithm for exponentiation.

Mathematica

In Mathematica there is an distinction between powering elements wise and as a matrix. So m^2 will give m with each element squared. To do matrix exponentation we use the function MatrixPower. It can handle all types of numbers for the power (integers, floats, rationals, complex) but also symbols for the power, and all types for the matrix (numbers, symbols et cetera), and will always keep the result exact if the matrix and the exponent is exact. <lang Mathematica>

a = {{3, 2}, {4, 1}};
MatrixPower[a, 0]
MatrixPower[a, 1]
MatrixPower[a, -1]
MatrixPower[a, 4]
MatrixPower[a, 1/2]
MatrixPower[a, Pi]

</lang> gives back:

$\left({\begin{array}{cc}1&0\\0&1\end{array}}\right)$

$\left({\begin{array}{cc}3&2\\4&1\end{array}}\right)$

$\left({\begin{array}{cc}-{\frac {1}{5}}&{\frac {2}{5}}\\{\frac {4}{5}}&-{\frac {3}{5}}\end{array}}\right)$

$\left({\begin{array}{cc}417&208\\416&209\end{array}}\right)$

$\left({\begin{array}{cc}{\frac {2{\sqrt {5}}}{3}}+{\frac {i}{3}}&{\frac {\sqrt {5}}{3}}-{\frac {i}{3}}\\{\frac {2{\sqrt {5}}}{3}}-{\frac {2i}{3}}&{\frac {\sqrt {5}}{3}}+{\frac {2i}{3}}\end{array}}\right)$

$\left({\begin{array}{cc}{\frac {(-1)^{\pi }}{3}}+2{\frac {5^{\pi }}{3}}&{\frac {5^{\pi }}{3}}-{\frac {1}{3}}(-1)^{\pi }\\2{\frac {5^{\pi }}{3}}-{\frac {2}{3}}(-1)^{\pi }&{\frac {2(-1)^{\pi }}{3}}+{\frac {5^{\pi }}{3}}\end{array}}\right)$

Symbolic matrices like {{i,j},{k,l}} to the power m give general solutions for all possible i,j,k,l, and m: <lang Mathematica>

MatrixPower[{{i, j}, {k, l}}, m] // Simplify

</lang> gives back (note that the simplification is not necessary for the evaluation, it just gives a shorter output):

$\left({\begin{array}{cc}{\frac {2^{-m-1}\left(\left({\sqrt {i^{2}-2il+4jk+l^{2}}}-i+l\right)\left(-{\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}+\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i-l\right)\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}\right)}{\sqrt {i^{2}-2il+4jk+l^{2}}}}&{\frac {j2^{-m}\left(\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}-\left(-{\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}\right)}{\sqrt {i^{2}-2il+4jk+l^{2}}}}\\{\frac {k2^{-m}\left(\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}-\left(-{\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}\right)}{\sqrt {i^{2}-2il+4jk+l^{2}}}}&{\frac {2^{-m-1}\left(\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i-l\right)\left(-{\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}+\left({\sqrt {i^{2}-2il+4jk+l^{2}}}-i+l\right)\left({\sqrt {i^{2}-2il+4jk+l^{2}}}+i+l\right)^{m}\right)}{\sqrt {i^{2}-2il+4jk+l^{2}}}}\end{array}}\right)$

Final note: Do not confuse MatrixPower with MatrixExp; the former is for matrix exponentiation, and the latter for the matrix exponential (E^m).

MATLAB

For exponents in the form of A*A*A*A*...*A, A must be a square matrix: <lang Matlab>function [output] = matrixexponentiation(matrixA, exponent)

  output = matrixA^(exponent);</lang>

Otherwise, to take the individual array elements to the power of an exponent (the matrix need not be square): <lang Matlab>function [output] = matrixexponentiation(matrixA, exponent)

  output = matrixA.^(exponent);</lang>

Octave

Of course GNU Octave handles matrix and operations on matrix "naturally".

Output:

ans =

   1   0
   0   1

ans =

   3   2
   2   1

ans =

   13    8
    8    5

ans =

  -1.0000   2.0000
   2.0000  -3.0000

ans =

   1.48931 + 0.13429i   0.92044 - 0.21729i
   0.92044 - 0.21729i   0.56886 + 0.35158i

(Of course this is not an implementation, but it can be used as reference for the results)

Python

Using matrixMul from Matrix multiplication#Python <lang python>>>> from operator import mul >>> def matrixMul(m1, m2):

 return map(
   lambda row:
     map(
       lambda *column:
         sum(map(mul, row, column)),
       *m2),
   m1)

>>> def identity(size): size = range(size) return [[(i==j)*1 for i in size] for j in size]

>>> def matrixExp(m, pow): assert pow>=0 and int(pow)==pow, "Only non-negative, integer powers allowed" accumulator = identity(len(m)) for i in range(pow): accumulator = matrixMul(accumulator, m) return accumulator

>>> def printtable(data): for row in data: print ' '.join('%-5s' % ('%s' % cell) for cell in row)

>>> m = [[3,2], [2,1]] >>> for i in range(5): print '\n%i:' % i printtable( matrixExp(m, i) )

0: 1 0 0 1

1: 3 2 2 1

2: 13 8 8 5

3: 55 34 34 21

4: 233 144 144 89 >>> printtable( matrixExp(m, 10) ) 1346269 832040 832040 514229 >>> </lang>

R

<lang R> library(Biodem) m <- matrix(c(3,2,2,1), nrow=2) mtx.exp(m, 0)

[,1] [,2]
[1,] 1 0
[2,] 0 1

mtx.exp(m, 1)

[,1] [,2]
[1,] 3 2
[2,] 2 1

mtx.exp(m, 2)

[,1] [,2]
[1,] 13 8
[2,] 8 5

mtx.exp(m, 3)

[,1] [,2]
[1,] 55 34
[2,] 34 21

mtx.exp(m, 10)

[,1] [,2]
[1,] 1346269 832040
[2,] 832040 514229

</lang> Note that non-integer powers are not supported with this function.

Ruby

The

Library: matrix.rb

class defines the exponentiation operator "**" like this [1]:

<lang ruby>class Matrix

 def ** (other)
   if other.kind_of?(Integer)
     x = self
     if other <= 0
       x = self.inverse
       return Matrix.identity(self.column_size) if other == 0
       other = -other
     end
     z = x
     n = other  - 1
     while n != 0
       while (div, mod = n.divmod(2)
              mod == 0)
         x = x * x
         n = div
       end
       z *= x
       n -= 1
     end
     z
   elsif other.kind_of?(Float) || defined?(Rational) && other.kind_of?(Rational)
     Matrix.Raise ErrOperationNotDefined, "**"
   else
     Matrix.Raise ErrOperationNotDefined, "**"
   end
 end

end</lang>

$ irb
irb(main):001:0> require 'matrix'
=> true
irb(main):002:0> m=Matrix[[3,2],[2,1]]
=> Matrix[[3, 2], [2, 1]]
irb(main):003:0> m**0
=> Matrix[[1, 0], [0, 1]]
irb(main):004:0> m ** 1
=> Matrix[[3, 2], [2, 1]]
irb(main):005:0> m ** 2
=> Matrix[[13, 8], [8, 5]]
irb(main):006:0> m ** 5
=> Matrix[[987, 610], [610, 377]]
irb(main):007:0> m ** 10
=> Matrix[[1346269, 832040], [832040, 514229]]
irb(main):008:0> m ** 1.5
ExceptionForMatrix::ErrOperationNotDefined: This operation(**) can't defined
        from /usr/lib/ruby/1.8/matrix.rb:665:in `**'
        from (irb):8

Tcl

Using code at Matrix multiplication#Tcl and Matrix Transpose#Tcl <lang tcl>package require Tcl 8.5 namespace path {::tcl::mathop ::tcl::mathfunc}

proc matrix_exp {m pow} {

   if { ! [string is int -strict $pow]} {
       error "non-integer exponents not implemented"
   }
   if {$pow < 0} {
       error "negative exponents not implemented"
   }
   lassign [size $m] rows cols
   # assume square matrix
   set temp [identity $rows]
   for {set n 1} {$n <= $pow} {incr n} {
       set temp [matrix_multiply $temp $m]
   }
   return $temp

}

proc identity {size} {

   set i [lrepeat $size [lrepeat $size 0]]
   for {set n 0} {$n < $size} {incr n} {lset i $n $n 1}
   return $i

}</lang>

% print_matrix [matrix_exp {{3 2} {2 1}} 1]
3 2 
2 1 
% print_matrix [matrix_exp {{3 2} {2 1}} 0]
1 0 
0 1 
% print_matrix [matrix_exp {{3 2} {2 1}} 2]
13 8 
 8 5 
% print_matrix [matrix_exp {{3 2} {2 1}} 3]
55 34 
34 21 
% print_matrix [matrix_exp {{3 2} {2 1}} 4]
233 144 
144  89 
% print_matrix [matrix_exp {{3 2} {2 1}} 10]
1346269 832040 
 832040 514229