Matrix-exponentiation operator
You are encouraged to solve this task according to the task description, using any language you may know.
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 A2n: <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**4);
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 = 233 144 144 89 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
);
<lang algol68># 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
);</lang>
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
C doesn't support classes or allow operator overloading. The following is code that defines a function, SquareMtxPower that will raise a matrix to a positive integer power. <lang c>#include <math.h>
- include <stdio.h>
- include <stdlib.h>
typedef struct squareMtxStruct {
int dim; double *cells; double **m;
} *SquareMtx;
/* function for initializing row r of a new matrix */ typedef void (*FillFunc)( double *cells, int r, int dim, void *ff_data);
SquareMtx NewSquareMtx( int dim, FillFunc fillFunc, void *ff_data ) {
SquareMtx sm = malloc(sizeof(struct squareMtxStruct));
if (sm) {
int rw;
sm->dim = dim;
sm->cells = malloc(dim*dim * sizeof(double));
sm->m = malloc( dim * sizeof(double *));
if ((sm->cells != NULL) && (sm->m != NULL)) {
for (rw=0; rw<dim; rw++) {
sm->m[rw] = sm->cells + dim*rw;
fillFunc( sm->m[rw], rw, dim, ff_data );
}
}
else {
if (sm->m) free(sm->m);
if (sm->cells) free(sm->cells);
free(sm);
printf("Square Matrix allocation failure\n");
return NULL;
}
}
else {
printf("Malloc failed for square matrix\n");
}
return sm;
}
void ffMatxSquare( double *cells, int rw, int dim, SquareMtx m0 ) {
int col, ix;
double sum;
double *m0rw = m0->m[rw];
for (col = 0; col < dim; col++) {
sum = 0.0;
for (ix=0; ix<dim; ix++)
sum += m0rw[ix] * m0->m[ix][col];
cells[col] = sum;
}
}
void ffMatxMulply( double *cells, int rw, int dim, SquareMtx mplcnds[] ) {
SquareMtx mleft = mplcnds[0]; SquareMtx mrigt = mplcnds[1]; double sum; double *m0rw = mleft->m[rw]; int col, ix;
for (col = 0; col < dim; col++) {
sum = 0.0;
for (ix=0; ix<dim; ix++)
sum += m0rw[ix] * mrigt->m[ix][col];
cells[col] = sum;
}
}
void MatxMul( SquareMtx mr, SquareMtx left, SquareMtx rigt) {
int rw; SquareMtx mplcnds[2]; mplcnds[0] = left; mplcnds[1] = rigt;
for (rw = 0; rw < left->dim; rw++)
ffMatxMulply( mr->m[rw], rw, left->dim, mplcnds);
}
void ffIdentity( double *cells, int rw, int dim, void *v ) {
int col; for (col=0; col<dim; col++) cells[col] = 0.0; cells[rw] = 1.0;
} void ffCopy(double *cells, int rw, int dim, SquareMtx m1) {
int col; for (col=0; col<dim; col++) cells[col] = m1->m[rw][col];
}
void FreeSquareMtx( SquareMtx m ) {
free(m->m); free(m->cells); free(m);
}
SquareMtx SquareMtxPow( SquareMtx m0, int exp ) {
SquareMtx v0 = NewSquareMtx(m0->dim, ffIdentity, NULL); SquareMtx v1 = NULL; SquareMtx base0 = NewSquareMtx( m0->dim, ffCopy, m0); SquareMtx base1 = NULL; SquareMtx mplcnds[2], t;
while (exp) {
if (exp % 2) {
if (v1)
MatxMul( v1, v0, base0);
else {
mplcnds[0] = v0; mplcnds[1] = base0;
v1 = NewSquareMtx(m0->dim, ffMatxMulply, mplcnds);
}
{t = v0; v0=v1; v1 = t;}
}
if (base1)
MatxMul( base1, base0, base0);
else
base1 = NewSquareMtx( m0->dim, ffMatxSquare, base0);
t = base0; base0 = base1; base1 = t;
exp = exp/2;
}
if (base0) FreeSquareMtx(base0);
if (base1) FreeSquareMtx(base1);
if (v1) FreeSquareMtx(v1);
return v0;
}
FILE *fout; void SquareMtxPrint( SquareMtx mtx, const char *mn ) {
int rw, col; int d = mtx->dim;
fprintf(fout, "%s dim:%d =\n", mn, mtx->dim);
for (rw=0; rw<d; rw++) {
fprintf(fout, " |");
for(col=0; col<d; col++)
fprintf(fout, "%8.5f ",mtx->m[rw][col] );
fprintf(fout, " |\n");
}
fprintf(fout, "\n");
}
void fillInit( double *cells, int rw, int dim, void *data) {
double theta = 3.1415926536/6.0; double c1 = cos( theta); double s1 = sin( theta);
switch(rw) {
case 0:
cells[0]=c1; cells[1]=s1; cells[2]=0.0;
break;
case 1:
cells[0]=-s1; cells[1]=c1; cells[2]=0;
break;
case 2:
cells[0]=0.0; cells[1]=0.0; cells[2]=1.0;
break;
}
}
int main() {
SquareMtx m0 = NewSquareMtx( 3, fillInit, NULL); SquareMtx m1 = SquareMtxPow( m0, 5); SquareMtx m2 = SquareMtxPow( m0, 9); SquareMtx m3 = SquareMtxPow( m0, 2);
// fout = stdout;
fout = fopen("matrx_exp.txt", "w");
SquareMtxPrint(m0, "m0"); FreeSquareMtx(m0);
SquareMtxPrint(m1, "m0^5"); FreeSquareMtx(m1);
SquareMtxPrint(m2, "m0^9"); FreeSquareMtx(m2);
SquareMtxPrint(m3, "m0^2"); FreeSquareMtx(m3);
fclose(fout);
return 0;
}</lang> Output:
m0 dim:3 = | 0.86603 0.50000 0.00000 | |-0.50000 0.86603 0.00000 | | 0.00000 0.00000 1.00000 | m0^5 dim:3 = |-0.86603 0.50000 0.00000 | |-0.50000 -0.86603 0.00000 | | 0.00000 0.00000 1.00000 | m0^9 dim:3 = | 0.00000 -1.00000 0.00000 | | 1.00000 0.00000 0.00000 | | 0.00000 0.00000 1.00000 | m0^2 dim:3 = | 0.50000 0.86603 0.00000 | |-0.86603 0.50000 0.00000 | | 0.00000 0.00000 1.00000 |
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.
Common Lisp
This Common Lisp implementation uses 2D Arrays to represent matrices, and checks to make sure that the arrays are the right dimensions for multiplication and square for exponentiation. <lang lisp>(defun multiply-matrices (matrix-0 matrix-1)
"Takes two 2D arrays and returns their product, or an error if they cannot be multiplied"
(let* ((m0-dims (array-dimensions matrix-0))
(m1-dims (array-dimensions matrix-1))
(m0-dim (length m0-dims))
(m1-dim (length m1-dims)))
(if (or (/= 2 m0-dim) (/= 2 m1-dim))
(error "Array given not a matrix")
(let ((m0-rows (car m0-dims))
(m0-cols (cadr m0-dims))
(m1-rows (car m1-dims))
(m1-cols (cadr m1-dims)))
(if (/= m0-cols m1-rows)
(error "Incompatible dimensions")
(do ((rarr (make-array (list m0-rows m1-cols)
:initial-element 0) rarr)
(n 0 (if (= n (1- m0-cols)) 0 (1+ n)))
(cc 0 (if (= n (1- m0-cols))
(if (/= cc (1- m1-cols))
(1+ cc) 0) cc))
(cr 0 (if (and (= (1- m0-cols) n)
(= (1- m1-cols) cc))
(1+ cr)
cr)))
((= cr m0-rows) rarr)
(setf (aref rarr cr cc)
(+ (aref rarr cr cc)
(* (aref matrix-0 cr n)
(aref matrix-1 n cc))))))))))
(defun matrix-identity (dim)
"Creates a new identity matrix of size dim*dim"
(do ((rarr (make-array (list dim dim)
:initial-element 0) rarr)
(n 0 (1+ n)))
((= n dim) rarr)
(setf (aref rarr n n) 1)))
(defun matrix-expt (matrix exp)
"Takes the first argument (a matrix) and multiplies it by itself exp times"
(let* ((m-dims (array-dimensions matrix))
(m-rows (car m-dims))
(m-cols (cadr m-dims)))
(cond
((/= m-rows m-cols) (error "Non-square matrix"))
((zerop exp) (matrix-identity m-rows))
((= 1 exp) (do ((rarr (make-array (list m-rows m-cols)) rarr)
(cc 0 (if (= cc (1- m-cols))
0
(1+ cc)))
(cr 0 (if (= cc (1- m-cols))
(1+ cr)
cr)))
((= cr m-rows) rarr)
(setf (aref rarr cr cc) (aref matrix cr cc))))
((zerop (mod exp 2)) (let ((me2 (matrix-expt matrix (/ exp 2))))
(multiply-matrices me2 me2)))
(t (let ((me2 (matrix-expt matrix (/ (1- exp) 2))))
(multiply-matrices matrix (multiply-matrices me2 me2)))))))</lang>
Output (note that this lisp implementation uses single-precision floats for decimals by default). We can also use rationals:
CL-USER> (setf 5x5-matrix
(make-array '(5 5)
:initial-contents
'((0 1 -1 -2 2)
(0.4 4 3.2 -3 -10)
(4.5 -2 0.5 1 7)
(10 1 0 1.5 -2)
(4 5 -3 -2 1))))
#2A((0 1 -1 -2 2)
(0.4 4 3.2 -3 -10)
(4.5 -2 0.5 1 7)
(10 1 0 1.5 -2)
(4 5 -3 -2 1))
CL-USER> (matrix-expt 5x5-matrix 3)
#2A((-163.25 -19.5 92.25 -7.5999985 -184.3)
(156.6 -412.09998 0.7999954 331.45 597.4)
(-129.82501 401.25 -66.975 -302.55 -390.15)
(-148.9 39.25 -5.200001 -67.225006 -7.300003)
(-495.05 -231.5 310.85 33.0 -328.5))
CL-USER> (setf 4x4-matrix
(make-array '(4 4)
:initial-contents
'(( 1/2 -1/2 4 8)
(-3/4 7/3 8/5 -2)
(-5 17 20/3 -5/2)
( 3/2 -1 -7/3 6))))
#2A((1/2 -1/2 4 8) (-3/4 7/3 8/5 -2) (-5 17 20/3 -5/2) (3/2 -1 -7/3 6))
CL-USER> (matrix-expt 4x4-matrix 3)
#2A((-233/8 182723/720 757/30 353/6)
(-73517/480 838241/2160 77789/450 -67537/180)
(-5315/9 66493/45 90883/135 -54445/36)
(37033/144 -27374/45 -15515/54 12109/18))
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
Factor
There is already a built-in word (m^n) that implements exponentiation. Here is a simple and less efficient implementation.
<lang factor>USING: kernel math math.matrices sequences ;
- my-m^n ( m n -- m' )
dup 0 < [ "no negative exponents" throw ] [
[ drop length identity-matrix ]
[ swap '[ _ m. ] times ] 2bi
] if ;</lang>
( scratchpad ) { { 3 2 } { 2 1 } } 0 my-m^n .
{ { 1 0 } { 0 1 } }
( scratchpad ) { { 3 2 } { 2 1 } } 4 my-m^n .
{ { 233 144 } { 144 89 } }
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.
J
<lang j>mp=: +/ .* NB. Matrix multiplication pow=: 4 : 'mp&x^:y =i.#x'</lang>
or, from the j wiki, and faster for large exponents:
<lang j>pow=: 4 : 'mp/ mp~^:(I.|.#:y) x'</lang>
This implements an optimization where the exponent is represented in base 2, and repeated squaring is used to create a list of relevant powers of the base matrix, which are then combined using matrix multiplication. Note, however, that these two definitions treat a zero exponent differently.
Example use:
(3 2,:2 1) pow 3 55 34 34 21
JavaScript
for the print() and Array.forEach() functions.
Extends Matrix Transpose#JavaScript and Matrix multiplication#JavaScript <lang javascript>// IdentityMatrix is a "subclass" of Matrix function IdentityMatrix(n) {
this.height = n;
this.width = n;
this.mtx = [];
for (var i = 0; i < n; i++) {
this.mtx[i] = [];
for (var j = 0; j < n; j++) {
this.mtx[i][j] = (i == j ? 1 : 0);
}
}
} IdentityMatrix.prototype = Matrix.prototype;
// the Matrix exponentiation function // returns a new matrix Matrix.prototype.exp = function(n) {
var result = new IdentityMatrix(this.height);
for (var i = 1; i <= n; i++) {
result = result.mult(this);
}
return result;
}
var m = new Matrix([[3, 2], [2, 1]]); [0,1,2,3,4,10].forEach(function(e){print(m.exp(e)); print()})</lang> output
1,0 0,1 3,2 2,1 13,8 8,5 55,34 34,21 233,144 144,89 1346269,832040 832040,514229
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:
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):
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>
OCaml
We will use some auxiliary functions
<lang ocaml>(* identity matrix *) let eye n =
let a = Array.make_matrix n n 0.0 in for i=0 to n-1 do a.(i).(i) <- 1.0 done; (a)
(* matrix dimensions *) let dim a = Array.length a, Array.length a.(0);;
(* make matrix from list in row-major order *) let matrix p q v =
if (List.length v) <> (p * q)
then failwith "bad dimensions"
else
let a = Array.make_matrix p q (List.hd v) in
let rec g i j = function
| [] -> a
| x::v ->
a.(i).(j) <- x;
if j+1 < q
then g i (j+1) v
else g (i+1) 0 v
in
g 0 0 v
(* matrix product *) let matmul a b =
let n, p = dim a
and q, r = dim b in
if p <> q then failwith "bad dimensions" else
let c = Array.make_matrix n r 0.0 in
for i=0 to n-1 do
for j=0 to r-1 do
for k=0 to p-1 do
c.(i).(j) <- c.(i).(j) +. a.(i).(k) *. b.(k).(j)
done
done
done;
(c)
(* generic exponentiation, usual algorithm *) let pow one mul a n =
let rec g p x = function
| 0 -> x
| i ->
g (mul p p) (if i mod 2 = 1 then mul p x else x) (i/2)
in
g a one n
(* example with integers *) pow 1 ( * ) 2 16;; (* - : int = 65536 *)</lang>
Now matrix power is simply a special case of pow :
<lang ocaml>let matpow a n =
let p, q = dim a in if p <> q then failwith "bad dimensions" else pow (eye p) matmul a n;;
matpow (matrix 2 2 [ 1.0; 1.0; 1.0; 0.0 ]) 10;; (* - : float array array = [|[|89.; 55.|]; [|55.; 34.|]|] *) </lang>
Octave
Of course GNU Octave handles matrix and operations on matrix "naturally".
<lang octave>M = [ 3, 2; 2, 1 ]; M^0 M^1 M^2 M^(-1) M^0.5</lang>
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)
PicoLisp
Uses the 'matMul' function from Matrix multiplication#PicoLisp <lang PicoLisp>(de matIdent (N)
(let L (need N (1) 0)
(mapcar '(() (copy (rot L))) L) ) )
(de matExp (Mat N)
(let M (matIdent (length Mat))
(do N
(setq M (matMul M Mat)) )
M ) )
(matExp '((3 2) (2 1)) 3)</lang> Output:
-> ((55 34) (34 21))
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
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
TI-89 BASIC
Built-in exponentiation:
<lang ti89b>[3,2;4,1]^4</lang>
Output:
Ursala
For matrices of floating point numbers, the library function mmult can be used as shown. The user-defined id function takes a square matrix to the identity matrix of the same dimensions. The mex function takes a pair
representing a real matrix and a natural exponent to the exponentiation using the naive algorithm.
<lang Ursala>#import nat
- import lin
id = @h ^|CzyCK33/1.! 0.!* mex = ||id@l mmult:-0^|DlS/~& iota</lang> Alternatively, this version uses the fast binary algorithm. <lang Ursala>mex = ~&ar^?\id@al (~&lr?/mmult@llPrX ~&r)^/~&alrhPX mmult@falrtPXPRiiX</lang> This test program raises a 2 by 2 matrix to a selection of powers. <lang Ursala>#cast %eLLL
test = mex/*<<3.,2.>,<2.,1.>> <0,1,2,3,4,10></lang> output:
<
<
<1.000000e+00,0.000000e+00>,
<0.000000e+00,1.000000e+00>>,
<
<3.000000e+00,2.000000e+00>,
<2.000000e+00,1.000000e+00>>,
<
<1.300000e+01,8.000000e+00>,
<8.000000e+00,5.000000e+00>>,
<
<5.500000e+01,3.400000e+01>,
<3.400000e+01,2.100000e+01>>,
<
<2.330000e+02,1.440000e+02>,
<1.440000e+02,8.900000e+01>>,
<
<1.346269e+06,8.320400e+05>,
<8.320400e+05,5.142290e+05>>>