Uncategorized

# Programming Quantum Period Finding (Shor’s Algorithm)

For the last post of my quantum programming series, here I show you the most famous quantum algorithm “Shor’s algorithm” with Q# programming.
As I describe later, you can also solve integer factorization with polynomial time computation (in the input’s bit size) using this quantum algorithm, though classical one needs exponential time.

As well as other posts, first I describe the outline (ideas) of this algorithm, and later we see the sample program with Q# along with this algorithm.

Before starting this post, I recommend you to read my early post “Programming QFT and Quantum Phase Estimation” and “Programming Quantum Arithmetics“, if you’re not familiar with Quantum Fourier Transform and other related topics.

# Quantum Period-Finding (Shor’s Algorithm)

If positive integer $N$ and $a$ is coprime, number theory tells us there exists $r > 1$ such as $a^r = 1\:mod\:N$. We now assume that $r$ is positive least number such that $a^r = 1\:mod\:N$.
As you can see :

$a^{r+b}\:mod\:N = (a^r \cdot a^b)\:mod\:N = a^b\:mod\:N$

Therefore $a^x\:mod\:N$ is periodic (cyclic) for input x by r-cycle. (If there exists more least integer $r_0 < r$ such that $a^{x_0} = a^{x_0+r_0}\:mod\:N$, this implies $a^{r_0} = 1\:mod\:N$ and $r$ is not the least integer.)
For instance, if N = 11 and a = 5, $a^x\:mod\:N\;(x=1, 2, \ldots)$ equals to 5, 3, 4, 9, 1, 5, 3, 4, 9, 1, 5, 3, … The period is 5.

$5^1\:mod\:11 = 5$
$5^2\:mod\:11 = 3$
$5^3\:mod\:11 = 4$
$5^4\:mod\:11 = 9$
$5^5\:mod\:11 = 1$
$5^6\:mod\:11 = 5$
$5^7\:mod\:11 = 3$
$5^8\:mod\:11 = 4$

This period $r$ is called “order” of $a$ modulus $N$.

Our concern is how to find the period (order) for given coprime integers $N > 1$ and $a > 1$.
Of course, you can iteratively search by calculating with 1 to N – 1, but it needs $O(2^n)$ computation time (i.e, exponential time), where the integer N has n bit’s size.
Here we want more efficient algorithms, and Shor showed us polynomial-time algorithm with quantum computer as follows.

The outline of Shor’s algorithm is : transform to all possible values using quantum superpositions, apply some N modulus operations for these qubits, and then finally estimate the period $r$ using measured values.
Now let’s see the details of this algorithm and see why this can solve this problem.

In this algorithm, we need the following 2 register :

• The first register $\left| x \right>$ has n qubits for integer $n$ such that $N^2 \leq 2^n < 2 N^2$ .
• The second register $\left| y \right>$ has at least $\log_2 N$ qubits for storing $\left|a^x\:mod\:N\right>$.

To simplify explanation, here we assume that $2^n/r$ is some integer $m$, i.e, $2^n/r=m$. Note that $r$ is unknown for us now. (Later we will discuss about the case in which $2^n/r$ is not an integer…)

First we generate the superposition for x, i.e, $\left| x \right> \to H^{\otimes n} \left| x \right>$.
Then the state of $\left| x \right>$ becomes :

$\displaystyle \frac{1}{\sqrt{2^n}} \sum_{i=0}^{2^n - 1} \left| i \right>$

Next we transform $\left| y \right>$ as follows.
For this implementation, please see my previous post “Programming Quantum Arithmetics“.

$\left| x \right> \left| y \right> = \left| x \right> \left| 0 \right> \to \left| x \right>\left|a^x\:mod\:N\right>$

Now we measure the register $\left| y \right>$, but we don’t need the measurement results. Once we’ve measured $\left| y \right>$, the register $\left| y \right>$ is garbage and we don’t need $\left| y \right>$ any more.

After you’ve measured the register $\left| y \right>$, $\left| y \right>$ results into $\left|a^{x_0}\:mod\:N\right>$ for some $x_0$.
Hence $\left| x \right>$ results into the following state, because $2^n/r=m$.

$\displaystyle \frac{1}{\sqrt{m}} \sum_{j=0}^{m-1} \left| x_0 + jr \right>$

Now we apply Quantum Fourier Transform (QFT) for register $\left| x \right>$. (See my early post for QFT.) Then $\left| x \right>$ results into the following state :

