diff --git a/CompactMatrix.py b/CompactMatrix.py
new file mode 100644
index 0000000..f40296c
--- /dev/null
+++ b/CompactMatrix.py
@@ -0,0 +1,197 @@
+from sympy import sqrt, I, pi, exp, N, simplify
+from math import log2
+
+class Matrix:
+  def __init__(self, m=None):
+    if m is None:
+      self.state = [{}]
+    else:
+      #print(matrix)
+      self.shape = (len(m), len(m[0]))
+      #self.state = [{}] * self.shape[1]
+      self.state = []
+      for i in range(self.shape[1]):
+        self.state.append({})
+      #print(self.state[0])
+      for i in range(self.shape[0]):
+        for j in range(self.shape[1]):
+          #print(self.state, m[i][j])
+          if m[i][j] != 0:
+            self.state[j][i] = m[i][j]
+  def __repr__(self):
+    str = '['
+    for i in range(self.shape[1]):
+      str += '{\n'
+      for j in sorted(self.state[i]):
+        str += f'  |{j:0{int(log2(self.shape[0]))}b}>: {N(abs(self.state[i][j]) ** 2) * 100:8.4f}% [{self.state[i][j]}],\n'
+      str = str.strip(',\n')
+      str += '\n}, \n'
+    str = str.strip(', \n')
+    str += ']'
+    return str
+  def get_matrix(self):
+    m = [[0 for j in range(self.shape[1])] for i in range(self.shape[0])]
+    for i in range(self.shape[0]):
+      for j in range(self.shape[1]):
+        if i in self.state[j]:
+          m[i][j] = self.state[j][i]
+    return m
+  def kron(self, m):
+    #print(m.state)
+    #print(n.state)
+    #print(m.shape[0], m.shape[1], n.shape[0], n.shape[1])
+    C = zeros(self.shape[0] * m.shape[0], self.shape[1] * m.shape[1])
+    for i in range(self.shape[0]):
+      for j in self.state[i]: # range(m.shape[1]):
+        for k in range(m.shape[0]):
+          for l in m.state[k]: # range(n.shape[1]):
+            #print(i,j,k,l)
+            #if i in m.state[j] and k in n.state[l]:
+            C.state[i * m.shape[0] + k][j * m.shape[1] + l] = self.state[i][j] * m.state[k][l]
+            #print(C.state)
+    return C
+  def dot(self, m):
+    #print(self.shape[0], self.shape[1], m.shape[0], m.shape[1])
+    C = zeros(self.shape[0], m.shape[1])
+    for j in range(m.shape[1]):
+      for k in m.state[j]: #range(m.shape[1]):
+        for i in self.state[k]: # range(m.shape[0]):
+          #print(i,j,k)
+          #if i in m.state[k] and k in n.state[j]:
+          if not i in C.state[j]:
+            C.state[j][i] = 0
+          C.state[j][i] += self.state[k][i] * m.state[j][k]
+          C.state[j][i] = simplify(C.state[j][i])
+          if C.state[j][i] == 0:
+            del C.state[j][i]
+    #print(C.state)
+    return C
+
+class GateMatrix:
+  def __init__(self, m=None):
+    if m is None:
+      self.state = [eye(2)]
+    else:
+      #print(matrix)
+      self.shape = (len(m), 2, 2)
+      #self.state = [{}] * self.shape[1]
+      self.state = m
+      #for i in range(self.shape[1]):
+      #  self.state.append({})
+      #print(self.state[0])
+      #for i in range(self.shape[0]):
+      #  for j in range(self.shape[1]):
+          #print(self.state, m[i][j])
+      #    if m[i][j] != 0:
+      #      self.state[j][i] = m[i][j]
+  def __repr__(self):
+    s = str(self.state)
+  #  str = '['
+  #  for i in range(self.shape[1]):
+  #    str += '{\n'
+  #    for j in self.state[i]:
+  #      str += f'  |{j:0{int(log2(self.shape[0]))}b}>: {self.state[i][j]},\n'
+  #    str = str.strip(',\n')
+  #    str += '\n}, \n'
+  #  str = str.strip(', \n')
+  #  str += ']'
+    return s
+  def get_matrix(self):
+    m = [[0 for j in range(self.shape[1])] for i in range(self.shape[0])]
+    for i in range(self.shape[0]):
+      for j in range(self.shape[1]):
+        if i in self.state[j]:
+          m[i][j] = self.state[j][i]
+    return m
+  def kron(self, m):
+    #print(m.state)
+    #print(n.state)
+    #print(m.shape[0], m.shape[1], n.shape[0], n.shape[1])
+    C = zeros(self.shape[0] * m.shape[0], self.shape[1] * m.shape[1])
+    for i in range(self.shape[0]):
+      for j in self.state[i]: # range(m.shape[1]):
+        for k in range(m.shape[0]):
+          for l in m.state[k]: # range(n.shape[1]):
+            #print(i,j,k,l)
+            #if i in m.state[j] and k in n.state[l]:
+            C.state[i * m.shape[0] + k][j * m.shape[1] + l] = self.state[i][j] * m.state[k][l]
+            #print(C.state)
+    return C
+  def dot(self, m):
+    #print(self.shape[0], self.shape[1], m.shape[0], m.shape[1])
+    C = zeros(self.shape[0], m.shape[1])
+    for j in range(m.shape[1]):
+      for k in m.state[j]: #range(m.shape[1]):
+        for i in self.state[k]: # range(m.shape[0]):
+          #print(i,j,k)
+          #if i in m.state[k] and k in n.state[j]:
+          if not i in C.state[j]:
+            C.state[j][i] = 0
+          C.state[j][i] += self.state[k][i] * m.state[j][k]
+          if C.state[j][i] == 0:
+            del C.state[j][i]
+    #print(C.state)
+    return C
+
+def zeros(x, y):
+  m = Matrix()
+  m.shape = (x, y)
+  m.state = []
+  for i in range(y):
+    m.state.append({})
+  return m
+
+def eye(x):
+  m = Matrix()
+  m.shape = (x, x)
+  m.state = []
+  for i in range(x):
+    m.state.append({})
+    m.state[i][i] = 1
+  return m
+
+#def kron(m, n):
+#  #print(m.state)
+#  #print(n.state)
+#  #print(m.shape[0], m.shape[1], n.shape[0], n.shape[1])
+#  C = zeros(m.shape[0] * n.shape[0], m.shape[1] * n.shape[1])
+#  for i in range(m.shape[0]):
+#    for j in m.state[i]: # range(m.shape[1]):
+#      for k in range(n.shape[0]):
+#        for l in n.state[k]: # range(n.shape[1]):
+#          #print(i,j,k,l)
+#          #if i in m.state[j] and k in n.state[l]:
+#          C.state[i * n.shape[0] + k][j * n.shape[1] + l] = m.state[i][j] * n.state[k][l]
+#          #print(C.state)
+#  return C
+
+#def dot(m, n):
+#  #print(m.shape[0], m.shape[1], n.shape[0], n.shape[1])
+#  C = zeros(m.shape[0], n.shape[1])
+#  for j in range(n.shape[1]):
+#    for k in n.state[j]: #range(m.shape[1]):
+#      for i in m.state[k]: # range(m.shape[0]):
+#        #print(i,j,k)
+#        #if i in m.state[k] and k in n.state[j]:
+#        if not i in C.state[j]:
+#          C.state[j][i] = 0
+#        C.state[j][i] += m.state[k][i] * n.state[j][k]
+#        if C.state[j][i] == 0:
+#          del C.state[j][i]
+#  #print(C.state)
+#  return C
+
+#I = eye(2)
+#X = Matrix([[0,1],[1,0]])
+
+GATES = {
+  'I':   Matrix([[1, 0], [0, 1]]),
+  'X':   Matrix([[0, 1], [1, 0]]),
+  'Y':   Matrix([[0, -I], [I, 0]]), 
+  'Z':   Matrix([[1, 0], [0, -1]]),
+  'S':   Matrix([[1, 0], [0, I]]),
+  'SDG': Matrix([[1, 0], [0, -I]]),
+  'T':   Matrix([[1, 0], [0, exp(I * pi / 4)]]),
+  'TDG': Matrix([[1, 0], [0, exp(-I * pi / 4)]]),
+  'H':   Matrix([[sqrt(2) / 2, sqrt(2) / 2], [sqrt(2) / 2, -sqrt(2) / 2]]),
+}
\ No newline at end of file
diff --git a/quantum.py b/quantum.py
index 07ef9cc..99c5c29 100644
--- a/quantum.py
+++ b/quantum.py
@@ -6,6 +6,7 @@
 import numpy as np
 import sympy
 import sympy.physics.quantum
