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

Demonstrate how to implement matrix exponentiation as an operator.

#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.

Content is available under GNU Free Documentation License 1.2.