from sympy import *
from sympy.physics.quantum.state import Ket, Bra
from sympy.physics.quantum import TensorProduct
import numpy as np
init_printing(use_unicode=True)

# Q bit gates
H=Matrix([[1/sqrt(2), 1/sqrt(2)],[1/sqrt(2), -1/sqrt(2)]])
S=Matrix([[1, 0],[0, I]])
X=Matrix([[0, 1],[1, 0]])
Y=Matrix([[0, -I],[I, 0]])
Z=Matrix([[1, 0],[0, -1]])
Identity=Matrix([[1, 0],[0, 1]])
H, S, X, Y, Z, Identity

#Q  bit gates
CCX=Matrix([[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]])
CNOT=Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
SWAP=Matrix([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]])
CIIX=(TensorProduct(Identity, Identity, SWAP)  #4 qbit CNOT, Cqbit=0 and acted qbit=3
      *TensorProduct(Identity, SWAP, Identity)
      *TensorProduct(CNOT, Identity, Identity)
      *TensorProduct(Identity, SWAP, Identity)
      *TensorProduct(Identity, Identity, SWAP))
CCX, SWAP,CIIX, CNOT

zero=Matrix([[1],[0]])
one=Matrix([[0],[1]])
plus=Matrix([[1/sqrt(2)], [1/sqrt(2)]])
minus=Matrix([[1/sqrt(2)], [-1/sqrt(2)]])
zero, one, plus, minus

circuit=(TensorProduct(Identity, Identity, CNOT)
        * TensorProduct(Identity, CCX)
        * CIIX
        * TensorProduct(H, Identity, H, H)
        * TensorProduct(zero, one, zero, one))
circuit

lhs=circuit
lhs

#option A
rhs=TensorProduct(minus, one, plus, minus)
print (lhs)
print (rhs)
if (lhs==rhs):
    print ("It is true")
else:
    print ("It is false")

#option B
rhs=(1/sqrt(2))*(TensorProduct(minus, one, zero, minus)+TensorProduct(minus, one, one, minus))
print (lhs)
print (rhs)
if (lhs==rhs):
    print ("It is true")
else:
    print ("It is false")

#option C
rhs=(-1)*(TensorProduct(plus, one, plus, minus))
print (lhs)
print (rhs)
if (lhs==rhs):
    print ("It is true")
else:
    print ("It is false")

#option D
rhs=(TensorProduct(plus, one, plus, minus))
print (lhs)
print (rhs)
if (lhs==rhs):
    print ("It is true")
else:
    print ("It is false")

#option A and B
opt_b=(1/sqrt(2))*(TensorProduct(minus, one, zero, minus)+TensorProduct(minus, one, one, minus))
opt_b

opt_a=TensorProduct(minus, one, plus, minus)
opt_a