$\displaystyle \frac{1}{\sqrt{2^n}} \frac{1}{\sqrt{m}} \sum_{k=0}^{2^n-1} \sum_{j=0}^{m-1} \exp(\frac{2\pi i (x_0+jr)k}{2^n}) \left| k \right>$

$\displaystyle = \frac{1}{\sqrt{2^n}} \frac{1}{\sqrt{m}} \sum_{k=0}^{2^n-1} \Biggl[ \sum_{j=0}^{m-1} \exp(\frac{2\pi ijrk}{2^n}) \Biggr] \exp(\frac{2\pi i x_0 k}{2^n}) \left| k \right>$

$\displaystyle = \frac{1}{\sqrt{2^n}} \frac{1}{\sqrt{m}} \sum_{k=0}^{2^n-1} \Biggl[ \sum_{j=0}^{m-1} \exp(\frac{2\pi ijk}{m}) \Biggr] \exp(\frac{2\pi i x_0 k}{2^n}) \left| k \right>$

Now let me fix some $k$ ( $k \in \{0, \ldots ,2^n-1\}$ ) such that $k/m$ is not integer. In this case, the terms in the square brackets in the above equation cancels each others and then results into 0, because $m = 2^p$ for some integer p. As a result, only terms with $k$, in which $k/m$ is integer, resides in the above equation.

Next let me fix some $k$ ( $k \in \{0, \ldots ,2^n-1\}$ ) such that $k/m$ is integer. In this case, $\exp(\frac{2\pi ijk}{m})$ equals to 1. And the number of such $k$ is r, such as $0,\:\frac{2^n}{r},\:\frac{2^n}{r} \times 2,\:\ldots\:,\:\:\frac{2^n}{r} \times (r - 1)$.
Therefore we can replace above equation with $k = \frac{2^n}{r} \times l \;\; (l = 0,\:\ldots\:, r-1)$, and we get :

$\displaystyle = \frac{1}{\sqrt{2^n}} \frac{1}{\sqrt{m}} \sum_{l=0}^{r-1} m \exp(\frac{2\pi i x_0 \times 2^n l / r}{2^n}) \left| \frac{2^n l}{r} \right>$

$\displaystyle = \frac{1}{\sqrt{r}} \sum_{l=0}^{r-1} \exp(\frac{2\pi i x_0 l}{r}) \left| \frac{2^n l}{r} \right>$

Now measure $\left| x \right>$ and we assume that the results is some integer $C$.

As you can easily see, we get $C = \frac{2^n l_0}{r}$ for some unknown integer $l_0$.
If $l_0$ and $r$ is coprime, you can get $r$ using greatest common divisor (by Euclidean algorithm) with known $C$ and $2^n$.
But it’s not so simple. Please remember our assumption : $2^n = mr$. For instance, when $r$ is odd, you cannot simply get $r$ by GCD.

When $2^n \neq mr$, it’s known that $r$ is approximated by the following $r^{\prime}$ using continued fraction expansion. If the following $r^{\prime}$ and $d$ is coprime, then you can get $r=r^{\prime}$. If it’s not coprime (check if the result $r^{\prime}$ is an order or not), you repeat this algorithm until you find an order $r$. (Here I don’t describe about this approximation, but please refer p.18 in the original paper “Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer” by Peter Shor.)

$\displaystyle \left| \frac{C}{2^n} - \frac{d}{r^{\prime}} \right| < \frac{1}{2 \cdot 2^n}$
where $r^{\prime} < N$ and $d \in \mathbb{N}$

Note : There exists another approach by phase estimation (see my early post for phase estimation) for approximation.
See “Applications – Shor’s Algorithm” in Q# document for details.

Number theory tells us that this takes place with probability at least $1 / \log \log r$ (see Hardy-Ramanujan theorem), and the entire period-fining (order-finding) has $O(n^2 \log n \log \log n)$ computation complexity. That is, this algorithm has polynomial computation complexity in the input size n.

# Programming Quantum Period-Finding

As I described in my previous post, I note that my code is written straight forward in order to help you understand the logic flow. (For instance, so many QFTs are called inside my function. When you have learned core essence in this post, please proceed to optimize your code.)

Now let’s start Q# programming along with above algorithm.

First we prepare the register $\left| x \right>$ and $\left| y \right>$, which states are all initialized as $\left| 0 \right>$.

// Get least integer n1 such as : num^2 <= 2^n1
let n1 = BitSizeI(num) * 2;
let n2 = BitSizeI(num);
...

