shor

name

shor

description

Use this skill when the user asks for Shor integer factorization, quantum order-finding, period estimation, or implementing/running/debugging ShorAlgorithm in this repository (especially matrix/operator methods, IQFT-based phase post-processing, and continued-fraction period recovery). Relevant keywords include shor, factor N, order finding, period finding, modular exponentiation, continued fraction, quantum factoring, and ShorAlgorithm.

Shor's Algorithm

Purpose

Shor's algorithm factors an integer $N$ in polynomial time, e.g. $O((\log N)^3)$ with textbook arithmetic; the exact exponent depends on the modular multiplication/arithmetic implementation and the chosen cost model. This improves over the best known classical sub-exponential factoring algorithms by reducing integer factorization to quantum period-finding via the Quantum Fourier Transform.

Use this skill when you need to:

• Factor a composite integer using a quantum simulator.
• Understand the quantum period-finding subroutine.

Overview

1. Classical pre-checks: Even numbers and perfect powers are handled classically.
2. Random base selection: Choose random $a$ with $1 < a < N$ and $\gcd(a, N) = 1$.
3. Quantum period-finding: Build a circuit that implements $f(x) = a^x \bmod N$ and applies the inverse QFT to find the period $r$ of $f$.
4. Classical post-processing: Use the continued-fraction algorithm to extract $r$ from the QPE measurement, then compute $\gcd(a^{r/2} \pm 1, N)$ to find factors.
5. Retry loop: Automatically retries with a new $a$ if the current attempt fails (up to max_retries times).

Two circuit methods are supported:

• 'matrix': Uses a matrix-based modular exponentiation circuit.
• 'operator': Uses a modular multiplication operator circuit.

Prerequisites

• Quantum Fourier Transform (QFT) and inverse QFT.
• Modular arithmetic and Euler's theorem.
• Continued fractions algorithm.
• Python: numpy, Circuit, unitarylab.library.QFT, unitarylab.library.IQFT.

Using the Provided Implementation

from unitarylab_algorithms import ShorAlgorithm

algo = ShorAlgorithm()
result = algo.run(
    N=15,              # Number to factor
    method='matrix',   # 'matrix' or 'operator'
    backend='torch',
    max_retries=15
)

print(result['status'])           # 'ok' on success, 'failed' if exhausted
print(result['factors'])          # List of factors, e.g. [3, 5]; None on failure
print(result['period'])           # Found period r, or None for classical path
print(result['Selected base'])    # Random base a used
print(result['circuit_path'])     # Path to SVG circuit diagram (None for classical path)
print(result['plot'])             # List of saved output files: [{'format': 'svg'/'txt', 'filename': '...'}]

Core Parameters Explained

Parameter	Type	Default	Description
N	int	required	The composite integer to factor.
method	str	'matrix'	Circuit method: 'matrix' or 'operator'.
max_retries	int	15	Maximum number of random base $a$ attempts.
backend	str	'torch'	Simulation backend ('torch' supported).
device	str	'cpu'	Torch device, e.g. 'cpu' or 'cuda'.
dtype	dtype	np.complex128	NumPy dtype for simulation.

algo_dir is set at construction time via ShorAlgorithm(algo_dir=...), not in run().

Common misunderstandings:

• Not every attempt succeeds; the loop retries automatically. The returned factors may come from a classical shortcut (even N, shared GCD).
• period is None when a classical path (e.g., even N) is taken.
• Circuit qubit counts: n_count = 2 * N.bit_length() (counting register), plus n_work qubits for the work register.

Return Fields

_build_return_dict(success, circuit_path, filepath, circuit) assembles the result. It merges the base fields with self.output (set by self.update_output() before each return).

Key	Type	Description
status	str	'ok' on success, 'failed' if max retries exhausted.
circuit_path	str or None	Path to saved SVG circuit diagram; None for classical path.
plot	list	Saved output files: [{'format': 'svg'|'txt', 'filename': '<path>'}].
circuit	Circuit or None	The Circuit object (only present on failure path; None otherwise).
factors	list or None	Non-trivial factors found, e.g. [3, 5]; None on failure.
period	int or None	Period $r$ found via quantum path; None for classical shortcut.
Selected base	int or None	Random base $a$ chosen; None for even-$N$ shortcut.
Computation time (s)	float	Wall-clock time for qc.execute() (quantum path only).
Measurement	int	Integer value of the counting-register measurement (quantum path only).
Total qubits	int	Total qubits in the circuit (quantum path only).

