Shor’s Algorithm

Shor’s algorithm, introduced by Peter Shor in 1994 (algorithms for quantum computation: discrete logarithms and factoring), is the breakthrough quantum algorithm that factors large integers in polynomial time, a task widely believed to be intractable for classical computers. Given an integer N, the algorithm finds its prime factors in O\left((\log N)^2 (\log \log N)(\log \log \log N)\right) time on a quantum computer. This is exponentially faster than the best classical algorithms, such as the general number field sieve, which runs in sub-exponential time O\left(e^{c (\log N)^{1/3} (\log \log N)^{2/3}}\right), where c = \frac{1}{3}(92 + 26\sqrt{13})^{1/3}.

Peter Shor, Photo from this link.

The discovery of Shor’s algorithm showed that a large quantum computer would be capable of breaking widely used public-key cryptosystems, most notably RSA, whose security relies on the assumed hardness of integer factorization. This result initiated both the global effort to build scalable quantum computers and the rapid development of post-quantum cryptography, which seeks cryptographic protocols secure even against quantum attackers.

RSA Cryptosystem

Since Shor’s algorithm was primarily designed to break the RSA public-key cryptosystem. For completeness, we first review the core principles of RSA. The RSA cryptosystem, along with its variants, is one of the most widely deployed methods for securely transmitting messages over public channels, such as the Internet, where it protects confidentiality and enables digital signatures. Its security relies fundamentally on a single, long-standing computational assumption: that factoring a large composite integer into its prime factors is extraordinarily difficult for classical computers.

It is worth emphasizing that testing whether a number is prime is computationally feasible, whereas factoring a large composite number into its prime factors remains extraordinarily difficult. RSA cryptography leverages precisely this asymmetry: while anyone can efficiently verify the primality of the factors used in key generation, recovering those factors from the public modulus is, for classical computers, practically intractable. This one-way computational difficulty forms the cornerstone of RSA’s security, enabling secure encryption and decryption of messages.

The RSA cryptosystem operates as follows:

  1. Alice prepares two very large prime numbers p and q, which she keeps secret. She computes N = pq and publishes N. She also chooses a number, called the public exponent (or encoding exponent), e < N, which is required to be relatively prime to (p-1)(q-1), that is, \mathbf{gcd}(e,(p-1)(q-1)) = 1. Alice then computes the multiplicative inverse of e modulo (p-1)(q-1), namely a number d satisfying

de \equiv 1 \pmod{(p-1)(q-1)}

She keeps d secret, thus called the secret exponent (or decoding exponent).

  1. When Bob wants to send a message to Alice, he first encodes the message as an integer m with m < N. He encrypts it as

M = m^{e} \pmod{N}

  1. Upon receiving M, Alice decrypts it using

m \equiv M^{d} \pmod{N}

where d is the multiplicative inverse of the public exponent e modulo (p-1)(q-1).

To see that the RSA cryptosystem correctly recovers the message m using m \equiv M^{d} \pmod{N}, we need to check

m \equiv (m^e)^{d} \pmod{N}.

From de \equiv 1 \pmod{(p-1)(q-1)}, there exists an integer k such that de = k (p-1)(q-1) + 1. Hence

(m^e)^{d} = m \bigl[ m^{k(q-1)} \bigr]^{p-1}.

Suppose m is not a multiple of p. Fermat’s little theorem gives

\bigl[ m^{k(q-1)} \bigr]^{p-1} \equiv 1 \pmod{p}.

If m is a multiple of p, then m^{de} \equiv 0 \pmod{p}.

Similarly,

\bigl[ m^{k(p-1)} \bigr]^{q-1} \equiv 1 \pmod{q}

when m is not a multiple of q, and m^{de} \equiv 0 \pmod{q} if m is a multiple of q.

We now verify m^{de} \equiv m \pmod{N} by considering several cases.

  1. If m is a multiple of both p and q, then it is automatically a multiple of N = pq, and m^{de} \equiv m \pmod{N} holds trivially.
  2. If m is a multiple of p but not of q, then

m \equiv 0 \pmod{p}, \quad m^{de} \equiv 0 \pmod{p},

and

\bigl[m^{k(p-1)}\bigr]^{q-1} \equiv 1 \pmod{q}.

Hence,

m^{de} \equiv m \pmod{q}, \quad m^{de} \equiv m \pmod{p}.

By the Chinese remainder theorem{}^1,

