# 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 randoly 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  + pe[qubitState] + pe[qubitState] + [0, 0] + pe[qubitState] + pe[qubitState] + 

:

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]+fockState[q] != 1:
return None
else:
L.append(fockState[q])
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]:

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

# H gates
for q in [f1, f2]:

:

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")
simulator = backend(circ)


### 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 = simulator.probampli(pcvl.BasicState(istate), 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> 1.734723475976807e-18j
|0,0,0,1> (0.05555555555555555+2.8912057932946782e-18j)
|0,0,1,0> 0j
|0,0,1,1> 0j
|0,1,0,0> (-0.05555555555555555+0j)
|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> (-6.938893903907228e-18+0j)
|1,0,1,0> 0j
|1,0,1,1> (-0.05555555555555555+0j)
|1,1,0,0> (2.0816681711721685e-17+0j)
|1,1,0,1> (2.7755575615628914e-17+0j)
|1,1,1,0> (0.05555555555555555+3.046145044787777e-18j)
|1,1,1,1> (6.938893903907228e-18+1.734723475976807e-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, obtainig 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.