Implementation Architecture

ShorAlgorithm in algorithm.py uses a retry loop around the five-stage QPE framework. Classical shortcuts and random base selection are handled before the quantum circuit pipeline.

run(N, method='matrix', max_retries=15, backend='torch', device='cpu', dtype=np.complex128) — Outer Loop + Five Stages:

Stage	Code Action	Algorithmic Role
Pre-check	Handles N % 2 == 0 and gcd(a, N) > 1 cases without any quantum circuit	Classical short-circuit for trivial cases
1 — Parameter Setup	Picks random a; computes n_work = N.bit_length(), n_count = 2 * n_work, n_work_actual (method-dependent), total_qubits	Determines qubit counts for the chosen method
2 — Circuit Construction	Creates Circuit(total_qubits, name=...); applies H to counting register via qc.h(range(n_count)); sets work register to $	1\rangle$ via qc.x(n_count); calls _build_modular_matrix_circuit or _build_modular_operator_circuit; appends IQFT(n_count) via qc.append(IQFT(n_count), range(n_count))
3 — Simulation	result = qc.execute(backend=backend, device=device, dtype=dtype); measure_bin = result.measure(range(n_count), endian='little'); measure_int = int(measure_bin, 2)	Single-shot measurement of counting register
4 — Classical Post-Processing	phase = measure_int / (2**n_count); Fraction(phase).limit_denominator(N).denominator gives r; checks r % 2 == 0; computes gcd(a^(r/2)±1, N)	Continued fractions + factor extraction
5 — Export	self.save_circuit(qc) and self.save_txt() (on success and failure)	Saves SVG circuit diagram and text results

Two Circuit Methods:

• _build_modular_matrix_circuit(qc, n_count, n_work, a, N) — For each counting qubit q, computes power_factor = a^(2^q) mod N, builds a permutation matrix via _get_modular_matrix(power_factor, N, n_work), and calls qc.unitary(matrix, work_qubits, control=q). Simple and direct, but requires explicit matrix construction.

• _build_modular_operator_circuit(qc, n_count, n_work, n_work_actual, a, N) — Uses an algebraic operator decomposition. Core sub-components:

  • _multiple_mod(n_qubits, a, N) — Controlled modular multiplier; built from _Add_constant_mod_opt sub-circuits.
  • _Add_constant_mod_opt(n_qubits, a, N) — Quantum modular adder using QFT-domain phase rotations (_Ph, _Controlled_Ph, QFT, IQFT).

Helper Methods:

• _get_modular_matrix(a, N, n_qubits) — Builds the permutation matrix for modular multiplication.
• _Ph(n, a, qc) / _Controlled_Ph(n, a, qc, control_qubit, data_qubits) — Apply phase rotations in QFT domain.
• self.update_output(output) / _build_return_dict(success, circuit_path, filepath, circuit) — Store runtime fields and package result dict.

Data flow: N → random a → circuit method dispatch → Circuit → execute() → result.measure() → continued fractions → gcd() → factors → _build_return().

Understanding the Key Quantum Components

The $n_{\text{count}} = 2 \cdot \lfloor\log_2 N\rfloor$ counting qubits are placed in uniform superposition via Hadamard: $$\frac{1}{\sqrt{2^{n_{\text{count}}}}}\sum_{x=0}^{2^{n_{\text{count}}}-1}|x\rangle|1\rangle_{\text{work}}$$

2. Modular Exponentiation Oracle

For each counting qubit $k$, a controlled operation applies $\times a^{2^k} \bmod N$ to the work register: $$\sum_x |x\rangle |a^x \bmod N\rangle$$ The work register starts at $|1\rangle$ so that $a^x \bmod N$ is correctly accumulated.

3. Inverse QFT (Period Extraction)

After applying the modular exponentiation, the IQFT on the counting register produces peaks at multiples of $2^{n_{\text{count}}}/r$, where $r$ is the period. Measuring the counting register gives a value $\approx k \cdot 2^{n_{\text{count}}}/r$.

4. Continued Fractions (Classical)