+import CompactMatrix as cm
 
 pi = sympy.pi
 I = sympy.I
@@ -16,11 +17,16 @@
 real = sympy.re
 imag = sympy.im
 
-eye = sympy.eye
-zeros = sympy.zeros
+#eye = sympy.eye
+#zeros = sympy.zeros
 #ones = sympy.ones
-kron = sympy.physics.quantum.TensorProduct
-Matrix = sympy.Matrix
+#kron = sympy.physics.quantum.TensorProduct
+#Matrix = sympy.Matrix
+eye = cm.eye
+zeros = cm.zeros
+#kron = cm.kron
+#dot = cm.dot
+Matrix = cm.Matrix
 
 start_time = time.time()
 last_time = start_time
@@ -45,7 +51,7 @@
   def __init__(self, num_qubits):
     self.num_qubits = num_qubits
     self.state = zeros(2**num_qubits, 1)
-    self.state[0] = 1 # Initialize to |0>
+    self.state.state[0][0] = 1 # Initialize to |0>
 
   def apply_gate(self, gate, target_qubits=None, control_qubits=None):
     if target_qubits is None:
@@ -64,11 +70,14 @@
       #  target_qubits = [target_qubits]
       #target_qubits = sorted(target_qubits)
     gate_matrix = eye(1)
