QCQI — Chapter 4: Quantum Circuits & Universality

# 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 Rx(a): return np.array([[np.cos(a/2), -1j*np.sin(a/2)],[ -1j*np.sin(a/2), np.cos(a/2)]], dtype=complex)
def Ry(a): return np.array([[np.cos(a/2), -np.sin(a/2)],[ np.sin(a/2),  np.cos(a/2)]], dtype=complex)
def Rz(a): return np.array([[np.exp(-1j*a/2), 0],[0, np.exp(1j*a/2)]], dtype=complex)
def global_phase_fix(U):
    det = np.linalg.det(U)
    return U / np.sqrt(det)  # normalize to det=1 (choose branch arbitrarily)
def zyz_decompose(U):
    # Ensure SU(2)
    U = global_phase_fix(U)
    # Handle edge cases
    if np.isclose(abs(U[0,0]), 0, atol=1e-12):
        beta = np.pi
        alpha = np.angle(-U[0,1]) - np.angle(U[1,0])
        gamma = 0.0
    else:
        beta = 2*np.arccos(np.clip(abs(U[0,0]), 0.0, 1.0))
        # compute phases
        a = U[0,0]/abs(U[0,0])
        b = U[0,1]/abs(U[0,0]) if not np.isclose(abs(U[0,0]),0) else 1.0
        alpha = np.angle(a) - np.angle(b)
        c = U[1,0]/abs(U[0,0])
        gamma = np.angle(a) - np.angle(c)
    # wrap to [-pi,pi)
    def wrap(x): 
        x = (x + np.pi) % (2*np.pi) - np.pi
        return x
    return wrap(alpha), wrap(beta), wrap(gamma)
def synth_from_zyz(alpha,beta,gamma):
    return Rz(alpha) @ Ry(beta) @ Rz(gamma)
def unitary_distance(U,V):
    # U,V in U(2). Use phase-insensitive distance via |Tr(U^\dagger V)|
    m = np.abs(np.trace(U.conj().T @ V))/2.0
    return float(np.clip(1-m, 0, 1))  # 0=identical, 1=orthogonal (heuristic)
np.set_printoptions(precision=4, suppress=True)

rng = np.random.default_rng(4)
# random angles
a0, b0, g0 = [float(x) for x in (rng.uniform(-np.pi, np.pi, size=3))]
U0 = synth_from_zyz(a0,b0,g0)
a,b,g = zyz_decompose(U0)
U1 = synth_from_zyz(a,b,g)
print("angles (true)     =", np.round([a0,b0,g0],4))
print("angles (recovered)=", np.round([a,b,g],4))
print("distance =", unitary_distance(U0,U1))

def rz_clifford_t_coarse(theta):
    # Approximate Rz(theta) ~ Rz(k * pi/4)
    step = np.pi/4
    k = int(np.round(theta/step))
    return Rz(k*step), k  # integer power of T up to phases
th = 0.63
U = Rz(th)
Uc, k = rz_clifford_t_coarse(th)
print("k ~", k, " distance =", unitary_distance(U,Uc))

# Starting state |+0>, apply controlled-Rz(theta) with control on qubit 0
def entanglement_from_cRz(theta):
    # Work in statevector space to keep framework-agnostic proof-of-principle
    H = (1/np.sqrt(2))*np.array([[1,1],[1,-1]], dtype=complex)
    Rz2 = Rz(theta)
    # Build 4x4 unitary for controlled-Rz on target qubit 1
    I = np.eye(2, dtype=complex)
    P0 = np.array([[1,0],[0,0]], dtype=complex)
    P1 = np.array([[0,0],[0,1]], dtype=complex)
    UcRz = np.kron(P0, I) + np.kron(P1, Rz2)
    # state |+0>
    psi = np.kron(H @ np.array([[1],[0]], dtype=complex), np.array([[1],[0]], dtype=complex))
    phi = UcRz @ psi
    rho = phi @ phi.conj().T
    # partial trace over qubit 1
    R00 = rho[0:2,0:2]; R11 = rho[2:4,2:4]
    rhoA = R00 + R11
    w = np.linalg.eigvalsh(rhoA)
    w = w[w>1e-14]
    S = float(-(w*np.log2(w)).sum())
    return S