m^{de} \equiv m \pmod{N}.

We are using a basic consequence of the Chinese remainder theorem: if two integers a and b satisfy a \equiv b \pmod{p} and a \equiv b \pmod{q}, with \gcd(p,q) = 1, then a \equiv b \pmod{pq}. In particular, if p and q are prime, it is easier to see p \mid (a-b) and q \mid (a-b) implies pq \mid (a-b).

  1. If m is a multiple of q but not of p, the argument is symmetric, giving

m^{de} \equiv m \pmod{N}.

  1. If m is not a multiple of either p or q, then Fermat’s little theorem gives

m^{de} \equiv m \pmod{p}, \quad m^{de} \equiv m \pmod{q}.

Again, by the Chinese remainder theorem,

m^{de} \equiv m \pmod{N}.

This completes the proof of m^{de} \equiv m \pmod{N}.

We see from the above analysis that if one could factor N efficiently, then the RSA system would be completely compromised. Indeed, factoring N reveals the primes p and q, from which one can compute (p-1)(q-1), determine the public exponent e, and then recover the secret exponent d. With d in hand, any encrypted message can be decrypted.

Shor’s Algorithm

Definition (The Factoring Problem)
The input to Shor’s algorithm is an odd composite integer N = p q, where p and q are distinct large primes. The goal is to find the prime factors p and q.

Shor’s algorithm is a hybrid quantum-classical algorithm that factors N in polynomial time on a quantum computer. It consists of two main phases:

  • Classical reduction: Transform factoring N into an instance of the order-finding problem (see this post).
  • Quantum order-finding: Solve the order-finding problem efficiently using a quantum computer.

The classical reduction (which succeeds with probability > 1/2 for random choices) is as follows:

1. Choose a random integer a \in {2, 3, \dots, N-1} and compute d = \gcd(a, N):

  • If 1 < d < N, then d is a non-trivial factor of N. Factoring complete.
  • If d = 1, continue.

2. Find the multiplicative order r of a modulo N, i.e., the smallest positive integer r such that
a^r \equiv 1 \pmod{N}.

3. If r is odd, return to the first step with a new a.

4. If r is even, compute b = a^{r/2} \mod N:

  • If b \equiv -1 \pmod{N}, return to the first step.
  • If b \not\equiv \pm 1 \pmod{N}, proceed.

5. Compute the factors:
p = \gcd(b - 1, N), q = \gcd(b + 1, N).
At least one of these is a non-trivial factor of N.

The only step that is computationally infeasible on a classical computer is step 2: finding the order r.
Shor’s breakthrough contribution is a quantum algorithm that performs this order-finding step in polynomial time using quantum phase estimation on the unitary operator x \mapsto a x \mod N.

Once r is obtained quantumly, the remaining classical post-processing immediately yields the prime factors p and q. Thus, on a sufficiently large fault-tolerant quantum computer, Shor’s algorithm solves the integer factorization problem efficiently, rendering RSA and related cryptosystems insecure.

Remark (Why the Reduction Works)

Assume that r is even and that a^{r/2} \not\equiv -1 \pmod{N} (the cases explicitly filtered out by the algorithm). Since r is the order of a modulo N, we have
a^r \equiv 1 \pmod{N}.
This implies
a^r - 1 = (a^{r/2} - 1)(a^{r/2} + 1) \equiv 0 \pmod{N},
so N \mid (a^{r/2} - 1)(a^{r/2} + 1), or equivalently,
pq \mid (a^{r/2} - 1)(a^{r/2} + 1).

Because p and q are distinct primes, the prime factors of the product (a^{r/2}-1)(a^{r/2}+1) can be distributed in only four logical ways:

  1. p divides a^{r/2}-1 and q divides a^{r/2}+1
  2. p divides a^{r/2}+1 and q divides a^{r/2}-1
  3. pq divides a^{r/2}-1 (i.e., N divides a^{r/2}-1)
  4. pq divides a^{r/2}+1 (i.e., N divides a^{r/2}+1)

Cases 1 and 2 are precisely what we want: they yield non-trivial factors via the gcd computations
p = \gcd(a^{r/2}-1, N) and q = \gcd(a^{r/2}+1, N) (or vice versa).
Cases 3 and 4, however, must be excluded:

  • Case 3 is impossible: If N divides a^{r/2}-1, then a^{r/2} \equiv 1 \pmod{N}.
    But r is the order of a modulo N, so r would have to divide r/2, which is impossible unless r = 0, a contradiction.
  • Case 4 is impossible under our assumption: If N divides a^{r/2}+1, then a^{r/2} \equiv -1 \pmod{N}.
    This situation is explicitly detected and rejected in the algorithm (step 5).

