#The low
e1 = matrix(8,8,{(0,1):-1, (1,0):1, (2,3):-1, (3,2):1, (4,5):-1, (5,4):1, (6,7):-1, (7,6):1})
e2 = matrix(8,8,{(0,2):-1, (2,0):1, (1,3):1, (3,1):-1, (4,6):1, (6,4):-1, (5,7):-1, (7,5):1})
e3 = matrix(8,8,{(0,3):-1, (1,2):-1, (2,1):1,(3,0):1, (4,7):-1, (5,6):-1, (6,5):1, (7,4):1})
e4 = matrix(8,8,{(0,4):-1, (1,5):1, (2,6):-1,(3,7):1, (4,0):1, (5,1):-1, (6,2):1, (7,3):-1})
e5 = matrix(8,8,{(0,5):-1, (1,4):-1, (2,7):1,(3,6):1}) 
e5 = e5  -e5.transpose() 
e6 = matrix(8,8,{(0,6):-1, (1,7):1, (2,4):1,(3,5):-1})
e6 = e6 - e6.transpose()
e7 = matrix(8,8,{(0,7):-1, (1,6):-1, (2,5):-1,(3,4):-1})
e7 = e7 - e7.transpose()
em = matrix(8,8,{(0,0):1, (1,1):-1, (2,2):-1,(3,3):-1,(4,4):-1,(5,5):-1,(6,6):-1,(7,7):-1})
unit = [e1^4, e1,e2,e3,e4,e5,e6,e7]

#Def prod of Cayley alg
def OCN(X,Y):
    e1 = matrix(8,8,{(0,1):-1, (1,0):1, (2,3):-1, (3,2):1, (4,5):-1, (5,4):1, (6,7):-1, (7,6):1})
    e2 = matrix(8,8,{(0,2):-1, (2,0):1, (1,3):1, (3,1):-1, (4,6):1, (6,4):-1, (5,7):-1, (7,5):1})
    e3 = matrix(8,8,{(0,3):-1, (1,2):-1, (2,1):1,(3,0):1, (4,7):-1, (5,6):-1, (6,5):1, (7,4):1})
    e4 = matrix(8,8,{(0,4):-1, (1,5):1, (2,6):-1,(3,7):1, (4,0):1, (5,1):-1, (6,2):1, (7,3):-1})
    e5 = matrix(8,8,{(0,5):-1, (1,4):-1, (2,7):1,(3,6):1}) 
    e5 = e5  -e5.transpose() 
    e6 = matrix(8,8,{(0,6):-1, (1,7):1, (2,4):1,(3,5):-1})
    e6 = e6 - e6.transpose()
    e7 = matrix(8,8,{(0,7):-1, (1,6):-1, (2,5):-1,(3,4):-1})
    e7 = e7 - e7.transpose()
    em = matrix(8,8,{(0,0):1, (1,1):-1, (2,2):-1,(3,3):-1,(4,4):-1,(5,5):-1,(6,6):-1,(7,7):-1})
    unit = [e1^4, e1,e2,e3,e4,e5,e6,e7]
    summ = matrix([[0,0,0,0,0,0,0,0]]).transpose()
    for n in range(8):
        summ += X[0][n] * unit[n]  * Y.transpose()
    return summ.transpose()

#Def: Variables
var("x0 x1 x2 x3 x4 x5 x6 x7 y0 y1 y2 y3 y4 y5 y6 y7 z0 z1 z2 z3 z4 z5 z6 z7 ζ1 ζ2 ζ3 η1 η2 η3")
var("u0 u1 u2 u3 u4 u5 u6 u7  v0 v1 v2 v3 v4 v5 v6 v7 w0 w1 w2 w3 w4 w5 w6 w7")
unit = [e1^4, e1,e2,e3,e4,e5,e6,e7]

X = matrix([[x0,x1,x2,x3,x4,x5,x6,x7]])
Y = matrix([[y0,y1,y2,y3,y4,y5,y6,y7]])
Z = matrix([[z0,z1,z2,z3,z4,z5,z6,z7]])
U = matrix([[u0,u1,u2,u3,u4,u5,u6,u7]])
V = matrix([[v0,v1,v2,v3,v4,v5,v6,v7]])
W = matrix([[w0,w1,w2,w3,w4,w5,w6,w7]])

b1 = matrix([[ζ1,0,0,0,0,0,0,0]])
b2 = matrix([[ζ2,0,0,0,0,0,0,0]])
b3 = matrix([[ζ3,0,0,0,0,0,0,0]])
d1 = matrix([[η1,0,0,0,0,0,0,0]])
d2 = matrix([[η2,0,0,0,0,0,0,0]])
d3 = matrix([[η3,0,0,0,0,0,0,0]])

J1 = [[b1, Z, Y*em],[Z*em, b2, X],[Y,X*em,b3]]
J2 = [[d1, W, V*em],[W*em, d2, U],[V,U*em,d3]]

