Implementing Shor’s Algorithm

We seek the prime factors p,q of a semi-prime number n= p q using Shor’s algorithm. The algorithm requires that we find the order of the function  f(r) = (x^r) Mod n , where x (the seed) is an arbitrary integer co-prime with n. In this demonstration we  find the prime factors of n=p q. We need to construct the following circuit.

We allow the upper (red) register to contain n q qubits in order to find a factor of n ∈ {((0,1})^nq). We also allow the second (blue) register to be constructed from the same number of qubits.  So this quantum circuit requires gate operators, or matrices, that contain (2^nq)× (2^nq)elements. Let’s analyze the above diagram for each time segment (Subscript[t, i]) Below we factor the semi-prime number 21, using registers with nq=7 qubits.

(* we define the following *)
ClearAll["Global`*"];
n=21;
nq=7;
ket[n_]:=Table[If[i==n+1,1,0],{i,1,2^nq}];
(*(Subscript[ ket[n] is a matrix representation for |n⟩, 4]) *)
Had=HadamardMatrix[2];
(* Had is the single qubit Hadamard gate *)
unit=IdentityMatrix[2^nq];
(* unit is a (2^nq) × (2^nq) identity matrix *)

Section  (Subscript[t, 1])

We input, in matrix form, the ket  |0..0⟩⊗|0...0⟩.

register1=ket[0];
register2=ket[0];
ψ[t1]=Flatten[KroneckerProduct[register1,register2]];

Section  (Subscript[t, 2])

The Hadamard gates operate on each qubit of the red register, so that ψ(t)
({{
({{(H^(⊗nq))}})⊗1(Subscript[|0⟩, nq])⊗(Subscript[|0⟩, nq])=(1/(Sqrt[(2^nq)]))(∑_(x=0)^((2^nq)-1))(Subscript[|x⟩, nq])⊗(Subscript[|0⟩, nq])}})

Section  (Subscript[t, 3])

We now construct the gate (unitary matrix)  that implements gate (Subscript[U, f]). This gate is represented by a (2^nq)× (2^nq) matrix which is costly, in time and memory resources, to construct. Instead, we recognize that the second (blue) register contains only state |00...⟩ so

({{({{ (Subscript[U, f]) (H^(⊗ nq))}})(Subscript[|0 ⟩, nq])⊗ (Subscript[|0⟩, nq]) = (1/(Sqrt[(2^nq)])) (∑_(x=0)^((2^nq)-1))(Subscript[|x⟩, nq])⊗(Subscript[| f(x) ⟩, nq])}})

and we can use this expression to evaluate ψ( (Subscript[t, 3])) without needing the full matrix representation for  (Subscript[U, f]).

(* define seed x, and fuction f(x), make sure that x is co-prime with n, i.e. GCD[x,n]==1 *)
x=8;GCD[x,n]==1
f[r_]:=Mod[x^r,n];

ψ[t3]=1/Sqrt[2^nq]Sum[Flatten[KroneckerProduct[ket[i],ket[f[i]]]],{i,0,2^nq-1}];

Section  (Subscript[t, 4])

We now construct the discrete Fourier transform gate. Using the method described in Notebook 4.2, the Fourier
gate for the upper (red) register is given by a (2^4) ×  (2^4) matrix.

Ufourier=Table[ket[i].Fourier[ket[j]],{i,0,2^nq-1},{j,0,2^nq-1}];

and so, according to the circuit diagram,

ψ[t4]=KroneckerProduct[Ufourier,unit].ψ[t3];

Measurement

We perform a measurement of the four qubits in the upper (red) register. Using the Born rule, the probability to obtain configuration {(Subscript[x, 3]) (Subscript[x, 2])(Subscript[x, 1]) (Subscript[x, 0])} is given by the expression


({{(∑_(a=0)^((2^nq)-1))|⟨a|⊗ ⟨({{(Subscript[x, 0])}})(Subscript[x, 1])  (Subscript[x, 2]) ...(Subscript[x, nq-1])|ψ((Subscript[t, 4]))⟩((|)^2)}})

p[i_]:=Sum[Abs[Flatten[KroneckerProduct[ket[i],ket[a]]].ψ[t4]]^2,{a,0,2^nq-1}]

probabilities=Table[p[i],{i,0,2^nq-1}];
binaryoutput=Table[IntegerDigits[i,2,nq],{i,0,2^nq-1}];
measurement:=RandomChoice[probabilities->binaryoutput];

Perform a measurement using the measurement command

measurement

If the measurement does not yield the trivial value {0,0,0,. . .}, save the result and convert the register reading to base-10 form

measurement
k=FromDigits[%,2]

Construct the ratio k/(2^nq), which is likely to be an integer multiple of 1/r

ratio=k/2^nq
Convergents[ratio]

Guess for the order r

r=2

If r is an even number, and (x^(r/2)) Mod n ≠ ±1, evaluate

If[Abs[Mod[x^(r/2),n]]==1,Print[" \!\(TraditionalForm\`\*SuperscriptBox[\(x\), \(r/2\)] Mod\\ n\\  = \\ \(±1\)\)  Try Again"]]

p1=GCD[x^(r/2)-1,n];p2=GCD[x^(r/2)+1,n];Print[" p = ",p1," q= ",p2]

Check if you guessed correctly (if either p,q ≠ 1)

n==p1 p2

If you did not guess correctly, try with another measurement or re-run the algorithm with a different seed.

A (2) x 12 qubit Shor Simulator

Below we link to a 24- qubit Shor box simulator. It is able to factor the first 40, odd, semi- prime numbers.

Hyperlink["http://www.physics.unlv.edu/%7Ebernard/MATH_book/Chap4/Shorbox.nb"]