using ((x, y) = (Qubit[n1], Qubit[n2])) {
...
}

We apply $\left| x \right> \to H^{\otimes n} \left| x \right>$ and generate superposition as follows.

...
using ((x, y) = (Qubit[n1], Qubit[n2])) {
Microsoft.Quantum.Canon.ApplyToEachCA(H, x);
...
}

Next we apply $\left| x \right> \left| y \right> = \left| x \right> \left| 0 \right> \to \left| x \right>\left|a^x\:mod\:N\right>$.

This is not simple, but we’ve already learned algorithms for quantum arithmetic operations and saw how to implement this. Here I don’t describe again, but please refer my previous post. (In the following code, QuantumExponentForPeriodFinding() is just a modified version of QuantumExponentByModulus() in my previous post, and the following QuantumMultiplyByModulus() is also the same one which has already been implemented in my previous post.)

...
using ((x, y) = (Qubit[n1], Qubit[n2])) {
Microsoft.Quantum.Canon.ApplyToEachCA(H, x);

// |x⟩ |0 (=y)⟩ -> |x⟩ |a^x mod N⟩
QuantumExponentForPeriodFinding(a, num, x, y);
...
}
...
// Implement : |x⟩ |0 (=y)⟩ -> |x⟩ |a^x mod N⟩ for some integer a
// (where y should be |0⟩)
// This is modified version of QuantumExponentByModulus() in my post.
operation QuantumExponentForPeriodFinding (a : Int, N : Int, x : Qubit[], y : Qubit[]) : Unit {
let n1 = Length(x);
let n2 = Length(y);

// set |y⟩ = |0...01⟩
X(y[n2 - 1]);

for(idx in 0 .. n1 - 1) {
// a^(2^((n1-1) - idx)) is too big, then we reduce beforehand
mutable a_mod = 1;
for(power in 1 .. 2^((n1-1) - idx)) {
set a_mod = (a_mod * a) % N;
}
// apply decomposition elements
Controlled QuantumMultiplyByModulus([x[idx]], (N, a_mod, y));
}
}

Next we measure $\left| y \right>$ register and apply Quantum Fourier Transform (QFT) to $\left| x \right>$ register.
See my early post “Programming QFT and Quantum Phase Estimation” for the following QFTImpl() operation. (You can also use Q# built-in operation for QFT.)

  ...

// measure y and reset
mutable tmpResult = new Result[n2];
for (idx in 0 .. n2 - 1) {
set tmpResult w/= idx <-MResetZ(y[idx]);
}

// QFT for x
QFTImpl(x);
...

Now we measure the register $\left| x \right>$, and set results into variable “realResult“.

  // Measure x and reset
mutable realResult = new Result[n1];
for (idx in 0 .. n1 - 1) {
set realResult w/= idx <-MResetZ(x[idx]);
}

As I mentioned above, we should get approximated fraction $d / r^{\prime}$ (where $r^{\prime} \leq N$) from the fraction $C / 2^n$ (where C is measured value, i.e, above “realResult”) by running continued fraction expansion.
Here we use Microsoft.Quantum.Math.ContinuedFractionConvergentL() in Q# library for continued fraction, but you can also implement Euclidean by yourself. (We have also used Microsoft.Quantum.Math.GreatestCommonDivisorL() in Q# library for GCD.)

  ...

// get integer's result from measured array (ex : |011⟩ -> 3)
let resultBool = [false] + Microsoft.Quantum.Convert.ResultArrayAsBoolArray(realResult); // for making unsigned positive integer, add first bit
let resultBool_R = Microsoft.Quantum.Arrays.Reversed(resultBool); // because BoolArrayAsBigInt() is Little Endian order
let resultIntL = Microsoft.Quantum.Convert.BoolArrayAsBigInt(resultBool_R);

// get period candidate by continued fraction expansion (thanks to Euclid !)
let gcdL = GreatestCommonDivisorL(resultIntL, 2L^n1);
let calculatedNumerator = resultIntL / gcdL;
let calculatedDenominator = 2L^n1 / gcdL;
let numL = Microsoft.Quantum.Convert.IntAsBigInt(num);
let approximatedFraction =
ContinuedFractionConvergentL(BigFraction(calculatedNumerator, calculatedDenominator), numL);
let (approximatedNumerator, approximatedDenominator) = approximatedFraction!;
mutable periodCandidateL = 0L;
if(approximatedDenominator < 0L) {
set periodCandidateL = approximatedDenominator * -1L;
}
else {
set periodCandidateL = approximatedDenominator;
}
set periodCandidate = ReduceBigIntToInt(periodCandidateL);

// output for debugging
Message($"Measured Fraction : {resultIntL} / {2L^n1}"); Message($"Approximated Fraction : {approximatedNumerator} / {approximatedDenominator}");
Message($"Period Candidate : {periodCandidate}"); ... // This is helper function to convert BigInt to Int ... operation ReduceBigIntToInt(numL : BigInt) : Int { // Check if numL is not large Microsoft.Quantum.Diagnostics.Fact(BitSizeL(numL) <= 32,$"Cannot convert to Int. Input is too large");

mutable resultInt = 0;
let numArray = Microsoft.Quantum.Convert.BigIntAsBoolArray(numL);
let numArray_R = Microsoft.Quantum.Arrays.Reversed(numArray); // because BigIntAsBoolArray() is Little Endian order
let nSize = Length(numArray_R);
for (idx in 0 .. nSize - 1) {
if(numArray_R[idx] and ((nSize - 1) - idx <= 31)) {
set resultInt = resultInt + (2 ^ ((nSize - 1) - idx));
}
}
return resultInt;
}

