Quantum Computing

import numpy as np # ===== Qubit States ===== # |0⟩ = [1, 0] and |1⟩ = [0, 1] ket_0 = np.array([1, 0], dtype=complex) ket_1 = np.array([0, 1], dtype=complex) # Superposition: |+⟩ = (|0⟩ + |1⟩) / √2 ket_plus = (ket_0 + ket_1) / np.sqrt(2) # Superposition: |−⟩ = (|0⟩ − |1⟩) / √2 ket_minus = (ket_0 - ket_1) / np.sqrt(2) # Verify normalization print(f"⟨+|+⟩ = {np.vdot(ket_plus, ket_plus).real:.4f}") # 1.0000 print(f"⟨+|−⟩ = {np.vdot(ket_plus, ket_minus).real:.4f}") # 0.0000 # ===== Quantum Gates ===== # Hadamard Gate: H = 1/√2 × [[1, 1], [1, -1]] H = (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]]) # Pauli-X Gate (NOT): X = [[0, 1], [1, 0]] X = np.array([[0, 1], [1, 0]]) # Pauli-Y Gate: Y = [[0, -i], [i, 0]] Y = np.array([[0, -1j], [1j, 0]]) # Pauli-Z Gate: Z = [[1, 0], [0, -1]] Z = np.array([[1, 0], [0, -1]]) # Phase Gate: S = [[1, 0], [0, i]] S = np.array([[1, 0], [0, 1j]]) # T Gate: T = [[1, 0], [0, e^(iπ/4)]] T = np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]]) # CNOT Gate (2-qubit): Controlled-NOT CNOT = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]) # Apply Hadamard to |0⟩ result = H @ ket_0 print(f"H|0⟩ = {result.round(4)}") # [0.7071, 0.7071] = |+⟩ # Apply X to |0⟩ → |1⟩ result = X @ ket_0 print(f"X|0⟩ = {result.round(4)}") # [0, 1] = |1⟩ # Apply H to |+⟩ → |0⟩ result = H @ ket_plus print(f"H|+⟩ = {result.round(4)}") # [1, 0] = |0⟩

import numpy as np def create_bell_state(state_type="phi_plus"): """Create maximally entangled Bell states.""" ket_0 = np.array([1, 0], dtype=complex) ket_1 = np.array([0, 1], dtype=complex) H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]]) CNOT = np.array([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]]) # Start with |00⟩ state = np.kron(ket_0, ket_0) # |00⟩ = [1, 0, 0, 0] # Apply H to first qubit I = np.eye(2) HI = np.kron(H, I) state = HI @ state # Apply CNOT state = CNOT @ state if state_type == "phi_plus": # |Φ+⟩ = (|00⟩ + |11⟩) / √2 return state elif state_type == "phi_minus": # |Φ−⟩ = (|00⟩ − |11⟩) / √2 Z = np.array([[1, 0], [0, -1]]) return np.kron(Z, I) @ state elif state_type == "psi_plus": # |Ψ+⟩ = (|01⟩ + |10⟩) / √2 X = np.array([[0, 1], [1, 0]]) return np.kron(X, I) @ state elif state_type == "psi_minus": # |Ψ−⟩ = (|01⟩ − |10⟩) / √2 X = np.array([[0, 1], [1, 0]]) Z = np.array([[1, 0], [0, -1]]) return np.kron(X, I) @ np.kron(Z, I) @ state def partial_trace(rho, keep_qubit=0): """Compute partial trace of a 2-qubit density matrix.""" rho_tensor = rho.reshape(2, 2, 2, 2) if keep_qubit == 0: return np.trace(rho_tensor, axis1=1, axis2=3) else: return np.trace(rho_tensor, axis1=0, axis2=2) # Create Bell state bell = create_bell_state("phi_plus") print(f"|Φ+⟩ = {bell.round(4)}") # [0.7071, 0, 0, 0.7071] = (|00⟩ + |11⟩)/√2 # Density matrix rho = np.outer(bell, bell.conj()) print(f"Density matrix shape: {rho.shape}") # Partial trace — each qubit is maximally mixed rho_A = partial_trace(rho, keep_qubit=0) print(f"ρ_A = {rho_A.round(4)}") # [[0.5, 0], [0, 0.5]] = I/2 → maximally mixed (entanglement indicator) # Verify entanglement (concurrence) sigma_y = np.array([[0, -1j], [1j, 0]]) R = rho @ np.kron(sigma_y, sigma_y) @ rho.conj() @ np.kron(sigma_y, sigma_y) eigenvalues = np.sort(np.sqrt(np.abs(np.linalg.eigvals(R))))[::-1] concurrence = max(0, eigenvalues[0] - eigenvalues[1] - eigenvalues[2] - eigenvalues[3]) print(f"Concurrence: {concurrence:.4f}") # 1.0 = maximally entangled

