# quantum.py3 computer simulation. This is meant to be mathematically
#             correct and obey the laws of quantum physics.
#             Yes, very picky about row and column, as Matlab would be.
#             Notation may vary, Open Source, modify to suit yourself.
#             Report typos or bugs to  squire@umbc.edu

import math # we implement our own KP Kronecker Product, TP Tensor Product,
            # DP Dot Product, etc,
import cmath # complex arithmetic is standard, this is for exp(i x)

# initialize constants, matrices, utility functions, etc.
#
def initialize():
  global ket0, ket1, ket00, ket01, ket10, ket11
  global ket000, ket001, ket010, ket011, ket100, ket101, ket110, ket111
  global bra0, bra1, bra00, bra01, bra10, bra11
  global H, PX, Not, PY, PZ, Cnot, CCnot, Swap, CSwap # matrices 
  #
  ket0 = [1,0] # |0>  technically a column vector, complex entries
  ket1 = [0,1] # |1>  technically a column vector
  bra0 = [1,0] # <0|  technically a row vector
  bra1 = [0,1] # <1|  technically a row vector
  ketx = [0,0] # temporaries
  ketxx = [0,0,0,0] # some may call these "register" 
  ketxxx = [0,0,0,0,0,0,0,0] # |x>|x>|x> = |xxx>
  ketxxxx = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
  reg2 = [0,0,0,0] # just concatenation of two qubits
  reg3 = [0,0,0,0,0,0,0,0] # just concatenation
  #
  ket00 = TP(ket0, ket0) # |0>|0>  =  |00> can be read as binary numbers
  ket01 = TP(ket0, ket1) # |0>|1>  =  |01>
  ket10 = TP(ket1, ket0) # |1>|0>  =  |10>
  ket11 = TP(ket1, ket1) # |1>|1>  =  |11>
  ket000 = TP3(ket0, ket0, ket0) # |0>|0>|0> = |000>
  ket001 = TP3(ket0, ket0, ket1) # |0>|0>|1> = |001>
  ket010 = TP3(ket0, ket1, ket0) # |0>|1>|0> = |010>
  ket011 = TP3(ket0, ket1, ket1) # |0>|1>|1> = |011>
  ket100 = TP3(ket1, ket0, ket0) # |1>|0>|0> = |100>
  ket101 = TP3(ket1, ket0, ket1) # |1>|0>|1> = |101>
  ket110 = TP3(ket1, ket1, ket0) # |1>|1>|0> = |110>
  ket111 = TP3(ket1, ket1, ket1) # |1>|1>|1> = |111>
  bra00 = TP(bra0, bra0) # <0|<0| = <00| same values, just row vector
  bra01 = TP(bra0, bra1) # <0|<1| = <01|
  bra10 = TP(bra1, bra0) # <1|<0| = <10|
  bra11 = TP(bra1, bra1) # <1|<1| = <11|
  # gates and unitary matrices
  He = 1.0/math.sqrt(2)
  H=[[He, He],[He, -He]] # Hadamard
  PX = [[0.0,1.0],[1.0,0.0]] # Pauli X "not" |0> to |1>, |1> to |0>
  Not = PX
  PY = [[0.0, 0.0-1j],[0.0+1j, 0.0]] # Pauli Y
  PZ = [[1, 0],[0, -1]] # Pauli Z
  Cnot = [[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]]
  CCnot = [[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,0,1,0,0,0,0,0],[0,0,0,1,0,0,0,0],[0,0,0,0,1,0,0,0],[0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,1],[0,0,0,0,0,0,1,0]]
  Swap = [[1,0,0,0],[0,0,1,0],[0,1,0,0],[0,0,0,1]]
  CSwap = [[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,0,1,0,0,0,0,0],[0,0,0,1,0,0,0,0],[0,0,0,0,1,0,0,0],[0,0,0,0,0,0,1,0],[0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,1]]
  #
  print("quantum function initialization complete")

# end initialize

def KP(A, B): # Kronecker Product of two matrix, returning matrix
  nar = len(A)     # n by m
  nac = len(A[0])
  nbr = len(B)     # p by q
  nbc = len(B[0])
  nyr = nar*nbr
  nyc = nac*nbc
  Y = [[0 for j in range(nyc)] for i in range(nyr)]
  for n in range(nar): # starts at zero
    for m in range(nac):
      for p in range(nbr):
        for q in range(nbc):
          Y[n*nbr+p][m*nbc+q] = A[n][m]*B[p][q]

        # end q

      # end p

    # end m

  # end n
  return Y

# end kp

def TP(A, B): # Tensor Product of two vectors returning vector, kets or bras
  na = len(A)
  nb = len(B)
  nc = na*nb
  Y = [0 for i in range(nc)]
  for i in range(na):
    for j in range(nb):
      Y[i*nb+j] = A[i]*B[j]

    # end j

  # end i
  return Y

