# C++: Matrix Exponentiation Operator

Posted in C++

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;
}
```

```// 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.

SOURCE

Content is available under GNU Free Documentation License 1.2.