Shors Algorithm

1. Shors Algorithm, RSA & Secrets
2. RSA
3. Quantum Fourier Transformation
4. Classical Factoring
5. Phase Kickback
6. Eigenvalue and Phase Estimation
7. Order and Period Finding
8. Shor’s Factoring Algorithm
9. Appendix A: Applying a Circuit Multiple Times
10. Appendix B: Geometric Series and Roots of Unity
11. Sources

1. Shors Algorithm, RSA & Secrets

‘Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer’ by Peter Shor, 1994, [1]

The described quantum algorithm can ‘efficiently’ factor large integers – something classical computer struggle with

Many of todays encryption systems, including RSA (Rivest–Shamir–Adleman), rely on the difficulty of integer factorization.

‘Harvest now, decrypt later’ (HNDL), see e.g. a whitepaper from Mastercard(!) from 2025: ‘Migration to post quantum cryptography – white paper’, [2]

2. RSA [3]

RSA is a public key cryptosystem based on factoring, using two distinct, but related keys

E.g. the 20^th Fermat number

is known to be composite [4], but non of its prime factors are known!

To generate both public and private keys for RSA cryptosystems, choose two very large prime numbers p and q:

i. n = pq

ii. Pick e, relatively prime to (p-1) (q-1),  i.e. gcd(e, (p-1)(q-1)) = 1)

iii. Then by Eulers theorem d = e^-1 (mod (p-1)(q-1)) exists, and d and n are relatively prime

iv. e and n build the public key, d the corresponding private key! (although p and q are not relevant anymore, they should never be revealed!)

To encrypt a message, break it into blocks and assign a number m < n.

Encryption: c = m^e (mod n) 

Decryption: m = c^d (mod n)

Find a detailed proof of the above in [3]

Example:

1. p = 3, q = 11, n = pq = 33
2. (p-1)(q-1) = 2 * 10 = 20 = Ф(33),  Ф is the Euler totient function: number of relative primes less then n. For    n = 33: 1, 2, 4, 5, 7, 8, 10, 13, 14, 16, 17, 19, 20, 23, 25, 26, 28, 29, 31, 32
3. Eulers Theorem:

Hence: choose 1 < e < Ф(33) and e and Ф(33) to be relatively prime: e.g. with e = 7 is d = e^-1 = 7¹⁹ mod 20 = 3 

Hence: public key e = 7 and n = 33; private key d = 3

Assume one coded text block is m = 2.

Encryption: c = 2⁷ mod 33 = 128 mod 33 = 29

Decryption: m = 29³ mod 33 = 24389 mod 33 = 2

3. Quantum Fourier Transformation [5], [6]

I hope you already have an idea of how useful Fourier Transformations are, because otherwise it could be difficult for now to explain why the following is really, really useful … 😊 (maybe because of waves in wave functions and the importance of phases …)

(please see ‚the getting started‘ blog for primitive roots)

Example:  N = 2, ω = -1

Hence:

A QFT is a representation of an input as a superposition of frequences (by a synchronized rotation of qubits), e.g.:

The QFT quantum circuit

Recall:

and

(which is identical to

)

Define  the rotation gate

and hence

With this, QFT4 is

4. Classical Factoring

Shors’s algorithm can factor large integers almost exponentially faster than classical methods! Largest prime known (Jan 26):

(GIMPS: Great Internet Mersenne Prime Search [7])

(Really) Classical Method for finding primes: Sieve of Eratosthenes

Method for finding factors (Fermat’s method)

n = p*q, p and q both primes. Then: p + q and p-q are even, hence p + q = 2r, p-q = 2k and therefore

p = r + k, q = r-k

  n = p q = r² –k²

  n + k² = r²

5. Phase Kickback

Hence:

To summarize:

1) Bring the control qubit q₀ in a superposition, q₁ into |1>, which is an eigenstate of

2) Then:

applied to q₁ gives a phase to q₀ (kickback!)

3) On measurement, q₁ is always in |1> and q₀ half of the time in |0> and in |1> otherwise. The state can be represented by a Tensorproduct, hence no entanglement!

6. Eigenvalue and Phase Estimation

Measuring a register of entangled qubits in superposition will give one single result with only some probability that it is the correct orthonormal basis ket (see ‘Inversion about the mean’ in the Getting Started blog)

Idea: transform the register of entangled qubits to another state, that singles-out the correct result (with the constraint that the sum of all probabilities equals one)

U be a quantum gate on n qubits with a 2ⁿx2ⁿ matrix representation. |ψ> be an eigenvector of U, i.e.

