#             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)

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> =",
  print ket00
  print "ket01 = TP(ket0, ket1) = |0>|1>  = |01> =",
  print ket01
  print "ket10 = TP(ket1, ket0) = |1>|0>  = |10> =",
  print ket10
  print "ket11 = TP(ket1, ket1) = |1>|1>  = |11> =",
  print ket11
  print " "
  print "ket000 = TP(ket0, TP(ket0, ket0)) = |0>|0>|0> = |000> =",
  print ket000
  print "ket000 = TP3(ket0, ket0, ket0) = |0>|0>|0> = |000> =",
  print ket000
  print "ket001 = TP3(ket0, ket0, ket1) = |0>|0>|1> = |001> =",
  print ket001
  print "ket010 = TP3(ket0, ket1, ket0) = |0>|1>|0> = |010> =",
  print ket010
  print "ket011 = TP3(ket0, ket1, ket1) = |0>|1>|1> = |011> =",
  print ket011
  print "ket100 = TP3(ket1, ket0, ket0) = |1>|0>|0> = |100> =",
  print ket100
  print "ket101 = TP3(ket1, ket0, ket1) = |1>|0>|1> = |101> =",
  print ket101
  print "ket110 = TP3(ket1, ket1, ket0) = |1>|1>|0> = |110> =",
  print ket110
  print "ket111 = TP3(ket1, ket1, ket1) = |1>|1>|1> = |111> =",
  print 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> = ",
  print c
  c = DP(bra0, ket1) # <0||1> = value 
  print "DP(bra0, ket1) = <0|1> = <0||1> = ",
  print c
  c = DP(bra1, ket0) # <1||0> = value 
  print "DP(bra1, ket0) = <1|0> = <1||0> = ",
  print c
  c = DP(bra1, ket1) # <1||1> = value 
  print "DP(bra1, ket1) = <1|1> = <1||1> = ",
  print c
  print " "

  print "demonstrate gates, common unitary matrices"

  print "Hadamard H=",
  print H
  print " "
  print "  H[0][0]=",
  print H[0][0]
  print "  H[1][1]=",
  print H[1][1]
  print " "

  ketx = mulvm(bra0, H)
  print " mulvm(bra0, H) = <0| H =",
  print ketx
  ketx = mulvm(bra1, H)
  print " mulvm(bra1, H) = <1| H =",
  print ketx
  ketx = mulmv(H, ket0)
  print " mulmv(H, ket0) = H |0> =",
  print ketx
  ketx = mulmv(H, ket1)
  print " mulmv(H, ket1) = H |1> =",
  print ketx
  print " "

  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 = ",
  print prob
  prob = measure(bra1, keth1)
  print "  probability = measure(bra1, H |1>) = ( <1| (H |1>) )^2 = ",
  print prob
  print " "

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

  ketx = mulmv(PX, ket0)
  print " PX |0> = |1> = ",
  print ketx
  ketx = mulmv(PX, ket1)
  print " PX |1> = |0> = ",
  print ketx
  print " "
  ketx = mulvm(bra0, PX)
  print " <0| PX = <1| = ",
  print ketx
  ketx = mulvm(bra1, PX)
  print " <1| PX = <0| = ",
  print ketx
  print " "

  print "Pauli Y PY = ", 
  print PY
  print " "

  ketx = mulmv(PY, ket0)
  print " PY |0> =",
  print ketx
  ketx = mulmv(PY, ket1)
  print " PY |1> =",
  print ketx
  print " "

  print "Pauli Z = PZ =",
  print PZ
  print " "

  ketx = mulmv(PZ, ket0)
  print " PZ |0> =",
  print ketx
  ketx = mulmv(PZ, ket1)
  print " PZ |1> =",
  print ketx
  print " "

  # test rotation 
  R8 = Rk(4)
  print "rotation 2/2^4 = 1/8 full circle = R8 = Rk(4) = ",
  print R8
  print " "
  print "  mulmv(R8, ket0) = R8 |0> = ",
  print mulmv(R8, ket0)
  print " "
  print "  mulvm(R8, ket1) = R8 |1> = ",
  print mulmv(R8, ket1)
  print " "
  prob = measure(bra1, mulmv(R8, ket1))
  print "  probability of |1> is ( <1| (R8 |1>) )^2 = ",
  print 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 = ",
  print R8
  print " "
  print " "


  CR4 = CRk(4)
  print "controlled rotation CRk(4) = CR8 = ",
  print CR4
  print "  mulmv(CR4, ket00) = CR4 |0>|0> = ",
  print mulmv(CR4, ket00)
  print "  mulmv(CR4, ket01) = CR4 |0>|1> = ",
  print mulmv(CR4, ket01)
  print "  mulmv(CR4, ket10) = CR4 |1>|0> = ",
  print 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 = ",
  print ketx
  print " "
  print " "


  print "two input controlled not gate Cnot =",
  print Cnot
  print " "
  print "  pure Cnot, first ket controls, |1> inverts second ket"
  print "  Cnot |00>  = |00> = ",
  ketxx = mulmv(Cnot, ket00)
  print ketxx
  print " "

  print "  Cnot |01> = |01> = ",
  ketxx = mulmv(Cnot, ket01)
  print ketxx
  print " "

  print "  Cnot |10> = |11> = ",
  ketxx = mulmv(Cnot, ket10)
  print ketxx
  print " "

  print "  Cnot |11> = |10> = ",
  ketxx = mulmv(Cnot, ket11)
  print ketxx
  print " "

  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
  print "  extract inverted ket, (Physics?) =",
  ketx = getket(1, ketxx) # kets numbered from zero
  print ketx
  print " "

  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
  print "  extract inverted ket =",
  ketx = getket(1, ketxx) # kets numbered from zero
  print ketx
  print " "


  # check out CCnot controlled Cnot
  print "three input controlled not gate CCnot = Toffoli ="
  print CCnot
  print " "
  print "CCnot controlled Cnot for 8 cases"
  ketxxx = mulvm(ket000, CCnot)
  print "  mulvm(ket000, CCnot) = |000> CCnot = ",
  print ketxxx
  print "  factor(previous) =",
  print factor(ketxxx)
  ketxxx = mulvm(ket001, CCnot)
  print "  mulvm(ket001, CCnot) = |001> CCnot = ",
  print ketxxx
  print "  factor(previous) =",
  print factor(ketxxx)
  ketxxx = mulvm(ket010, CCnot)
  print "mulvm(ket010, CCnot) = |010> CCnot = ",
  print ketxxx
  print "  factor(previous) =",
  print factor(ketxxx)
  ketxxx = mulvm(ket011, CCnot)
  print "  mulvm(ket011, CCnot) = |011> CCnot = ",
  print ketxxx
  print "  factor(previous) =",
  print factor(ketxxx)
  ketxxx = mulvm(ket100, CCnot)
  print "  mulvm(ket100, CCnot) = |100> CCnot = ",
  print ketxxx
  print "  factor(previous) =",
  print factor(ketxxx)
  ketxxx = mulvm(ket101, CCnot)
  print "  mulvm(ket101, CCnot) = |101> CCnot = ",
  print ketxxx
  print "  factor(previous) =",
  print factor(ketxxx)
  ketxxx = mulvm(ket110, CCnot)
  print "  mulvm(ket110, CCnot) = |110> CCnot = ",
  print ketxxx
  print "  factor(previous) =",
  print factor(ketxxx)
  ketxxx = mulvm(ket111, CCnot)
  print "  mulvm(ket111, CCnot) = |111> CCnot = ",
  print ketxxx
  print "  factor(previous) =",
  print factor(ketxxx)
  print " "
  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 =",
  print Swap
  print " "

  ketxx = mulvm(ket00, Swap)
  print "  mulvm(ket00, Swap) = |0>|0>Swap = |00> Swap =",
  print ketxx
  print "  factor(previous) =",
  print factor(ketxx)
  print factor_str(ketxx)
  ketxx = mulvm(ket01, Swap)
  print "  mulvm(ket01, Swap) = |0>|1>Swap = |01> Swap =",
  print ketxx
  print "  factor(previous) =",
  print factor(ketxx)
  print factor_str(ketxx)
  ketxx = mulvm(ket10, Swap)
  print "  mulvm(ket10, Swap) = |1>|0>Swap = |10> Swap =",
  print ketxx
  print "  factor(previous) =",
  print factor(ketxx)
  print factor_str(ketxx)
  ketxx = mulvm(ket11, Swap)
  print "  mulvm(ket11, Swap) = |1>|1>Swap = |11> Swap =",
  print ketxx
  print "  factor(previous) =",
  print factor(ketxx)
  print factor_str(ketxx)
  print " "


  # check out CSwap
  print "three input gate, controlled swap, CSwap = Fredkin = ",
  print CSwap
  print " "
  ketxxx = mulmv(CSwap, ket001)
  print "  mulmv(CSwap, ket001) = CSwap |001> = ",
  print ketxxx
  ketab, ketcd, ketef = factor3(ketxxx)
  print "  factor and last two kets not swapped = "
  print "  ",
  print ketcd
  print "  ",
  print ketef
  ketxxx = mulmv(CSwap, ket101)
  print "  mulmv(CSwap, ket101) = CSwap |101> = ",
  print ketxxx
  ketab, ketcd, ketef = factor3(ketxxx)
  print "  factor and last two kets swapped = "
  print "  ",
  print ketcd
  print "  ",
  print ketef
  print " "

  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 = "
  print S
  T = addm(PZ, negm(PX))
  T = divm(T, math.sqrt(2))
  print "T = PZ + PX = "
  print T
  QS = KP(Q, S)
  print "QS = KP(Q, S) = "
  print QS
  RS = KP(R, S)
  print "RS = KP(R, S) = "
  print RS
  RT = KP(R, T)
  print "RT = KP(R, T) = "
  print RT
  QT = KP(Q, T)
  print "QT = KP(Q, T) = "
  print QT
  s0110 = addv(ket01, negv(ket10))
  s0110 = divv(s0110, math.sqrt(2))
  print "s0110 = (ket01-ket10)/sqrt(2) = "
  print s0110
  junk = mulmv(QT, s0110)
  val = DP(s0110, junk)
  print "<s0110|QT|s0110> = ",
  print val
  s1001 = addv(ket10, negv(ket01))
  s1001 = divv(s1001, math.sqrt(2))
  print "s1001 = (ket10-ket01)/sqrt(2) = "
  print s1001
  junk = mulmv(QT, s1001)
  val = DP(s1001, junk)
  print "<s1001|QT|s1001> = ",
  print val
  print " "

  print "test utility functions"
  test()
  print "quantum.py finished"

