# Showing the Toffoli gate can be constructed from one-qubit gates using computer matrix multiplication

I need help coming up with the code that I can use in a computer algebra system like Sage, Python, or Mathematica that can evaluate the circuit.

I am not a programmer and used ChatGPT to try to help me create the code. However, I believe these are incorrect because they do not evaluate.

Here is what ChatGPT came up with in Mathematica:

TensorProduct[(1/Sqrt[2]){{1, 1}, {1, -1}}, IdentityMatrix[2], IdentityMatrix[2]] .

TensorProduct[{{1, 0}, {0, e^(-i*pi/4)}}, IdentityMatrix[2], IdentityMatrix[2]] .

TensorProduct[{{1, 0}, {0, e^(i*pi/4)}}, IdentityMatrix[2], IdentityMatrix[2]] .

TensorProduct[{{1, 0}, {0, e^(-i*pi/4)}}, IdentityMatrix[2], IdentityMatrix[2]] .

TensorProduct[{{1, 0}, {0, e^(i*pi/4)}}, {{1, 0}, {0, e^(-i*pi/4)}}, IdentityMatrix[2]] .

TensorProduct[(1/Sqrt[2]){{1, 1}, {1, -1}}, IdentityMatrix[2], IdentityMatrix[2]] .

TensorProduct[IdentityMatrix[2], {{1, 0}, {0, e^(-i*pi/4)}}, IdentityMatrix[2]] .

TensorProduct[IdentityMatrix[2], {{1, 0}, {0, -1}}, {{1, 0}, {0, e^(i*pi/4)}}]


Here is what ChatGPT came up with in Sage/Python:

from sage.matrix.constructor import matrix
from sage.symbolic.constants import pi
from sage.functions.trig import exp

# Define matrices
A = (1/sqrt(2)) * matrix([[1, 1], [1, -1]])
B = matrix.identity(2)

# Compute the result of the sequential matrix products
result = A.tensor_product(B).tensor_product(B) * \
matrix([[1, 0], [0, exp(-I*pi/4)]]) \
* matrix([[1, 0], [0, exp(I*pi/4)]]) \
* matrix([[1, 0], [0, exp(-I*pi/4)]]) \
* matrix([[1, 0], [0, exp(I*pi/4)]]) \
* matrix([[1, 0], [0, exp(-I*pi/4)]]) \
* A.tensor_product(B).tensor_product(B) \
* B.tensor_product(matrix([[1, 0], [0, exp(-I*pi/4)]])) \
* B.tensor_product(matrix([[1, 0], [0, -1]])) \
* matrix([[1, 0], [0, exp(I*pi/4)]])

# Display the result
print(result)


The main problem is that ChatGPT completely ignores CNOT gates in both cases (and $$S$$ gate as well). What is more, you need not only CNOT for two adjacent qubits but also CNOT for non-adjacent qubits. How to construct a matrix for controlled gates with non-adjacent qubits is described in my other answer here. In the answer, I also provided matrix for CNOT acting on first and last qubit in your circuit.

But back to the code you are looking for. Below is the code for MatLab (or Octave). You can easily rewrite it to any other language. You can inspire by answer by ChatGPT in your question but replace matrices with correct ones. Note that function kron stands for tensor (Kronecker) product.

Please note that my code is for upside-down Toffoli, i.e. a target qubit is the lowest one (see a schematic here). If you run it, in variable TOF, there is this matrix: $$\begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\ \end{pmatrix},$$ which is clearly Toffoli gate.

So, here is the code:

%define gates
ID = eye(2); %identity operator
H = (1/sqrt(2))*[1 1; 1 -1]; %Hadamard
S = [1 0; 0 exp(i*pi/2)]; %S gate
T = [1 0; 0 exp(i*pi/4)]; %T gate
Tdag = [1 0; 0 exp(-i*pi/4)]; %inverse of T gate
X = [0 1; 1 0]; %Pauli X - negation
CNOT_12 = [eye(2) zeros(2); zeros(2) X]; %CNOT on adjacent qubits
CNOT_13 = [eye(4) zeros(4); zeros(4) kron(eye(2), X)]; %CNOT between first and last qubits

%construct Toffoli gate
TOF = kron(kron(ID,ID),H);
TOF = kron(ID, CNOT_12) * TOF;
TOF = kron(kron(ID,ID),Tdag) * TOF;
TOF = CNOT_13 * TOF;

TOF = kron(kron(ID,ID),T) * TOF;
TOF = kron(ID, CNOT_12) * TOF;
TOF = kron(kron(ID,ID),Tdag) * TOF;
TOF = CNOT_13 * TOF;

TOF = kron(kron(ID,Tdag),T)*TOF;
TOF = kron(CNOT_12,H) * TOF;
TOF = kron(kron(ID,Tdag),ID)*TOF;
TOF = kron(CNOT_12,ID) * TOF;
TOF = kron(kron(T,S),ID)*TOF;

TOF %print matrix stored in variable TOF