# qfunct.py functions for use in quantum simulation
#             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, keti, 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 
  global Fadd, Fadds, Faddc, FA, FAs, FAc

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

  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|
  keti = [0+1j] # complex
  keti00 = TP3(keti, ket0, ket0)
  keti01 = TP3(keti, ket0, ket1)
  keti10 = TP3(keti, ket1, ket0)
  keti11 = TP3(keti, ket1, ket1)
  ket1010 = TP4(ket1, ket0, ket1, ket0)

  # typical variables
  ketx = [0,0] # temporaries
  ketxx = [0,0,0,0] # some may call these "register" or "state" 
  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

  # 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]]
  FA = [[1,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0],[0,0,1,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,0,1,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,0,0,1,0,0,0,0]]
  FAs = [[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[1,0,0,0,0,0,0,0],[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0]]
  FAc = [[1,0,0,0,0,0,0,0],[1,0,0,0,0,0,0,0],[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,1,0,0,0,0,0,0]]

  print "quantum function initialization complete"

# constant definitions
def def_ket0(): # vectors
  return [1,0]

def def_ket1():
  return [0,1]

def def_keti():
  return [0+1j,0]

def def_H():  # matrices 
  return H

def def_PX():
  return PX

def def_Not():
  return PX

def def_PY():
  return PY

def def_PZ():
  return PZ

def def_Cnot():
  return Cnot

def def_CCnot():
  return CCnot

def def_Swap():
  return Swap

def  def_CSwap():
  return CSwap

def def_Fadd():
  Fadd = transp(FA)
  return Fadd

def def_Fadds():
  Fadds = transp(FAs)
  return Fadds

def def_Faddc():
  Faddc = transp(FAc)
  return Faddc

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].conjugate())

  return Y

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

def TP4(A, B, C, D): # Tensor Product of four vectors returning vector
  abc = TP3(A, B, C)
  return TP(abc, D)
 
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 state2int(S):
  n = len(S)
  if n==2: # 1 ket
    y = 0
    if DP(bra1, k)>0.5:
      y = 1

    return y

  elif n==4: # 2 ket
    y = 0
    a, b = factor(S)
    if DP(bra1, b)>0.5:
      y = 1

    if DP(bra1, a)>0.5:
      y = y+2

    return y

  elif n==8: # 3 ket
    y = 0
    a, b, c = factor3(S)
    if DP(bra1, c)>0.5:
      y = 1

    if DP(bra1, b)>0.5:
      y = y+2

    if DP(bra1, a)>0.5:
      y = y+4

    return y

  elif n==16: # 4 ket
    y = 0
    a, b, c, d = factor4(S)
    if DP(bra1, d)>0.5:
      y = 1

    if DP(bra1, c)>0.5:
      y = y+2

    if DP(bra1, b)>0.5:
      y = y+4

    if DP(bra1, a)>0.5:
      y = y+8

    return y

  else:
    print "limitation: state more than 4 kets"
    return 0