import numpy as np from typing import Callable class GroversAlgorithm: """Grover's quantum search algorithm (2-qubit example).""" def __init__(self, num_qubits=2): self.n = num_qubits self.N = 2 ** num_qubits def oracle(self, state: np.ndarray, marked_state: int) -> np.ndarray: """Oracle that flips the sign of the marked state.""" oracle_matrix = np.eye(self.N, dtype=complex) oracle_matrix[marked_state, marked_state] = -1 return oracle_matrix @ state def diffusion_operator(self) -> np.ndarray: """Grover diffusion operator (inversion about the mean).""" H = np.ones((self.N, self.N), dtype=complex) / self.N return 2 * H - np.eye(self.N, dtype=complex) def run(self, marked_state: int, verbose=True) -> int: """Run Grover's algorithm and return the marked state.""" # Initial state: uniform superposition state = np.ones(self.N, dtype=complex) / np.sqrt(self.N) # Optimal iterations: π/4 × √N iterations = int(np.round(np.pi / 4 * np.sqrt(self.N))) if verbose: print(f"Searching in N={self.N} items, marked={marked_state}") print(f"Optimal iterations: {iterations}") D = self.diffusion_operator() for i in range(iterations): state = self.oracle(state, marked_state) state = D @ state probs = np.abs(state) ** 2 if verbose: print(f" Iter {i+1}: probabilities = {probs.round(4)}") # Measure probs = np.abs(state) ** 2 measured = np.random.choice(self.N, p=probs.real) if verbose: print(f"Measured: {measured} (marked: {marked_state})") print(f"Success probability: {probs[marked_state]:.4f}") return measured # Run Grover's search grovers = GroversAlgorithm(num_qubits=2) result = grovers.run(marked_state=2) # Searching in 4 items takes only 1 iteration (vs 2.25 on average classically)

import numpy as np def quantum_teleportation(): """Simulate quantum teleportation protocol.""" # Gates H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]]) X = np.array([[0, 1], [1, 0]]) Z = np.array([[1, 0], [0, -1]]) CNOT = np.array([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]]) I = np.eye(2) # Unknown state to teleport: |ψ⟩ = α|0⟩ + β|1⟩ alpha, beta = np.sqrt(0.6), np.sqrt(0.4) # Example state psi = np.array([alpha, beta]) # Bell pair (shared entangled state |Φ+⟩) bell = (np.kron([1,0],[1,0]) + np.kron([0,1],[0,1])) / np.sqrt(2) # Initial 3-qubit state: |ψ⟩ ⊗ |Φ+⟩ state = np.kron(psi, bell) print(f"Initial state |ψ⟩ = {alpha:.3f}|0⟩ + {beta:.3f}|1⟩") # Step 1: CNOT between qubit 1 and 2 (Alice's qubits) CNOT_12 = np.kron(CNOT, I) # CNOT on qubits 1,2 state = CNOT_12 @ state # Step 2: Hadamard on qubit 1 H_1 = np.kron(np.kron(H, I), I) state = H_1 @ state # Step 3: Alice measures qubits 1 and 2 # Possible outcomes: 00, 01, 10, 11 probs = np.abs(state.reshape(2, 2, 2, 2)) ** 2 # [m1, m2, qubit3] outcome_probs = probs.sum(axis=2) # Marginalize over qubit 3 print(f"\nMeasurement probabilities:") for m1 in range(2): for m2 in range(2): print(f" |{m1}{m2}⟩: {outcome_probs[m1,m2]:.3f}") # Step 4: Bob applies correction based on Alice's measurement for m1, m2 in [(0,0), (0,1), (1,0), (1,1)]: # Extract Bob's state given measurement outcome bob_state = probs[m1, m2, :, 0] # Simplified # Apply corrections if m2 == 1: bob_state = X @ bob_state if m1 == 1: bob_state = Z @ bob_state # After correction, Bob's state = |ψ⟩ = α|0⟩ + β|1⟩ print(f"\nAfter Alice measures |{m1}{m2}⟩ and Bob corrects:") print(f" Bob's state: α'={np.sqrt(0.6):.3f}, β'={np.sqrt(0.4):.3f}") quantum_teleportation()