# end TP

def TP3(A, B, C): # Tensor Product of three vectors returning vector, a state
  bc = TP(B, C)
  return TP(A, bc)

# end TP3

def TP4(A, B, C, D): # Tensor Product of four vectors returning vector
  bcd = TP3(B, C, D)
  return TP(A, bcd)

# end TP4

def DP(A, B): # Dot Product of two vectors returning a value
  na = len(A)
  nb = len(B)
  if na != nb :
    print("error, vectors different length in dot product",na," ",nb)

  # end if
  Y = 0
  for i in range(na):
    Y += A[i]*B[i]

  # end i
  return Y

# end DP

def negv(A): # negate all entries in a ket
  na = len(A)
  Y = [0 for i in range(na)]
  for i in range(na):
    Y[i] = -A[i]

  # end for
  return Y

# end negv

def negm(A): # negate all entries in a matrix
  nr = len(A)
  nc = len(A[0])
  Y = [[0 for j in range(nc)]for i in range(nr)]
  for i in range(nr):
    for j in range(nc):
      Y[i][j] = -A[i][j]

    # end j

  # end i
  return Y

# end negm

def addv(A, B): # add vector A to B
  na = len(A)
  nb = len(B)
  if na != nb :
    print("error, vectors different length in addv",na," ",nb)

  # end if
  Y = [0 for i in range(na)]
  for i in range(na):
    Y[i] = A[i] + B[i]

  # end i
  return Y

# end addv

def addm(A, B): # add matrix A to matrix B
  nar = len(A)
  nac = len(A[0])
  nbr = len(B)
  nbc = len(B[0])
  if nar != nbr :
    print("error, rows different in addm",nar," ", nbr)

  # end if
  if nac != nbc :
    print("error, columns different in addm",nac," ",nbc)

  # end if
  Y = [[0 for j in range(nac)]for i in range(nar)]
  for i in range(nar):
    for j in range(nac):
      Y[i][j] = A[i][j] + B[i][j]

    # end j

  # end i
  return Y

# end addm

def divv(A, d): # divide vector A by d
  na = len(A)
  Y = [0 for i in range(na)]
  for i in range(na):
    Y[i] = A[i]/d

  # end i
  return Y

# end divv

def divm(A, d): # divide  matrix A by d
  nr = len(A)
  nc = len(A[0])
  Y = [[0 for j in range(nc)]for i in range(nr)]
  for i in range(nr):
    for j in range(nc):
      Y[i][j] = A[i][j]/d

    # end j

  # end i
  return Y

# end divm

def entangle(A, B): # entangle ketA with KetB, caller must negate
  na = len(A)       # equivalent to add vectors with normalization
  nb = len(B)
  if na != nb :
    print("error, vectors different length in entangle",na," ",nb)

  # end if
  Y = [0 for i in range(na)]
  sum = 0
  for i in range(na):
    Y[i] = A[i]+B[i]
    sum += abs(Y[i]*Y[i])

  # end i
  for i in range(na):
    Y[i] = Y[i]/sum

  # end i
  return Y

# end entangle

def mulvm(V, M): # multiply row vector times matrix Y=V*M
                 # |# #| * |# # #| = |# # #|
                 #         |# # #|
  nv = len(V)
  nmr = len(M)
  nmc = len(M[0])
  Y = [0 for i in range(nmc)]
  if nv != nmr:
    print("error, mulvm vector, matrix size incompatible",nv," ",nmr)

  # end if
  for j in range(nmc):
    for i in range(nv):
      Y[j] += V[i]*M[i][j]

    # end j

  # end i
  return Y

# end mulvm

def mulmv(M, V): # multiply matrix times column vector Y=M*V
                 # | # # | * |#| = |#| 
                 # | # # |   |#|   |#|
                 # | # # |         |#|
  nv = len(V)
  nmr = len(M)     # number matrix rows
  nmc = len(M[0])  # number matrix columns
  Y = [0 for i in range(nmr)]
  if nv != nmc:
    print("error, mulmv matrix, vector size incompatible",nmc," ", nv)

  # end if
  for i in range(nmr):
    for j in range(nv):
      Y[i] += M[i][j]*V[j]

    # end j

  # end i
  return Y

# end mulmv

def mulmm(M1, M2): # multiply matrix times matrix Y=V*M
  nm1r = len(M1)
  nm1c = len(M1[0])
  nm2r = len(M2)
  nm2c = len(M2[0])
  Y = [[0 for j in range(nm2c)] for i in range(nm1r)]
  if nm2r != nm1c:
    print("error, mulmm matrix, matrix size incompatible")
    print(nm1r," ",nm1c," ",nm2r," ",nm2c)

  # end if
  for i in range(nm1r):
    for j in range(nm2c):
      for k in range(nm1c):
        Y[i][j] += M1[i][k]*M2[k][j]

      # end k

    # end j

  # end i
  return Y