Controlled U^k-gate:                                                                  

(n: |ψ> is a n-qubit quantum register, U^k is only applied if the above control bit is 1)

Recap: General phase kickback circuit:

With that, consider the following circuit:

Note the phase kickback: |ψ> is not affected, but the phase information of the eigenvalues (of the eigenvectors forming |ψ> by superposition) are propagated to the two qubits by U!

The input to the inverse QFT: (essentially a superposition of q₀ and q₁ with a phase shift)

with

(primitive 4^th root of unity), and because U is unitary:     

Still ‘only’ a phase shift, the physical (measurable) state represented by the orthonormal basis is still unchanged!

Generalize to

with QFT_m^-1, at the end of the circuit:  

Objective is to estimate φ by a measurement, i.e. the eigenvalue of one specific eigenvector |c>

Now:

φ in [0,1] (QFT is a rotation!) unknown, but there is a c in Z within

Why 1/2? Applied to the phase kickback, the QFT amplitude peaks at

exactly, which is in general not an integer. But the measurement gives an integer as some eigenvalue, hence: because j is an integer it holds as for every integer:

Hence define

with ⌊x⌋ the largest integer ≤ x and find c in a way that

is a good approximation to φ. Define in addition

|d| measuring the difference between φ and its rational approximation

The measurement of the {q₀, …, q_m-1}-register should be the eigenvalue of |c> as an approximation to

From d:

hence

Measurement of |c>: 

with j = c, hence 

Assume d=0, i.e. the approximation is accurate and hence

is a rational number (showing for d=0 the circuit before the QFT^-1 creates a QFT expression!):

and it holds the measurement of the upper qubit register yields with 100% probability |c>!

d ≠ 0 (which will be the case most of the time)?  As with Simon’s algorithm, run the circuit multiple times to reduce the probability of not getting the correct answer below a certain threshold (see Appendix A):

^{is a geometric series}

With

^{(the last step from the double angle formula), hence}

^{Sinus is a concave function}

hence

^{and therefore with}

^{and by symmetry it holds also for}

^{it follows with}

^{and in addition}

^{which is the probability to get it right. Running the circuit 28 times, p of never getting the right answer is less than}

7. Order and Period Finding

When start functions (if …) to repeat themselves? Define f_a(x) = a^x mod M, a coprime to M (i.e. gcd = 1): the smallest r for which f_a(x) = f_a(x+r) is the period of f_a

Another way of looking to this: What is the smallest r in N, such that a^r = 1 mod M? If it is existing, r is the order of a mod M. Since a^x+r = a^x mod M, r is also the period of f_a, i.e. period finding equals order finding

Why is r of interest? Assume r ≥ 0:  a^r = 1 mod M, hence a^r-1 = 0 mod M, which might give a factorization of M!

How? Let’s assume we find r via a quantum algorithm. Then with

M divides

but most probably neither of the factors individually, hence 

are likely non-trivial factors, factoring M! Using Euclid’s algorithm with

calculate gcd(x-1,M) with x-1 and M positive integers, hence there exist integers q and r, so that x-1 = Mq + r with 0 ≤ r < M; r is called the remainder, q the quotient. If n divides x-1 and M, then n divides (x-1) – Mq = r. Especially if n = gcd(a,b):  gcd(a,b) = gcd(b,r) = gcd(b, a mod b), i.e. some smaller numbers

Repeat this until r ≥ 0 eventually stops the process, and the last r > 0 is the gcd. Analogously for x+1

Example:

M = 15, a = 7, then for 7 mod 15: r = 4 (7⁴ = 2401 = 160*15+1), hence x = 7² = 49, and so x-1 = 48, x+1 = 50 Now Euclid: gcd(48, 15) = 3, gcd(50, 15) = 5! So, M factors into 3*5 (surprise 😊)

Now, how to find r (with the help of phase estimation …) in order to find a factorization of M?

A register of l qubits can represent 2^l states. To represent an integer M, one needs enough states to cover the range [0,M-1], therefore 2^l ≥ M. Define

(the ceiling, so l is an integer) as the minimum number of qubits needed to represent M. For y in {0,1}^l define

with a coprime to M (i.e gcd(a,M) = 1). Hence there are b, c in Z (Bézout’s identity or lemma) giving 1 = ab + Mc = ab mod M  (Mc is an exact multiple of M!). Apply U multiple times

U²|y> = UU |y> = U |ay mod M> = |a²y mod M>

This can be generalized to Uⁿ