for t in [0.0, np.pi/8, np.pi/4, np.pi/2, np.pi]:
    print("theta=", round(t,3), "  entanglement entropy ~", round(entanglement_from_cRz(t), 4))

if backend == "cirq":
    import cirq, numpy as np
    q = cirq.LineQubit(0)
    alpha, beta, gamma = 0.7, -1.1, 0.3
    circ = cirq.Circuit(cirq.rz(alpha)(q), cirq.ry(beta)(q), cirq.rz(gamma)(q), cirq.measure(q, key='m'))
    print("ZYZ circuit:"); print(circ)
    # Coarse Clifford+T approx for Rz(alpha)
    step = np.pi/4; k = int(np.round(alpha/step))
    approx = cirq.Circuit(cirq.rz(k*step)(q), cirq.ry(beta)(q), cirq.rz(gamma)(q), cirq.measure(q, key='m2'))
    print("Coarse C+T approx:"); print(approx)
    sim = cirq.Simulator()
    res = sim.simulate(circ[:-1])
    res2 = sim.simulate(approx[:-1])
    U_target = np.array(res.final_state_vector[:2]).reshape(2,1)  # visualize
    U_approx = np.array(res2.final_state_vector[:2]).reshape(2,1)
    print("State overlap ~", float(np.abs(np.vdot(U_target, U_approx))))

if backend == "pennylane":
    import pennylane as qml, numpy as np
    dev = qml.device("default.qubit", wires=1, shots=None)
    alpha, beta, gamma = 0.7, -1.1, 0.3
    @qml.qnode(dev)
    def synth():
        qml.RZ(alpha, 0); qml.RY(beta, 0); qml.RZ(gamma, 0)
        return qml.state()
    psi = synth()
    print("state norm =", np.linalg.norm(psi))

if backend == "braket":
    from braket.circuits import Circuit
    from braket.devices import LocalSimulator
    import numpy as np
    theta = np.pi/4
    circ = Circuit().h(0).cnot(0,1).rz(1, float(theta)).cnot(0,1).probability()
    dev = LocalSimulator()
    res = dev.run(circ, shots=1000).result()
    print("Probabilities (|00>,|01>,|10>,|11|):", np.round(res.values[0], 3))

if backend == "pyquil":
    from pyquil import Program
    from pyquil.gates import H, RZ, RX, CZ, CNOT, CCNOT, MEASURE
    from pyquil.api import get_qc
    import numpy as np
    p = Program()
    ro = p.declare('ro', 'BIT', 1)
    p += H(0)
    # emulate controlled-Rz via CNOT and single-qubit rotations:
    # cRz(theta) on target t with control c: CNOT(c,t); Rz(theta) on t; CNOT(c,t)
    theta = np.pi/4
    p += CNOT(0,1); p += RZ(theta, 1); p += CNOT(0,1)
    p += MEASURE(0, ro[0])
    qc = get_qc('2q-qvm')
    res = qc.run(p.wrap_in_numshots_loop(200))
    bits = res.get_register_map()['ro']
    print("Sample bit 0 mean:", 1 - 2*bits.mean())

if backend == "qiskit":
    from qiskit import QuantumCircuit
    from qiskit_aer import AerSimulator
    import numpy as np
    sim = AerSimulator()
    alpha, beta, gamma = 0.7, -1.1, 0.3
    qc = QuantumCircuit(2,1)
    qc.ry(np.pi/2, 0); qc.ry(-np.pi/2, 0)  # no-op just to show structure
    qc.rz(alpha, 1); qc.ry(beta, 1); qc.rz(gamma, 1)
    qc.h(0); qc.crz(np.pi/4, 0, 1); qc.measure(0,0)
    res = sim.run(qc, shots=1000).result().get_counts()
    print(res)