# 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"

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]

  return Y

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]

  return Y

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

def TP4(A, B, C, D): # Tensor Product of four vectors returning vector
  bcd = TP3(B, C, D)
  return TP(A, bcd)
 
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",
    print na, nb

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

  return Y

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]

  return Y
  
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]

  return Y

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",
    print na, nb

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

  return Y
  
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",
    print nar, nbr

  if nac != nbc :
    print "error, columns different in addm",
    print nac, nbc

  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]

  return Y

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

  return Y

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

  return Y
  
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",
    print na, nb

  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])

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

  return Y

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",
    print nv, nmr

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

  return Y
 
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",
    print nmc, nv

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

  return Y

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

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

  return Y

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]

  return Y

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

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

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

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

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

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

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


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

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

    j = j/2
    n = n/2

  return Y

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

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

    j = j/2
    n = n/2

  return Y

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

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

    j = j/2
    n = n/2

  return Y

# 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

    return [a, b], [c, d]

  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

    return [a, b], [c, d]

  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

    return [a, b], [c, d]

  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

    return [a, b], [c, d]


# 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]]

  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]

def test():
  # test
  m1 = [[1,2,3,4],[5,6,7,8]]
  print "m1 ="
  print m1
  m2 = [[10,11,12],[13,14,15],[16,17,18],[19,20,21]]
  print "m2 ="
  print m2
  v1 = [10,11]
  print "v1 ="
  print v1
  print "mulvm(v1,m1) = v1*m1 ="
  vx = mulvm(v1,m1)
  print vx
  v2 = [10,11,12,13]
  print "v2="
  print v2
  print "mulmv(m1,v2) = m1*v2 ="
  vx = mulmv(m1,v2)
  print vx
  print "m3 = m1*m2 ="
  m3 = mulmm(m1,m2)
  print m3
  print " "
  print "test getket by number, 0, 1, 2, ..."
  print "getket(0,ket000) =",
  print getket(0,ket000)
  print "getket(1,ket010) =",
  print getket(1,ket010)
  print "getket(2,ket100) =",
  print getket(2,ket100)
  print "getket(0,ket001) =",
  print getket(0,ket001)
  print "getket(1,ket011) =",
  print getket(1,ket011)
  print "getket(2,ket101) =",
  print getket(2,ket101)
  print "getbra(1,bra01) =",
  print getbra(1,bra01)
  print " "
  print "test factor on a pure kets, bras same results, added string printout"
  print "factor(ket0) =",
  print factor(ket0)
  print "factor_str(ket0) =",
  print factor_str(ket0)
  print "factor(ket00) =",
  print factor(ket00)
  print factor_str(ket00)
  print "factor(ket000) =",
  print factor(ket000)
  print factor_str(ket000)
  print "factor(ket001) =",
  print factor(ket001)
  print factor_str(ket001)
  print "factor(ket010) =",
  print factor(ket010)
  print factor_str(ket010)
  print "factor(ket011) =",
  print factor(ket011)
  print factor_str(ket011)
  print "factor(ket100) =",
  print factor(ket100)
  print factor_str(ket100)
  print "factor(ket101) =",
  print factor(ket101)
  print factor_str(ket101)
  print "factor(ket110) =",
  print factor(ket110)
  print factor_str(ket110)
  print "factor(ket111) =",
  print factor(ket111)
  print factor_str(ket111)
  print " "
  print "factor_strb(bra01) = ",
  print 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> = ",
  print a,
  print b
  print "|cd> = ",
  print c,
  print d
  s = [a*c, a*d, b*c, b*d] # tensor product
  ketab, ketcd = factor2(s)
  print "factor2(",
  print s,
  print ") = "
  print "ketab |ab> = ",
  print ketab
  print "ketcd |cd> = ",
  print 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> = ",
  print kab
  print "|cd> = ",
  print kcd
  print "|ef> = ",
  print 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(",
  print s,
  print ") = "
  print "ketab |ab> = ",
  print ketab
  print "ketcd |cd> = ",
  print ketcd
  print "ketef |ef> = ",
  print 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 " "

if __name__ == '__main__':
  main()