# Let O1 be an octonion algebra with base {1,i,j,ij,k,ik,(ij)k} and O2 be an octonion algebra with base {1,I,J,IJ,K,IK,JK,(IJ)K}. O1xO2 has dimension 64, we will represent its elements with a list of length 64 which corresponds to the coefficients taken in the base of O1xO2 given by gen_OxO:
gen_OxO=['1X1','1XI','1XJ','1XIJ','1XK','1XIK','1XJK','1X(IJ)K'
,'iX1','iXI','iXJ','iXIJ','iXK','iXIK','iXJK','iX(IJ)K',
'jX1','jXI','jXJ','jXIJ','jXK','jXIK','jXJK','jX(IJ)K',
'ijX1','ijXI','ijXJ','ijXIJ','ijXK','ijXIK','ijXJK','ijX(IJ)K',
'kX1','kXI','kXJ','kXIJ','kXK','kXIK','kXJK','kX(IJ)K',
'ikX1','ikXI','ikXJ','ikXIJ','ikXK','ikXIK','ikXJK','ikX(IJ)K',
'jkX1','jkXI','jkXJ','jkXIJ','jkXK','jkXIK','jkXJK','jkX(IJ)K',
'(ij)kX1','(ij)kXI','(ij)kXJ','(ij)kXIJ','(ij)kXK','(ij)kXIK','(ij)kXJK','(ij)kX(IJ)K']