Finally you check if $r^{\prime}$ (above “periodCandidate”) is the period or not. If not, please repeat once again.

As a result, the complete Q# code is as follows.

open Microsoft.Quantum.Intrinsic;
open Microsoft.Quantum.Math;
open Microsoft.Quantum.Measurement;

operation QuantumPeriodFinding (num : Int, a : Int) : Unit {
// Get least integer n1 such as : num^2 <= 2^n1
let n1 = BitSizeI(num) * 2;
let n2 = BitSizeI(num);
mutable periodCandidate = 1;
repeat {
using ((x, y) = (Qubit[n1], Qubit[n2])) {
Microsoft.Quantum.Canon.ApplyToEachCA(H, x);

// |x⟩ |0 (=y)⟩ -> |x⟩ |a^x mod N⟩
QuantumExponentForPeriodFinding(a, num, x, y);

// measure y and reset
mutable tmpResult = new Result[n2];
for (idx in 0 .. n2 - 1) {
set tmpResult w/= idx <-MResetZ(y[idx]);
}

// QFT for x
QFTImpl(x);

// Measure x and reset
mutable realResult = new Result[n1];
for (idx in 0 .. n1 - 1) {
set realResult w/= idx <-MResetZ(x[idx]);
}

// get integer's result from measured array (ex : |011⟩ -> 3)
let resultBool = [false] + Microsoft.Quantum.Convert.ResultArrayAsBoolArray(realResult); // for making unsigned positive integer, add first bit
let resultBool_R = Microsoft.Quantum.Arrays.Reversed(resultBool); // because BoolArrayAsBigInt() is Little Endian order
let resultIntL = Microsoft.Quantum.Convert.BoolArrayAsBigInt(resultBool_R);

// get period candidate by continued fraction expansion (thanks to Euclid !)
let gcdL = GreatestCommonDivisorL(resultIntL, 2L^n1);
let calculatedNumerator = resultIntL / gcdL;
let calculatedDenominator = 2L^n1 / gcdL;
let numL = Microsoft.Quantum.Convert.IntAsBigInt(num);
let approximatedFraction =
ContinuedFractionConvergentL(BigFraction(calculatedNumerator, calculatedDenominator), numL);
let (approximatedNumerator, approximatedDenominator) = approximatedFraction!;
mutable periodCandidateL = 0L;
if(approximatedDenominator < 0L) {
set periodCandidateL = approximatedDenominator * -1L;
}
else {
set periodCandidateL = approximatedDenominator;
}
set periodCandidate = ReduceBigIntToInt(periodCandidateL);

// output for debugging
Message($"Measured Fraction : {resultIntL} / {2L^n1}"); Message($"Approximated Fraction : {approximatedNumerator} / {approximatedDenominator}");
Message($"Period Candidate : {periodCandidate}"); } } until ((periodCandidate != 0) and (ExpModI(a, periodCandidate, num) == 1)) fixup { } // output for debugging Message("Found period " + Microsoft.Quantum.Convert.IntAsString(periodCandidate)); Message(""); } // Implement : |x⟩ |0 (=y)⟩ -> |x⟩ |a^x mod N⟩ for some integer a // (where y should be |0⟩) // This is modified version of QuantumExponentByModulus() in my post. // See https://tsmatz.wordpress.com/2019/05/22/quantum-computing-modulus-add-subtract-multiply-exponent/ operation QuantumExponentForPeriodFinding (a : Int, N : Int, x : Qubit[], y : Qubit[]) : Unit { let n1 = Length(x); let n2 = Length(y); // set |y⟩ = |0...01⟩ X(y[n2 - 1]); for(idx in 0 .. n1 - 1) { // a^(2^((n1-1) - idx)) is too big, then we reduce beforehand mutable a_mod = 1; for(power in 1 .. 2^((n1-1) - idx)) { set a_mod = (a_mod * a) % N; } // apply decomposition elements Controlled QuantumMultiplyByModulus([x[idx]], (N, a_mod, y)); } } // This is helper function to convert BigInt to Int ... operation ReduceBigIntToInt(numL : BigInt) : Int { // Check if numL is not large Microsoft.Quantum.Diagnostics.Fact(BitSizeL(numL) <= 32,$"Cannot convert to Int. Input is too large");

mutable resultInt = 0;
let numArray = Microsoft.Quantum.Convert.BigIntAsBoolArray(numL);
let numArray_R = Microsoft.Quantum.Arrays.Reversed(numArray); // because BigIntAsBoolArray() is Little Endian order
let nSize = Length(numArray_R);
for (idx in 0 .. nSize - 1) {
if(numArray_R[idx] and ((nSize - 1) - idx <= 31)) {
set resultInt = resultInt + (2 ^ ((nSize - 1) - idx));
}
}
return resultInt;
}