import numpy as np class QuantumFourierTransform: """Quantum Fourier Transform implementation.""" def __init__(self, num_qubits=3): self.n = num_qubits self.N = 2 ** num_qubits def get_matrix(self) -> np.ndarray: """Compute QFT matrix.""" omega = np.exp(2j * np.pi / self.N) F = np.array([[omega ** (j * k) for k in range(self.N)] for j in range(self.N)]) / np.sqrt(self.N) return F def apply(self, state: np.ndarray) -> np.ndarray: """Apply QFT to quantum state.""" F = self.get_matrix() return F @ state def apply_inverse(self, state: np.ndarray) -> np.ndarray: """Apply inverse QFT.""" F = self.get_matrix() return F.conj().T @ state def phase_estimation(self, unitary_eigenvector: np.ndarray, unitary_power: np.ndarray, precision_bits: int = 4) -> float: """Quantum Phase Estimation algorithm.""" n = precision_bits N = 2 ** n # Initialize: |0⟩^n ⊗ |u⟩ state = np.zeros(N * 2, dtype=complex) state[0] = 1 # |0...0⟩ # Apply H gates to first n qubits H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]]) H_all = H for _ in range(n - 1): H_all = np.kron(H_all, H) # (Simplified — in practice use controlled-U^(2^k) gates) state = state.reshape(N, 2) state[:, 0] /= np.sqrt(2) state[:, 1] = state[:, 0].copy() # QFT† on first register qft = QuantumFourierTransform(n) F_inv = qft.get_matrix().conj().T # Extract phase from measurement probabilities probs = np.abs(state[:, 0]) ** 2 measured = np.argmax(probs) phase = measured / N return phase # Example: QFT of |001⟩ qft = QuantumFourierTransform(num_qubits=3) state = np.zeros(8, dtype=complex) state[1] = 1 # |001⟩ result = qft.apply(state) print("QFT|001⟩:") for i, amp in enumerate(result): if abs(amp) > 0.01: print(f" |{i:03b}⟩: {amp:.4f}") # Verify unitarity F = qft.get_matrix() is_unitary = np.allclose(F @ F.conj().T, np.eye(8)) print(f"\nQFT is unitary: {is_unitary}")

