Hands-on with Shor’s Algorithm: Cracking Secrets the Quantum Way🔓

Today felt like peeking behind the universe’s lockbox 🔐 — where Shor’s algorithm whispered how even the toughest numbers can crack when quantum steps in. Not with brute force, but with rhythm, interference, and a sprinkle of weirdness. ✨

It’s not every day you break big numbers into prime factors (classically HARD!) using the Quantum magic✨.

What’s Shor’s Algorithm? 🧮

At its quantum heart, Shor’s algorithm does what classic computers struggle with: factorizing large numbers into primes — a problem that underpins modern cryptography. The trick? Quantum computers can find a “period” (order) of modular exponential functions crazy-fast, letting us extract nontrivial divisors easily. Imagine undoing the locks on RSA encryption using quantum computing — yikes (and wow)! 🔓

My Hands-On: Completing and Testing the Code

I tackled a notebook exercise where much of Shor’s quantum core was left for me to finish. Here’s what I did:

• Built the modular multiplication unitary (Uₐ) so the quantum computer can “compute in superposition.”
• Created controlled-powered unitaries to encode the modular exponentiation for ALL quantum states… at once!
• Implemented a hand-made QFT (Quantum Fourier Transform), skipping the comfort of Qiskit’s built-in magic for the raw thrill of building it from scratch.
• Pieced it all together into a full quantum order-finding circuit and glued it to the classical post-processing needed to turn quantum insights into actual factors.

import math
import numpy as np
from fractions import Fraction
from collections import Counter
from qiskit import QuantumCircuit, transpile
from qiskit_aer import AerSimulator
from qiskit.circuit.library import UnitaryGate
from qiskit import Aer

# Helper math functions
def iscoprime(a, b):
    return math.gcd(a, b) == 1

def continuedfractionexpansion(x, maxden=1000000):
    frac = Fraction(x).limit_denominator(maxden)
    return frac.numerator, frac.denominator

def fractionfromphase(phase, qcount):
    y_over_2q = phase / (2**qcount)
    num, den = continuedfractionexpansion(y_over_2q, maxden=2**qcount)
    return num, den

def getmostlikelyresult(counts):
    measuredstr, _ = Counter(counts).most_common(1)[0]
    return int(measuredstr, 2), measuredstr

# Quantum building blocks
def buildmultiplicationmodmatrix(a, N):
    n = int(np.ceil(np.log2(N)))
    dim = 2**n
    U = np.zeros((dim, dim))
    for x in range(dim):
        if x < N:
            y = (a * x) % N
            U[y, x] = 1
        else:
            U[x, x] = 1
    return U

def controlledpoweredunitarygate(Umatrix, power):
    Upow = np.linalg.matrix_power(Umatrix, power)
    gate = UnitaryGate(Upow)
    return gate

