We will need a few functions to do necessary mathematical operations on qubits (or more specifically on the state vector or probability vector of a quantum system). Note that I follow what I believe is the most common numbering of inputs with qubit 0 being the first input from above in the circuit and the first bit from the left in the state vector index, i.e the opposite of how we usually number bits in binary, from right to left (e.g if probability[1] in a three qubit system is 100% then the bits in the binary representation are 001 with a 1 in bit number 0, which is qubit number 2). You can see this reversal when we are doing actual bit operations (bit shifts) in the following functions.
We need a matrix multiplication routine (dot product):
def mult(a,b):
c=[[0 for x in range(len(b[0]))] for x in range(len(a))]
for i in range(len(a)):
for j in range(len(b[i])):
for k in range(len(b)):
c[i][j]+=a[i][k]*b[k][j]
return c
We need a matrix tensor product routine (cross product or more specifically the Kronecker product):
def kron(a,b):
c=[[0 for x in range(len(a[0])*len(b[0]))]\
for x in range(len(a)*len(b))]
for i in range(len(a)):
for k in range(len(b)):
for j in range(len(a[i])):
for l in range(len(b[k])):
c[i*len(b)+k][j*len(b[k])+l]=a[i][j]*b[k][l]
return c
We need a routine for measuring qubits (or more specifically making a measurement on the probability vector, also note that this function for simplicity measures ALL qubits in the system, but that is ok as the state is not stored back, instead we can set wanted measurements back with a qbitset function, see later):
def M(a):
c=[[0] for x in range(len(a))]
r=random()
j=len(a)-1
s=0.0
for i in range(j):
s+=abs(a[i][0])**2
if s>=r:
j=i
break
c[j][0]=1
return c
We need a routine for storing the measured state of a qubit (0 or 1) back into the probability vector (by zeroing out states we now know have probability zero, after which a normalization of the probabilities are needed, see later):
def qbitset(p,v,a):
for i in range(len(a)):
if (i>>round(log(len(a),2)-p-1))&1!=v:
a[i]=[0]
return norm(a)
We need a routine for normalizing a probability vector (the sum of the probabilities must be 1):
def norm(a):
s=0.0
for i in range(len(a)):
s+=abs(a[i][0])**2
t=sqrt(s)
for i in range(len(a)):
a[i][0]/=t
return a
We need a function for returning the value of a qubit (0 or 1) after a measurement (by finding the state that now has probability 1, otherwise if the qubit is still in superposition just return -1) :
def qbit(p,a):
for i in range(len(a)):
if a[i][0]==1:
if (i>>round(log(len(a),2)-p-1))&1==1:
return 1
else:
return 0
return -1