Thus, whenever the algorithm reaches step 6, we are guaranteed to be in Case 1 or Case 2, and the two gcd operations recover the prime factors p and q of N.

Example

Let us illustrate Shor’s algorithm with the integer N = 15.

  1. Choose a random number: Let a = 7. Compute \gcd(7, 15) = 1.
    Since \gcd(7, 15) = 1, we proceed to the next step.
  2. Find the period r: Using a quantum computer, we determine the period of the function
    f(x) = 7^x \bmod 15.
    The result is r = 4, since 7^4 \equiv 1 \pmod{15}.
  3. Use the period to find the factors:
    • Since r = 4 is even, compute 7^{r/2} - 1 = 7^2 - 1 = 48.
    • Compute the gcd: \gcd(48, 15) = 3, which is one factor of 15.
    • Similarly, 7^{r/2} + 1 = 7^2 + 1 = 50, and \gcd(50, 15) = 5, which is the other factor.

Thus, we have successfully factored 15 = 3 \times 5.

Quantum Order-Finding Algorithm

The order-finding problem is a fundamental problem that plays a central role in several quantum algorithms; in particular, it underlies Shor’s celebrated factoring algorithm. Classically, order-finding is hard—no efficient classical algorithm is known—but quantum computers can solve it efficiently. Efficiently determining the order of an integer modulo N allows integers to be factored in polynomial time, which would compromise RSA encryption, as we will discuss in the subsequent sections.

The order-finding problem is defined rigorously as follows:

Definition (Order-Finding Problem) Let N \in \mathbf{Z}_+ be a positive integer, and let a \in \mathbf{Z}_N^{\times} be an integer coprime to N (that is, the greatest common divisor \gcd(a,N) of a and N equals one, \gcd(a,N) = 1). The notation \mathbf{Z}_N^\times denotes the set of invertible elements in \mathbf{Z}_N under multiplication—those integers that are coprime to N. The order of a modulo N, denoted \mathrm{ord}_N(a) or simply r, is the smallest positive integer r such that
a^r \equiv 1 (\mathrm{mod}\, N)
The order-finding problem is: given a and N, compute r = \mathrm{ord}_N(a).

At first glance, this definition may appear somewhat abstract.
Do not worry—shortly, we will unpack each of the mathematical ingredients and illustrate their meaning with concrete examples.

The quantum algorithm for order-finding exploits the periodic structure of the function

f_{a,N}(x) \equiv a^x (\mathrm{mod}\, N) \ \ \ \ \ (1)

which is periodic with period exactly r. By preparing a quantum superposition over all possible exponents and applying the quantum phase estimation algorithm to the unitary operator

U |y\rangle = |a y (\mathrm{mod}\, N)\rangle,

we obtain an efficient quantum procedure that succeeds with high probability.

Preliminaries of Modular Arithmetic

Before presenting the complete quantum algorithm, let us first review a few essential results from number theory. This section may appear somewhat technical compared to others; if you prefer not to delve into the details, simply remember that for coprime integers a and N, there exists a positive order r, and the function in Eq.(1) is periodic and turn to the next section.

An integer N > 1 is called a prime number if it has no nontrivial divisors; for example, N = 2, 3, 5, 7, 11, \dots are prime numbers. Otherwise, N is called a composite number. We use a \mid b to denote that b is divisible by a, for instance 3 \mid 12, 2 \mid 6. Similarly, a \nmid b means that a does not divide b, for example 3 \nmid 14.

Every integer N > 1 can be uniquely expressed as a product of prime powers:

N = p_1^{e_1} \times p_2^{e_2} \times \cdots \times p_m^{e_m},

where p_1, \ldots, p_m are distinct prime numbers and e_1, \ldots, e_m are positive integers. This result is known as the Fundamental Theorem of Arithmetic.