# end mulmm

def transp(A): # transpose a matrix
  nr = len(A)
  nc = len(A[0])
  Y = [[0 for j in range(nr)]for i in range(nc)]
  for i in range(nr):
    for j in range(nc):
      Y[j][i] = A[i][j]

    # end j

  # end i
  return Y

# end transp

def Rk(k): # rotation gate [1,0],[0,e^i*2pi/(2^k)]
  Y = [[1,0],[0,0]]
  z = (1j * 2.0*math.pi) / pow(2,k)
  Y[1][1] = cmath.exp(z)
  return Y

# end RK

def CRk(k): # controlled rotation matrix |1 0 0 0 |
            #                            |0 1 0 0 |
            #                            |0 0 1 0 |
            #                            |0 0 0 e^i*2pi/(2^k)|
  Y = [[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,0]]
  z = (1j * 2.0*math.pi) / pow(2,k)
  Y[3][3] = cmath.exp(z)
  return Y

# end CRK

def measure(V1, V2):
  z = DP(V1,V2)
  # yuk! no cmath.pow(z, 2)
  return abs(cmath.exp(2.0*cmath.log(z)))

# end measure

def getket(i, K): # i counts 0, 1, 2, ...
  Y = [0,0]
  if i == 0:
    Y[0] = K[0]
    Y[1] = K[1]
    return Y

  # end id
  j = pow(2,i) # need int(math.pow(2.i))
  Y[0] = K[j]
  Y[1] = K[j+1]
  return Y

# end getket

def getbra(i, B): # i counts 0, 1, 2, ...
  Y = [0,0]
  if i == 0:
    Y[0] = B[0]
    Y[1] = B[1]
    return Y

  # end if
  j = pow(2,i) # 2^i
  Y[0] = B[j]
  Y[1] = B[j+1]
  return Y

# end getbra

def factor(K): # factor pure kets, assumes exactly one 1
  n = len(K)
  Y = []
  for i in range(n):
    if K[i] == 1:
      j = i

    # end if

  # end i
  while n>1.01:
    if math.fmod(j,2) == 0:
      Y.insert(0, [1,0])

    else:
      Y.insert(0, [0,1])

    # end if
    j = j/2
    n = n/2

  # end while
  return Y

# end factor

def factor_str(K): # factor pure kets, assumes exactly one 1
  n = len(K)
  Y = []
  for i in range(n):
    if K[i] == 1:
      j = i

    # end i

  # end i
  while n>1.01:
    if math.fmod(j,2) == 0:
      Y.insert(0, "|0> ")

    else:
      Y.insert(0, "|1> ")

    # end if
    j = j/2
    n = n/2

  # end while
  return Y

# end factor_str

def factor_strb(B): # factor pure bras, assumes exactly one 1
  n = len(B)
  Y = []
  for i in range(n):
    if B[i] == 1:
      j = i

    # enf if

  # end i
  while n>1.01:
    if math.fmod(j,2) == 0:
      Y.insert(0, "<0| ")

    else:
      Y.insert(0, "<1| ")

    # end if
    j = j/2
    n = n/2

  # end while
  return Y

# end factor_strb

# factor2  given [ac, ad, bc, bd] a^2+b^2=1, c^2+d^2=1, 
#  find [a,b,c,d] return kets |ab> and |ac>
#  if ad != 0, let v = ac/ad = c/d, let ac^2/ad^2 = v^2 = vv = c^2/d^2
#  d^2 = 1-c^2, thus vv = c^2/(1-c^2), solve for c^2, then c = sqrt(c^2)
#  vv*(1-c^2) = c^2, vv = c^2 + vv*c^2 = vv = (1+vv)*c^2
#  thus c^2 = vv/(1+vv), compute: v = ac/ad  c = sqrt((v*v)/(1.0+v*v))
#  d = sqrt(1-c*c) , a = ac/c , b = bc/c
def factor2(s):
  ac = s[0]
  ad = s[1]
  bc = s[2]
  bd = s[3]
  if abs(ad) > 0.01:
    v = ac/ad
    c = cmath.sqrt((v*v)/(1.0+v*v))
    d = cmath.sqrt(1-c*c)
    if abs(c) < 0.001:
      a = 1
      b = 0
    else:
      a = ac/c
      b = bc/c

    # end if
    return [a, b], [c, d]

  # end if
  if abs(bd) > 0.01:
    v = bc/bd
    c = cmath.sqrt((v*v)/(1.0+v*v))
    d = cmath.sqrt(1-c*c)
    if abs(c) < 0.001:
      a = 0
      b = 1
    else:
      a = ac/c
      b = bc/c

    # end if
    return [a, b], [c, d]

  # end if
  if abs(ac) > 0.01:
    v = ad/ac
    d = cmath.sqrt((v*v)/(1.0+v*v))
    c = cmath.sqrt(1-d*d)
    if abs(c) < 0.001:
      a = 0
      b = 1
    else:
      a = ac/c
      b = bc/c

    # end if
    return [a, b], [c, d]

  # end if
  if abs(bc) > 0.01:
    v = bd/bc
    d = cmath.sqrt((v*v)/(1.0+v*v))
    c = cmath.sqrt(1-d*d)
    if abs(c) < 0.001:
      a = 0
      b = 1
    else:
      a = ac/c
      b = bc/c

    # end if
    return [a, b], [c, d]

  # end if
  
