Matrix-exponentiation operator: Difference between revisions
(lang tags, grammar, omit from Java) |
(add Tcl) |
||
Line 506: | Line 506: | ||
(Of course this is not an implementation, but it can be used as reference for the results) |
(Of course this is not an implementation, but it can be used as reference for the results) |
||
=={{header|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 exponenets 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> |
|||
<pre>% 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 |
|||
</pre> |
Revision as of 11:31, 29 April 2009
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**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.
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)
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 exponenets 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