SA Shor's Algorithm Simulation

Theory

Table of Contents

Intro

A semiprime $N$ is the product of two primes. For an integer $a < N$ coprime to $N$, the function

f(x) = a^x mod N

is periodic. Its period $r$ satisfies:

a^r == 1 (mod N)

If $r$ is even, then:

a^r - 1 == 0 (mod N)
(a^(r/2) - 1)(a^(r/2) + 1) == 0 (mod N)

This means non-trivial factors of $N$ may be recovered with:

gcd(a^(r/2) - 1, N)
gcd(a^(r/2) + 1, N)

The quantum part of Shor's algorithm is a period-finding routine. The classical post-processing then turns a useful period into factors.

Classical Pre-Processing

Before simulating the quantum period-finding step, the code checks for easy classical cases:

If any of these checks finds factors, no quantum simulation is needed. Otherwise, the code proceeds with a coprime base $a$.

Quantum Period Finding

Registers

Let:

n = ceil(log2(N))
Q = 2^(2n)
M = 2^n

The simulator uses:

Using Q ~= N^2 gives enough resolution for continued-fraction period recovery.

Hadamard Gates

The first register starts in state $|0⟩$. Hadamard gates create an equal superposition over all first-register states:

|ψ1⟩ = (1/sqrt(Q)) sum_x |x⟩|0⟩

In mode="matrix", this is represented by an explicit Hadamard matrix on the first register tensor an identity matrix on the second register.

Modular Oracle

The oracle encodes the periodic function in the second register:

|x⟩|0⟩ -> |x⟩|a^x mod N⟩

The explicit matrix implementation uses a reversible XOR oracle:

|x⟩|y⟩ -> |x⟩|y xor (a^x mod N)⟩

This makes the oracle a permutation matrix and therefore unitary. Starting from $|y⟩ = |0⟩$, it produces the same encoded function values needed for period finding.

After the oracle, the state has the form:

|ψ2⟩ = (1/sqrt(Q)) sum_x |x⟩|a^x mod N⟩

Inverse Quantum Fourier Transform

The inverse quantum Fourier transform is applied to the first register. It concentrates probability near values c satisfying:

c / Q ~= s / r

where:

The resulting first-register probabilities are plotted by probability_plot.py.

Continued Fractions

find_period.py sorts measurement outcomes by probability. For each high-probability measured value c, it forms:

c / Q

and uses continued fractions to recover denominator candidates. Candidate denominators and their multiples are accepted only if they pass validation:

pow(a, r, N) == 1
gcd(a^(r/2) - 1, N) and gcd(a^(r/2) + 1, N) are non-trivial factors

This replaces the earlier peak-spacing heuristic. It also correctly rejects cases where a period exists but cannot factor $N$, such as $N=33, a=2$.

Classical Post-Processing

Once a useful period r is found, the code computes:

gcd(a^(r/2) - 1, N)
gcd(a^(r/2) + 1, N)

If both are non-trivial and multiply to $N$, factorization succeeds. Otherwise, Shor's algorithm must be retried with a different base $a$.

Simulation Modes

The repository supports two period-finding modes:

For N=15, a=2, both modes produce matching first-register probabilities up to floating-point error.

Example

For $N = 15$ and $a = 2$:

2^0 == 1  (mod 15)
2^1 == 2  (mod 15)
2^2 == 4  (mod 15)
2^3 == 8  (mod 15)
2^4 == 1  (mod 15)

The period is $r = 4$. Since $r$ is even:

2^(r/2) = 2^2 = 4
gcd(4 - 1, 15) = gcd(3, 15) = 3
gcd(4 + 1, 15) = gcd(5, 15) = 5

So the factors of $15$ are $3$ and $5$.

The visualization helpers can show:

For register-flow diagrams and Qiskit-generated circuit drawings, see CIRCUITS.md.

References

Primary Sources

Algorithm Background

Implementation References

Author

Sid Richards

License

MIT. See LICENSE.