# Custom inverse QFT
def custom_inverse_qft_circuit(num_qubits):
    qc = QuantumCircuit(num_qubits)
    for j in range(num_qubits // 2):
        qc.swap(j, num_qubits - 1 - j)
    for j in range(num_qubits):
        for k in range(j):
            qc.cp(-np.pi / float(2 ** (j - k)), k, j)
        qc.h(j)
    return qc

# Quantum order finding
def shororderfinding(N, a, qcount=None, shots=1024, seedsim=42, verbose=True):
    if not (1 < a < N and math.gcd(a, N) == 1):
        raise ValueError("a must be between 1 and N and coprime to N")

    n = math.ceil(math.log2(N))
    if qcount is None:
        qcount = 2 * n

    U = buildmultiplicationmodmatrix(a, N)
    totalqubits = qcount + n
    qc = QuantumCircuit(totalqubits, qcount)

    countingqubits = list(range(qcount))
    targetqubits = list(range(qcount, qcount + n))

    qc.x(targetqubits[0])  
    qc.h(countingqubits) 

    for j in range(qcount):
        power = 2 ** j
        Upowgate = controlledpoweredunitarygate(U, power).control(1)
        qc.append(Upowgate, [countingqubits[j]] + targetqubits)

    invqft = custom_inverse_qft_circuit(qcount)  # Custom IQFT
    qc.compose(invqft, qubits=countingqubits, inplace=True)

    qc.measure(countingqubits, list(range(qcount)))

    sim = AerSimulator(seed_simulator=seedsim)
    transpiled = transpile(qc, sim)
    result = sim.run(transpiled, shots=shots).result()
    counts = result.get_counts()

    yint, ystr = getmostlikelyresult(counts)
    if verbose:
        print(f"Most frequent measurement counting register: {ystr} (int: {yint})")

    s, r_cand = fractionfromphase(yint, qcount)
    if verbose:
        print(f"Continued fraction approx: s={s}, r_candidate={r_cand}")

    for mult in range(1, 11):
        r = r_cand * mult
        if r != 0 and pow(a, r, N) == 1:
            if verbose:
                print(f"Found order r: {r}")
            return r
    return None

# Classical postprocessing and Shor's factorization
def shorfactor(N, shots=1024, tries=5, verbose=True):
    if N % 2 == 0:
        return 2, N // 2

    for attempt in range(tries):
        a = np.random.randint(2, N-1)
        if math.gcd(a, N) != 1:
            d = math.gcd(a, N)
            if verbose:
                print(f"Non-trivial factor found by gcd: {d}")
            return d, N // d
        if verbose:
            print(f"Attempt {attempt+1}: trying base a={a}")
        r = shororderfinding(N, a, shots=shots, verbose=verbose)
        if r is None or r % 2 != 0:
            continue

        apow = pow(a, r // 2, N)
        candidate1 = math.gcd(apow - 1, N)
        candidate2 = math.gcd(apow + 1, N)

        if candidate1 != 1 and candidate1 != N:
            return candidate1, N // candidate1
        if candidate2 != 1 and candidate2 != N:
            return candidate2, N // candidate2
    return None

# Usage example
print(shorfactor(15))
print(shorfactor(21))

PS: This code snippet is adapted from my hands-on notebook.

The (Hand-made!) Inverse QFT Circuit 🌀

def custom_inverse_qft_circuit(num_qubits):
    qc = QuantumCircuit(num_qubits)
    for j in range(num_qubits//2):
        qc.swap(j, num_qubits - 1 - j)
    for j in range(num_qubits):
        for k in range(j):
            qc.cp(-np.pi / float(2 ** (j - k)), k, j)
        qc.h(j)
    return qc

GHZ State: Quantum Entanglement Playground! 😍

• Built and visualized a nice 3-qubit GHZ state circuit;
• Confirmed measurement outcomes show near-perfect “000” and “111” states — pure quantum!

def build_ghz_circuit(num_qubits):
    qc = QuantumCircuit(num_qubits, num_qubits)
    qc.h(0)
    for i in range(1, num_qubits):
        qc.cx(0, i)
    qc.measure(range(num_qubits), range(num_qubits))
    return qc

def simulate_circuit(qc, shots=1024):
    sim = Aer.get_backend("qasm_simulator")
    transpiled = transpile(qc, sim)
    result = sim.run(transpiled, shots=shots).result()
    counts = result.get_counts()
    return counts

num_qubits = 3 
qc = build_ghz_circuit(num_qubits)
print("Generated GHZ circuit:")
print(qc.draw())

counts = simulate_circuit(qc, shots=1024)
print("Measurement counts:", counts)

Key Takeaways..?

→ Quantum programming can be quirky, lots of debugging and reading circuit outputs.

→ The magic happens at the intersection of clever math (continued fractions for order extraction!) and raw quantum weirdness.

→ Writing your own QFT by hand solidifies why quantum computers can do things the classical world just can’t.

Reflections..!

More than just getting code to work, this challenge was about intuition. Visualizing superpositions, debugging probability distributions, and that “aha!” moment when factors appear from quantum output. ❤️🔍 It’s hands-on proof that even the quirkiest, most abstract quantum algorithms are real.. and can actually run on simulators today.

If you’re pondering a dive into quantum algorithms, DO IT — write it, break it, and celebrate every quantum oddity you find!

*“What feels impossible today may just be waiting for a quantum twist of perspective.” 🔄✨*