#the class of bioctonionalgebras. When an instance is created, automatically the matrix to calculate the multiplication effeciently is determined.
class OxOparent():
    def __init__(self,O1,O2):
        if isinstance(O1,octonion_parent) and isinstance(O1,octonion_parent) and O1.K==O2.K:
            self.O1=O1
            self.O2=O2
            #We have that gen_OxO[i]*gen_OxO[j]=a*gen_OxO[k] for some a\in K. We will define mult_mat[i][j]=[k,a]
            mult_mat=[[0 for r in [0..63]] for r in [0..63]]
            #gen_OxO[i]=x1\times y1
            for i in [0..63]:
                x1=[0 for k in [0..63]];x1[i//8]=1;x1=octonion(x1,self.O1)
                y1=[0 for k in [0..63]];y1[i%8]=1;y1=octonion(y1,self.O2)
                #gen_OxO[j]=x2\times y2
                for j in [0..63]:
                    x2=[0 for k in [0..63]];x2[j//8]=1;x2=octonion(x2,self.O1)
                    y2=[0 for k in [0..63]];y2[j%8]=1;y2=octonion(y2,self.O2)
                    tensorproduct=self.embed_in_OxO(x1*x2,y1*y2)
                    #Determination of index k
                    k=0
                    while tensorproduct.L[k]==0 and k<63:
                        k+=1
                    mult_mat[i][j]=[k,tensorproduct.L[k]]
            self.mult_mat=mult_mat
        else:
            raise TypeError("O1 and O2 should be octonion algebras over the same field.")

    def __repr__(self):
        return "The algebra O1XO2 over K, O1 has parameters %s %s %s, O2 has parameters %s %s %s" %(self.O1.a, self.O1.b, self.O1.c,self.O2.a, self.O2.b, self.O2.c)
        
    def __eq__(self, y):
        return self.O1==y.O1 and self.O2==y.O2
   
    #generates a random element of the octonion with coefficients in Q
    def random(self):
        return OxO( [QQ.random_element() for i in [0..63]],self)
    #generates a random skew element of the octonion with coefficients in Q
    def random_skew(self):
    	return (self.embed_in_OxO(self.O1.random_skew(),self.O2.one())+self.embed_in_OxO(self.O1.one(),self.O2.random_skew()))

    def zero(self):
        return OxO([0 for i in [0..63]],self)  

    #input: x\in O1, y\in O2 returns x\otimes y as an instance of OxO   
    def embed_in_OxO(self,x,y):
        if x.par==self.O1 and y.par==self.O2:
            ret=[]
            for i in [0..7]:
                for j in [0..7]:
                    ret.append((x.L[i])*(y.L[j]))
            return OxO(ret,self) 
        raise TypeError("x\in O1, y\in O2 has to hold!")

    def __call__(self,L):
        if L in self.O1.K:
            ret=[0 for i in [0..63]]; ret[0]=L
            return OxO(ret,self)
        elif len(L)==64:
            return OxO(L,self)
        else:
           raise TypeError("input should either be in K or a list of 64 elements") 

#instances of OxO are elements of the bioctonion algebra
class OxO():
    def __init__(self,L,bioct):
        if len(L)==64:
            self.L=L
            self.par=bioct
        else:
            raise TypeError("input should be a list of 64 coordinates")
            
    def __repr__(self):
        ret=""        
        written=false
        for i in [0..63]:
            if self.L[i]!=0:
                if written:
                    ret+=' + '
                if self.L[i]==1 and i!=0:
                    ret+=gen_OxO[i]
                else:        
                    ret+='(%s) %s' %(self.L[i],gen_OxO[i])
                written=true
        if not written:
            ret='0'
        return ret
    #addition    
    def __add__(self,y):
        if self.par==y.par:
           return OxO([self.L[i]+y.L[i] for i in [0..63]],self.par)
        raise TypeError("Elements should be in the same algebra.")
        
    def __sub__(self,y):
        if self.par==y.par:
            return OxO([self.L[i]-y.L[i] for i in [0..63]],self.par)
        raise TypeError("Elements should be in the same algebra.")
          
    def __neg__(self):
           return OxO([-self.L[i] for i in [0..63]],self.par)
           
    def __mul__(self,y):
       if y in self.par.O1.K:
           return OxO([(y*self.L[i]) for i in [0..63]],self.par)
       if self.par==y.par:
           #we use mult_mat and use the distributivity of the multiplication.
           ret=[0 for i in [0..63]]
           for i in [0..63]:
               if self.L[i]!=0:
                   for j in [0..63]:
                       if y.L[j]!=0:
                           ret[self.par.mult_mat[i][j][0]]+=self.par.mult_mat[i][j][1]*self.L[i]*y.L[j]
           return OxO(ret,self.par)
       raise TypeError(" Elements should be in the same algebra or one should be in K.")
    
    def __rmul__(self,y):
        return self*y
    
    #division by elements in the basefield
    def __div__(self,y):
        if y in K:
           return OxO([((1/y)*self.L[i]) for i in [0..63]],self.par)
        raise TypeError("One element should be in K.")           
    
    #taking powers should only be done for low coefficients
    def __pow__(self,n):
        if n==0:
            ret=[1]
            ret.extend([0 for i in [0..62]])
            return OxO(ret,self.par)
        elif n>0:
            power=self
            k=1
            while(k<n):
                power=power*self
                k+=1
            return power
        
        raise TypeError("n should be >=0")
            
    def __eq__(self,y):
        return self.par==y.par and self.L==y.L
            
    def __getitem__(self,i):
        return self.L[i]
                
    def __setitem__(self, i, value):
        self.L[i]=value 

    #the standard involution on a bioctonion   
    def si(self):
        ret=[self.L[0]]
        ret.extend([-self.L[i] for i in [1..7]])
        for i in [8..63]:
            if i%8==0:                
                ret.append(-self.L[i])
            else:
                ret.append(self.L[i])
        return OxO(ret,self.par)
    
    #if s=s1x1+1xs2 is skew, split_skew(s) returns [s1,s2]
    def split_skew(self):
        if self.L[0]!=0:
            raise TypeError("given element should be skew")
        o1=[0]
        for i in [8..63]:
            if i%8==0: 
                o1.append(self.L[i])
            else: 
                if self.L[i]!=0:
                    raise TypeError("given element should be skew")
        o1=octonion(o1,self.par.O1)
        o2=[0];o2.extend(self.L[1:8])
        o2=octonion(o2,self.par.O2)
        return [o1,o2]
    #the albert quadratic form on skew elements
    def q_albert(self):
    	L=self.split_skew()
        return L[0].norm()-L[1].norm()
    #the linearization of the albert quadratic form
    def f_albert(self,o2):
	    return (self+o2).q_albert()-self.q_albert()-o2.q_albert()

    #the inverse of a skew element if it is invertible
    def inverse(self):
    	sk=self.split_skew()
        alb=sk[0].norm()-sk[1].norm()
        if alb==0:
		    raise TypeError("%s is not invertible!" %self)
        return (self.par.embed_in_OxO(-sk[0]/alb,self.par.O2.one())+self.par.embed_in_OxO(self.par.O1.one(),sk[1]/alb))
    
    #tests if a given element is skew
    def is_skew(self):
    	skew=True
    	if self.L[0]!=0:
		    skew=False
        for i in [8..63]:
		    if i%8!=0 and self.L[i]!=0: 
			    skew=False
        return skew    
                       
    #the V-operator of a structurable algebra
    def V(self,Y,Z):
        Ybar=Y.si()
        return (self*Ybar)*Z+(Z*Ybar)*self-(Z*self.si())*Y