## 6 Implementation
### The Unitary Class  
The _Unitary_ is instantiated with a string and an array. The content of this array is the unitary matrix, which is of size \(2^{n}\times 2^{n}\) for an \(n\)-qubit gate. The complete Quantum Shannon Decomposition is computed only when \"decompose()\" is called, and the calculated angles for the resulting rotation gates are added to a list. This is done so that the _Unitary_ can be used multiple times in a program without the recalculation of the whole decomposition.  
However, before the decomposition is started, it is first checked whether the input matrix is unitary. If this is the case, all of the intermediate matrices will also be unitary [19], so this check is only necessary once. Furthermore, all of the Gray code matrices needed for the multi-controlled rotation gates are added to a lookup table so they do not need to be calculated anew at each decomposition step.  
To make certain that the decomposition is correct, each single intermediate decomposition is checked. For each step, only three matrices need to be multiplied, and this saves any calculations that might be performed on an incorrect matrix. If any step of the decomposition is not correct, an exception is thrown and the decomposition is stopped.  
The Eigen [23] library is used to perform singular value decomposition (SVD), eigenvalue decomposition and matrix multiplication. The recursion is centered on a main function, which takes as parameters a unitary matrix and the number of qubits. The latter is to keep track of the level of recursion.  
Computation of the CSD is done using the method from [19], which uses SVD. The demultiplexing function uses Schur matrix decomposition for (sub)matrices smaller than \(2^{6}\times 2^{6}\) and eigenvalue decomposition for bigger matrices. This is done because Schur matrix decomposition is faster for small matrices [23].  
The algorithm is recursive, and the demultiplexing step calls on the main function again for the decomposition of the smaller unitary matrices. If the matrices are of size \(2\times 2\), the rotation angles for the one-qubit rotation gates are calculated using _ZYZ_ decomposition as in Section 5.3.  
Because the _Unitary_ does not have access to the qubit numbers of the circuit, only the angles for the multi-controlled \(R_{y}\) and \(R_{z}\) are calculated at this point. This is done as in Section 3.6 by solving the following matrix equalities:\[M^{k}\begin{bmatrix}\theta_{1}\\ \vdots\\ \theta_{2^{k}}\end{bmatrix}=\begin{bmatrix}\alpha_{1}\\ \vdots\\ \alpha_{2^{k}}\end{bmatrix} \tag{19}\]  
where \(M^{k}\) is a square matrix where all the entries are either \"\(+1\)\" or \"\(-1\)\", which are calculated using Gray code using Equation (20).  
\[M^{k}_{ij}=(-1)^{b_{i}(-1)\cdot\gamma(j-1)} \tag{20}\]  
where the exponent is the bit-wise inner product of two binary vectors: \(b_{i}\) and \(\gamma_{j}\). \(b_{i}\) is the integer \(i\), and \(\gamma_{j}\) is the \(j\)th value of the Gray code.  
For the multi-controlled \(R_{y}\) gate, the values of \(\alpha_{i}\) are calculated by taking the arc sine of the diagonal entries of the \(S\)-matrix from the CSD.  
\[\begin{bmatrix}\alpha_{1}\\ \vdots\\ \alpha_{2^{k}}\end{bmatrix}=2\cdot arcsin(S_{i,i}) \tag{21}\]  
For the multi-controlled \(R_{z}\) gates, the values of \(\alpha_{i}\) is calculated by taking the natural logarithm of the \(D\)-matrix from the demultiplexing.  
\[\begin{bmatrix}\alpha_{1}\\ \vdots\\ \alpha_{2^{k}}\end{bmatrix}=-2\cdot\sqrt{-1}\cdot ln\Big{(}D_{(i,i)}\Big{)} \tag{22}\]  
All the angles for all rotation gates are added to a list, which is used to generate the correct gates when the _Unitary_ is added to a circuit.