diff --git a/quantum.py b/quantum.py index 65cdbf8..1c34477 100644 --- a/quantum.py +++ b/quantum.py @@ -2,13 +2,26 @@ import sys import time import re -import math import psutil import numpy as np -#import sympy +import sympy +import sympy.physics.quantum -#pi = sympy.pi -#e = sympy.e +pi = sympy.pi +#E = sympy.E +I = sympy.I +exp = sympy.exp +sqrt = sympy.sqrt +log = sympy.log +eye = sympy.eye +zeros = sympy.zeros +#ones = sympy.ones +abs = sympy.Abs +kron = sympy.physics.quantum.TensorProduct +real = sympy.re +imag = sympy.im + +Matrix = sympy.Matrix start_time = time.time() @@ -31,11 +44,11 @@ class QuantumRegister: def __init__(self, num_qubits): self.num_qubits = num_qubits - self.state = np.zeros(2**num_qubits, dtype=np.complex128) + self.state = zeros(2**num_qubits, 1) self.state[0] = 1.0 # Initialize to |0> def apply_gate(self, gate, target_qubit=None, control_qubits=None): - #print_gate("gate", gate) + #print("gate", gate) if control_qubits is None: # Uncontrolled gate if target_qubit is None: @@ -43,16 +56,16 @@ #if isinstance(target_qubits, int): # target_qubits = [target_qubits] #target_qubits = sorted(target_qubits) - gate_matrix = np.kron(np.eye(int(2**(self.num_qubits - target_qubit - 1)), dtype=np.complex128), np.kron(gate, np.eye(2**target_qubit, dtype=np.complex128))) + gate_matrix = kron(eye(int(2**(self.num_qubits - target_qubit - 1))), kron(gate, eye(2**target_qubit))) else: # Controlled gate if target_qubit is None or control_qubits is None: raise ValueError("Both target and control qubits must be specified for controlled gates.") if isinstance(control_qubits, int): control_qubits = [control_qubits] - gate_matrix = np.kron(np.eye(int(2**(self.num_qubits - target_qubit - 1)), dtype=np.complex128), np.kron(gate, np.eye(2**target_qubit, dtype=np.complex128))) + gate_matrix = kron(eye(int(2**(self.num_qubits - target_qubit - 1))), kron(gate, eye(2**target_qubit))) #print_gate("init gate_matrix", gate_matrix) - controlled_gate = np.zeros((2**self.num_qubits, 2**self.num_qubits), dtype=np.complex128) + controlled_gate = zeros(2**self.num_qubits, 2**self.num_qubits) for i in range(2**(self.num_qubits)): if all((i // 2**q) % 2 == 1 for q in control_qubits): for j, a in enumerate(gate_matrix[i,:]): @@ -63,100 +76,18 @@ #print_gate("controlled_gate", controlled_gate) #print_gate("gate_matrix", gate_matrix) self.state = gate_matrix @ self.state - def apple_gate(self, gate, target_qubit=None, control_qubit=None): - print_gate("gate", gate) - if gate.shape[0] == 2: #single qubit gate - if target_qubit is None: - raise ValueError("Target qubit must be specified for single qubit gate.") - gate_matrix = np.kron(np.eye(int(2**(self.num_qubits - target_qubit - 1)), dtype=np.complex128), np.kron(gate, np.eye(2**target_qubit, dtype=np.complex128))) - else: #two+ qubit gate - if target_qubit is None or control_qubit is None: - raise ValueError("Both target and control qubits must be specified for two qubit gate.") - # Apply controlled gate - #control_qubit, target_qubit = min(control_qubit, target_qubit), max(control_qubit, target_qubit) - #print("control: ", control_qubit) - #print("target: ", target_qubit) - gate_matrix = np.eye(2**(self.num_qubits - 2), dtype=np.complex128) - #gate_matrix = np.kron(gate, gate_matrix) - #gate_matrix = np.kron(gate_matrix, gate) - #gate_matrix = np.kron(gate, gate_matrix) if target_qubit > 0 else np.kron(gate_matrix, gate) - gate_matrix = np.kron(gate_matrix, gate) if target_qubit == self.num_qubits - 1 else np.kron(gate, gate_matrix) - #gate_matrix = np.kron(np.eye(int(2**(self.num_qubits - target_qubit - 1)), dtype=np.complex128), np.kron(gate[2:,2:], np.eye(2**target_qubit, dtype=np.complex128))) - print_gate("init gate_matrix", gate_matrix) - controlled_gate = np.zeros((2**self.num_qubits, 2**self.num_qubits), dtype=np.complex128) - for i in range(2**(self.num_qubits)): - if (i // 2**control_qubit) % 2 == 1: - controlled_gate[i, i ^ (1 << target_qubit)] = 1 - else: - controlled_gate[i, i] = 1 - gate_matrix = controlled_gate @ gate_matrix - #print_gate("controlled_gate", controlled_gate) - if self.num_qubits < 5: - print_gate("gate_matrix", gate_matrix) - self.state = gate_matrix @ self.state - #self.state /= np.linalg.norm(self.state) - def appled_gate(self, gate, target_qubits, control_qubits=None): - # Ensure target_qubits and gate dimensions match - if isinstance(target_qubits, int): - target_qubits = [target_qubits] - - if control_qubits is None: - if gate.shape[0] != 2**len(target_qubits): - raise ValueError("Gate dimensions do not match the number of target qubits.") - - # Ensure control_qubits and gate dimensions match - if control_qubits is not None: - if isinstance(control_qubits, int): - control_qubits = [control_qubits] - if gate.shape[0] != 2**len(control_qubits) + 2**len(target_qubits): - raise ValueError("Gate dimensions do not match the number of control qubits.") - - # Apply gate directly to target qubits if no control qubits are specified - if control_qubits is None: - # Construct the gate matrix for the target qubits - gate_matrix = gate - for target_qubit in reversed(target_qubits): - gate_matrix = np.kron(gate_matrix, np.eye(2**(target_qubit), dtype=np.complex128)) - gate_matrix = np.kron(np.eye(int(2**(self.num_qubits - target_qubit - 1)), dtype=np.complex128), gate_matrix) - - # Apply gate with control qubits - else: - # Ensure control_qubits are sorted in ascending order - control_qubits = sorted(control_qubits) - - # Construct the controlled gate matrix - control_mask = sum(1 << qubit for qubit in control_qubits) - control_states = [(state, state ^ control_mask) for state in range(2**self.num_qubits)] - controlled_gate_matrix = np.zeros((2**self.num_qubits, 2**self.num_qubits), dtype=np.complex128) - for control_state, target_state in control_states: - controlled_gate_matrix[control_state, target_state] = 1 - - # Construct the gate matrix for the target qubits - gate_matrix = gate - for target_qubit in reversed(target_qubits): - gate_matrix = np.kron(gate_matrix, np.eye(2**(target_qubit), dtype=np.complex128)) - gate_matrix = np.kron(np.eye(2**(self.num_qubits - target_qubit - 2), dtype=np.complex128), gate_matrix) - - # Apply the controlled gate - print(gate_matrix.shape, controlled_gate_matrix.shape, self.num_qubits) - gate_matrix = controlled_gate_matrix.dot(gate_matrix) - - # Apply the gate to the state vector - self.state = gate_matrix.dot(self.state) - - # Normalize the state vector - self.state /= np.linalg.norm(self.state) def measure(self): - probabilities = np.abs(self.state)**2 - probabilities /= np.sum(probabilities) + probabilities = abs(self.state).applyfunc(lambda x: x ** 2) + probabilities = np.array(probabilities, dtype=np.float64).flatten() + #probabilities /= sum(probabilities) outcome = np.random.choice(range(len(probabilities)), p=probabilities) - self.state = np.zeros_like(self.state) + self.state = zeros(*self.state.shape) self.state[outcome] = 1.0 #print(self.state) return bin(outcome)[2:].zfill(self.num_qubits) def measure_qubit(self, target_qubit): num_states = len(self.state) - probabilities = np.abs(self.state)**2 + probabilities = abs(self.state)**2 # Calculate the range of states that correspond to the target qubit being 0 or 1 states_0 = [i for i in range(num_states) if (i >> target_qubit) % 2 == 0] @@ -173,11 +104,11 @@ if outcome == 0: for state in states_1: self.state[state] = 0 - normalization_factor = np.sqrt(prob_0) + normalization_factor = sqrt(prob_0) else: for state in states_0: self.state[state] = 0 - normalization_factor = np.sqrt(prob_1) + normalization_factor = sqrt(prob_1) # Normalize the state vector self.state /= normalization_factor @@ -191,9 +122,9 @@ binary_state = bin(i)[2:].zfill(self.num_qubits) # Calculate the probability of the current state probability = self.state[i] - #probability = self.state[i] + #print(i, binary_state, probability) # Store the state and its probability in the dictionary - if probability.real > 1e-10 or probability.real < -1e-10 or probability.imag > 1e-10 or probability.imag < -1e-10: + if probability != 0 or real(probability) > 1e-10 or real(probability) < -1e-10 or imag(probability) > 1e-10 or imag(probability) < -1e-10: if percentages: probability = abs(probability)**2 states_with_probabilities[binary_state] = probability @@ -201,9 +132,9 @@ def probabilities(self, percentages = True): num_states = len(self.state) if percentages: - probabilities = np.zeros((self.num_qubits, 2)) # Initialize probabilities array + probabilities = zeros(self.num_qubits, 2) # Initialize probabilities array else: - probabilities = np.zeros((self.num_qubits, 2), np.complex128) + probabilities = zeros(self.num_qubits, 2) for i in range(num_states): # Convert the state index to binary representation binary_state = bin(i)[2:].zfill(self.num_qubits) @@ -226,51 +157,51 @@ for i, ket in enumerate(ga): if i > 0: print(', ', end='') - if abs(ket.real) > 0.01 or abs(ket.imag) > 0.01: + if abs(real(ket)) > 0.01 or abs(imag(ket)) > 0.01: print('\033[31m', end='') - print("{0.real:+2.0f}{0.imag:+2.0f}".format(ket), end='') - if abs(ket.real) > 0.01 or abs(ket.imag) > 0.01: + print("{0:+2.0f}{1:+2.0f}".format(real(ket), imag(ket)), end='') + if abs(real(ket)) > 0.01 or abs(imag(ket)) > 0.01: print('\033[0m', end='') print(']') -# Identity gate -I = np.array([[1, 0], - [0, 1]]) - # Pauli-X gate (bit-flip gate) -X = np.array([[0,1], - [1,0]], dtype=np.complex128) +X = Matrix([[0, 1], + [1, 0]]) # Pauli-Y gate -Y = np.array([[0,-1j], - [1j,0]], dtype=np.complex128) +Y = Matrix([[0, -I], + [I, 0]]) # Pauli-Z gate -Z = np.array([[1,0], - [0,-1]], dtype=np.complex128) +Z = Matrix([[1,0], + [0,-1]]) # Hadamard gate -H = np.array([[math.sqrt(0.5),math.sqrt(0.5)], - [math.sqrt(0.5),-math.sqrt(0.5)]], dtype=np.complex128) +H = Matrix([[sqrt(2) / 2,sqrt(2) / 2], + [sqrt(2) / 2,-sqrt(2) / 2]]) # Phase gate (S, P) -S = np.array([[1, 0], - [0, 1j]], dtype=np.complex128) +S = Matrix([[1, 0], + [0, I]]) -SD = np.array([[1, 0], - [0, -1j]], dtype=np.complex128) +SD = Matrix([[1, 0], + [0, -I]]) # pi / 8 gate (T) -T = np.array([[1, 0], - [0, np.exp(1j * np.pi / 4)]], dtype=np.complex128) +T = Matrix([[1, 0], + [0, exp(I * pi / 4)]]) -TD = np.array([[1, 0], - [0, np.exp(-1j * np.pi / 4)]], dtype=np.complex128) +TD = Matrix([[1, 0], + [0, exp(-I * pi / 4)]]) -SWAP = np.array([[1, 0, 0, 0], +SWAP = Matrix([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], - [0, 0, 0, 1]], dtype=np.complex128) + [0, 0, 0, 1]]) + +# Identity gate +EYE = Matrix([[1, 0], + [0, 1]]) GATES = { 'X': X, @@ -281,20 +212,10 @@ 'SD': SD, 'T': T, 'TD': TD, + 'EYE': EYE, # 'SWAP': SWAP, } -def build_controlled_gate(gate, num_control_gates): - # Calculate the size of the controlled gate matrix - num_qubits = num_control_gates + int(math.log(gate.shape[0], 2)) - gate_size = 2 ** num_qubits - controlled_gate = np.eye(gate_size, dtype=np.complex128) - # Calculate the position to insert the gate - target_position = 2 ** num_qubits - gate.shape[0] - # Insert the gate into the controlled gate matrix - controlled_gate[target_position:target_position + gate.shape[0], target_position:target_position + gate.shape[0]] = gate - return controlled_gate - def quantum_parse(quantum_code): instructions = quantum_code.strip().split("\n") parsed_instructions = [] @@ -305,15 +226,15 @@ parsed_instructions.append(parsed_instruction) return parsed_instructions -#@profile def quantum_interpreter(quantum_code, qreg=None): parsed_instructions = quantum_parse(quantum_code) if qreg is None: num_qubits = int(parsed_instructions[0][1]) qreg = QuantumRegister(num_qubits) + print(f"initialized qreg with {num_qubits} qubits") parsed_instructions = parsed_instructions[1:] for instruction in parsed_instructions: - #print(*instruction) + #print(f"{time.time() - start_time:8.3f}s:", *instruction) cmd = re.findall(r'(C*)(\d*)(\w+)', instruction[0])[0] if cmd[2] in GATES.keys(): control_qubits = [] @@ -344,30 +265,27 @@ probabilities_p = qreg.probabilities() print("qubit probabilities:") for q, q_p in zip(probabilities, probabilities_p): - print("[{0.real:+.2f}{0.imag:+.2f}j, {1.real:+.2f}{1.imag:+.2f}j] ({2:3.0f}%, {3:3.0f}%)".format(q[0], q[1], q_p[0] * 100, q_p[1] * 100)) + print("[{0:+.2f}{1:+.2f}j, {2:+.2f}{3:+.2f}j] ({4:3.0f}%, {5:3.0f}%)".format(real(q[0]), imag(q[0]), real(q[1]), imag(q[1]), q_p[0] * 100, q_p[1] * 100)) elif instruction[0] == 'PEEK_S': states = qreg.list_states(percentages=False) states_p = qreg.list_states() print("qreg probabilities:") for s in states: - print("|{0}>: [{1.real:+.2f}{1.imag:+.2f}] ({2:6.2f}%)".format(s, states[s], states_p[s] * 100)) + print("|{0}>: [{1:+.2f}{2:+.2f}] ({3:6.2f}%)".format(s, real(states[s]), imag(states[s]), states_p[s] * 100)) elif instruction[0] == 'STATE': - for i in range(2**qreg.num_qubits): - print('|{:08b}>'.format(i), end='') - print() - print('[', end='') - for i, s in enumerate(qreg.state): - #print(s) - if i > 0: - print(', ', end='') - print("{0.real:+4.1f}{0.imag:+4.1f}".format(s), end='') - print(']') + print(qreg.state) + #for i in range(2**qreg.num_qubits): + # print('|{:08b}>'.format(i), end='') + #print() + #print('[', end='') + #for i, s in enumerate(qreg.state): + # #print(s) + # if i > 0: + # print(', ', end='') + # print("{0:+4.1f}{1:+4.1f}".format(real(s), imag(s)), end='') + #print(']') elif instruction[0] == 'MEASURE': - #probabilities = qreg.probabilities() - #states = qreg.list_states() result = qreg.measure() - #print("qubit probabilities:", *probabilities) - #print("qreg probabilities:", states) print("measured state: |{}>".format(result)) elif instruction[0] == 'MEASURE_Q': result = qreg.measure_qubit(int(instruction[1]))