from numpy import *
from scipy.linalg import hadamard
import math
from random import randint
from tabulate import tabulate
The nontrival factors of $C$ can be found classically via Euclid's algorithm:
# Pick two large prime numbers
P1 = 3
P2 = 5
# Compute their product
C = P1*P2
print(f"C is {C}")
# Create a boolean
w = True
# Step 1
# Check that C is odd and is not a power of some smaller integer
if C%2 == 0:
w = False
r = int(math.log(C,3))
k = linspace(2,r,r-1)
for i in range(len(k)):
c1 = int(C**(1/k[i]))
c2 = c1**k[i]
if c2 == C:
w = False
print(f'{C} is the power {int(k[i])} of the smaller integer {c1}.')
while w is True:
# Step 2
# Pick any integer in the range 1 < a < C
a = randint(2,C-1)
# Step 3
# Find the greatest commmon demoninator
g = math.gcd(a,C)
if g > 1:
f1 = g
f2 = int(C/g)
print(f"factors are f1 = {f1} and f2 = {f2}.")
break
else:
# Step 4
# Find the smallest integer p > 1
p = 2
while True:
if (a**p)%C == 1:
break
p += 1
# Step 5
# If p is odd or has modular congruence -1, go back to step 2
if (p%2)!=0:
continue
else:
if (a**(p/2))%C != -1:
# Step 6
# Find the nontrivial factors of C
f1 = math.gcd(a**(int(p/2))+1,C)
f2 = math.gcd(a**(int(p/2))-1,C)
print(f"factors are f1 = {f1} and f2 = {f2}.")
break
else:
continue
print(f"a is {a}")
The only step in which a quantum computer is needed is step 4, which signifies that $a^{p}-1$ is an integer multiple of $C$. This step is called "period finding", because the function $f(x)=a^{x} (mod C)$ is periodic with period $p$.
L = 3 # qubits in x-register
M = 4 # qubits in f-register
N = L+M # total number of qubits
s = 2**N # total number of states
I = eye(2) # 2d identity input in had()
H = hadamard(2) # 2d hadamard input in had()
factor = 1/sqrt(2) # hadamard factor
def had(i,N,H,I):
"""Compute the hadamard acting on the i-th qubit"""
if i == N:
product = H
else:
product = I
for j in range(N-1,0,-1):
if j!=i:
product = kron(I,product)
else:
product = kron(H,product)
return array(product).astype(int)
def allstates(s):
"""Represent all possible state in binary"""
a = linspace(0,s-1,s,dtype=int) # numerical possible states
b = list() # holds binary possible states
for i in range(len(a)):
qstr = '{:07b}'.format(a[i])
b.append(qstr)
return b
def kcol(kbin, qcontrol, A):
"""Given a column and controlling-qubit, calculate
the corresponding row that has the value 1 for a
controlled phase-shift gate"""
q = qcontrol - 1
l0 = kbin[q]
k = int(kbin,2) # column number
j = int
f = int
fp = int
fpbin = str
jbin = str
Amatrix = int
if qcontrol == 3:
Amatrix = A[0]
elif qcontrol == 2:
Amatrix = A[1]
elif qcontrol == 1:
Amatrix = A[2]
if int(l0) == 0:
j = k
elif int(l0) == 1:
f = int(kbin[3:7],2) # m3m2m1m0 as decimal
if f >= C:
j = k
elif f < C:
fp = Amatrix*f%C
fpbin = '{:04b}'.format(fp)
jbin = kbin[0:3] + fpbin
j = int(jbin,2)
return j
def fx(b,qcontrol,A):
"""Create the controled-phase shift gate for a controlling-qubit"""
m = zeros((s,s))
for i in range(len(b)):
j = kcol(b[i],qcontrol,A)
m[j,i] = 1
return m
def unitary(m,s):
"""Determine whether a matrix is unitary (true/false)"""
z = m.dot(m.conj().T)
return array_equal(z,eye(s))
# For qubit 1, 2, or 3 controlling qubit 1, 2, or 3
def CNOT(Cgate,row,col,b,qbitA,qbitB):
"""Construct the full 128x128 matrix given a phase shift matrix"""
row = array(list(b[row])).astype(int)
col = array(list(b[col])).astype(int)
qA = qbitA-1
qB = qbitB-1
p,pp = row[qA],col[qA]
q,qp = row[qB],col[qB]
boo = True
# kronecker delta (T/F) corresponds to (1,0)
for i in range(len(row)):
if (i!=qA) or (i!=qB):
if row[i] != col[i]:
boo = False
if boo == True:
rs = str(q) + str(p)
ri = int(rs,2)
cs = str(qp) + str(pp)
ci = int(cs,2)
return Cgate[ri][ci].real
else:
return 0
The quantum computation requires Hadamard gates, controlled phase-shift gates, and permutation matricies. The result of this algorithm is the period.
# Create initial state of the system
qb = zeros(s,dtype=int)
qb[1] = 1 # Turn on state '0000001'
qb = array([qb]).T
# Create binary of all possible states
b = allstates(s)
# Compute the Hadamard gates on qubits 1, 2, 3
H1, H2, H3 = factor*had(1,N,H,I), factor*had(2,N,H,I), factor*had(3,N,H,I)
# Values of 'a' to iterate through
avals = [2,4,7,8,11,13,14]
# Initialize phase-shift matricies
Cgate1 = array([[1,0,0,0],
[0,1,0,0],
[0,0,1,0],
[0,0,0,exp(1j*pi/2)]])
Cgate2 = array([[1,0,0,0],
[0,1,0,0],
[0,0,1,0],
[0,0,0,exp(1j*pi/4)]])
# Compute controlled phase-shift gates
Cnot12 = zeros((s,s))
Cnot13 = zeros((s,s))
Cnot23 = zeros((s,s))
for i in range(0,s):
for j in range(0,s):
Cnot12[i][j] = CNOT(Cgate1,i,j,b,1,2)
Cnot13[i][j] = CNOT(Cgate2,i,j,b,1,3)
Cnot23[i][j] = CNOT(Cgate1,i,j,b,2,3)
def Shor(a,qstates):
"""Shor's Algorithm for 7-qubits"""
A = [(a)%C, (a**2)%C, (a**4)%C]
M1, M2, M3 = fx(b,1,A), fx(b,2,A), fx(b,3,A)
qstates = H3@Cnot23@H2@Cnot13@Cnot12@H1@M3@M2@M1@H3@H2@H1@qstates
poss = list()
xtilde = list()
for i in range(len(qstates)):
if qstates[i] > 0:
poss.append(int(b[i][3:7],2))
xstr = b[i][2] + b[i][1] + b[i][0]
xtilde.append(int(xstr,2))
poss = unique(poss)
xtilde = unique(xtilde)
xx = xtilde/(2**L)
return (poss,xtilde,xx)
# Create the table
table = list()
for a in avals:
poss, xtilde, xx = Shor(a,qb)
table.append([a,poss,xtilde,xx])
print(tabulate(table, headers=["a","f(X) values","x-tilde", "results"]))
Because the measured $\frac{\tilde{x}}{2^{L}}$ values were multiples of $\frac{1}{4}$, they give a guess $p=4$ for the period (the correct answer).
table = list()
avals = linspace(2,C-1,C-2,dtype=int)
for a in avals:
A = [(a)%C, (a**2)%C, (a**4)%C]
M1, M2, M3 = fx(b,1,A), fx(b,2,A), fx(b,3,A)
u1, u2, u3 = unitary(M1,s), unitary(M2,s), unitary(M3,s)
table.append([a,u1,u2,u3])
print("Permutation matricies marked as True are definied as unitary.")
print("This is why the previous values of 'a' were chosen as [2,4,7,8,11,13,14].")
print(tabulate(table,headers=["a","M1","M2","M3"]))