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 and is coprime, number theory tells us there exists such as . We now assume that is positive least number such that .

As you can see :

Therefore is periodic (cyclic) for input x by r-cycle. (If there exists more least integer such that , this implies and is not the least integer.)

For instance, if N = 11 and a = 5, equals to 5, 3, 4, 9, 1, 5, 3, 4, 9, 1, 5, 3, … The period is 5.

…

This period is called “order” of modulus .

Our concern is how to find the period (order) for given coprime integers and .

Of course, you can iteratively search by calculating with 1 to N – 1, but it needs 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 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 has n qubits for integer such that .
- The second register has at least qubits for storing .

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

First we generate the superposition for x, i.e, .

Then the state of becomes :

Next we transform as follows.

For this implementation, please see my previous post “Programming Quantum Arithmetics“.

Now we measure the register , but we don’t need the measurement results. Once we’ve measured , the register is garbage and we don’t need any more.

After you’ve measured the register , results into for some .

Hence results into the following state, because .

Now we apply Quantum Fourier Transform (QFT) for register . (See my early post for QFT.) Then results into the following state :

Now let me fix some ( ) such that 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 for some integer p. As a result, only terms with , in which is integer, resides in the above equation.

Next let me fix some ( ) such that is integer. In this case, equals to 1. And the number of such is r, such as .

Therefore we can replace above equation with , and we get :

Now measure and we assume that the results is some integer .

As you can easily see, we get for some unknown integer .

If and is coprime, you can get using greatest common divisor (by Euclidean algorithm) with known and .

But it’s not so simple. Please remember our assumption : . For instance, when is odd, you cannot simply get by GCD.

When , it’s known that is approximated by the following using continued fraction expansion. If the following and is coprime, then you can get . If it’s not coprime (check if the result is an order or not), you repeat this algorithm until you find an order . (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.)

where and

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 (see Hardy-Ramanujan theorem), and the entire period-fining (order-finding) has 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 and , which states are all initialized as .

```
// 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 and generate superposition as follows.

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

Next we apply .

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.
// 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));
}
}
```

Next we measure register and apply Quantum Fourier Transform (QFT) to 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 , 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 (where ) from the fraction (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 (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;
// answer should be 5
N = 11;
a = 5;
QuantumPeriodFinding.Run(qsim, N, a).Wait();
// answer should be 4
N = 15;
a = 7;
QuantumPeriodFinding.Run(qsim, N, a).Wait();
}
}
```

Please download my source code : https://github.com/tsmatz/Quantum-Shor-Algorithm

# 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 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 , it’s known that at least one of is a non-trivial factor of N. Then overall algorithm of integer factorization for integer is as follows.

The following step 5 is quantum step, and others are classical parts.

- If is even, return 2 as a factor.

If not, proceed to the next step. - Check if for some integer and . (Check for all . This is n-polynomial computation for n-bit’s integer N.)

If so, return as a factor. If not, proceed to the next step. - Pick a random integer and get the great common divisor with Euclidean algorithm.
- If gcd(N,a) > 1, return gcd(N,a) as a factor.

If not, proceed to the next step. - If gcd(N,a) = 1, and is coprime and then run above quantum period-finding algorithm for and . (i.e, find the least integer such as .)
- If the period is odd, go back to step 3.
- If the period is even and , go back to step 3.
- Otherwise, check if is a nontrivial factor of . 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

Categories: Uncategorized

## 2 replies »