# end factor2

# factor3  given [ace, acf, ade, adf, bce, bcf, bde, bdf] a state
#  a^2+b^2=1, c^2+d^2=1, e^2+f^2=1, find [a,b,c,d,e,f] to get
#  |ab> |cd> |ef>  three ket output
#  if acf != 0, let v = ace/acf = e/f, let ace^2/acf^2 = v^2 = vv = e^2/f^2
#  f^2 = 1-e^2, thus vv = e^2/(1-e^2), solve for e^2, then e = sqrt(e^2)
#  vv*(1-e^2) = e^2, vv = e^2 + vv*e^2 = vv = (1+vv)*e^2
#  thus e^2 = vv/(1+vv), compute: v = ace/acf  e = sqrt((v*v)/(1.0+v*v))
#  f = sqrt(1-e*e), solve for c and d similarly, or
#  now reduce problem to [ace/e, ade/e, bce/e, bde/e] use factor2 
#  This could be generalized to a recursive general routine.
def factor3(s):
  ace = s[0]
  acf = s[1]
  ade = s[2]
  adf = s[3]
  bce = s[4]
  bcf = s[5]
  bde = s[6]
  bdf = s[7]
  e   = 0
  if abs(acf) > 0.01:
    v = ace/acf
    e = cmath.sqrt((v*v)/(1.0+v*v))
    f = cmath.sqrt(1-e*e)

  elif abs(adf) > 0.01:
    v = ade/adf
    e = cmath.sqrt((v*v)/(1.0+v*v))
    f = cmath.sqrt(1-e*e)

  elif abs(bcf) > 0.01:
    v = bce/bcf
    e = cmath.sqrt((v*v)/(1.0+v*v))
    f = cmath.sqrt(1-e*e)

  elif abs(bdf) > 0.01:
    v = bde/bdf
    e = cmath.sqrt((v*v)/(1.0+v*v))
    f = cmath.sqrt(1-e*e)

  elif abs(ace) > 0.01:
    v = acf/ace
    f = cmath.sqrt((v*v)/(1.0+v*v))
    e = cmath.sqrt(1-f*f)

  elif abs(ade) > 0.01:
    v = adf/ade
    f = cmath.sqrt((v*v)/(1.0+v*v))
    e = cmath.sqrt(1-f*f)

  elif abs(bce) > 0.01:
    v = bcf/bce
    f = cmath.sqrt((v*v)/(1.0+v*v))
    e = cmath.sqrt(1-f*f)

  elif abs(bde) > 0.01:
    v = bdf/bde
    f = cmath.sqrt((v*v)/(1.0+v*v))
    e = cmath.sqrt(1-f*f)

  else:
    print("can not factor")
    return [[0,0],[0,0],[0,0]]

  # end if
  if abs(e) < 0.001:
    ketab, ketcd = factor2([acf/f, adf/f, bcf/f, bdf/f])
    return ketab, ketcd, [e, f]
  else:
    ketab, ketcd = factor2([ace/e, ade/e, bce/e, bde/e])
    return ketab, ketcd, [e, f]

  # end if

# end factor3

