Matrix-exponentiation operator: Difference between revisions
No edit summary |
m (→{{header|C++}}: lang tag) |
||
Line 258: | Line 258: | ||
=={{header|C++}}== |
=={{header|C++}}== |
||
This is an implementation in C++. |
This is an implementation in C++. |
||
⚫ | |||
<pre> |
|||
⚫ | |||
#include <cmath> |
#include <cmath> |
||
#include <iostream> |
#include <iostream> |
||
Line 307: | Line 306: | ||
} |
} |
||
return d; |
return d; |
||
} |
}</lang> |
||
</pre> |
|||
This is the task part. |
This is the task part. |
||
<lang cpp> |
|||
<pre style="background-color:#ffe"> |
|||
// C++ does not have a ** operator, instead, ^ (bitwise Xor) is used. |
// C++ does not have a ** operator, instead, ^ (bitwise Xor) is used. |
||
Mx operator^(int n) { |
Mx operator^(int n) { |
||
Line 322: | Line 320: | ||
return d; |
return d; |
||
} |
} |
||
</pre> |
|||
<pre> |
|||
}; |
}; |
||
Line 341: | Line 337: | ||
return 0; |
return 0; |
||
}</lang> |
|||
} |
|||
</pre> |
|||
Output: |
Output: |
||
<pre> |
<pre> |
Revision as of 00:20, 4 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 integer and real only.
The following programs demonstrates how to implement a "complex number" 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.
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
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.