QCQI — Chapter 7: Shor’s Algorithm & Phase Estimation

# Backend selector
backend = "backend"  # change: "pennylane", "braket", "pyquil", "qiskit"
print("Selected backend:", backend)
# Optional installs (if needed)
# !pip install cirq pennylane amazon-braket-sdk pyquil qiskit qiskit-aer

import numpy as np
def phase_unitary(phi):
    # One-qubit U with eigenstate |1> and eigenphase exp(2πi phi)
    return np.array([[1,0],[0,np.exp(2j*np.pi*phi)]], dtype=complex)

if backend == "cirq":
    import cirq, numpy as np
    def qpe_cirq(phi, t=6, shots=2048):
        ctrl = [cirq.LineQubit(i) for i in range(t)]
        tgt = cirq.LineQubit(t)
        c = cirq.Circuit()
        # prepare eigenstate |1> on target
        c.append(cirq.X(tgt))
        # put controls in |+>
        c.append(cirq.H.on_each(*ctrl))
        # controlled-U^{2^k}
        for k, q in enumerate(ctrl):
            angle = 2*np.pi*phi*(2**k)
            # U^{2^k} = Rz(2π phi 2^k) on |1>, implement via controlled phase (Z-rotation)
            c.append(cirq.CZPowGate(exponent=angle/np.pi).on(q, tgt))
        # inverse QFT on controls
        # build manually
        for j in range(t//2):
            c.append(cirq.SWAP(ctrl[j], ctrl[t-1-j]))
        for j in range(t):
            for k2 in range(j):
                c.append(cirq.CZPowGate(exponent=-1/2**(j-k2)).on(ctrl[j], ctrl[k2]))
            c.append(cirq.H(ctrl[j]))
        c.append(cirq.measure(*ctrl, key='m'))
        sim = cirq.Simulator()
        res = sim.run(c, repetitions=shots).histogram(key='m')
        return res
    res = qpe_cirq(phi=0.34375, t=6, shots=1024)  # 0.34375 = 22/64
    print("Top outcomes:", sorted(res.items(), key=lambda x: x[1], reverse=True)[:5])

if backend == "pennylane":
    import pennylane as qml, numpy as np
    def qpe_pl(phi=0.3, t=5, shots=2048):
        dev = qml.device("default.qubit", wires=t+1, shots=shots)
        @qml.qnode(dev)
        def circuit():
            # target eigenstate |1>
            qml.PauliX(t)
            # put controls in |+>
            for j in range(t): qml.Hadamard(j)
            # controlled U^{2^k}: phase on |11> branch equals 2π phi 2^k
            for k in range(t):
                qml.ControlledPhaseShift(2*np.pi*phi*(2**k), wires=[k, t])
            # inverse QFT (manual, MSB-first; with swaps)
            for j in range(t//2):
                qml.SWAP(wires=[j, t-1-j])
            for j in range(t):
                for k2 in range(j):
                    qml.ControlledPhaseShift(-np.pi/2**(j-k2-1), wires=[j, k2])
                qml.Hadamard(j)
            return [qml.sample(qml.PauliZ(j)) for j in range(t)]
        bits = circuit()
        # Convert +/-1 samples to 0/1 by (1-s)/2 and aggregate majority
        maj = []
        for s in bits:
            b = (1 - s)/2
            maj.append(int(np.round(b.mean())))
        return maj
    print("Estimated bits (MSB..LSB):", qpe_pl(phi=0.3125, t=5))

if backend == "braket":
    from braket.circuits import Circuit
    from braket.devices import LocalSimulator
    import numpy as np
    def qpe_braket(phi=0.3, t=5, shots=2048):
        c = Circuit()
        # controls: 0..t-1, target: t
        c.x(t)
        for j in range(t): c.h(j)
        for k in range(t):
            c.cphaseshift(k, t, 2*np.pi*phi*(2**k))
        # inverse QFT (simple, not optimized)
        for j in range(t//2):
            c.swap(j, t-1-j)
        for j in range(t):
            for k2 in range(j):
                c.cphaseshift(j, k2, -np.pi/2**(j-k2-1))
            c.h(j)
        for j in range(t): c.measure(j)
        dev = LocalSimulator()
        res = dev.run(c, shots=shots).result().measurement_counts
        return res
    print("Counts (top few):", list(qpe_braket(0.28125, 6).items())[:5])

if backend == "pyquil":
    from pyquil import Program
    from pyquil.gates import H, X, CPHASE, SWAP, MEASURE
    from pyquil.api import get_qc
    import numpy as np
    def qpe_pyquil(phi=0.375, t=5, shots=1024):
        p = Program(); ro = p.declare('ro','BIT',t)
        tgt = t
        p += X(tgt)
        for j in range(t): p += H(j)
        for k in range(t):
            p += CPHASE(2*np.pi*phi*(2**k), k, tgt)
        for j in range(t//2): p += SWAP(j, t-1-j)
        for j in range(t):
            for k2 in range(j):
                p += CPHASE(-np.pi/2**(j-k2-1), j, k2)
            p += H(j)
            p += MEASURE(j, ro[j])
        qc = get_qc(f'{t+1}q-qvm')
        res = qc.run(p.wrap_in_numshots_loop(shots))
        bits = res.get_register_map()['ro']
        return bits
    print("Sample (first 5 runs):", qpe_pyquil(0.375, 5, 10)[:5])

if backend == "qiskit":
    from qiskit import QuantumCircuit
    from qiskit_aer import AerSimulator
    import numpy as np
    def qpe_qiskit(phi=0.34375, t=6, shots=2048):
        qc = QuantumCircuit(t+1, t)
        qc.x(t)  # target eigenstate |1>
        for j in range(t): qc.h(j)
        for k in range(t):
            qc.cp(2*np.pi*phi*(2**k), k, t)
        # inverse QFT on first t qubits
        for j in range(t//2):
            qc.swap(j, t-1-j)
        for j in range(t):
            for k2 in range(j):
                qc.cp(-np.pi/2**(j-k2-1), j, k2)
            qc.h(j)
        qc.measure(range(t), range(t))
        sim = AerSimulator()
        counts = sim.run(qc, shots=shots).result().get_counts()
        return counts
    print("Top counts:", qpe_qiskit(0.34375,6))

from fractions import Fraction
def recover_r_from_fraction(y, t, N, max_den=None):
    if max_den is None:
        max_den = N
    frac = Fraction(y, 2**t).limit_denominator(max_den)
    return frac.numerator, frac.denominator
print("Example CF recovery:", recover_r_from_fraction(22, 6, N=2**6))

import numpy as np
import math
def order_of(a, N):
    x=1
    for r in range(1, 10*N):
        x = (x*a) % N
        if x==1: return r
    return None
def shor_toy(N=15, a=2, t=8):
    # Classical order for reference
    r_true = order_of(a,N)
    # Simulate the QPE measurement outcome y (ideal noiseless)
    # Pick s uniformly from {0,...,r-1}
    s = np.random.randint(0, r_true)
    y = int(np.round((s/(r_true))*(2**t))) % (2**t)
    s_hat, r_hat = recover_r_from_fraction(y, t, N)
    from math import gcd
    if r_hat % 2 == 1 or pow(a, r_hat//2, N) == N-1:
        success = False
        factors = None
    else:
        f1 = math.gcd(pow(a, r_hat//2, N)-1, N)
        f2 = math.gcd(pow(a, r_hat//2, N)+1, N)
        success = (f1 not in [1,N]) or (f2 not in [1,N])
        factors = (f1, f2)
    return {"y": y, "r_true": r_true, "r_hat": r_hat, "success": success, "factors": factors}
print("Toy Shor run:", shor_toy(15,2,8))

def ipea_bits(phi, t=6):
    # LSB-first version; returns list of bits [b_t ... b_1] (MSB..LSB)
    est = 0.0
    bits = []
    for j in range(t,0,-1):
        # simulate outcome as nearest bit considering residual phase
        resid = (phi - est) % 1.0
        bit = int(resid >= 0.5)
        bits.append(bit)
        est += bit * (1/2**j)
    return bits[::-1]
print("IPEA bits for phi=11/32:", ipea_bits(11/32, t=5))