Now: introduce repeated squaring

x⁹ = xxxxxxxxx, i.e. 9 multiplications. But it can be reduced to 4: x⁹ = x((x²)²)², with 9 = 1001₂, the right 1 being the least significant bit (since 2⁰=1, but it determines whether the number is even or odd). This procedure is called ‘repeated squaring’.

In general: if n bits represent an exponent b, one can compute x^b with 2n multiplications at most. Let’s see the procedure:

Compute x¹³, 13 = 1101₂ , so n = 4

Hence 3 squarings and 3 multiplications = 6 multiplications < O(n³)

With the binary form

one gets

To calculate a^z mod M, the values 

 can be built by successive squaring: a¹, a² = (a)², a⁴ = (a²)², …

With that define

This is important for Shor’s algorithm: repeated squaring is of O(k), multiplying two k-bit numbers is of O(k²), hence the |z>|y> expression scales O(k³), e.g.

K = 4: a) exponent has 4 bits ~ 4 squaring; b) each multiplication: 4*4 = 16; c) sum up: 4*16 = 64 = 4³ operations

This translates to the number of gates needed for realization, i.e. Shor’s algorithm quantum circuit scales (polynomial) with O(l³) gates!

Now let’s proceed with U, with r the order of a mod M: 

are eigenvectors of U for 0 ≤ j ≤ r. Let’s see:

The last step holds, because a^r = a⁰ mod M, i.e. the sum is the same. So, each U applied to the target state kickbacks a phase to the corresponding control qubit, adding to the binary representation |z> |a^z y mod M> above. As we will see, applying one inverse QFT maps the superposition of the control qubits with order r to multiples of 

The Circuit

Use the phase estimation to determine the eigenvalues of |w_j> and derive r. How to create |w_j>? It can’t be created directly! Create a superposition instead:

1. Apply H to the control qubits, generating a uniform superposition
2. Regarding the target qubit: see below
3. Each controlled U-gate multiplies the target register with  mod M. Note: because of the superposition in the control state, |x>|a^x mod M> is computed in one step!
4. Inverse QFT: converts the periodic structure of the state into phases, their amplitudes peaking around multiples of  
5. Measurement yields ~ 

The initial target register’s state is a superposition of the eigenvectors of U:

(in line 3, the summation formula is used, see the ‘Getting Started’ blog)

That’s the |1> of the target qubit in the circuit!

So far so good, but still some way to go for r 😊

Now let’s analyze the circuit.

After all the U gates the joint state is

(note that |1> is also n (= ⌈ log₂M ⌉) qubits representing M)

The target register |a^x mod M> is periodic in r, i.e.

x = lr + s, s= 0, 1, …, r-1 (repeating |a⁰>, |a¹>, …, |a^r-1>), and

[ to see l: for s, s+r, s+2r, …  x and s share the same remainder: x = lr + s for l = 0, 1, 2, … ; the largest x is 2²ⁿ-1:  

since 0 ≤ s ≤ r-1 extended to 1 ≤ 1+s ≤ r and hence    is between    and 1. Therefore the error is at most 1, which can be ignored against    for large 2²ⁿ)      ]

Hence, we have  

Apply

to the control register:

where the term in brackets is a geometric series: 

when

is in Z and 0 otherwise (see Appendix B). Hence it is non-zero for

With that the inverse QFT expression yields

Putting it together for the full state after the inverse QFT:

And with

i.e.

one gets

Finally ignoring the target register (not measured):

Hence the relevant states are those for

and therefore the measurement result m gives for some j

and hence

is a rational number close to the result, recall the

from the beginning!

Why ~ ? r does not necessarily divide

hence the inverse QFT peaks not at the exact value, but at the nearest integer:

So, closer, but still not r 😊 How to determine r from the approximation to j/r? The next bit needed: continuous fraction

Continuous Fraction [8]

Express any real number by integers alone. For a_i integers:

a_i is called denominators, the 1’s are called numerators. Example of a simple (because all numerators are 1s) continued fraction:

‘Best approximation property’: The nth convergent is the ratio of two succession terms in a sequence of numerators and denominators, defined recursively for n ≥ 1:

P_-1 = 1, P₀ = a₀,  P_n = a_nP_n-1 + P_n-2

Q_-1 = 0, Q₀ = 1, Q_n = a_nQ_n-1 + Q_n-2

The sequence of convergents

provides the best rational approximations to the real number in question! In general: the procedure only terminates for rational numbers. The convergents oscillate around x:

r a real number be the value of the continued fraction [x₀; x₁, x₂, …. ] It holds:

a) Each convergent is a reduced fraction, i.e. gcd(P_n, Q_n) = 1

b) If k > j, then Q_k > Q_j

c) The denominators are increasing exponentially:

d) r can be approximated as closely as needed:

f) If

then 

must appear as a convergent in the continued fraction expansion of x!

Number of Qubits & Accuracy of the Phase Measurement

Recall: The work register has n = ⌈log₂M⌉ qubits and hence for any measured phase φ_measured (remember from the beginning: 

 as the rational approximation to the real phase φ): 

hence:

n qubits represent 

 equally spread values on [0,1):

The true phase φ does in general not fall exactly on one of those values, but the inverse QFT measurement collapses the state to the nearest value to φ in the list! The width between two points is

so the worst case is φ exactly between two points and hence 

With that, let the control register have m = 2n+1 qubits. Then a measurement of the phase φ_m will necessarily be accurate as

and hence

and therefore 

must appear as a convergent in the continued fraction expansion of 

and the exact phase will be of the form 

There are more conditions for Shor’s algorithm to work, see e.g. [9] :

i) r needs to be even, because of

ii) It has to hold

in order that

is not a trivial root of unity:

then is

and N divides

and therefore the result is trivial:

iii) And of course,

8. Shor’s Factoring Algorithm

Now putting everything together:

1. Choose a random number a in {2, 3, …, M-1}. If gcd(a, M) = 1 go to step 2. (because otherwise this is already a non-trivial factor … )

2. With the quantum phase circuit measure the phase of the modular exponentiation operator U(a, M): U|y> = |a y mod M>_l  with

 and for an eigenvector

with j in {0, 1, 2, …, r-1}. Map the superposition by an inverse QFT to multiples of 

3. Use continued fraction to derive the exact period r from the measured approximation 

 of the phase by analyzing all the convergents such that

which is achieved by the number of control qubits being m = 2n+1 (there must be even more qubits related to the probability calculated for d not zero, but let’s skip this here, since it doesn’t really add to the schema). Check the convergents from the smallest to the largest: if r is odd, return to step 1; also if

If

the algorithm found a solution for r and one can proceed with

4. Factors of N are now given by

and the gcd can be derived using Euclid’s algorithm. If both gcds are trivial, start again with step 1 with a different a. If either gcd is not 1, the algorithm found a factor and is done!

In 2001, the factoring of 15 was realized with Shor’s algorithm experimentally, [10]

9. Appendix A: Applying a Circuit Multiple Times

Let p = 0.6 be the probability to get the wrong answer to some question by running a circuit once. Then 1-p = 0.4 is the probability of getting the right one

Probability to never get the correct answer to be less then 10^-6 :

pⁿ < 10^-6  := ε

n log p ≤ log ε

0 < p < 1:  log p < 0, hence  n > log⁡(ϵ) / log(⁡p)

Hence n > 27

10. Appendix B: Geometric Series and Roots of Unity

Be S equal to

Geometric series formula:

Now:

rk is an integer,

and hence z^N – 1 = 0. Two cases for the denominator:

i. rk/N is not an integer: then the numerator z-1 ≠ 0, so S = 0

ii. rk/N is an integer:

and the geometric series formula does not apply. But in this case:

In the context of Quantum Fourier Transformations: after the QFT (mapping periodicity to frequency), only amplitudes for k = 0, N/R. 2N/R … will be different from 0! In other words, for k = jN/R,  with j = 0, 1, …

11. Sources

[1] ‚Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer‘ by Shor, 1995 (www.arxiv.org/abs/quant-ph/9508027)

[2] Migration-to-post-quantum-cryptography-WhitePaper by Mastercard, 2025

[3] ‚Introduction To Cryptography With Open-Source Software‘ by McAndrew, CRC Press 2011

[4] The Twentieth Fermat Number is Composite by Young et al., 1988

[5] ‚Quantum Computation and Quantum Information‘ by Nielsen et al., Camebridge University Press, 2010

[6] ‚Dancing with Qubits‘ by Sutor, <packt>, second edition, 2024

[7] Great Internet Mersenne Prime Search

[8] ‚Simple Continued Fractions an Approach for High School Students‘, by Paraskevopoulos, 2023

[9] ‚Shor’s Factoring Algorithm and Modular Exponentiation Operators‚ by Singleton, 2023

[10] ‚Experimental realization of Shor’s quantum factoring algorithm using nuclear magnetic resonance‚ by Vandersypen et al., 2001