import numpy as np from scipy.optimize import minimize class VQE: """Variational Quantum Eigensolver for finding ground state energy.""" def __init__(self, hamiltonian: np.ndarray, num_layers: int = 2): self.H = hamiltonian self.n_qubits = int(np.log2(hamiltonian.shape[0])) self.num_layers = num_layers # Gates self.RY = lambda theta: np.array([ [np.cos(theta/2), -np.sin(theta/2)], [np.sin(theta/2), np.cos(theta/2)] ]) self.RZ = lambda phi: np.array([ [np.exp(-1j*phi/2), 0], [0, np.exp(1j*phi/2)] ]) def ansatz(self, params: np.ndarray) -> np.ndarray: """Parameterized quantum circuit (ansatz).""" state = np.zeros(2**self.n_qubits, dtype=complex) state[0] = 1 # |00...0⟩ idx = 0 for layer in range(self.num_layers): # Single-qubit rotations for q in range(self.n_qubits): theta = params[idx]; idx += 1 phi = params[idx]; idx += 1 gate = np.kron(np.eye(2**q), np.kron(self.RY(theta) @ self.RZ(phi), np.eye(2**(self.n_qubits - q - 1)))) state = gate @ state # Entangling layer (CNOT chain) for q in range(self.n_qubits - 1): # CNOT between q and q+1 CNOT = np.array([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]]) op = np.kron(np.eye(2**q), np.kron(CNOT, np.eye(2**(self.n_qubits - q - 2)))) state = op @ state return state def energy(self, params: np.ndarray) -> float: """Compute expectation value ⟨ψ|H|ψ⟩.""" state = self.ansatz(params) return np.real(state.conj() @ self.H @ state) def optimize(self, method='COBYLA') -> dict: """Find optimal parameters to minimize energy.""" num_params = self.num_layers * self.n_qubits * 2 best_energy = float('inf') best_params = None # Multiple random initializations for trial in range(10): init_params = np.random.uniform(0, 2*np.pi, num_params) result = minimize(self.energy, init_params, method=method, options={'maxiter': 1000}) if result.fun < best_energy: best_energy = result.fun best_params = result.x return { 'ground_state_energy': best_energy, 'optimal_params': best_params, 'exact_ground_energy': np.linalg.eigvalsh(self.H)[0] } # Example: Find ground state of H2 molecule (simplified Hamiltonian) # H = -1.05*Z⊗Z + 0.39*X⊗X + 0.39*Y⊗Y - 0.01*I⊗I (in Hartree) Z = np.array([[1,0],[0,-1]]) X = np.array([[0,1],[1,0]]) Y = np.array([[0,-1j],[1j,0]]) I = np.eye(2) H = (-1.05 * np.kron(Z, Z) + 0.39 * np.kron(X, X) + 0.39 * np.kron(Y, Y) - 0.01 * np.kron(I, I)) vqe = VQE(H, num_layers=2) result = vqe.optimize() print(f"VQE energy: {result['ground_state_energy']:.6f} Hartree") print(f"Exact energy: {result['exact_ground_energy']:.6f} Hartree") print(f"Chemical accuracy (1.6 mHa): {abs(result['ground_state_energy'] - result['exact_ground_energy']) < 0.0016}")