+    #print("kron start")
     for i in range(self.num_qubits):
       if a[i]:
-        gate_matrix = kron(gate, gate_matrix)
+        gate_matrix = gate.kron(gate_matrix)
       else:
-        gate_matrix = kron(eye(2), gate_matrix)
+        gate_matrix = eye(2).kron(gate_matrix)
+      #print(f"round {i} done")
+    #print("kron done")
       #gate_matrix_b = kron(eye(int(2**(self.num_qubits - target_qubits[0] - 1))), kron(gate, eye(2**target_qubits[0])))
       #print(gate_matrix)
       #print(gate_matrix_b)
@@ -79,38 +88,44 @@
       #gate_matrix = kron(eye(int(2**(self.num_qubits - target_qubits[0] - 1))), kron(gate, eye(2**target_qubits[0])))
       #print_gate("init gate_matrix", gate_matrix)
       controlled_gate = zeros(2**self.num_qubits, 2**self.num_qubits)
+      #print("begin build gate")
       for i in range(2**(self.num_qubits)):
+        #print("q1 ", end='')
         if all((i // 2**q) % 2 == 1 for q in control_qubits):
-          for j, a in enumerate(gate_matrix[i,:]):
-            controlled_gate[i,j] = a
+          for j in gate_matrix.state[i]:
+            controlled_gate.state[i][j] = gate_matrix.state[i][j]
         else:
-          controlled_gate[i, i] = 1
+          controlled_gate.state[i][i] = 1
+      #print("gate complete")
       gate_matrix = controlled_gate #@ gate_matrix
       #print_gate("controlled_gate", controlled_gate)
-    #print_gate("gate_matrix", gate_matrix)
-    self.state = gate_matrix @ self.state
+    #print_gate("gate_matrix", gate_matrix.get_matrix())
+    self.state = gate_matrix.dot(self.state)
+    #print(self.state.state)
+    #print_gate("state", self.state.get_matrix())
 
   def measure(self):
-    probabilities = abs(self.state).applyfunc(lambda x: x ** 2)
+    probabilities = np.abs(self.state.get_matrix()) ** 2
+    #print(probabilities)
     probabilities = np.array(probabilities, dtype=np.float64).flatten()
     #probabilities /= sum(probabilities)
     outcome = np.random.choice(range(len(probabilities)), p=probabilities)
     self.state = zeros(*self.state.shape)
-    self.state[outcome] = 1.0
+    self.state.state[0][outcome] = 1
     #print(self.state)
     return bin(outcome)[2:].zfill(self.num_qubits)
 
   def measure_qubit(self, target_qubit):
-    num_states = len(self.state)
-    probabilities = abs(self.state)**2
+    num_states = 2 ** self.num_qubits #len(self.state)
+    probabilities = np.abs(self.state.get_matrix())**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]
     states_1 = [i for i in range(num_states) if (i >> target_qubit) % 2 == 1]
 
     # Calculate the total probabilities of measuring 0 or 1 for the target qubit
-    prob_0 = sum(probabilities[state] for state in states_0)
-    prob_1 = sum(probabilities[state] for state in states_1)
+    prob_0 = sum(probabilities[state] for state in states_0)[0]
+    prob_1 = sum(probabilities[state] for state in states_1)[0]
 
     # Randomly choose the measurement outcome based on the probabilities
     outcome = np.random.choice([0, 1], p=[prob_0, prob_1])
@@ -118,18 +133,24 @@
     # Update the state vector based on the measurement outcome
     if outcome == 0:
       for state in states_1:
-        self.state[state] = 0
+        self.state.state[0][state] = 0
       normalization_factor = sqrt(prob_0)
     else:
       for state in states_0:
-        self.state[state] = 0
+        self.state.state[0][state] = 0
       normalization_factor = sqrt(prob_1)
 
     # Normalize the state vector
-    self.state /= normalization_factor
+    #self.state /= normalization_factor
 
     return outcome
 
+  def set_qubit(self, target_qubit, state):
+    #state = state.strip('|><')
+    if state != '0' or state != '1' or state != '+' or state != '-' or state != 'i' or state != '-i':
+      raise ValueError(f"state {state} not recognized")
+    print(f"set qubit {target_qubit} to state |{state}> (not implemented)")
+
   def list_states(self, percentages = True):
     num_states = len(self.state)
     states_with_probabilities = {}
@@ -233,6 +254,8 @@
 #  'SWAP': SWAP,
 }
 