Now let’s invoke from your C# program (or Python).

using Microsoft.Quantum.Simulation.Simulators;

static void Main(string[] args)
{
// For Period Finding Only
using (var qsim = new QuantumSimulator())
{
int N, a;

N = 11;
a = 5;
QuantumPeriodFinding.Run(qsim, N, a).Wait();

N = 15;
a = 7;
QuantumPeriodFinding.Run(qsim, N, a).Wait();
}
}

# Integer Factorization

Using above quantum period-finding, we can extend this algorithm to integer factorization, which is also difficult to be solved by classical efficient algorithms.

Note : As you can easily see, it needs $O(2^n)$ when you apply the trivial iteration search for integer factorization, where the integer has n bit’s size. There exist several challenges (studies) for reducing computation complexity with classical algorithms, but it is known that no existence with n-polynomial complexity for a n-bits number. (See “Integer factorization” in Wikipedia.)

When r is even and $a^{r/2}\:mod\:N \neq N - 1$, it’s known that at least one of $gcd(N, a^{r/2} \pm 1)$ is a non-trivial factor of N. Then overall algorithm of integer factorization for integer $N$ is as follows.
The following step 5 is quantum step, and others are classical parts.

1. If $N$ is even, return 2 as a factor.
If not, proceed to the next step.
2. Check if $N = j^k$ for some integer $j > 1$ and $k > 1$. (Check $\sqrt[k]{N}$ for all $k \leq \log_2 N$. This is n-polynomial computation for n-bit’s integer N.)
If so, return $j = \sqrt[k]{N}$ as a factor. If not, proceed to the next step.
3. Pick a random integer $1 < a < N$ and get the great common divisor $gcd(N,a)$ with Euclidean algorithm.
4. If gcd(N,a) > 1, return gcd(N,a) as a factor.
If not, proceed to the next step.
5. If gcd(N,a) = 1, $N$ and $a$ is coprime and then run above quantum period-finding algorithm for $N$ and $a$. (i.e, find the least integer $r$ such as $a^r\:mod\:N = 1$.)
6. If the period $r$ is odd, go back to step 3.
7. If the period $r$ is even and $a^{r/2}\:mod\:N = N - 1$, go back to step 3.
8. Otherwise, check if $gcd(a^{r/2} \pm 1,N)$ is a nontrivial factor of $N$. If not, go back to step 3.

# Some Note

With current quantum apparatus, unfortunately it’s difficult to run this Shor’s algorithm for a large number with scalability. (Remember that quantum system has some amount of errors due to several noises in the current apparatus.)
But in the near future, the quantum technology’s adoption for the production will come closer to reality.

At that time, this algorithm (Shor’s algorithm) would have significant impacts for current technologies such as cryptography, because these current technologies depend on the difficulty of factoring large composite integers by classical approaches.
For the coming future, there are also several studies for the quantum key distribution’s techniques. (See here.)

Source Code :

https://github.com/tsmatz/Quantum-Shor-Algorithm

Reference :

Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer

Principles of Quantum Computation and Information