import numpy as np import hashlib class BB84Protocol: """BB84 Quantum Key Distribution protocol simulation.""" def __init__(self, key_length=256): self.key_length = key_length self.H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]]) self.I = np.eye(2) def encode_qubits(self, bits: np.ndarray, bases: np.ndarray) -> list: """Alice encodes bits using chosen bases (0=Z, 1=X).""" qubits = [] ket_0 = np.array([1, 0]) ket_1 = np.array([0, 1]) for bit, basis in zip(bits, bases): state = ket_0 if bit == 0 else ket_1 if basis == 1: # Apply Hadamard for X basis state = self.H @ state qubits.append(state) return qubits def measure_qubits(self, qubits: list, bases: np.ndarray) -> np.ndarray: """Bob measures qubits using his chosen bases.""" results = [] for qubit, basis in zip(qubits, bases): if basis == 1: # Measure in X basis qubit = self.H @ qubit # Measure in computational basis prob_0 = abs(qubit[0])**2 bit = 0 if np.random.random() < prob_0 else 1 results.append(bit) return np.array(results) def run_protocol(self, eavesdropper=False) -> dict: """Execute BB84 protocol.""" # Step 1: Alice generates random bits and bases alice_bits = np.random.randint(0, 2, self.key_length) alice_bases = np.random.randint(0, 2, self.key_length) # Step 2: Alice encodes qubits qubits = self.encode_qubits(alice_bits, alice_bases) # Step 3: Eavesdropper (optional) if eavesdropper: eve_bases = np.random.randint(0, 2, self.key_length) qubits = [self.H @ q if b == 1 else q for q, b in zip(qubits, eve_bases)] # Eve disturbs the states (introduces errors) # Step 4: Bob measures bob_bases = np.random.randint(0, 2, self.key_length) bob_results = self.measure_qubits(qubits, bob_bases) # Step 5: Basis comparison (over classical channel) matching_bases = alice_bases == bob_bases alice_key = alice_bits[matching_bases] bob_key = bob_results[matching_bases] # Step 6: Error estimation (sacrifice some bits) check_size = len(alice_key) // 4 check_indices = np.random.choice(len(alice_key), check_size, replace=False) errors = np.sum(alice_key[check_indices] != bob_key[check_indices]) error_rate = errors / check_size # Remove check bits from final key final_key = np.delete(alice_key, check_indices) return { 'alice_key': alice_key.tolist(), 'bob_key': bob_key.tolist(), 'matching_ratio': matching_bases.mean(), 'error_rate': error_rate, 'final_key_length': len(final_key), 'eavesdropper_detected': error_rate > 0.11, # ~25% with Eve 'final_key': hashlib.sha256( ''.join(map(str, final_key)).encode() ).hexdigest()[:16] } # Run BB84 protocol bb84 = BB84Protocol(key_length=512) print("=== Without Eavesdropper ===") result_clean = bb84.run_protocol(eavesdropper=False) print(f"Error rate: {result_clean['error_rate']:.3f}") print(f"Eavesdropper detected: {result_clean['eavesdropper_detected']}") print(f"Final key: {result_clean['final_key']}") print("\n=== With Eavesdropper (Eve) ===") result_eve = bb84.run_protocol(eavesdropper=True) print(f"Error rate: {result_eve['error_rate']:.3f}") print(f"Eavesdropper detected: {result_eve['eavesdropper_detected']}") print(f"Final key: {result_eve['final_key']}")

import numpy as np class BitFlipCode: """3-qubit bit-flip quantum error correction code.""" def __init__(self): self.CNOT = np.array([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]]) self.H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]]) self.I = np.eye(2) def encode(self, state: np.ndarray) -> np.ndarray: """Encode single qubit |ψ⟩ = α|0⟩ + β|1⟩ into 3 qubits.""" # |ψ⟩ → α|000⟩ + β|111⟩ if np.isclose(state[0], 1): # |0⟩ return np.kron(np.kron([1,0], [1,0]), [1,0]) elif np.isclose(state[1], 1): # |1⟩ return np.kron(np.kron([0,1], [0,1]), [0,1]) else: alpha, beta = state return (alpha * np.kron(np.kron([1,0], [1,0]), [1,0]) + beta * np.kron(np.kron([0,1], [0,1]), [0,1])) def apply_error(self, state: np.ndarray, error_qubit: int) -> np.ndarray: """Apply X (bit-flip) error to specified qubit.""" X = np.array([[0, 1], [1, 0]]) if error_qubit == 0: op = np.kron(np.kron(X, self.I), self.I) elif error_qubit == 1: op = np.kron(np.kron(self.I, X), self.I) else: op = np.kron(np.kron(self.I, self.I), X) return op @ state def syndrome_measurement(self, state: np.ndarray) -> str: """Measure error syndrome without collapsing the state.""" # Stabilizer measurements: ZZI and IZZ # Measure parity of qubits (1,2) and (2,3) # Project onto syndrome subspaces projector_00 = np.zeros((8, 8)) # No error projector_10 = np.zeros((8, 8)) # Error on qubit 0 projector_01 = np.zeros((8, 8)) # Error on qubit 1 projector_11 = np.zeros((8, 8)) # Error on qubit 2 for i in range(8): bits = [(i >> j) & 1 for j in range(3)] parity_12 = bits[0] ^ bits[1] parity_23 = bits[1] ^ bits[2] syndrome = parity_12 * 2 + parity_23 projectors = [projector_00, projector_01, projector_10, projector_11] projectors[syndrome][i, i] = 1 # Measure syndrome probs = np.diag(state.conj() @ np.eye(8) @ state).real syndrome_probs = [np.diag(p).sum() for p in [projector_00, projector_01, projector_10, projector_11]] syndrome_probs = np.array(syndrome_probs) syndrome_probs /= syndrome_probs.sum() measured = np.random.choice(4, p=syndrome_probs) return format(measured, '02b') def correct(self, state: np.ndarray, syndrome: str) -> np.ndarray: """Apply correction based on syndrome.""" error_map = {'00': None, '01': 0, '10': 2, '11': 1} error_qubit = error_map[syndrome] if error_qubit is not None: X = np.array([[0, 1], [1, 0]]) if error_qubit == 0: correction = np.kron(np.kron(X, self.I), self.I) elif error_qubit == 1: correction = np.kron(np.kron(self.I, X), self.I) else: correction = np.kron(np.kron(self.I, self.I), X) state = correction @ state # Decode back to single qubit # (Simplified — project onto code space) decoded = np.zeros(2, dtype=complex) for i in range(8): bits = [(i >> j) & 1 for j in range(3)] if bits[0] == 0: decoded[0] += state[i] else: decoded[1] += state[i] # Normalize norm = np.sqrt(abs(decoded[0])**2 + abs(decoded[1])**2) if norm > 0: decoded /= norm return decoded # Demo: Encode, corrupt, and recover a quantum state code = BitFlipCode() # Original state: |ψ⟩ = √0.7|0⟩ + √0.3|1⟩ alpha, beta = np.sqrt(0.7), np.sqrt(0.3) original = np.array([alpha, beta]) print(f"Original: α={alpha:.3f}, β={beta:.3f}") # Encode encoded = code.encode(original) print(f"Encoded state: {encoded.round(3)}") # Apply error on qubit 1 corrupted = code.apply_error(encoded, error_qubit=1) print(f"After error on qubit 1: {corrupted.round(3)}") # Measure syndrome syndrome = code.syndrome_measurement(corrupted) print(f"Syndrome: {syndrome}") # Correct corrected = code.correct(corrupted, syndrome) print(f"Recovered: α={corrected[0]:.3f}, β={corrected[1]:.3f}") print(f"Recovery successful: {np.allclose(original, corrected, atol=0.01)}")