#Def: Cul of M(Cayley,3)
def JCN(J1,J2):
    c11 = OCN(J1[0][0],J2[0][0]) + OCN(J1[0][1],J2[1][0]) + OCN(J1[0][2],J2[2][0])
    c12 = OCN(J1[0][0],J2[0][1]) + OCN(J1[0][1],J2[1][1]) + OCN(J1[0][2],J2[2][1])
    c13 = OCN(J1[0][0],J2[0][2]) + OCN(J1[0][1],J2[1][2]) + OCN(J1[0][2],J2[2][2])
    c21 = OCN(J1[1][0],J2[0][0]) + OCN(J1[1][1],J2[1][0]) + OCN(J1[1][2],J2[2][0])
    c22 = OCN(J1[1][0],J2[0][1]) + OCN(J1[1][1],J2[1][1]) + OCN(J1[1][2],J2[2][1])
    c23 = OCN(J1[1][0],J2[0][2]) + OCN(J1[1][1],J2[1][2]) + OCN(J1[1][2],J2[2][2])
    c31 = OCN(J1[2][0],J2[0][0]) + OCN(J1[2][1],J2[1][0]) + OCN(J1[2][2],J2[2][0])
    c32 = OCN(J1[2][0],J2[0][1]) + OCN(J1[2][1],J2[1][1]) + OCN(J1[2][2],J2[2][1])
    c33 = OCN(J1[2][0],J2[0][2]) + OCN(J1[2][1],J2[1][2]) + OCN(J1[2][2],J2[2][2])
    return [[c11,c12,c13],[c21,c22,c23],[c31,c32,c33]]

#Def: Trace
def O_trace(J1,J2):
    return JCN(J1,J2)[0][0] + JCN(J1,J2)[1][1] + JCN(J1,J2)[2][2]

#Def: Jordan Product
def J_prod(J1,J2):
    jordan = [[JCN(J1,J2)[0][0]/2 + JCN(J2,J1)[0][0]/2,JCN(J1,J2)[0][1]/2 + JCN(J2,J1)[0][1]/2,JCN(J1,J2)[0][2]/2 + JCN(J2,J1)[0][2]/2],
              [JCN(J1,J2)[1][0]/2 + JCN(J2,J1)[1][0]/2,JCN(J1,J2)[1][1]/2 + JCN(J2,J1)[1][1]/2,JCN(J1,J2)[1][2]/2 + JCN(J2,J1)[1][2]/2],
              [JCN(J1,J2)[2][0]/2 + JCN(J2,J1)[2][0]/2,JCN(J1,J2)[2][1]/2 + JCN(J2,J1)[2][1]/2,JCN(J1,J2)[2][2]/2+ JCN(J2,J1)[2][2]/2]]
    return jordan[0][0] + jordan[1][1] + jordan[2][2]

def J_prod_o(J1,J2):
    jordan = [[JCN(J1,J2)[0][0]/2 + JCN(J2,J1)[0][0]/2,
               JCN(J1,J2)[0][1]/2 + JCN(J2,J1)[0][1]/2,
               JCN(J1,J2)[0][2]/2 + JCN(J2,J1)[0][2]/2],
              [JCN(J1,J2)[1][0]/2 + JCN(J2,J1)[1][0]/2,
               JCN(J1,J2)[1][1]/2 + JCN(J2,J1)[1][1]/2,
               JCN(J1,J2)[1][2]/2 + JCN(J2,J1)[1][2]/2],
              [JCN(J1,J2)[2][0]/2 + JCN(J2,J1)[2][0]/2,
               JCN(J1,J2)[2][1]/2 + JCN(J2,J1)[2][1]/2,
               JCN(J1,J2)[2][2]/2+ JCN(J2,J1)[2][2]/2]]
    return jordan

#修正中
def Cross(J1,J2):
    
    #Tr J1 in Real
    trX = J1[0][0] + J1[1][1] + J1[2][2]
    trX = trX[0][0]
    
    #J3 is TrJ1 * J2
    J3 = [[trX * J2[0][0],trX * J2[0][1],trX * J2[0][2]],
          [trX * J2[1][0], trX * J2[1][1], trX * J2[1][2]],
          [trX * J2[2][0],trX * J2[2][1],trX * J2[2][2]]]
    
    #Tr J2 = Tr Y
    trY = J2[0][0] + J2[1][1] + J2[2][2]
    trY = trY[0][0]
    
    #Tr J2 * J1
    J4 = [[trY * J1[0][0],trY * J1[0][1],trY * J1[0][2]],
          [trY * J1[1][0], trY * J1[1][1], trY * J1[1][2]],
          [trY * J1[2][0],trY * J1[2][1],trY * J1[2][2]]]
    
    #KKK is TrJ1 * TrJ2 - X dot Y
    KKK = matrix([[trX*trY,0,0,0,0,0,0,0]]) - J_prod(J1,J2)
    
    #UUU is Mat C,3
    UUU = [[KKK,matrix(1,8),matrix(1,8)],
           [matrix(1,8),KKK,matrix(1,8)],
           [matrix(1,8),matrix(1,8),KKK]]

    return [[1/2*expand((((2*J_prod_o(J1,J2)[i][j])-J3[i][j]-J4[i][j] + UUU[i][j]))) for j in range(3)] for i in range(3)]