def test():
  # test
  m1 = [[1,2,3,4],[5,6,7,8]]
  print("m1 =",m1)
  m2 = [[10,11,12],[13,14,15],[16,17,18],[19,20,21]]
  print("m2 =",m2)
  v1 = [10,11]
  print("v1 =",v1)
  vx = mulvm(v1,m1)
  print("mulvm(v1,m1) = v1*m1 =",vx)
  v2 = [10,11,12,13]
  print("v2=",v2)
  vx = mulmv(m1,v2)
  print("mulmv(m1,v2) = m1*v2 =",vx)
  m3 = mulmm(m1,m2)
  print("m3 = m1*m2 =",m3)
  print(" ")
  print("test getket by number, 0, 1, 2, ...")
  print("getket(0,ket000) =",getket(0,ket000))
  print("getket(1,ket010) =",getket(1,ket010))
  print("getket(2,ket100) =",getket(2,ket100))
  print("getket(0,ket001) =",getket(0,ket001))
  print("getket(1,ket011) =",getket(1,ket011))
  print("getket(2,ket101) =",getket(2,ket101))
  print("getbra(1,bra01) =",getbra(1,bra01))
  print(" ")
  print("test factor on a pure kets, bras same results, added string printout")
  print("factor(ket0) =",factor(ket0))
  print("factor_str(ket0) =",factor_str(ket0))
  print("factor(ket00) =",factor(ket00))
  print(factor_str(ket00))
  print("factor(ket000) =",factor(ket000))
  print(factor_str(ket000))
  print("factor(ket001) =",factor(ket001))
  print(factor_str(ket001))
  print("factor(ket010) =",factor(ket010))
  print(factor_str(ket010))
  print("factor(ket011) =",factor(ket011))
  print(factor_str(ket011))
  print("factor(ket100) =",factor(ket100))
  print(factor_str(ket100))
  print("factor(ket101) =",factor(ket101))
  print(factor_str(ket101))
  print("factor(ket110) =",factor(ket110))
  print(factor_str(ket110))
  print("factor(ket111) =",factor(ket111))
  print(factor_str(ket111))
  print("factor_strb(bra01) = ",factor_strb(bra01))
  print(" ")
  # test factor2
  a = 0.1+0.2j
  b = cmath.sqrt(1.0-a*a)
  c = 0.4+0.3j
  d = cmath.sqrt(1.0-c*c)
  print("|ab> = ",a," ",b)
  print("|cd> = ",c," ",d)
  s = [a*c, a*d, b*c, b*d] # tensor product
  ketab, ketcd = factor2(s)
  print("factor2(",s,") = ")
  print("ketab |ab> = ",ketab)
  print("ketcd |cd> = ",ketcd)
  print(" ")
  # test factor3
  a = 0.1+0.2j
  b = cmath.sqrt(1.0-a*a)
  kab = [a, b]
  c = 0.4+0.3j
  d = cmath.sqrt(1.0-c*c)
  kcd = [c, d]
  e = 0.5+0.7j
  f = cmath.sqrt(1.0-e*e)
  kef = [e, f]
  print("|ab> = ",kab)
  print("|cd> = ",kcd)
  print("|ef> = ",kef)
  s = TP3(kab, kcd, kef);
  # s = [a*c*e, a*c*f, a*d*e, a*d*f, b*c*e, b*c*f, b*d*e, b*d*f] 
  ketab, ketcd, ketef = factor3(s)
  print("factor3(",s,") = ")
  print("ketab |ab> = ",ketab)
  print("ketcd |cd> = ",ketcd)
  print("ketef |ef> = ",ketef)
  print(" ")
  print("test factor3  [0.5, -0.5, -0.5, 0.0, -0.5, 0.0, 0.0, 0.0]")
  print("should be 0.5(|000> - |001> - |010> - |100>)")
  f3t = [0.5, -0.5, -0.5, 0.0, -0.5, 0.0, 0.0, 0.0]
  ketab, ketcd, ketef = factor3(f3t)
  print(ketab)
  print(ketcd)
  print(ketef)
  ketxxx = TP(ketab,TP(ketcd,ketef))
  print("rebuild ketab tp ketcd tp ketef")
  print(ketxxx)
  print(" ")
  print("KP(A,B) Kronecker Product 2,3 * 4,5 = big 8,15")
  A=[[1,2,3],[4,5,6]]
  print(A)
  B=[[1,1,1,1,2],[2,2,2,2,3],[3,3,3,3,4],[4,4,4,4,5]]
  print(B)
  Y = KP(A,B)
  print(Y)
  print(" ")