The rational approximation $k/r \approx m/2^{n_{\text{count}}}$ is extracted using the continued fractions algorithm (Fraction.limit_denominator(N)), yielding the period $r$.

5. Factor Extraction

If $r$ is even and $a^{r/2} \not\equiv -1 \pmod{N}$: $$\gcd(a^{r/2} - 1, N) \quad \text{and} \quad \gcd(a^{r/2} + 1, N)$$ yields non-trivial factors of $N$.

Theory-to-Code Mapping

README / Theory Concept	Code Object or Location
Counting register $H^{\otimes n}	0\rangle$
Work register initialized to $	1\rangle$
Controlled $a^{2^k} \bmod N$ (matrix method)	_build_modular_matrix_circuit() — qc.unitary(matrix, work, control=q)
Controlled $a^{2^k} \bmod N$ (operator method)	_build_modular_operator_circuit() via _multiple_mod() using QFT-domain adders
Inverse QFT (iQFT)	qc.append(IQFT(n_count), range(n_count)) from unitarylab.library
Measurement of counting register	result = qc.execute(...); result.measure(range(n_count), endian='little') → binary string → integer
Phase $\phi \approx k/r$	phase = measure_int / 2^n_count
Continued fractions algorithm	Fraction(phase).limit_denominator(N).denominator → r
Factor extraction $\gcd(a^{r/2}\pm1, N)$	math.gcd(guess - 1, N) and math.gcd(guess + 1, N)
Classical shortcut (even N or gcd>1)	Pre-loop checks before any quantum circuit construction

Notes on method differences: The 'matrix' method builds an explicit permutation matrix for each controlled power, making it circuit-efficient but memory-intensive for large $N$. The 'operator' method uses QFT-domain phase additions and is more gate-intensive but architecturally general. Both append the same IQFT from unitarylab.library. The retry loop is necessary because the continued fractions step can fail for unlucky choices of a or measurement outcomes.

Mathematical Deep Dive (the multiplicative order of $a$ modulo $N$), i.e., $a^r \equiv 1 \pmod{N}$.

Factoring step: Since $a^r - 1 = (a^{r/2}-1)(a^{r/2}+1) \equiv 0 \pmod{N}$, and $N$ divides the product but not (with high probability) either factor individually, $\gcd(a^{r/2} \pm 1, N)$ yields non-trivial factors.

Success probability: At least $50%$ of bases $a$ produce a useful period. Multiple attempts are needed with probability exponentially small to fail.

Complexity: Polynomial in $\log N$; a common textbook estimate is $O((\log N)^3)$ quantum gates, but the exact exponent depends on the modular multiplication/arithmetic implementation and cost model. Classical best: sub-exponential $e^{O((\log N)^{1/3}(\log\log N)^{2/3})}$.

Hands-On Example

from unitarylab_algorithms import ShorAlgorithm

algo = ShorAlgorithm()
result = algo.run(N=15, method='matrix', backend='torch', max_retries=20)
print(f"Status: {result['status']}")                  # 'ok' or 'failed'
print(f"Factors of 15: {result['factors']}")          # [3, 5] or [5, 3]
print(f"Period r: {result['period']}")
print(f"Selected base a: {result['Selected base']}")
print(f"Circuit path: {result['circuit_path']}")
for f in result['plot']:                               # saved output files
    print(f["format"], f["filename"])

Reference Implementation (Classiq)

Classiq can be used as a reference implementation for Shor-style period finding. It describes the modular multiplication based QPE workflow through high-level QMOD functions, automatic synthesis, and built-in modular arithmetic primitives.

Example A: Minimal Classiq Shor Period-Finding Run

from classiq import *


@qfunc
def period_finding(n: CInt, a: CInt, x: QNum, phase_var: QNum):
    x ^= 1
    qpe_flexible(
        lambda p: modular_multiply_constant_inplace(n, a**p, x),
        phase_var,
    )


modulo_num = 21
a_num = 11
x_len = modulo_num.bit_length()
phase_len = 2 * x_len


@qfunc
def main(phase_var: Output[QNum[phase_len, UNSIGNED, phase_len]]):
    x = QNum()
    allocate(x_len, x)
    allocate(phase_var)

    period_finding(modulo_num, a_num, x, phase_var)

    drop(x)


