Implementing Shor’s algorithm in Perceval

This notebook presents a simulation in Perceval of a 4-qubit 12-modes optical circuit performing Shor’s algorithm, based on Alberto Politi, Jonathan C.F. Matthews, and Jeremy L. O’brien. “Shor’s quantum factoring algorithm on a photonic chip.” Science 325.5945 (2009): 1221-1221.

Shor’s algorithm

The purpose of Shor’s algorithm is to find nontrivial factors of a given number \(N\) in polynomial time, yielding an near-exponential speedup compared to state of the art classical algortihms.

The main routine of Shor’s algorithm consists in finding the order \(r\) of a number \(a \in \mathbb{Z}_N\), i.e. the smallest integer \(r\) such that \(a^r = 1 \pmod N\).

If the order of a randomly chosen \(a\) which is coprime with \(N\) is even, then \((a^{r/2} - 1)(a^{r/2} + 1) = k N\). If none of these factors are multiples of \(N\), then \(\gcd(N, a^{r/2} - 1)\) and \(\gcd(N, a^{r/2} + 1)\) are nontrivial factors of \(N\).

Preliminaries

import sympy as sp
import perceval as pcvl

Path encoding functions

The following functions allow for conversion between the qubit and Fock state representations.

def toFockState(qubitState):
    # path encoding
    pe = {0:[1,0],  1:[0,1]}
    return [0] + pe[qubitState[0]] + pe[qubitState[2]] + [0, 0] + pe[qubitState[1]] + pe[qubitState[3]] + [0]

def toQubitState(fockState):
    # qubit modes
    x1 = [1, 2]
    f1 = [3, 4]
    x2 = [7, 8]
    f2 = [9, 10]
    # auxiliary modes
    am1 = [0, 5]
    am2 = [6, 11]

    # auxiliary modes
    for i in am1 + am2:
        if fockState[i]!=0:
            return None
    L=[]
    # qubit modes
    for q in [x1, x2, f1, f2]:
        if fockState[q[0]]+fockState[q[1]] != 1:
            return None
        else:
            L.append(fockState[q[1]])
    return L

def strState(state):
    return str(pcvl.BasicState(state))

The circuit

Quantum circuit

The quantum circuit has been optimized after choosing parameters \(N = 15\) and \(a = 2\), and aims to calculate \(r=4\). It features 5 qubits labelled \(x_0, x_1, x_2\) and \(f_1, f_2\). Qubits \(x_i\) and \(f_1\) are initially in state \(|0\rangle\), and \(f_2\) in state \(|1\rangle\). In the non-optimised Shor algorithm, qubits \(x_i\) encode a binary number representing a pre-image of the Modular Exponentiation Function (MEF) \(x \mapsto a^x \pmod N\), while qubits \(f_i\) hold the image obtained after applying the MEF to qubits \(x_i\). Applying the MEF when qubits \(x_i\) hold a superposition of different pre-images (obtained with H gates on qubits \(x_i\)) allows to efficiently compute the order \(r\) of parameter \(a\) modulo \(N\).

The circuit consists of \(\mathsf{H}\) gates being first applied to each \(x_i\) qubit, followed by \(\mathsf{CNOT}\) gates applied on \(x_1, f_1\) and \(x_2, f_2\) pairs, where \(x_i\) are control qubits; finally the inverse \(\mathsf{QFT}\) algorithm is applied on qubits \(x_i\).

\(\mathsf{CNOT}\) gates on \(x_i, f_i\) pairs (\(x_i\) being the control) are implemented using \(\mathsf{H}\) and \(\mathsf{CZ}\) gates: the \(\mathsf{CZ}\) gate is sandwiched between two applications of \(\mathsf{H}\) on \(f_i\).

The input state of the circuit after optimisation is \(|0\rangle_{x_0}|0\rangle_{x_1}|0\rangle_{x_2}|0\rangle_{f_1}|1\rangle_{f_2}\).

The expected output state is then \(\frac{1}{2} |0\rangle_{x_0} \otimes \left ( |0\rangle_{x_1}|0\rangle_{f_1} + |1\rangle_{x_1}|1\rangle_{f_1} \right ) \otimes \left ( |0\rangle_{x_2}|1\rangle_{f_2} + |1\rangle_{x_2}|0\rangle_{f_2} \right )\).

Photonic circuit

The optical circuit from the result by Politi et al features twelve modes (ordered from 0 to 11 from top to bottom).

During the execution, qubit \(x_0\) remains unentangled from the other qubits. It can therefore be removed from the optical implementation.

The qubits \(x_1, x_2, f_1, f_2\) are path encoded as modes \((1, 2)\), \((3, 4)\), \((7, 8)\), \((9, 10)\) respectively. The four remaining modes are used as auxiliary modes to implement the \(\mathsf{CZ}\) gates.

