An implementation of Shor’s algorithm using Quantum.NET

Shor’s algorithm, named after mathematician Peter Shor, is the most commonly cited example of quantum algorithm. It solves the integer factorization problem in polynomial time, substantially faster than the most efficient known classical factoring algorithm, the general number field sieve, which works in sub-exponential time.

It owes its fame to its forecasted ability to break public-key cryptography schemes, which are based on the assumption that factoring large numbers is computationally infeasible. However, the greatest number so far factorized using Shor’s algorithm, 21, demonstrates how the technical difficulty of building a quantum computer is a limiting factor to real-life applications of quantum computing.

The different steps of Shor’s algorithm to factor an integer \(n\) are as follows:

  1. Pick \(a < n\) randomly using a uniform distribution.
  2. Calculate the greatest common divisor of \(a\) and \(n\) using Euclid’s algorithm. If \(\gcd(a, n) \neq 1\), this is a non-trivial factor of \(n\) and stop.
  3. Pick \(m\) such that \(m = 2^l\) for some integer \(l\) and \(n^2 \leq m < 2n^2\).
  4. Prepare two quantum registers in compound state \(\left|0\right\rangle \left|0\right\rangle\). Let them be large enough to represent numbers up to \(m-1\), i.e. of length \(l\).
  5. Apply the Hadamard transform to the first register.
    The resulting state is \(\frac{1}{\sqrt{m}} \sum_{x=0}^{m-1} \left|x\right\rangle \left|0\right\rangle\).
  6. Apply a map defined by \(\left|x\right\rangle \left|0\right\rangle \rightarrow \left|x\right\rangle \left|a^x \bmod n\right\rangle\).
    This results in state \(\frac{1}{\sqrt{m}} \sum_{x=0}^{m-1} \left|x\right\rangle \left|a^x \bmod n\right\rangle\).
  7. Apply the quantum Fourier transform to the first register, which becomes \(\frac{1}{\sqrt{m}} \sum_{y=0}^{m-1} e^\frac{2 \pi x y i}{m} \left|y\right\rangle\).
    This results in global state \(\frac{1}{m} \sum_{x=0}^{m-1} \sum_{y=0}^{m-1}\) \(e^\frac{2 \pi x y i}{m} \left|y\right\rangle \left|a^x \bmod n\right\rangle\).
    The goal here is to find the period \(p\) of \(x \rightarrow a^x \bmod n\). Let \(x = p q + r\) with \(r < p\) and \(q < \lfloor\frac{m – r}{p}\rfloor\). The state can be rewritten as the following: \(\frac{1}{m} \sum_{y=0}^{m-1} \sum_{r=0}^{p-1} \sum_{q=0}^{\lfloor\frac{m – r – 1}{p}\rfloor}\) \(e^\frac{2 \pi (p q + r) y i}{m} \left|y\right\rangle \left|a^r \bmod n\right\rangle\).
  8. Observe the first register. The probability of observing a certain state \(\left|y\right\rangle \left|a^r \bmod n\right\rangle\) is \(\left| \frac{1}{m} \sum_{q=0}^{\lfloor\frac{m – r – 1}{p}\rfloor} e^\frac{2 \pi (p q + r) y i}{m} \right|^2\) \(= \frac{1}{m^2} \left| \sum_{q=0}^{\lfloor\frac{m – r – 1}{p}\rfloor} e^\frac{2 \pi p q y i}{m} \right|^2\).
    Summing the geometric series gives \(\frac{1}{m^2} \left| \frac{1 – e^\frac{2 \pi p \lfloor\frac{m – r}{p}\rfloor y i}{m}}{1 – e^\frac{2 \pi p y i}{m}} \right|^2\) \(= \frac{1 – cos(\frac{2 \pi p \lfloor\frac{m – r}{p}\rfloor y}{m})}{m^2 (1 – cos(\frac{2 \pi p y}{m}))}\), which attains its maximal values when \(\frac{p y}{m}\) is the closest to an integer. Therefore, there is a high probability of observing \(y\) close to a multiple of \(\frac{m}{p}\) in the first register.
  9. Perform continued fraction expansion on \(\frac{y}{m}\) to approximate it to \(\frac{k}{p’}\) with \(p’ < n\). The denominator \(p’\) is therefore likely to be the period \(p\) or a factor of it.
  10. Verify that \(a^{p’} \equiv 1 \bmod n\). If so, \(p = p’\): go to step 13.
  11. Otherwise, try the same thing with multiples of \(p’\). If it works, this gives \(p\): go to step 13.
  12. Otherwise, go back to step 4.
  13. If \(p\) is odd, go back to step 1.
  14. If \(a^{\frac{p}{2}} \equiv −1 \bmod n\), go back to step 1.
  15. Compute \(\gcd(a^{\frac{p}{2}} + 1, n)\) and \(\gcd(a^{\frac{p}{2}} – 1, n)\). They are both nontrivial factors of \(n\).