+GATES = cm.GATES
+
 def quantum_parse(quantum_code):
   instructions = quantum_code.strip().split("\n")
   parsed_instructions = []
@@ -276,6 +299,14 @@
         target_qubits.append(int(instruction[q + num_control_qubits + 1]))
       #target_qubit = int(instruction[num_control_qubits + 1])
       qreg.apply_gate(GATES[cmd[2]], target_qubits, control_qubits)
+    elif instruction[0] == 'QREG':
+      num_qubits = int(instruction[1])
+      qreg = QuantumRegister(num_qubits)
+      print(f"initialized qreg with {num_qubits} qubits")
+    elif instruction[0] == 'SET':
+      q = int(instruction[1])
+      s = instruction[2]
+      qreg.set_qubit(q, s)
     elif instruction[0] == 'Q':
       filename = instruction[1].lower()
       if not os.path.isfile(filename):
@@ -285,23 +316,23 @@
           code = file.read()
           code = code.format(*instruction[2:])
         qreg = quantum_interpreter(code, qreg)
-    elif instruction[0] == 'PEEK':
-      probabilities = qreg.probabilities()
-      states = qreg.list_states()
-      print("qubit probabilities:", *probabilities)
-      print("qreg probabilities:", states)
-    elif instruction[0] == 'PEEK_Q':
-      probabilities = qreg.probabilities(percentages=False)
-      probabilities_p = qreg.probabilities()
-      print("qubit probabilities:")
-      for q, q_p in zip(probabilities, probabilities_p):
-        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:+.2f}{2:+.2f}] ({3:6.2f}%)".format(s, real(states[s]), imag(states[s]), states_p[s] * 100))
+    #elif instruction[0] == 'PEEK':
+    #  probabilities = qreg.probabilities()
+    #  states = qreg.list_states()
+    #  print("qubit probabilities:", *probabilities)
+    #  print("qreg probabilities:", states)
+    #elif instruction[0] == 'PEEK_Q':
+    #  probabilities = qreg.probabilities(percentages=False)
+    #  probabilities_p = qreg.probabilities()
+    #  print("qubit probabilities:")
+    #  for q, q_p in zip(probabilities, probabilities_p):
+    #    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:+.2f}{2:+.2f}] ({3:6.2f}%)".format(s, real(states[s]), imag(states[s]), states_p[s] * 100))
     elif instruction[0] == 'STATE':
       print(qreg.state)
       #for i in range(2**qreg.num_qubits):
@@ -328,10 +359,12 @@
       print(*instruction[1:])
     elif instruction[0] == 'CONSOLE':
       while True:
-        print(f' Q{qreg.num_qubits}> ', end='')
+        print(f'{time.time() - start_time:8.3f}s: Q{qreg.num_qubits}> ', end='')
         instruction = input().upper()
         if instruction.startswith('BREAK'):
           break
+        elif instruction.startswith('EXIT'):
+          exit()
         qreg = quantum_interpreter(f'QREG 8\n{instruction}', qreg)
     elif os.path.isfile(instruction[0].lower() + ".q"):
       with open(instruction[0].lower() + ".q") as file:
@@ -342,7 +375,7 @@
       print("No command", instruction[0])
   return qreg
 
-quantum_code = "QREG 4\nCONSOLE"
+quantum_code = "QREG 4\nCONSOLE\nMEASURE"
 if len(sys.argv) > 1:
   filename = sys.argv[1]
   if os.path.isfile(filename):