With path encoding each \(\mathsf{H}\) gate in the quantum circuit is implemented with a beam splitter with reflectivity \(R=1/2\) between the two pathes corresponding to the qubit. In our implementation in Perceval, phase shifters are added to properly tune the phase between each path.

\(\mathsf{CZ}\) gates are implemented with three beam splitters with reflectivity \(R=2/3\) acting on six modes: one inner BS creates interference between the two qubits, and two outer BS balance detection probability using auxiliary modes. This optical implementation succesfully yields the output state produced by a \(\mathsf{CZ}\) gate with probability 1/9; otherwise it creates a dummy state, which can be removed by post-selection.

In the case \(r=4\) the QFT can be performed classically and doesn’t need to be implemented in the photonic circuit.

In perceval

Implementing the circuit

circ = pcvl.Circuit(12)

# qubit modes
# for qubit states 0, 1
x1 = [1, 2]
f1 = [3, 4]
x2 = [7, 8]
f2 = [9, 10]
# auxiliary modes
am1 = [0, 5]
am2 = [6, 11]


# H gates
for q in [x1, f1, x2, f2]:
    circ.add(q, pcvl.BS.H())

# CZ gates
for x, f, am in [(x1, f1, am1), (x2, f2, am2)]:
    circ.add((am[0], x[0]), pcvl.BS(pcvl.BS.r_to_theta(1/3))) # R = 1/3
    circ.add((x[1],  f[0]), pcvl.BS(pcvl.BS.r_to_theta(1/3)))
    circ.add((f[1], am[1]), pcvl.BS(pcvl.BS.r_to_theta(1/3)))

# H gates
for q in [f1, f2]:
    circ.add(q, pcvl.BS.H())

print("Optical circuit for Shor's algorithm")
pcvl.pdisplay(circ)

Optical circuit for Shor's algorithm

print("The associated matrix")
pcvl.pdisplay(circ.compute_unitary(use_symbolic=False))

The associated matrix

$\left[\begin{array}{cccccccccccc}\frac{\sqrt{3}}{3} & \frac{\sqrt{3} i}{3} & \frac{\sqrt{3} i}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\frac{\sqrt{6} i}{3} & \frac{\sqrt{6}}{6} & \frac{\sqrt{6}}{6} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & \frac{\sqrt{6}}{6} & - \frac{\sqrt{6}}{6} & \frac{\sqrt{3} i}{3} & \frac{\sqrt{3} i}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & \frac{\sqrt{6} i}{6} & - \frac{\sqrt{6} i}{6} & \frac{\sqrt{3}}{3} & 0 & \frac{\sqrt{3} i}{3} & 0 & 0 & 0 & 0 & 0 & 0\\0 & \frac{\sqrt{6} i}{6} & - \frac{\sqrt{6} i}{6} & 0 & \frac{\sqrt{3}}{3} & - \frac{\sqrt{3} i}{3} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & \frac{\sqrt{3} i}{3} & - \frac{\sqrt{3} i}{3} & \frac{\sqrt{3}}{3} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & \frac{\sqrt{3}}{3} & \frac{\sqrt{3} i}{3} & \frac{\sqrt{3} i}{3} & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & \frac{\sqrt{6} i}{3} & \frac{\sqrt{6}}{6} & \frac{\sqrt{6}}{6} & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{\sqrt{6}}{6} & - \frac{\sqrt{6}}{6} & \frac{\sqrt{3} i}{3} & \frac{\sqrt{3} i}{3} & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{\sqrt{6} i}{6} & - \frac{\sqrt{6} i}{6} & \frac{\sqrt{3}}{3} & 0 & \frac{\sqrt{3} i}{3}\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{\sqrt{6} i}{6} & - \frac{\sqrt{6} i}{6} & 0 & \frac{\sqrt{3}}{3} & - \frac{\sqrt{3} i}{3}\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{\sqrt{3} i}{3} & - \frac{\sqrt{3} i}{3} & \frac{\sqrt{3}}{3}\end{array}\right]$

Input state

In the preliminaries we provided functions to map qubit states to the corresponding Fock states and vice-versa.

A computational basis qubit state on qubits \(x_1, f_1, x_2, f_2\) is represented with a list of 4 boolean values.

A Fock state is represented with a list of twelve integer values. As described above, for Fock states, the modes are enumerated as follows: * mode pairs \((1,2), (3,4), (7,8), (9,10)\) respectively correspond to states \(0,1\) for qubits \(x_1, f_1, x_2, f_2\) * modes \(0,5,6,11\) are auxiliary modes used for CZ gates

