|
| 1 | +import numpy as np |
| 2 | +import matplotlib.pyplot as plt |
| 3 | +from matplotlib.animation import FuncAnimation |
| 4 | +from matplotlib import cm |
| 5 | + |
| 6 | + |
| 7 | + |
| 8 | +def unitary_operator(): |
| 9 | + # Controlled Phase Gate (CP) with a phase of pi/4 |
| 10 | + return np.array([[1, 0, 0, 0], |
| 11 | + [0, 1, 0, 0], |
| 12 | + [0, 0, np.exp(1j * np.pi / 4), 0], |
| 13 | + [0, 0, 0, np.exp(1j * np.pi / 4)]], dtype=complex) |
| 14 | + |
| 15 | +def eigenstate(): |
| 16 | + # Superposition state |ψ⟩ = (|0⟩ + |1⟩)/√2 |
| 17 | + return np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex) |
| 18 | + |
| 19 | +class QuantumPhaseEstimation: |
| 20 | + def __init__(self, n_qubits, unitary, eigenstate, inverse=False): |
| 21 | + self.n_qubits = n_qubits |
| 22 | + self.N = 2 ** n_qubits |
| 23 | + self.unitary = unitary # Unitary operator U |
| 24 | + self.eigenstate = eigenstate # Eigenstate |ψ⟩ |
| 25 | + self.inverse = inverse # Whether to reverse the QFT |
| 26 | + self.state = np.zeros(self.N, dtype=complex) |
| 27 | + self.state[0] = 1.0 # Initialize ancillary qubits to |0...0⟩ |
| 28 | + |
| 29 | + def apply_single_qubit_gate(self, gate, target): |
| 30 | + I = np.eye(2) |
| 31 | + ops = [I] * self.n_qubits |
| 32 | + ops[target] = gate |
| 33 | + U = ops[0] |
| 34 | + for op in ops[1:]: |
| 35 | + U = np.kron(U, op) |
| 36 | + self.state = U @ self.state |
| 37 | + |
| 38 | + def hadamard(self): |
| 39 | + return np.array([[1, 1], [1, -1]]) / np.sqrt(2) |
| 40 | + |
| 41 | + def apply_controlled_unitary(self, control, target, exp_power): |
| 42 | + """Apply controlled unitary operation U^(2^j).""" |
| 43 | + U = np.eye(self.N, dtype=complex) |
| 44 | + for i in range(self.N): |
| 45 | + b = format(i, f'0{self.n_qubits}b') |
| 46 | + if b[self.n_qubits - 1 - control] == '1': |
| 47 | + U[i, i] = np.exp(1j * (np.pi / (2 ** exp_power))) |
| 48 | + self.state = U @ self.state |
| 49 | + |
| 50 | + def apply_qft(self, inverse=False): |
| 51 | + """Apply Quantum Fourier Transform on ancillary qubits.""" |
| 52 | + for target in range(self.n_qubits): |
| 53 | + idx = self.n_qubits - 1 - target |
| 54 | + for offset in range(1, self.n_qubits - target): |
| 55 | + control = self.n_qubits - 1 - (target + offset) |
| 56 | + angle = np.pi / (2 ** offset) |
| 57 | + if inverse: |
| 58 | + angle *= -1 |
| 59 | + self.apply_controlled_unitary(control, idx, angle) |
| 60 | + self.apply_single_qubit_gate(self.hadamard(), idx) |
| 61 | + |
| 62 | + self.swap_registers() |
| 63 | + |
| 64 | + def swap_registers(self): |
| 65 | + perm = [int(format(i, f'0{self.n_qubits}b')[::-1], 2) for i in range(self.N)] |
| 66 | + self.state = self.state[perm] |
| 67 | + |
| 68 | + def measure(self): |
| 69 | + """Measure the ancillary qubits.""" |
| 70 | + probs = np.abs(self.state) ** 2 |
| 71 | + outcomes = np.random.choice(self.N, size=1024, p=probs) |
| 72 | + counts = {} |
| 73 | + for o in outcomes: |
| 74 | + bitstring = format(o, f'0{self.n_qubits}b') |
| 75 | + counts[bitstring] = counts.get(bitstring, 0) + 1 |
| 76 | + return counts |
| 77 | + |
| 78 | + def plot_probability_distribution(self, results): |
| 79 | + bitstrings = sorted(results.keys()) |
| 80 | + counts = [results[b] for b in bitstrings] |
| 81 | + plt.bar(bitstrings, counts) |
| 82 | + plt.xlabel("Bitstring") |
| 83 | + plt.ylabel("Counts") |
| 84 | + plt.title("QPE Measurement Results") |
| 85 | + plt.xticks(rotation=45) |
| 86 | + plt.tight_layout() |
| 87 | + plt.show() |
| 88 | + |
| 89 | + def apply_phase_estimation(self): |
| 90 | + """Perform Quantum Phase Estimation (QPE).""" |
| 91 | + # Apply Hadamard to all ancillary qubits |
| 92 | + for target in range(self.n_qubits): |
| 93 | + self.apply_single_qubit_gate(self.hadamard(), target) |
| 94 | + |
| 95 | + # Apply controlled unitaries (U^(2^j)) |
| 96 | + for j in range(self.n_qubits): |
| 97 | + self.apply_controlled_unitary(j, self.n_qubits, j) |
| 98 | + |
| 99 | + # Apply QFT |
| 100 | + self.apply_qft(inverse=True) |
| 101 | + |
| 102 | + # Measure the ancillary qubits |
| 103 | + results = self.measure() |
| 104 | + self.plot_probability_distribution(results) |
| 105 | + |
| 106 | +# --------------------------- |
| 107 | +# Example usage |
| 108 | +# --------------------------- |
| 109 | +if __name__ == "__main__": |
| 110 | + n_qubits = 3 # Number of ancillary qubits |
| 111 | + qpe = QuantumPhaseEstimation(n_qubits, unitary_operator(), eigenstate()) |
| 112 | + qpe.apply_phase_estimation() |
0 commit comments