Now the process has been explained, we’ll start with the interface of our implementation. Our algorithm will be run in the following way:

Note the use of dependency injection to avoid issues with pseudorandom number generator determinism.

The corresponding structure of the algorithm class will be as follows:

Shor’s algorithm can be divided in two main sections: the classical part and the quantum period-finding subroutine. The implementation of the classical part is rather straightforward and can be done as follows:

The greatest common divisor is calculated recursively using Euclid’s algorithm:

And we use exponentiation by squaring to calculate integer powers:

Exponentiation by squaring uses the fact that \(x^n= x (x^{2})^{\frac{n – 1}{2}}\) if \(n\) is odd and \((x^{2})^{\frac{n}{2}}\) if \(n\) is even, allowing powers to be calculated using multiplication only. A quick benchmark shows, without much surprise, that this approach is more than three times as fast as using .NET’s native double Math.Pow(double, double) implementation and casting it to an int.

Next is the period-finding subroutine, the quantum part of Shor’s algorithm. We will use Quantum.NET to simulate the required quantum circuit, shown below:

The quantum circuit used in Shor's algorithm

Initializing the input and output quantum registers is as simple as it gets:

And so is applying the Hadamard transform to the first register:

We use the following method to calculate \(m\):

Another quick benchmark shows it to be about 3.6 times as fast as the equivalent approach using .NET’s logarithm, (int) Math.Log(n * n - 1, 2) + 1. It works by removing chunks of 16, 8, 4, 2 and 1 bit(s) (we’re working on 32-bit integers) when possible and counting the number of bits removed.

Now to the crux of Shor’s algorithm, but also the bottleneck in terms of runtime: modular exponentiation. We have to build a gate (more like a sub-circuit) that applies the map defined by \(\left|x\right\rangle \left|0\right\rangle \rightarrow \left|x\right\rangle \left|a^x \bmod n\right\rangle\) to the full register.

Building such a circuit is a very complex matter and there are many approaches, but all start with a decomposition into modular multipliers. We use the fact that \(a^x \bmod n\) \(= \prod_{i=0}^{l-1} (a^{2^{i}} \bmod n)^{x_i} \bmod n\), where the \(x_i\) are the bits of \(x\) such that \(x = \sum_{i=0}^{l-1} x_i 2^i\), to split the circuit into a series of controlled modular multiplication gates, each using a different qubit from the input register as a control qubit.

However, for the sake of simplicity, we will represent this circuit with a unique modular exponentiation gate. The map leaves the input registry unchanged, therefore it can be represented by a block diagonal matrix of total order \(m^2\):

\(\mathbf{M} = \begin{bmatrix}
\mathbf{M}_{0} & 0 & \cdots & 0 \\ 0 & \mathbf{M}_{1} & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & \mathbf{M}_{m – 1}

It also overrides the output registry regardless of its initial value (with the result of the modular exponentiation of the input registry), therefore the elements of \(\mathbf{M}_{k}\) are \(\mathbf{{m}_{k}}_{i, j} = 1\) if \(i = a^k \bmod n\) and \(0\) otherwise, regardless of \(j\).

The resulting permutation matrix is built using the following:

Modular exponentiation is performed using a slightly modified version of exponentiation by squaring:

We use the newly created permutation matrix to create a quantum gate we apply to the global register, formed from the input and output registers.

The registers are now entangled. The next step consists in applying the quantum Fourier transform to the input register, but it needs to be processed along with the output register in order to preserve entanglement. The solution consists in creating a quantum gate from two smaller ones: a quantum Fourier transform for the input register and an identity gate for the output register, which will be left unchanged.

It is now time to observe the first register so that it collapses into a pure state, and read its value. fullRegister is of length 2 * registerLength, so we use the portionStart and portionLength parameters of QuantumRegister.GetValue to only read the first registerLength bits, returning a value \(y\):

The idea is now to perform continued fraction expansion on \(y / m\). As explained above, the denominator of the fraction approximation is likely to be the period \(p\) or a factor of it.

Every irrational number can be represented in precisely one way as an infinite continued fraction, which is an expression of the form:

\(a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{a_3 + {_\ddots}}}}\)

