Let us first consider the generating element for $C_2$ : $$M_1 = \left[\begin{array}{cc}0&1\\1&0\end{array}\right], \text{ and } P_1 = ({M_1})^2 = I_2 = \left[\begin{array}{cc}1&0\\0&1\end{array}\right]$$ If considering $I_2$ to be a number representing $+1$, it behaves just as the integer $-1$ does in the sense that $(-1)^N \approx {(M_1)}^N$ ( $\approx$ "corresponds" to ). Now strange things can happen when we do addition. For example $1+(-1)$ which should equal to 0 becomes the matrix: $\left[\begin{array}{cc}1&1\\1&1\end{array}\right]$. We can easily identify it as the sum of +1 and -1, but we also realize that we can not have uniqueness as any matrix $\left[\begin{array}{cc}N&N\\N&N\end{array}\right]$ will represent $N-N = 0$.
However we could just project any matrix we have ended up with onto $P_1$ and $M_1$ and plain remove the minimum from both. If we end up with the matrix:
$$\left[\begin{array}{cc}N&M\\M&N\end{array}\right] \text{ we remove } \min(M,N) \cdot \left[\begin{array}{cc}1&1\\1&1\end{array}\right]$$
This should give us $\left[\begin{array}{cc}N-M&0\\0&N-M\end{array}\right]$ or $\left[\begin{array}{cc}0&M-N\\M-N&0\end{array}\right]$ depending on which was largest. And then we have found a way to resolve the non-uniqueness.
Now we assume this works and then we plant such 2x2 blocks of elements in ${\mathbb Z_+}$ into $2\times 2$ blocks using the famous representation of complex numbers:$$a+bi \approx \left[\begin{array}{rr}a&b\\-b&a\end{array}\right]$$ This can easily be done / expressed with the help of Kronecker products or as block matrices:
$$\left[\begin{array}{cc}|a|P_1&|b|P_1\\|b|M_1&|a|P_1\end{array}\right]$$
And that any correction described above is applied on each of these 2x2 blocks.
For example we can calculate $(i+1)^8$ which we all know should be $(\sqrt{2})^8 = 2^4=16$
First out the $i+1$ matrix will be represented by:
$$\left[\begin{array}{cc|cc}0&0&1&0\\0&0&0&1\\\hline 0&1&0&0\\1&0&0&0\end{array}\right] + \left[\begin{array}{cc|cc}1&0&0&0\\0&1&0&0\\\hline 0&0&1&0\\0&0&0&1\end{array}\right] = \left[\begin{array}{cc|cc}1&0&1&0\\0&1&0&1\\\hline 0&1&1&0\\1&0&0&1\end{array}\right]$$
Carrying out matrix multiplications, we do get the matrix : $$\left[\begin{array}{cccc}72&56&64&64\\ 56&72&64&64\\64&64&72&56\\64&64&56&72\end{array}\right] = \left/\cases{64-64=0\\ 72-56=16}\right/ = \left[\begin{array}{cccc}16&0&0&0\\0&16&0&0\\0&0&16&0\\0&0&0&16\end{array}\right]$$
Which represents the complex number $+16 + 0i$ as it should. But now to the question...
Is this approach theoretically sound?
For a free headache I can toss in the following lines of octave
re = eye(4); % positive real one unit matrix
im = [0,0,1,0;0,0,0,1;0,1,0,0;1,0,0,0]; % positive imaginary one unit matrix
N_ = kron(eye(2),eye(2)(end:-1:1,:)); % real negative one unit matrix
for r_ = -8:8
for i_ = -8:8
abs(sum(conv2((N_^(r_<0)*re.*abs(r_) + N_^(i_<0)*im.*abs(i_))^3,-[1,-1],'valid')(1:2:end,1:2:end)(1:2).*[1,-i]) - (r_+1i*i_)^3)
endfor
endfor
They supposedly verify that the matrix implementation cube of each gaussian integer with $a,b \in [-8,8]$ equals the complex number implementation in the language.
This is just a tiny verification but of course no proof.