To obtain the prime factorization of a given number N, one typically follows these steps:

  1. Divide N by 2 repeatedly until the quotient is no longer divisible by 2.
    Let e_2 be the number of times we divide by 2, and denote the resulting quotient by N_2 = N / 2^{e_2}, for which 2 \nmid N_2.
  2. For the quotient N_2, proceed with the next prime, 3, and divide repeatedly until the quotient is no longer divisible by 3.
    Let e_3 be the number of divisions, and denote the resulting quotient by N_3 = N_2 / 3^{e_3}.
  3. Continue in this manner with successive primes 5, 7, 11, \dots until the quotient becomes 1.

Exercise. Decompose N = 84 into its prime factors: 84 = 2^2 \times 3 \times 7.

The key operation we will use is modular arithmetic.
In everyday life we count linearly: 1, 2, 3, …, but in modular arithmetic we count on a circle.
Take a clock that shows only 12 hours. When the hand reaches 12, it jumps back to 0.
Mathematically, we say time is measured modulo 12.
This simple idea—wrapping around after a fixed modulus N—is crucial in number theory and its applications in modern public-key cryptography.

Definition (Congruence) Fix an integer N \ge 2, called the modulus.
Two integers a and b are congruent modulo N if their difference is a multiple of N:

a \equiv b (\mathrm{mod}\, N) \quad \Longleftrightarrow \quad N \mid (a - b) \quad \Longleftrightarrow \quad \exists\, k \in \mathbf{Z} \text{ such that } a = b + kN.

We read a \equiv b (\mathrm{mod}\, N) as “a is congruent to b modulo N”.

The congruence class of a is usually represented by its smallest non-negative residue:

a (\mathrm{mod}\, N) \in \mathbf{Z}_N := \{0, 1, 2, \dots, N-1\}

It is also convenient to use [a] to denote the congruence class whenever there is no ambiguity about the modulus.

Example. Let N = 7. We have \mathbf{Z}_7 :=\{0, 1, 2, \dots, 6\}. Then 1 \equiv 8 (\mathrm{mod}\, 7) and 2 \equiv 9 (\mathrm{mod}\, 7). The congruence class [0] consists of \{0, 7, 14, 21, \dots\}, the congruence class [1] consists of \{1, 8, 15, 22, \dots\}, and [2] consists of \{2, 9, 16, 23, \dots\}, etc.

The order-finding problem considers the modular exponentiation for given integers a and N:

a^x \equiv b (\mathrm{mod}\, N)

The remainder b can take values in 0, 1, \dots, N-1:

a \equiv 0 (\mathrm{mod}\, N)
a \equiv 1 (\mathrm{mod}\, N)
\vdots
a \equiv N-1 (\mathrm{mod}\, N)

However, not all values of b = 0, 1, \dots, N-1 are necessarily attainable.

For example, let a = 2 and N = 6. We find that

2^0 \equiv 1 (\mathrm{mod}\, 6)
2^1 \equiv 2 (\mathrm{mod}\, 6)
2^2 \equiv 4 (\mathrm{mod}\, 6)
2^3 \equiv 2 (\mathrm{mod}\, 6)
2^4 \equiv 4 (\mathrm{mod}\, 6)
2^5 \equiv 2 (\mathrm{mod}\, 6)
\vdots

We observe a repeating pattern. In this case, the residues 0, 3, and 5 never appear.

Note also that a^0 \equiv 1 (\mathrm{mod}\, N) for any integers a and N. Hence, the order of a modulo N is defined as the smallest positive integer r such that a^r \equiv 1 (\mathrm{mod}\, N). However, we must still ensure that such an order actually exists. In fact, if \gcd(a, N) = d > 1, then no positive integer r satisfies a^r \equiv 1 (\mathrm{mod}\, N). This can be seen in the previous example with a = 2 and N = 6, where the sequence never returns to 1 modulo 6.

Exercise. Let N > 1 and \gcd(a, N) = d > 1. Then the only integer r satisfying

a^r \equiv 1 (\mathrm{mod}\, N)

is r = 0.

Proof. Suppose there exists an integer r such that a^r \equiv 1 (\mathrm{mod}\, N).
Then N divides a^r - 1, so in particular

d \mid (a^r - 1).

Since d \mid a, we also have d \mid a^r, and therefore

d \mid (a^r - 1) \implies d \mid 1,

which contradicts d > 1. Hence no such r \neq 0 exists.

For r = 0, the convention a^0 = 1 gives

a^0 = 1 \equiv 1 (\mathrm{mod}\, N)

trivially. Thus r = 0 is the only solution ∎.

