text: from __future__ import annotations
  import numpy as np
  from dataclasses import dataclass
  import util
  import time
  #import myModule
  @dataclass
  class CHState:
      N : int # number of qubits
      A : np.ndarray # NxN matrix of bytes (we are using as bits) partly determines U_C
      B : np.ndarray # NxN matrix of bytes (we are using as bits) partly determines U_C
      C : np.ndarray # NxN matrix of bytes (we are using as bits) partly determines U_C
      g : np.ndarray # gamma is in (Z / Z^4)^N
      v : np.ndarray # array of N bytes (which we are using as bits) determining U_H 
      s : np.ndarray #array of N bytes (which we are using as bits) - the initial state 
      phase : complex #initial phase
  
      @classmethod
      def basis(cls, N:int = None, s=None) -> CHState:
          """
          Return a computational basis state defined by the bitstring s
          """
          if N == None and s is None:
              #given no input we assume a single qubit in state |0>
              return cls(N=1,
                         A=np.eye(1, dtype=np.uint8),
                         B=np.eye(1, dtype=np.uint8),
                         C=np.zeros((1,1), dtype=np.uint8),
                         g=np.zeros(1, dtype=np.uint8),
                         v=np.zeros(1, dtype=np.uint8),
                         s=np.zeros(1, dtype=np.uint8),
                         phase = complex(1,0)
              )
          elif N != None and s is None:
              #we get given a number of qubits but no other information so return the state |0,0....0>
              return cls(N=N,
                         A=np.eye(N, dtype=np.uint8),
                         B=np.eye(N, dtype=np.uint8),
                         C=np.zeros((N,N), dtype=np.uint8),
                         g=np.zeros(N, dtype=np.uint8),
                         v=np.zeros(N, dtype=np.uint8),
                         s=np.zeros(N, dtype=np.uint8),
                         phase = complex(1,0)
              )
          elif N == None and not s is None:
              #we get given some bitstring so we return that computational basis state
              N = len(s)
              s = np.array(s, dtype=np.uint8) #we accept lists etc, but convert them to np arrays
              return cls(N=N,
                         A=np.eye(N, dtype=np.uint8),
                         B=np.eye(N, dtype=np.uint8),
                         C=np.zeros((N,N), dtype=np.uint8),
                         g=np.zeros(N, dtype=np.uint8),
                         v=np.zeros(N, dtype=np.uint8),
                         s=s,
                         phase = complex(1,0)
              )
          else:
              #both N and s are not none
              # if N <= len(s) we truncate s to length s and proceed as before
              # if N > len(s) we extend s by adding zeros at the end and proceed as before
              s = np.array(s, dtype=np.uint8)
              if N <= len(s):
                  return CHState.basis(N=None, s = s[:N])
              else:
                  return CHState.basis(N=None, s = np.concatenate((s, np.zeros(N-len(s), dtype=np.uint8))))  
  
      @property
      def F(self):
          return self.A
      @F.setter
      def F(self, mat):
          self.A = mat
      @property
      def G(self):
          return self.B
      @G.setter
      def G(self, mat):
          self.B = mat
      @property
      def M(self):
          return self.C
      @M.setter
      def M(self, mat):
          self.C = mat
      @property
      def gamma(self):
          return self.g
      @gamma.setter
      def gamma(self, mat):
          self.g = mat
      @property
      def w(self):
          return self.phase
      @w.setter
      def w(self, c):
          self.phase = c
  
      def __or__(self, other : CliffordGate):
          return other.applyCH(self)
  
      def _rowToStr(row):
          return "".join(map(str,row))
  
      def tab(self):
          """
          pretty "to string" method for small qubit numbers
          prints blocks F G M gamma v s
          with headings to indicate which is which
          """
          
          s  = str(self.N) + " "
          qubitNumberStrLen = len(s)
          matrix_width = self.N
          half_matrix_width = self.N//2
          s = "N" + " "*(qubitNumberStrLen -1 + half_matrix_width) + "F" + " "*matrix_width + "G" + " "*matrix_width + "M" + " "*(matrix_width-half_matrix_width) + "g v s w\n" + s
          
          for i, (Fr, Gr, Mr, gr, vr, sr) in enumerate(zip(self.F, self.G, self.M, self.g, self.v, self.s)):
              if i != 0:
                  s += " "*qubitNumberStrLen
              s += CHState._rowToStr(Fr) + " " + CHState._rowToStr(Gr) + " " + CHState._rowToStr(Mr) + " " + str(gr) + " " + str(vr) + " " + str(sr)
  
              if i == 0:
                  s += " " + str(self.phase)
              s += "\n"
          return s
  
          
      def __str__(self):
          """
          pretty "to string" method for small qubit numbers
          prints blocks F G M gamma v s
          """
          qubitNumberStrLen = None
          s = ""
          
          for i, (Fr, Gr, Mr, gr, vr, sr) in enumerate(zip(self.F, self.G, self.M, self.g, self.v, self.s)):
              if i == 0:
                  s = str(self.N) + " "
                  qubitNumberStrLen = len(s)
              if i != 0:
                  s += " "*qubitNumberStrLen
              s += CHState._rowToStr(Fr) + " " + CHState._rowToStr(Gr) + " " + CHState._rowToStr(Mr) + " " + str(gr) + " " + str(vr) + " " + str(sr)
  
              if i == 0:
                  s += " " + str(self.phase)
              s += "\n"
          return s
  
      def delete_qubit(self, k):
          mask = np.ones(self.N,dtype=bool)
          mask[k] = False
          mask2d = np.outer(mask,mask)
          
          return CHState(self.N-1, self.A[mask2d].reshape(self.N-1, self.N-1), self.B[mask2d].reshape(self.N-1, self.N-1), self.C[mask2d].reshape(self.N-1, self.N-1), self.g[mask], self.v[mask],self.s[mask], self.phase)
          
  
  
      def __sub__(self,other):
          return CHState(self.N,
                           (self.A - other.A)%np.uint8(2),
                           (self.B - other.B)%np.uint8(2),
                           (self.C - other.C)%np.uint8(2),
                           (self.g - other.g)%np.uint8(4),
                           (self.v - other.v)%np.uint8(2),
                           (self.s - other.s)%np.uint8(2),
                           (self.phase / other.phase))
  
      def __add__(self,other):
          return CHState(self.N,
                           (self.A + other.A)%np.uint8(2),
                           (self.B + other.B)%np.uint8(2),
                           (self.C + other.C)%np.uint8(2),
                           (self.g + other.g)%np.uint8(4),
                           (self.v + other.v)%np.uint8(2),
                           (self.s + other.s)%np.uint8(2),
                           (self.phase * other.phase))
  
      def equatorial_inner_product(self, A):
          """
          Given an equatorial state |phi_A> defined by a symmetric binary matrix A
          compute
          <|phi_A | self >
          """
          J = np.int64((self.M @ self.F.T) % np.uint8(2))
          J[np.diag_indices_from(J)] = self.g
          K = (self.G.T @ (A + J) @ self.G)
          prefactor = (2**(-(self.N + self.v.sum())/2)) * ((1j)**(self.s @ K @ self.s)) * ((-1)**(self.s @ self.v))
          B = (K + 2*np.diag(self.s + self.s @ K))[self.v == 1][:,self.v == 1] % np.uint8(4)
  
          #M = np.triu(B) % np.uint8(2) #upper triangular part including diagonal
          #M[np.diag_indices_from(M)] = np.uint8(0)
          #print("M = ", M)
          K = B[np.diag_indices_from(B)] % np.uint8(2)
  
          L = ((B[np.diag_indices_from(B)] - K) // np.uint8(2))  # the // forces integer division and makes sure the dtype remains uint8
  
          newL = np.append(L,0)
  
          newM = np.triu(B +np.outer(K,K)) %np.uint8(2)
          newM[np.diag_indices_from(newM)] = np.uint8(0)
          
          newM = np.concatenate((newM, np.array([K],dtype=np.uint8)), axis=0)
          newM = np.concatenate((newM, np.array([[0]*newM.shape[0]],dtype=np.uint8).T) ,  axis=1)
          newM = np.uint8(newM)
          
          newL = np.uint8(newL)
  
          # m2 = np.copy(newM)
          # l2 = np.copy(newL)
          # m3 = np.copy(newM)
          # l3 = np.copy(newL)
          
          #re = util.slowZ2ExponentialSum(newM, newL) / 2
          #newL[-1] = 1
          #im = util.slowZ2ExponentialSum(newM, newL) / 2
          
          re  = util.z2ExponentialSum(m2, l2)
          l2[-1] = 1
          im  = util.z2ExponentialSum(m2, l2)   
          re /= 2
          im /= 2
  
          #delta = time.monotonic()
          #re, im = util.z2DoubleExponentialSum2(newM,newL)
          
          
          #re, im = myModule.exponential_sum(newM,newL)
          #re /= 2
          #im /= 2
  
          #delta = time.monotonic() - delta
          # re2, im2 = util.z2DoubleExponentialSum2(m2,l2)
          # re2 /=2
          # im2 /=2
  
          # re3  = util.z2ExponentialSum(m3, l3)
          # l3[-1] = 1
          # im3  = util.z2ExponentialSum(m3, l3)
          
          # re3 /= 2
          # im3 /= 2
  
          # good = abs(re-re2) + abs(im-im2) < 1e-6
          # if not good:
          #     print(re, im)
          #     print(re2, im2)
          #     print(re3, im3)
          #     print()
  
          
          return self.phase.conjugate()*prefactor*complex(re, im )
  
          
          
      def __eq__(self, other):
          if not isinstance(other, CHState):
              return False
  
          if (self.F == other.F).all() and \
             (self.G == other.G).all() and \
             (self.M == other.M).all() and \
             (self.g == other.g).all() and \
             (self.v == other.v).all() and \
             (self.s == other.s).all() and abs(self.phase - other.phase) < 1e-10:
              return True
             
              
  