def main():
  print("quantum.py basic demonstration of many quantum operations") 

  initialize()

  print("Dirac notation and corresponding computer program names.")
  print(" ")
  print("A quantum bit, qubit, is represented in bra and ket notation")
  print("bra <x| is a two complex element row vector with length 1.0")
  print("    <x| = a<0| + b<1|  a and b complex numbers, abs(a^2+b^2)=1")
  print("ket |x> is a two complex element column vector with length 1.0")
  print("    |x> = a|0> + b|1>  a and b complex numbers, abs(a^2+b^2)=1")
  print(" ")
  print("A two qubit system will be on surface of Bloch unit sphere")
  print("  |y> = c|0> + d|1> then |x>|y> = |xy> = four entry vector =")
  print("  ac|00> + ad|01> + bc|10> + bd|11> written as a vector =")
  print("  | ac  ad  bc  bd | just four complex numbers, a 'state'")
  print("  mathematically, tensor product of  |a b| with |c d|")
  print("  The complex conjugate of |y> is used, thus d is -d in ad and bd")  
  print(" ")
  print("quantum physics computation must be reversible,")
  print("qubits can not be copied, measurement causes qubit to collapse,")
  print("Schrodinger's equation and Heisenberg Uncertainty Principle apply.")
  print("This will work on quantum superposition,")
  print("  vector |ac ad bc bd| can be factored, when |a b| is known")
  print("  into |a b| and |c d|")
  print("yet, will not work on quantum entanglement. Bell's theorem")
  print("a quantum entangled state can not be factored")
  print(" ")
  print("notation z^2 is  z raised to power 2")
  print("notation U^-1 is inverse of U, U^* is complex conjugate transpose of U")
  print("Unitary transforms, U: U^* U = U U^* = I, U^-1 = U ")
  print("                       eigenvalues of U are on complex unit circle")
  print(" a state times a unitary transform is a valid state")
  print(" ")

  print("Python definitions of values of ket0, ket1, bra0, bra1")
  print("ket0 = [1,0] = |0> =",)
  print(ket0)
  print("ket1 = [0,1] = |1> =",)
  print(ket1)
  print("bra0 = [1,0] = <0| =",)
  print(bra0)
  print("bra1 = [0,1] = <1| =",)
  print(bra1)
  print(" ")

  print("demonstrate use of functions and available constants")
  print("TP for tensor product")
  print("ket00 = TP(ket0, ket0) = |0>|0>  = |00> =",ket00)
  print("ket01 = TP(ket0, ket1) = |0>|1>  = |01> =",ket01)
  print("ket10 = TP(ket1, ket0) = |1>|0>  = |10> =",ket10)
  print("ket11 = TP(ket1, ket1) = |1>|1>  = |11> =",ket11)
  print(" ")
  print("ket000 = TP(ket0, TP(ket0, ket0)) = |0>|0>|0> = |000> =",ket000)
  print("ket000 = TP3(ket0, ket0, ket0) = |0>|0>|0> = |000> =",ket000)
  print("ket001 = TP3(ket0, ket0, ket1) = |0>|0>|1> = |001> =",ket001)
  print("ket010 = TP3(ket0, ket1, ket0) = |0>|1>|0> = |010> =",ket010)
  print("ket011 = TP3(ket0, ket1, ket1) = |0>|1>|1> = |011> =",ket011)
  print("ket100 = TP3(ket1, ket0, ket0) = |1>|0>|0> = |100> =",ket100)
  print("ket101 = TP3(ket1, ket0, ket1) = |1>|0>|1> = |101> =",ket101)
  print("ket110 = TP3(ket1, ket1, ket0) = |1>|1>|0> = |110> =",ket110)
  print("ket111 = TP3(ket1, ket1, ket1) = |1>|1>|1> = |111> =",ket111)
  print(" ")
  print("kets can be complex, j in Python, keti = |i> =",)
  keti = [0+1j,0]
  print(keti)
  print("complex keti01 = |i>|0>|1> =")
  keti01 = TP3(keti, ket0, ket1)
  print(keti01)
  keta, ketb, ketc = factor3(keti01)
  print("factor3(keti01)")
  print(keta,)
  print(ketb,)
  print(ketc)
  print("complex keti10 = |i>|1>|0> =")
  keti10 = TP3(keti, ket1, ket0)
  print(keti10)
  keta, ketb, ketc = factor3(keti10)
  print("factor3(keti10)")
  print(keta,)
  print(ketb,)
  print(ketc)
  print(" ")

  print("DP for dot product (actually a measurement returning a probability)")
  c = DP(bra0, ket0) # <0||0> = value 
  print("DP(bra0, ket0) = <0|0> = <0||0> = ",c)
  c = DP(bra0, ket1) # <0||1> = value 
  print("DP(bra0, ket1) = <0|1> = <0||1> = ",c)
  c = DP(bra1, ket0) # <1||0> = value 
  print("DP(bra1, ket0) = <1|0> = <1||0> = ",c)
  c = DP(bra1, ket1) # <1||1> = value 
  print("DP(bra1, ket1) = <1|1> = <1||1> = ",c)
  print(" ")
  print("demonstrate gates, common unitary matrices")
  print("Hadamard H=",H)
  print(" ")
  print("  H[0][0]=",H[0][0])
  print("  H[1][1]=",H[1][1])
  print(" ")

  ketx = mulvm(bra0, H)
  print(" mulvm(bra0, H) = <0| H =",ketx)
  ketx = mulvm(bra1, H)
  print(" mulvm(bra1, H) = <1| H =",ketx)
  ketx = mulmv(H, ket0)
  print(" mulmv(H, ket0) = H |0> =",ketx)
  ketx = mulmv(H, ket1)
  print(" mulmv(H, ket1) = H |1> =",ketx)
  print("  measure H |1>  in standard basis <1| bra1,")
  keth1 = mulmv(H, ket1)
  prob = math.pow(DP(bra1, ketx),2) # in general use measure function
  print("  probability of |1> is ( <1| (H |1>) )^2 = ",prob)
  prob = measure(bra1, keth1)
  print("  probability = measure(bra1, H |1>) = ( <1| (H |1>) )^2 = ",prob)
  print(" ")

  print("Pauli X = Not = PX = ",Not)
  print(" ")

  ketx = mulmv(PX, ket0)
  print(" PX |0> = |1> = ",ketx)
  ketx = mulmv(PX, ket1)
  print(" PX |1> = |0> = ",ketx)
  print(" ")
  ketx = mulvm(bra0, PX)
  print(" <0| PX = <1| = ",ketx)
  ketx = mulvm(bra1, PX)
  print(" <1| PX = <0| = ",ketx)
  print(" ")
  print("Pauli Y PY = ",PY)
  print(" ")
  #
  ketx = mulmv(PY, ket0)
  print(" PY |0> =",ketx)
  ketx = mulmv(PY, ket1)
  print(" PY |1> =",ketx)
  print(" ")
  print("Pauli Z = PZ =",PZ)
  print(" ")
  ketx = mulmv(PZ, ket0)
  print(" PZ |0> =",ketx)
  ketx = mulmv(PZ, ket1)
  print(" PZ |1> =",ketx)
  print(" ")
  #
  # test rotation 
  R8 = Rk(4)
  print("rotation 2/2^4 = 1/8 full circle = R8 = Rk(4) = ",R8)
  print("  mulmv(R8, ket0) = R8 |0> = ",mulmv(R8, ket0))
  print("  mulvm(R8, ket1) = R8 |1> = ",mulmv(R8, ket1))
  prob = measure(bra1, mulmv(R8, ket1))
  print("  probability of |1> is ( <1| (R8 |1>) )^2 = ",prob)
  print(" ")
  #
  # Euler's formula e^ix = cos(x)+i sin(x) "sqrt(-1) is j in Python"
  eix = cmath.exp(1j * math.pi/8.0) # e^i*pi*x
  R8 = [[1.0, 0.0],[0.0, eix]]
  print("  Euler's  pi/8  R8 = ",R8)
  print(" ")


  CR4 = CRk(4)
  print("controlled rotation CRk(4) = CR8 = ",CR4)
  print("  mulmv(CR4, ket00) = CR4 |0>|0> = ",mulmv(CR4, ket00))
  print("  mulmv(CR4, ket01) = CR4 |0>|1> = ",mulmv(CR4, ket01))
  print("  mulmv(CR4, ket10) = CR4 |1>|0> = ",mulmv(CR4, ket10))
  print("  mulmv(CR4, ket11) = CR4 |1>|1> = ",)
  ketxx =  mulmv(CR4, ket11)
  print(ketxx)
  print("  extract rotated ket, (Physics?) =",)
  ketx = getket(1, ketxx) # kets numbered from zero
  print(ketx)
  junk, ketx = factor2(ketxx)
  print("  must use 'getket' rather than 'factor2'")
  print("  junk, ketx = factor2(mulmv(CR4, Ket11))   ketx = ",ketx)
  print(" ")
  #
  print("two input controlled not gate Cnot =",Cnot)
  print("  pure Cnot, first ket controls, |1> inverts second ket")
  ketxx = mulmv(Cnot, ket00)
  print("  Cnot |00>  = |00> = ",ketxx)
  ketxx = mulmv(Cnot, ket01)
  print("  Cnot |01> = |01> = ",ketxx)
  ketxx = mulmv(Cnot, ket10)
  print("  Cnot |10> = |11> = ",ketxx)
  ketxx = mulmv(Cnot, ket11)
  print("  Cnot |11> = |10> = ",ketxx)
  print("  Cnot on a Hadamard ")
  print("  ( |1> ( |0>H )) Cnot =",)
  ketx = mulvm(ket0, H)
  ketxx = TP(ket1, ketx) # needed for input to raw Cnot
  ketxx = mulvm(ketxx, Cnot)
  print(ketxx)
  ketx = getket(1, ketxx) # kets numbered from zero
  print("  extract inverted ket, (Physics?) =",ketx)
  print("  ( |1> ( |1>H )) Cnot =",)
  ketx = mulvm(ket1, H)
  ketxx = TP(ket1, ketx) # needed for input to raw Cnot
  ketxx = mulvm(ketxx, Cnot)
  print(ketxx)
  ketx = getket(1, ketxx) # kets numbered from zero
  print("  extract inverted ket =",ketx)
  #
  # check out CCnot controlled Cnot
  print("three input controlled not gate CCnot = Toffoli =",CCnot)
  print("CCnot controlled Cnot for 8 cases")
  ketxxx = mulvm(ket000, CCnot)
  print("  mulvm(ket000, CCnot) = |000> CCnot = ",ketxxx)
  print("  factor(previous) =",factor(ketxxx))
  ketxxx = mulvm(ket001, CCnot)
  print("  mulvm(ket001, CCnot) = |001> CCnot = ",ketxxx)
  print("  factor(previous) =",factor(ketxxx))
  ketxxx = mulvm(ket010, CCnot)
  print("mulvm(ket010, CCnot) = |010> CCnot = ",ketxxx)
  print("  factor(previous) =",factor(ketxxx))
  ketxxx = mulvm(ket011, CCnot)
  print("  mulvm(ket011, CCnot) = |011> CCnot = ",ketxxx)
  print("  factor(previous) =",factor(ketxxx))
  ketxxx = mulvm(ket100, CCnot)
  print("  mulvm(ket100, CCnot) = |100> CCnot = ",ketxxx)
  print("  factor(previous) =",factor(ketxxx))
  ketxxx = mulvm(ket101, CCnot)
  print("  mulvm(ket101, CCnot) = |101> CCnot = ",ketxxx)
  print("  factor(previous) =",factor(ketxxx))
  ketxxx = mulvm(ket110, CCnot)
  print("  mulvm(ket110, CCnot) = |110> CCnot = ",ketxxx)
  print("  factor(previous) =",factor(ketxxx))
  ketxxx = mulvm(ket111, CCnot)
  print("  mulvm(ket111, CCnot) = |111> CCnot = ",ketxxx)
  print("  factor(previous) =",factor(ketxxx))
  print("  unknown how physics would factor this to get output")
  print("  last two bits are correct ket only if first two |1>")
  print("  I used an artificial factor function")
  print(" ")
  # check swap on pure and Hadamard
  print("two input gate Swap =",Swap)
  ketxx = mulvm(ket00, Swap)
  print("  mulvm(ket00, Swap) = |0>|0>Swap = |00> Swap =",ketxx)
  print("  factor(previous) =",factor(ketxx))
  print(factor_str(ketxx))
  ketxx = mulvm(ket01, Swap)
  print("  mulvm(ket01, Swap) = |0>|1>Swap = |01> Swap =",ketxx)
  print("  factor(previous) =",factor(ketxx))
  print(factor_str(ketxx))
  ketxx = mulvm(ket10, Swap)
  print("  mulvm(ket10, Swap) = |1>|0>Swap = |10> Swap =",ketxx)
  print("  factor(previous) =",factor(ketxx))
  print(factor_str(ketxx))
  ketxx = mulvm(ket11, Swap)
  print("  mulvm(ket11, Swap) = |1>|1>Swap = |11> Swap =",ketxx)
  print("  factor(previous) =",factor(ketxx))
  print(factor_str(ketxx))
  # check out CSwap
  print("three input gate, controlled swap, CSwap = Fredkin = ",CSwap)
  ketxxx = mulmv(CSwap, ket001)
  print("  mulmv(CSwap, ket001) = CSwap |001> = ",ketxxx)
  ketab, ketcd, ketef = factor3(ketxxx)
  print("  factor and last two kets not swapped = ",ketcd)
  print(ketef)
  ketxxx = mulmv(CSwap, ket101)
  print("  mulmv(CSwap, ket101) = CSwap |101> = ",ketxxx," ",ketab," ",ketcd)
  ketef = factor3(ketxxx)
  print("  factor and last two kets swapped = ",ketcd,"  ",ketef)
  print("test <s|Q T|s>")
  Q = PZ
  R = PX
  S = addm(negm(PZ), negm(PX))
  S = divm(S, math.sqrt(2))
  print("S = -PZ - PX = ",S)
  T = addm(PZ, negm(PX))
  T = divm(T, math.sqrt(2))
  print("T = PZ + PX = ",T)
  QS = KP(Q, S)
  print("QS = KP(Q, S) = ",QS)
  RS = KP(R, S)
  print("RS = KP(R, S) = ",RS)
  RT = KP(R, T)
  print("RT = KP(R, T) = ",RT)
  QT = KP(Q, T)
  print("QT = KP(Q, T) = ",QT)
  s0110 = addv(ket01, negv(ket10))
  s0110 = divv(s0110, math.sqrt(2))
  print("s0110 = (ket01-ket10)/sqrt(2) = ",s0110)
  junk = mulmv(QT, s0110)
  val = DP(s0110, junk)
  print("<s0110|QT|s0110> = ",val)
  s1001 = addv(ket10, negv(ket01))
  s1001 = divv(s1001, math.sqrt(2))
  print("s1001 = (ket10-ket01)/sqrt(2) = ",s1001)
  junk = mulmv(QT, s1001)
  val = DP(s1001, junk)
  print("<s1001|QT|s1001> = ",val)
  print("test utility functions")
  test()
  print("quantum.py3 finished")

# end main

if __name__ == '__main__':
  main()