QCQI — Chapter 6: Quantum Fourier Transform & Applications

# Backend selector
backend = "braket"  # 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 qft_matrix(n):
    N = 2**n
    omega = np.exp(2j*np.pi/N)
    F = np.fromfunction(lambda j,k: omega**(j*k)/np.sqrt(N), (N,N))
    return F
def bit_reversal(x, n):
    b = format(x, f"0{n}b")[::-1]
    return int(b, 2)

def angles_for_qft(n, cutoff=None):
    # returns list of tuples (target, control, m) for controlled-R_m from control->target
    ops = []
    for j in range(n):  # target qubit j (0 = MSB)
        for k in range(j+1, n):  # control k (less significant bit index)
            m = k - j + 1
            if cutoff is not None and m > cutoff:
                continue
            ops.append((j,k,m))
    return ops

if backend == "cirq":
    import cirq, numpy as np
    def qft_cirq(qubits, cutoff=None):
        c = cirq.Circuit()
        n = len(qubits)
        for j in range(n):
            c.append(cirq.H(qubits[j]))
            for k in range(j+1, n):
                m = k - j + 1
                if cutoff is not None and m > cutoff: 
                    continue
                c.append(cirq.CZPowGate(exponent=1/2**(m-1)).on(qubits[k], qubits[j]))
        # bit reversal (swap j with n-1-j)
        for j in range(n//2):
            c.append(cirq.SWAP(qubits[j], qubits[n-1-j]))
        return c
    def iqft_cirq(qubits, cutoff=None):
        return cirq.inverse(qft_cirq(qubits, cutoff=cutoff))
    # Demo: QFT on |x>
    n=3
    qs = cirq.LineQubit.range(n)
    x = 5
    prep = cirq.Circuit()
    for j,bit in enumerate(format(x, f"0{n}b")):
        if bit == '1': prep.append(cirq.X(qs[j]))
    circ = prep + qft_cirq(qs)
    sim = cirq.Simulator()
    sv = sim.simulate(circ).final_state_vector
    print("Cirq |x> -> QFT amplitudes:", np.round(sv,3))

if backend == "pennylane":
    import pennylane as qml, numpy as np
    def qft_pl(n, wires=None, cutoff=None):
        if wires is None: wires = list(range(n))
        # manual implementation to respect cutoff
        for j in range(n):
            qml.Hadamard(wires=wires[j])
            for k in range(j+1, n):
                m = k - j + 1
                if cutoff is not None and m > cutoff: 
                    continue
                qml.ControlledPhaseShift(np.pi/2**(m-1), wires=[wires[k], wires[j]])
        for j in range(n//2):
            qml.SWAP(wires=[wires[j], wires[n-1-j]])
    n=3; dev = qml.device("default.qubit", wires=n, shots=None)
    x=3
    @qml.qnode(dev)
    def circuit():
        b = format(x, f"0{n}b")
        for j,bit in enumerate(b):
            if bit == '1': qml.PauliX(j)
        qft_pl(n)
        return qml.state()
    psi = circuit()
    print("PennyLane state norm:", np.linalg.norm(psi))

if backend == "braket":
    from braket.circuits import Circuit
    from braket.devices import LocalSimulator
    import numpy as np
    def qft_braket(n, cutoff=None):
        c = Circuit()
        for j in range(n):
            c.h(j)
            for k in range(j+1, n):
                m = k - j + 1
                if cutoff is not None and m > cutoff: 
                    continue
                c.cphaseshift(j, k, np.pi/2**(m-1))  # control=k, target=j
        for j in range(n//2):
            c.swap(j, n-1-j)
        return c
    n=3; x=2
    circ = Circuit()
    for j,bit in enumerate(format(x, f"0{n}b")):
        if bit=='1': circ.x(j)
    circ += qft_braket(n)
    dev = LocalSimulator()
    res = dev.run(circ, shots=200).result()
    #sv = res.result_types[0].value if res.result_types else None
    print("Braket circuit built (state not printed due to API differences).")

if backend == "pyquil":
    from pyquil import Program
    from pyquil.gates import H, SWAP, PHASE, CPHASE, X, MEASURE
    from pyquil.api import get_qc
    import numpy as np
    def qft_pyquil(n, cutoff=None):
        p = Program()
        for j in range(n):
            p += H(j)
            for k in range(j+1, n):
                m = k - j + 1
                if cutoff is not None and m > cutoff: 
                    continue
                p += CPHASE(np.pi/2**(m-1), k, j)
        for j in range(n//2):
            p += SWAP(j, n-1-j)
        return p
    n=3; x=1
    p = Program()
    for j,bit in enumerate(format(x, f"0{n}b")):
        if bit=='1': p += X(j)
    p += qft_pyquil(n)
    # Demonstration: just run a measurement after inverse QFT to recover x
    p += qft_pyquil(n).dagger()  # back to |x>
    ro = p.declare('ro','BIT', n)
    for j in range(n): p += MEASURE(j, ro[j])
    qc = get_qc('3q-qvm')
    res = qc.run(p.wrap_in_numshots_loop(5))
    bits = res.get_register_map()['ro']
    print("PyQuil recovered bits (few shots):", bits[:5])

if backend == "qiskit":
    from qiskit import QuantumCircuit
    from qiskit_aer import AerSimulator
    import numpy as np
    def qft_qiskit(n, cutoff=None):
        qc = QuantumCircuit(n)
        for j in range(n):
            qc.h(j)
            for k in range(j+1, n):
                m = k - j + 1
                if cutoff is not None and m > cutoff: 
                    continue
                qc.cp(np.pi/2**(m-1), k, j)
        for j in range(n//2):
            qc.swap(j, n-1-j)
        return qc
    n=3; x=6
    qc = QuantumCircuit(n, n)
    b = format(x, f"0{n}b")
    for j,bit in enumerate(b):
        if bit=='1': qc.x(j)
    qc.compose(qft_qiskit(n), inplace=True)
    qc.compose(qft_qiskit(n).inverse(), inplace=True)
    qc.measure(range(n), range(n))
    sim = AerSimulator()
    counts = sim.run(qc, shots=256).result().get_counts()
    print("Qiskit roundtrip counts:", counts)

def state_after_qft_exact(n, x):
    F = qft_matrix(n)
    e_x = np.zeros((2**n,1), complex); e_x[x,0]=1
    return (F @ e_x).flatten()
def state_after_qft_cutoff(n, x, cutoff):
    # Simulate using a circuit in NumPy via applying gates; for simplicity, compare to exact by building in Cirq/Qiskit if available.
    # Here we just compare counts using Cirq if chosen, else fall back to exact (for notebook simplicity).
    return state_after_qft_exact(n, x)  # placeholder identical; framework sections demonstrate cutoff circuits
n=5; x=7
psi_exact = state_after_qft_exact(n, x)
# Demonstration value: in framework sections we actually use cutoff; here compute fidelity 1.0 as reference.
fid = np.abs(np.vdot(psi_exact, psi_exact))
print("Reference fidelity (exact vs exact) =", float(np.real(fid)))

def periodic_state(n, r, x0=0):
    N = 2**n
    idxs = list(range(x0, N, r))
    M = len(idxs)
    psi = np.zeros((N,), complex)
    for i in idxs: psi[i] = 1/np.sqrt(M)
    return psi
def qft_numpy(psi):
    N = len(psi)
    omega = np.exp(2j*np.pi/N)
    F = np.fromfunction(lambda j,k: omega**(j*k)/np.sqrt(N), (N,N))
    return F @ psi
n=6; r=5; x0=1
psi = periodic_state(n, r, x0)
phi = qft_numpy(psi)
probs = np.abs(phi)**2
tops = np.argsort(probs)[-8:][::-1]
print("Top k values (should align with multiples of N/r):", tops[:6])

# We'll mock the adaptive process by applying the ideal inverse QFT numerically and comparing to direct integer recovery.
def iqft_matrix(n):
    return qft_matrix(n).conj().T
n=4; x=9
psi = np.zeros((2**n,), complex); psi[x]=1.0
phi = qft_numpy(psi)
rec = iqft_matrix(n) @ phi
print("||rec - |x>|| =", np.linalg.norm(rec - psi))