The input state of the circuit is \(|0\rangle_{x_1}|0\rangle_{x_2}|0\rangle_{f_1}|1\rangle_{f_2}\). With path encoding this corresponds to sending 4 photons in total in the optical circuit, in the input modes corresponding to the inital state of each qubit.

qubit_istate = [0,0,0,1]
istate = toFockState(qubit_istate)

print("Input qubit state:", strState(qubit_istate))
print("Corresponding input Fock state:", strState(istate))

Input qubit state: |0,0,0,1>
Corresponding input Fock state: |0,1,0,1,0,0,0,1,0,0,1,0>

Simulation

backend = pcvl.BackendFactory().get_backend("Naive")
backend.set_circuit(circ)
backend.set_input_state(pcvl.BasicState(istate))

Computing the output state

Using perceval we compute the output state obtained with the optical circuit.

Amplitudes

We first decompose the output state in the computational basis and plot the corresponding amplitudes.

output_qubit_states = [
    [x1,x2,f1,f2]
    for x1 in [0,1] for x2 in [0,1] for f1 in [0,1] for f2 in [0,1]
]

print("Output state amplitudes: (post-selected on qubit states, not renormalized)")
print("|x1,x2,f1,f2>")
for oqstate in output_qubit_states:
    ostate = toFockState(oqstate)
    a = backend.prob_amplitude(pcvl.BasicState(ostate))
    print(strState(oqstate), a)

Output state amplitudes: (post-selected on qubit states, not renormalized)
|x1,x2,f1,f2>
|0,0,0,0> 0j
|0,0,0,1> (0.055555555555555566+0j)
|0,0,1,0> 0j
|0,0,1,1> 0j
|0,1,0,0> (-0.05555555555555557-1.734723475976807e-18j)
|0,1,0,1> (-1.3877787807814457e-17-1.734723475976807e-18j)
|0,1,1,0> 0j
|0,1,1,1> 0j
|1,0,0,0> 0j
|1,0,0,1> (-1.3877787807814457e-17+0j)
|1,0,1,0> 0j
|1,0,1,1> (-0.05555555555555557-1.734723475976807e-18j)
|1,1,0,0> (1.3877787807814457e-17+0j)
|1,1,0,1> 0j
|1,1,1,0> (0.05555555555555559+0j)
|1,1,1,1> (1.3877787807814457e-17+5.204170427930421e-18j)

The amplitudes obtained with perceval correspond to the expected output state

\[\frac{1}{2} \left ( |0\rangle_{x_1}|0\rangle_{f_1} + |1\rangle_{x_1}|1\rangle_{f_1} \right ) \otimes \left ( |0\rangle_{x_2}|1\rangle_{f_2} + |1\rangle_{x_2}|0\rangle_{f_2} \right )\]

up to numerical precision, without renormalization.

Distribution

We now compute the probabilities to obtain each outcome corresponding to measuring the expected output state in the computational basis.

input_states = {
    pcvl.BasicState(pcvl.BasicState(istate)): strState(qubit_istate)
}

expected_output_states = {
    pcvl.BasicState(toFockState([x1,x2,x1,1-x2])): strState([x1,x2,x1,1-x2])
    for x1 in [0,1] for x2 in [0,1]
}

p = pcvl.Processor("Naive", circ)

ca = pcvl.algorithm.Analyzer(p, input_states, expected_output_states)
ca.compute()

print("Output state distribution: (post-selected on expected qubit states, not renormalized)")
print("|x1,x2,f1,f2>")
pcvl.pdisplay(ca)

Output state distribution: (post-selected on expected qubit states, not renormalized)
|x1,x2,f1,f2>

	|0,0,0,1>	|0,1,0,0>	|1,0,1,1>	|1,1,1,0>
|0,0,0,1>	0.003086	0.003086	0.003086	0.003086

The distribution computed with Perceval is uniform over each outcome, which corresponds to the expected distribution obtained in the paper when \(x_0 = 0\).

Interpretation of the outcomes

For each outcome we consider the values of qubits \(x_2, x_1, x_0\) (where \(x_0\) is 0) which represent a binary number between 0 and 7, here corresponding to 0, 4, 2 and 6 in the order of the table above. After sampling the circuit, obtaining outcomes 2 or 6 allows to successfully compute the order \(r=4\). Obtaining outcome 0 is an expected failure of the quantum circuit, inherent to Shor’s algorithm. Outcome 4 is an expected failure as well, as it only allows to compute the trivial factors 1 and 15.

Since the distribution is uniform the circuit successfully yields a successful outcome with probability 1/2. This probability can be amplified exponentially close to \(1\) by sampling the circuit multiple times.

Reference

Alberto Politi, Jonathan C.F. Matthews, and Jeremy L. O’brien. “Shor’s quantum factoring algorithm on a photonic chip.” Science 325.5945 (2009): 1221-1221.