In Shor’s algorithm, we either assume \gcd(a, N) = 1, or first compute the greatest common divisor classically, which already yields a nontrivial factor of N if d > 1.

Another useful concept we will use is Euler’s Totient function: \varphi(N) counts the number of integers k \in {1, 2, \dots, N-1} coprime to N.
For example, \varphi(1) = 1, \varphi(2) = 1, \varphi(3) = 2, \varphi(4) = 2, etc.

For fixed N, consider \mathbf{Z}_N.
We call a \in \mathbf{Z}_N invertible if there exists an integer c such that a c \equiv 1 (\mathrm{mod}\, N).

Proposition (Bézout’s Identity). Let a and b be integers with greatest common divisor d = \gcd(a,b). Then there exist integers x and y such that

a x + b y = d.

Moreover, the integers of the form a z + b t are exactly the multiples of d.

The set \mathbf{Z}_N is not a group under multiplication modulo N in general. We denote by (\mathbf{Z}_N^{\times}, \times) the set of invertible elements, which indeed forms a group. By Bézout’s identity, \mathbf{Z}_N^{\times} consists of all integers coprime to N, and this group has order \phi(N), where \phi denotes Euler’s totient function.

From a group-theoretic perspective, for any a \in \mathbf{Z}_N^{\times}, there exists a positive integer r such that

a^r \equiv 1 (\mathrm{mod}\, N).

Since a must be invertible, it follows that \gcd(a, N) = 1. In the multiplicative group \mathbf{Z}_N^{\times}, the order of any element divides the order of the group, that is,

r \mid |\mathbf{Z}_N^{\times}| = \phi(N).

Theorem (Euler’s Theorem). If \gcd(a, N) = 1, then

a^{\phi(N)} \equiv 1 (\mathrm{mod}\, N)

and hence the order of a modulo N divides \phi(N).

By Euler’s theorem, for integers a and N with \gcd(a, N) = 1, there always exists a positive integer r such that

a^r \equiv 1 (\mathrm{mod}\, N).

The order r = \mathrm{ord}_N(a) always satisfies 1 \le r \le \phi(N).

Since in the cyclic subgroup generated by [a] \in \mathbf{Z}_N^{\times} we have

{[a^0], [a^1], [a^2], \dots, [a^{r-1}]}, \quad \text{with } [a^r] = [a^0] = [1],

the sequence a^n (\mathrm{mod}\, N) is periodic with period equal to the order r:

a^0 \equiv 1 (\mathrm{mod}\, N)
a^1 \equiv a (\mathrm{mod}\, N)
\dots
a^{r-1} \not\equiv 1 (\mathrm{mod}\, N)
a^r \equiv 1 (\mathrm{mod}\, N)
a^{r+1} \equiv a (\mathrm{mod}\, N)
a^{r+2} \equiv a^2 (\mathrm{mod}\, N)
a^{r+3} \equiv a^3 (\mathrm{mod}\, N)

Thus, the sequence repeats with period r.

Example. Set a = 2 and N = 7. We have f_{a,N}(x) = 2^x (\mathrm{mod}\, 7) as:

2^0 \equiv 1 (\mathrm{mod}\, 7)
2^1 \equiv 2 (\mathrm{mod}\, 7)
2^2 \equiv 4 (\mathrm{mod}\, 7)
2^3 \equiv 8 \equiv 1 (\mathrm{mod}\, 7)
2^4 \equiv 16 \equiv 2 (\mathrm{mod}\, 7)
2^5 \equiv 32 \equiv 4 (\mathrm{mod}\, 7)
2^6 \equiv 64 \equiv 1 (\mathrm{mod}\, 7)
2^7 \equiv 128 \equiv 2 (\mathrm{mod}\, 7)
2^8 \equiv 256 \equiv 4 (\mathrm{mod}\, 7)

We see that the order is r = \mathrm{ord}_7(2) = 3.

Quantum Order-Finding Algorithm

The quantum order-finding algorithm can be divided into two main steps:

  1. Encode the order into the phase of an eigenvalue of a suitable unitary operator, thereby reducing the order-finding problem to a quantum phase estimation problem.
  2. Recover the order r from the estimated phase \varphi using the method of continued fractions.

Let us first examine how to encode the order into the phase of a unitary operator.

Fix integers a and N such that \gcd(a, N) = 1.
Define the unitary operator U_{a,N} acting on computational basis states |y\rangle by

