Matrix exponentiation operator
From Rosetta Code
Programming Task
This is a programming task. It lays out a problem which Rosetta Code users are encouraged to solve, using languages they know.
The following programs demonstrates how to implement a "complex number" matrix exponentiation (**) as an operator.
Contents |
[edit] 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:
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;
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.
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;
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.
[edit] 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));
[edit] C++
This is an implementation in C++.
#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;
}
This is the task part.
// 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;
}
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.
[edit] D
This is a implementation by 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 ;
}
This is the task part.
// 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) ;
}
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 preseted. For instance, arithemic + * , 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-comutative operator should be chose, or implement a corresponding opXXX_r overloading that may throw static assert/error.
Details see Operator Loading in D
[edit] 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).
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]]
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.
Categories: Programming Tasks | Matrices | Ada | ALGOL 68 | C++ | D | Haskell