def int2state(n, ns):
  if n==0:
    if ns==1:
      return ket0
    else:
      return ket1 

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 measure3(V1, M, V2):
  Vt = mulmv(M, V2)
  z = DP(V1, Vt)
  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 or i
  n = len(K)
  Y = []
  j = n-1
  for i in range(n):
    if abs(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 or i
  n = len(B)
  Y = []
  for i in range(n):
    if abs(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]

# end factor3

# factor4  given [aceg, aceh, acfg, acfh,
#                 adeg, adeh, adfg, adfh,
#                 bceg, bceh, bcfg, bcfh,
#                 bdeg, bdeh, bdfg, bdfh] a state
#  a^2+b^2=1, c^2+d^2=1, e^2+f^2=1, g*2+h^2=1,
#  find [a,b,c,d,e,f,g,h] to get
#  |ab> |cd> |ef>  |gh> four 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 factor4(s):
  aceg = s[0]
  aceh = s[1]
  acfg = s[2]
  acfh = s[3]
  adeg = s[4]
  adeh = s[5]
  adfg = s[6]
  adfh = s[7]
  bceg = s[8]
  bceh = s[9]
  bcfg = s[10]
  bcfh = s[11]
  bdeg = s[12]
  bdeh = s[13]
  bdfg = s[14]
  bdfh = s[15]
  e   = 0
  if abs(aceh) > 0.01:
    v = aceg/aceh
    g = cmath.sqrt((v*v)/(1.0+v*v))
    h = cmath.sqrt(1-g*g)

  if abs(acfh) > 0.01:
    v = acfg/acfh
    g = cmath.sqrt((v*v)/(1.0+v*v))
    h = cmath.sqrt(1-g*g)

  elif abs(adeh) > 0.01:
    v = adeg/adeh
    g = cmath.sqrt((v*v)/(1.0+v*v))
    h = cmath.sqrt(1-g*g)

  elif abs(adfh) > 0.01:
    v = adfg/adfh
    g = cmath.sqrt((v*v)/(1.0+v*v))
    h = cmath.sqrt(1-g*g)

  elif abs(bceh) > 0.01:
    v = bceg/bceh
    g = cmath.sqrt((v*v)/(1.0+v*v))
    h = cmath.sqrt(1-g*g)

  elif abs(bcfh) > 0.01:
    v = bcfg/bcfh
    g = cmath.sqrt((v*v)/(1.0+v*v))
    h = cmath.sqrt(1-g*g)

  elif abs(bdeh) > 0.01:
    v = bdeg/bdeh
    g = cmath.sqrt((v*v)/(1.0+v*v))
    h = cmath.sqrt(1-g*g)

  elif abs(bdfh) > 0.01:
    v = bdfg/bdfh
    g = cmath.sqrt((v*v)/(1.0+v*v))
    h = cmath.sqrt(1-g*g)

  elif abs(aceg) > 0.01:
    v = aceh/aceg
    h = cmath.sqrt((v*v)/(1.0+v*v))
    g = cmath.sqrt(1-h*h)

  elif abs(acfg) > 0.01:
    v = acfh/acfg
    h = cmath.sqrt((v*v)/(1.0+v*v))
    g = cmath.sqrt(1-h*h)

  elif abs(adeg) > 0.01:
    v = adeh/adeg
    h = cmath.sqrt((v*v)/(1.0+v*v))
    g = cmath.sqrt(1-h*h)

  elif abs(adfg) > 0.01:
    v = adfh/adfg
    h = cmath.sqrt((v*v)/(1.0+v*v))
    g = cmath.sqrt(1-h*h)

  elif abs(bceg) > 0.01:
    v = bceh/bceg
    h = cmath.sqrt((v*v)/(1.0+v*v))
    g = cmath.sqrt(1-h*h)

  elif abs(bcfg) > 0.01:
    v = bcfh/bcfg
    h = cmath.sqrt((v*v)/(1.0+v*v))
    g = cmath.sqrt(1-h*h)

  elif abs(bdeg) > 0.01:
    v = bdeh/bdeg
    h = cmath.sqrt((v*v)/(1.0+v*v))
    g = cmath.sqrt(1-h*h)

  elif abs(bdfg) > 0.01:
    v = bdfh/bdfg
    h = cmath.sqrt((v*v)/(1.0+v*v))
    g = cmath.sqrt(1-h*h)

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

  if abs(g) < 0.001:
    ketab, ketcd, ketef = factor3([aceh/h, acfh/h, adeh/h, adfh/h,
                                   bceh/h, bcfh/h, bdeh/h, bdfh/h])
    return ketab, ketcd, ketef, [g, h]
  else:
    ketab, ketcd, ketef = factor3([aceg/g, acfg/g, adeg/g, adfg/g, 
                                  bceg/g, bcfg/g, bdeg/g, bdfg/g])
    return ketab, ketcd, ketef, [g, h]

# end factor4