qprog = synthesize(
    main,
    preferences=Preferences(qasm3=True, optimization_level=1),
    constraints=Constraints(optimization_parameter="width"),
)

sample_result = execute(qprog).get_sample_result()
df = sample_result.dataframe
print(df)

Minimal Manual Implementation (UnitaryLab)

Below is a skeleton that reconstructs Shor's algorithm at the component level, matching the ShorAlgorithm structure in algorithm.py.

# Simplified reconstruction — mirrors ShorAlgorithm.run() (matrix method)
import math, random
from fractions import Fraction
import numpy as np
from unitarylab.core import Circuit
from unitarylab.library import IQFT

def get_modular_matrix(mult: int, N: int, n_work: int) -> np.ndarray:
    """Build the 2^n_work x 2^n_work unitary for |x> -> |x*mult mod N>."""
    dim = 2**n_work
    U = np.zeros((dim, dim))
    for x in range(dim):
        y = (x * mult) % N if x < N else x
        U[y, x] = 1.0
    return U.astype(complex)

def build_shor_circuit(a: int, N: int, n_count: int, n_work: int) -> Circuit:
    """Construct the quantum phase-estimation circuit for order finding."""
    total = n_count + n_work
    qc = Circuit(total)

    # 1. Superpose counting register; set work register to |1>
    qc.h(range(n_count))
    qc.x(n_count)

    # 2. Controlled-U^{2^j} gates (matrix method): |j>|x> -> |j>|x * a^{2^j} mod N>
    for j in range(n_count):
        mult = pow(a, 2**j, N)
        U = get_modular_matrix(mult, N, n_work)
        qc.unitary(U, range(n_count, total), j, '1')  # controlled on qubit j

    # 3. Inverse QFT on counting register
    qc.append(IQFT(n_count), range(n_count))
    return qc

def shor_factor(N: int, max_retries: int = 15, backend: str = 'torch'):
    """Full Shor factoring loop."""
    if N % 2 == 0:
        return [2, N // 2]

    n_work  = N.bit_length()
    n_count = 2 * n_work

    for _ in range(max_retries):
        a = random.randint(2, N - 1)
        g = math.gcd(a, N)
        if g > 1:
            return [g, N // g]   # classical shortcut

        qc = build_shor_circuit(a, N, n_count, n_work)

        # Measure counting register
        result = qc.execute(backend=backend)
        meas = result.measure(range(n_count), endian='little')
        phase_int = int(meas, 2)

        # Classical post-processing: continued fractions -> period r
        phase = phase_int / (2**n_count)
        r = Fraction(phase).limit_denominator(N).denominator

        if r % 2 != 0 or r == 0:
            continue
        guess = pow(a, r // 2, N)
        if guess in (1, N - 1):
            continue

        p = math.gcd(guess - 1, N)
        q = math.gcd(guess + 1, N)
        if p * q == N:
            return [p, q]
    return None  # failed

Component roles:

• get_modular_matrix — builds the $2^n \times 2^n$ permutation matrix for $|x\rangle \to |x \cdot \text{mult} \bmod N\rangle$. This is the _build_modular_matrix_circuit bypass used in the 'matrix' method.
• build_shor_circuit — assembles $H^{\otimes n_{\text{count}}}$ + controlled-$U^{2^j}$ gates + IQFT, the three-stage quantum phase estimation structure.
• shor_factor — the outer retry loop: picks random $a$, builds + measures the circuit, applies continued fractions to extract period $r$, then computes $\gcd(a^{r/2}\pm 1, N)$.

Debugging Tips

1. Simulation is slow for large N: Qubit count scales as $O(\log N)$, but the state vector grows as $2^{O(\log N)} = N^O$. Keep $N$ small (e.g., 15, 21, 35) for feasible simulation.
2. status='failed' after max_retries: Increase max_retries. Some choices of $N$ and $a$ lead to trivial periods more often.
3. period is odd: Odd periods are discarded; the loop retries. This is expected for some $a$.
4. factors = [1, N]: Trivial factors — the period-finding step gave an $a^{r/2} \equiv \pm 1 \pmod N$ case; retries automatically.
5. method='operator' vs 'matrix': 'matrix' uses dense matrix multiplication; 'operator' uses a modular multiplication operator circuit. Both are correct; 'matrix' is simpler.