U_{a,N} |y\rangle = | a y (\mathrm{mod}\, N) \rangle, \quad y = 0, 1, \dots, N-1

We will see that the order can be mapped to the phase of an eigenvalue of U_{a,N}.

Proposition. Let N \ge 2 and let a be an integer coprime to N.
The operator U_{a,N} on \mathbf{C}^N defined by

U_{a,N} |y\rangle = | a y (\mathrm{mod}\, N) \rangle, \quad y = 0, 1, \dots, N-1

is unitary.

Proof. Since \gcd(a, N) = 1, multiplication by a modulo N defines a bijection on the set {0, 1, \dots, N-1}.
Thus U_{a,N} merely permutes the computational basis { |0\rangle, |1\rangle, \dots, |N-1\rangle }. Any permutation of an orthonormal basis extends to a unitary operator, so U_{a,N} is unitary.

Alternatively, suppose U_{a,N} |y\rangle = U_{a,N} |y'\rangle. Then

a y \equiv a y' (\mathrm{mod}\, N) \implies a (y - y') \equiv 0 (\mathrm{mod}\, N) \implies N \mid a (y - y').

Because \gcd(a, N) = 1, we may cancel a (more precisely, multiply both sides by the modular inverse of a modulo N) to obtain

N \mid (y - y').

Since 0 \le y, y' < N, the only possibility is y - y' = 0, so y = y'.
Thus U_{a,N} is injective on a basis of \mathbf{C}^N.
Being an isometry that maps a basis to a basis, it is unitary ∎.

Proposition. Let r be the order of a modulo N, i.e., the smallest positive integer such that

a^r \equiv 1 (\mathrm{mod}\, N).

For each s = 0, 1, \dots, r-1, define the state

|u_s\rangle := \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{-2\pi i s k / r} \, | a^k (\mathrm{mod}\, N) \rangle.

Then, for each s = 0, 1, \dots, r-1, |u_s\rangle is an eigenstate of U_{a,N} with eigenvalue

e^{2\pi i s / r}.

Proof. Apply U_{a,N} to |u_s\rangle:

U_{a,N} |u_s\rangle = \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{-2\pi i s k / r} \, U_{a,N} | a^k (\mathrm{mod}\, N) \rangle

\quad\; = \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{-2\pi i s k / r} \, | a^{k+1} (\mathrm{mod}\, N) \rangle

Reindex the sum by letting k \mapsto k+1:

\quad\; = \frac{1}{\sqrt{r}} \sum_{k=1}^{r} e^{-2\pi i s (k-1) / r} \, | a^k (\mathrm{mod}\, N) \rangle

\quad\; = \frac{1}{\sqrt{r}} \sum_{k=1}^{r} e^{-2\pi i s k / r} \, e^{2\pi i s / r} \, | a^k (\mathrm{mod}\, N) \rangle

\quad\; = e^{2\pi i s / r} \frac{1}{\sqrt{r}} \sum_{k=1}^{r} e^{-2\pi i s k / r} \, | a^k (\mathrm{mod}\, N) \rangle

The sum now runs from k = 1 to r, but the k = r term is

e^{-2\pi i s r / r} \, | a^r (\mathrm{mod}\, N) \rangle = e^{-2\pi i s} \, | 1 (\mathrm{mod}\, N) \rangle = | 1 (\mathrm{mod}\, N) \rangle,

while the k = 0 term in the original definition of |u_s\rangle is

e^0 \, | a^0 (\mathrm{mod}\, N) \rangle = | 1 (\mathrm{mod}\, N) \rangle.

Thus the missing k = 0 term is exactly replaced by the k = r term, and we obtain

U_{a,N} |u_s\rangle = e^{2\pi i s / r} \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{-2\pi i s k / r} \, | a^k (\mathrm{mod}\, N) \rangle = e^{2\pi i s / r} |u_s\rangle.

This holds for every s = 0, 1, \dots, r-1 ∎.

Now the procedure becomes “straightforward”. We apply quantum phase estimation (see this post) to the operators U_a^{2^j} with the eigenstate |u_s\rangle, which yields the phases

\theta_s = \frac{s}{r}, \quad s = 0, 1, \dots, r-1

The output of the phase estimation algorithm is thus

0.\theta_0 \theta_1 \dots \theta_{n}

represented as binary numbers. A classical algorithm called continued fractions can be applied to obtain the order r.