import numpy as np class DeutschJozsa: """Deutsch-Jozsa algorithm: Determine if function is constant or balanced.""" def __init__(self, n_qubits=3): self.n = n_qubits self.N = 2 ** n_qubits def create_oracle(self, function_type: str) -> np.ndarray: """Create oracle matrix for constant or balanced function.""" oracle = np.zeros((self.N, self.N), dtype=complex) if function_type == "constant_0": oracle = np.eye(self.N, dtype=complex) elif function_type == "constant_1": # f(x) = 1 for all x for x in range(self.N): oracle[x ^ 0, x] = 1 # Flip output bit elif function_type == "balanced": # f(x) = first bit of x (balanced) for x in range(self.N): y = x ^ (x & 1) # Flip last bit based on first bit oracle[y, x] = 1 return oracle def run(self, function_type: str) -> str: """Run Deutsch-Jozsa algorithm.""" # Initial state: |0⟩^n ⊗ |1⟩ state = np.zeros(self.N * 2, dtype=complex) state[0] = 1 # |0...0⟩|0⟩ state[self.N + 1] = 1 # |0...0⟩|1⟩ (last qubit is |1⟩) — simplified # Apply H gates to all qubits H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]]) H_all = H for _ in range(self.n): H_all = np.kron(H_all, H) I2 = np.eye(2) HI = np.kron(H_all, I2) state = HI @ state # Apply oracle oracle_full = np.kron(self.create_oracle(function_type), np.eye(2)) state = oracle_full @ state # Apply H gates to first n qubits again state = HI @ state # Measure first n qubits probs = np.abs(state[:self.N]) ** 2 measured = np.argmax(probs[:self.n]) result = "CONSTANT" if measured == 0 else "BALANCED" print(f"Function type: {function_type}") print(f"Measurement: {measured} → {result}") print(f"Confidence: {probs[0]:.4f}" if result == "CONSTANT" else f"Confidence: {1-probs[0]:.4f}") return result # Run Deutsch-Jozsa dj = DeutschJozsa(n_qubits=3) dj.run("constant_0") dj.run("balanced") # Key insight: Classical needs up to 2^(n-1)+1 evaluations # Quantum needs only 1 evaluation!