Limiting the expression to a certain number of terms gives rational approximations called the convergents of the continued fraction. They are \(a_0 = \frac{a_0}{1}\), \(a_0 + \frac{1}{a_1} = \frac{a_1 a_0 + 1}{a_1 }\), \(a_0 + \frac{1}{a_1 + \frac{1}{a_2}} = \frac{a_2 (a_1 a_0 + 1) + a_0}{a_2 a_1 + 1}\), and can generally be expressed as \(\frac{h_n}{k_n}\), with the numerators and denominators recursively defined by \(h_n = a_n h_{n − 1} + h_{n − 2}\) and \(k_n = a_n k_{n − 1} + k_{n − 2}\).

This property allows us to find the greatest denominator \(p’\) such that\(p’ < n\), using the following implementation:

Now we have a candidate for \(p\), we verify if \(p’\) or one of its multiples is the period we’re looking for.

If that’s the case, then we are done with the period-finding subroutine, and Shor’s algorithm as a whole. If not, we start the subroutine over:

This article now comes to an end. We have successfully modeled Shor’s algorithm, and while it won’t help us factor anything (the classical alternatives being incomparably faster than simulating a system whose complexity grows exponentially with regard to the size of the input), it provides us with great insight as to how the algorithm can be implemented in practice.

RSA is not broken yet.

Quantum.NET v1.1.0 is now available

Version 1.1.0 of Quantum.NET has been released today. It introduces several new features aimed at facilitating the building of more complex quantum circuits and the obtaining of the results they yield.

First is the introduction of identity gates:

  • QuantumGate.IdentityGate
  • QuantumGate.IdentityGateOfLength(int registerLength)

While useless in themselves, they become very powerful when combined with another new feature, quantum gate stacking, or the creation of quantum gates from other, smaller quantum gates (variadic constructor, also works with QuantumGate[] and IEnumerable<QuantumGate>):

This allows to apply a gate to a subsection of a quantum register, leaving the rest unchanged while preserving entanglement.

Another new feature is the reading of the values contained in pure-state quantum registers. Once an algorithm is run, we will be left with a register to observe which, once collapsed, will contain the data we need. Optional offset and length parameters can be used to read a subsection of the quantum register only, as often required.

Quantum.NET is available as a NuGet package under Lachesis.QuantumComputing, and the source code can be found on GitHub at phbaudin/quantum-computing.

Quantum.NET: a new quantum computing library in C#

I’ve just published the first version of my latest project, a quantum computing library in C# called Quantum.NET that allows the manipulation of qubits and the modeling of quantum circuits.

It is available as a NuGet package under Lachesis.QuantumComputing, and the source code can be found on GitHub at phbaudin/quantum-computing.

The way it works is pretty straightforward. A qubit can be created from its probability amplitudes:

Or from its amplitudes’ real and imaginary parts:

It can also be created from its colatitude and longitude on the Bloch sphere:

Shortcuts are available for notable qubits:

A qubit is a quantum register of length 1 and can be manipulated as such. The Qubit class is merely a subclass of QuantumRegister designed for ease of use:

Shortcuts are also available for notable quantum registers:

A quantum register can be created from other quantum registers (variadic constructor, also works with QuantumRegister[] and IEnumerable<QuantumRegister>):

Or from the 2n complex probability amplitudes of each of its pure states (variadic constructor, also works with Complex[] and IEnumerable<Complex>):

Quantum registers are mostly used to represent numbers and can therefore be created from integers (this will naturally generate pure states):

A quantum register can be observed and collapse into a pure state (note: use your own Random instance to avoid issues with pseudorandom number generator determinism):

Quantum gates are required to operate on quantum registers. Shortcuts are also available for notable quantum gates:

  • QuantumGate.HadamardGate
  • QuantumGate.HadamardGateOfLength(int registerLength)
  • QuantumGate.NotGate
  • QuantumGate.PauliYGate
  • QuantumGate.PauliZGate
  • QuantumGate.SquareRootNotGate
  • QuantumGate.PhaseShiftGate(double phase)
  • QuantumGate.SwapGate
  • QuantumGate.SquareRootSwapGate
  • QuantumGate.ControlledNotGate
  • QuantumGate.ControlledGate(QuantumGate gate)
  • QuantumGate.ToffoliGate
  • QuantumGate.FredkinGate
  • QuantumGate.QuantumFourierTransform(int registerLength)

Quantum gates can also be created from a bidimensional array of complex numbers:

Applying a quantum gate to a quantum register is as simple as using the multiplication operator on them:

Unary gates only operate on one qubit, binary gates on two, etc.:

Please note that this is not a literal arithmetic library. While measures have been taken to circumvent a range of errors caused by floating-point precision, the use of QuantumRegister.AlmostEquals might be required in places:

Hope you find this useful!