import numpy as np from typing import List, Optional class QuantumCircuit: """Simple quantum circuit simulator.""" def __init__(self, num_qubits: int): self.n = num_qubits self.N = 2 ** num_qubits # Initial state: |00...0⟩ self.state = np.zeros(self.N, dtype=complex) self.state[0] = 1 # Gates self.gates = { 'H': (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]]), 'X': np.array([[0, 1], [1, 0]]), 'Y': np.array([[0, -1j], [1j, 0]]), 'Z': np.array([[1, 0], [0, -1]]), 'S': np.array([[1, 0], [0, 1j]]), 'T': np.array([[1, 0], [0, np.exp(1j*np.pi/4)]]), 'I': np.eye(2), } self.CNOT = np.array([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]]) def apply_gate(self, gate: np.ndarray, target: int): """Apply single-qubit gate to target qubit.""" op = np.array([[1]], dtype=complex) for i in range(self.n): if i == target: op = np.kron(op, gate) else: op = np.kron(op, np.eye(2)) self.state = op @ self.state def apply_cnot(self, control: int, target: int): """Apply CNOT gate.""" # Build CNOT for arbitrary control/target pair op = np.zeros((self.N, self.N), dtype=complex) for i in range(self.N): bits = [(i >> j) & 1 for j in range(self.n)] if bits[control] == 1: bits[target] ^= 1 j = sum(b << idx for idx, b in enumerate(bits)) op[j, i] = 1 self.state = op @ self.state def h(self, target: int): self.apply_gate(self.gates['H'], target) def x(self, target: int): self.apply_gate(self.gates['X'], target) def z(self, target: int): self.apply_gate(self.gates['Z'], target) def s(self, target: int): self.apply_gate(self.gates['S'], target) def t(self, target: int): self.apply_gate(self.gates['T'], target) def cnot(self, control: int, target: int): self.apply_cnot(control, target) def rx(self, theta: float, target: int): """Rotation around X axis.""" gate = np.array([ [np.cos(theta/2), -1j*np.sin(theta/2)], [-1j*np.sin(theta/2), np.cos(theta/2)] ]) self.apply_gate(gate, target) def ry(self, theta: float, target: int): """Rotation around Y axis.""" gate = np.array([ [np.cos(theta/2), -np.sin(theta/2)], [np.sin(theta/2), np.cos(theta/2)] ]) self.apply_gate(gate, target) def rz(self, phi: float, target: int): """Rotation around Z axis.""" gate = np.array([ [np.exp(-1j*phi/2), 0], [0, np.exp(1j*phi/2)] ]) self.apply_gate(gate, target) def measure(self, shots: int = 1024) -> dict: """Measure all qubits.""" probs = np.abs(self.state) ** 2 outcomes = np.random.choice(self.N, size=shots, p=probs.real) results = {} for outcome in outcomes: bitstring = format(outcome, f'0{self.n}b') results[bitstring] = results.get(bitstring, 0) + 1 return {k: v/shots for k, v in sorted(results.items(), key=lambda x: -x[1])} def get_statevector(self) -> np.ndarray: return self.state.copy() def print_state(self): print("Quantum State:") for i, amp in enumerate(self.state): if abs(amp) > 0.001: print(f" |{format(i, f'0{self.n}b')}⟩: {amp:.4f} " f"(|amp|² = {abs(amp)**2:.4f})") # Example: Create Bell state with circuit qc = QuantumCircuit(2) qc.h(0) qc.cnot(0, 1) qc.print_state() # |00⟩: 0.7071 (|amp|² = 0.5000) # |11⟩: 0.7071 (|amp|² = 0.5000) print("\nMeasurements:") results = qc.measure(shots=1000) for state, prob in results.items(): print(f" |{state}⟩: {prob:.3f}")