CAlgRank4SFromList:=function(L)
local B0,B1,B2,B3,T0,T1,T2,T3,T4,T5,T6,T7,T8,T9,T;

B0:=IdentityMat(4);

B1:=[[0,L[10],0,0],[1,L[10]-L[1]-L[4]-1,L[1],L[4],],[0,L[10]-L[2]-L[5],L[2],L[5]],[0,L[10]-L[3]-L[6],L[3],L[6]]];

B2:=[[0,0,L[11],0],[0,L[1],L[11]-L[1]-L[7],L[7]],[1,L[2],L[11]-L[2]-L[8]-1,L[8]],[0,L[3],L[11]-L[3]-L[9],L[9]]];

B3:=[[0,0,0,L[12]],[0,L[4],L[7],L[12]-L[7]-L[4]],[0,L[5],L[8],L[12]-L[8]-L[5]],[1,L[6],L[9],L[12]-L[9]-L[6]-1]];

T1:=B1^2-(L[10]*B0+(L[10]-L[1]-L[4]-1)*B1+(L[10]-L[2]-L[5])*B2+(L[10]-L[3]-L[6])*B3);
T1:=Concatenation(T1);
T2:=B1*B2-(L[1]*B1+L[2]*B2+L[3]*B3);
T2:=Concatenation(T2);
T0:=Concatenation(T1,T2);
T3:=B1*B3-(L[4]*B1+L[5]*B2+L[6]*B3);
T3:=Concatenation(T3);
T0:=Concatenation(T0,T3);
T4:=B2^2-(L[11]*B0+(L[11]-L[1]-L[7])*B1+(L[11]-L[2]-L[8]-1)*B2+(L[11]-L[3]-L[9])*B3);
T4:=Concatenation(T4);
T0:=Concatenation(T0,T4);
T5:=B2*B3-(L[7]*B1+L[8]*B2+L[9]*B3); 
T5:=Concatenation(T5);
T0:=Concatenation(T0,T5);
T6:=B3^2-(L[12]*B0+(L[12]-L[4]-L[7])*B1+(L[12]-L[5]-L[8])*B2+(L[12]-L[6]-L[9]-1)*B3);
T6:=Concatenation(T6);
T0:=Concatenation(T0,T6);
T7:=B1*B2-B2*B1; 
T7:=Concatenation(T7);
T0:=Concatenation(T0,T7);
T8:=B1*B3-B3*B1;
T8:=Concatenation(T8);
T0:=Concatenation(T0,T8);
T9:=B2*B3-B3*B2;
T9:=Concatenation(T9);
T0:=Concatenation(T0,T9);

T:=AsSet(T0); 

return T; 
end; 
########################################
###########################################
EnumerationFromListRank4S:=function(n) 
local L,k,k1,k2,k3,i1,i2,i3,i4,i5,i6,i7,i8,i9,S,T,t,V,j; 

L:=[];
S:=[];
for k1 in [1..(n-3)] do 
k:=Minimum(k1,n-k1-2); 
for k2 in [1..k] do 
k3:=n-k1-k2-1;
if 0<k3 and k3<k1+1 then 
for i1 in [0..(k1-1)] do 
for i2 in [0..(k1-1)] do 
for i4 in [0..(k1-1)] do
for i8 in [0..(k2-1)] do 
T:=[];
i3:=(1/k3)*(k1*k2-i1*k1-i2*k2); 
if IsInt(i3) and not(i3 < 0) then 
i5:=(1/k2)*(i3*k3); 
if IsInt(i5) and not(i5 < 0) then 
i7:=(1/k1)*(i3*k3); 
if IsInt(i7) and not(i7 < 0) then 
i6:=(k1/k3)*(k3-i4-i7); 
if IsInt(i6) and not(i6 < 0) then 
i9:=(k2/k3)*(k3-i5-i8); 
if IsInt(i9) and not(i9 < 0) then 
  T:=[i1,i2,i3,i4,i5,i6,i7,i8,i9,k1,k2,k3]; 
  Add(S,T);
fi; 
fi; 
fi;
fi; 
fi; 
od; 
od; 
od; 
od; 
fi; 
od; 
od; 

for T in S do 
  t:=CAlgRank4SFromList(T); 
  if Size(t)=1 then 
    Add(L,T); 
  fi; 
od; 

return L; 
end; 
##############################################
#######################################
RegularMatricesRank4S:=function(V)
local B0,B1,B2,B3,B;

B0:=IdentityMat(4);
B1:=[[0,V[10],0,0],[1,V[10]-V[1]-V[4]-1,V[1],V[4],],[0,V[10]-V[2]-V[5],V[2],V[5]],[0,V[10]-V[3]-V[6],V[3],V[6]]];

B2:=[[0,0,V[11],0],[0,V[1],V[11]-V[1]-V[7],V[7]],[1,V[2],V[11]-V[2]-V[8]-1,V[8]],[0,V[3],V[11]-V[3]-V[9],V[9]]];

B3:=[[0,0,0,V[12]],[0,V[4],V[7],V[12]-V[7]-V[4]],[0,V[5],V[8],V[12]-V[8]-V[5]],[1,V[6],V[9],V[12]-V[9]-V[6]-1]];

B:=[B0,B1,B2,B3];

return B;
end; 
######################################
ITAsRnk4SOrder:=function(n)
local L,l,i,V,t,j,k,S,T,B,b,G,b1,B1,g12,g1234,g; 

S:=[]; 
T:=[];
L:=EnumerationFromListRank4S(n);
l:=Length(L);
for i in [1..l] do 
   B:=RegularMatricesRank4S(L[i]); 
   t:=0;
   for b in B do 
    for j in [1..4] do 
     for k in [1..4] do 
       if b[j][k]<0 then 
         t:=1; 
       fi; 
     od; 
    od; 
   od;
   if t=0 then Add(S,B); fi; 
 od; 
g12:=PermutationMat((1,2),4); 
g1234:=PermutationMat((1,2,3,4),4); 
G:=Group([g12,g1234]); 
for B in S do 
 t:=0;
 for g in G do 
   B1:=[];
   for b1 in B do 
    Add(B1,b1^g); 
   od; 
   for b in T do 
   if IsEqualSet(B1,b) then   
      t:=1; 
   fi; 
   od; 
 od; 
 if t=0 then    
  Add(T,B); 
 fi;
od;

return T; 
end;  
###################################

###################################
CharTableOfTA:=function(b) 
local r,n,i,m,K,K1,k,c,q,i1,F,D,f,d,L,Q,G,a,D1,L1,l1,j,F1,f1,d1,v,d2,J,J1,L3,L4,L5,m1,m2,m3,t; 

r:=Size(b);
n:=0; 
for i in [1..r] do 
  m:=Sum(b[i][1]); 
  n:=n+m; 
od; 
Print(n,"\n");
K:=[];
c:=0; 
for i in [2..r] do 
 q:=Degree(MinimalPolynomial(b[i]));
 if q=r then 
   c:=1; 
   i1:=i;
   Print("polynomial",Factors(MinimalPolynomial(b[i])),"\n");
 fi;
od; 
# Print(c,"\n"); 

if c=1 then 
 F:=Factors(MinimalPolynomial(b[i1])); 
 if Length(F)=r then 
   K:=Eigenvectors(Rationals,b[i1]); 
 else 
   D:=[1]; 
   for f in F do 
    d:=Degree(f);
    if d=2 then 
     Add(D,ER(Discriminant(f)));
    fi;
   od; 
   L:=Field(D); 
   Q:=Eigenvectors(L,b[i1]);
   if Length(Q)=r then 
     K:=Q; 
   else
    f:=F[Length(F)];
    G:=TransitiveGroup(Degree(f),GaloisType(f)); 
    if not(IsAbelian(G)) then 
      K:=["Not Cyclotomic",f,G]; 
    else
      a:=PrimitiveElement(L); 
      d:=Degree(f); 
      D1:=Discriminant(f); 
      Print("cyclotomic","\r");
      L1:=[];
      for m in [4..100] do
        if IsInt(Phi(m)/d) and IsInt(D1/m) then  
        Add(L1,Field([a,E(m)])); 
        fi;
      od; 
      L1:=AsSet(L1);
    #  Print(L1,"\n"); 
      for L in L1 do 
        l1:=Length(Factors(PolynomialRing(L),f));
    #    Print(L," ",l1,"\n");  
        if l1=d then 
          break; 
        fi; 
      od;  
    # for L in L1 do 
      Q:=Eigenvectors(L,b[i1]);
    # Print(L,Length(Q),"\r"); 
      if Length(Q)=r then 
          K:=Q; 
    #    break;
      fi; 
    # od; 
    fi; 
   fi; 
 fi; 
fi;

if c=1 and K=[] then 
  Print("cyclotomic, requires larger field","\n"); 
fi; 

if c=0 then 
 F:=[];
 for j in [2..r] do 
  F1:=AsSet(Factors(MinimalPolynomial(b[j]))); 
  F:=Concatenation(F,F1); 
  F:=AsSet(F); 
 od; 
 F:=SortedList(F); 
 Print(F,"\n");
 D:=[E(2)];
 d:=1; 
 for f in F do 
  d:=Maximum(Degree(f),d);
 od;
 Print(d,"\n"); 
 if d=1 then 
   L:=Field(D); 
 fi; 
 if d=2 then 
   for f in F do 
    if Degree(f)=2 then 
     d1:=Discriminant(f); 
     Add(D,ER(d1)); 
    fi; 
   od; 
   L:=Field(D); 
 fi;
 if d>2 then 
  G:=TransitiveGroup(Degree(f),GaloisType(f)); 
   if not(IsAbelian(G)) then 
     K:=["Not Cyclotomic",f,G];  
   fi;
 fi;
 if K=[] then   
  D:=[E(2)];  
  for f in F do 
   d:=Degree(f); 
   if d=2 then 
    d1:=Discriminant(f); 
    Add(D,ER(d1)); 
   fi; 
   if d>2 then 
    v:=Discriminant(f); 
    d1:=AbsoluteValue(v); 
    for d2 in Factors(d1) do 
    Add(D,E(d2)); 
    od; 
   fi; 
  od; 
  L:=Field(D);
  Print(L,"\n"); 
  d1:=1;
  for f in F do 
   F1:=Factors(PolynomialRing(L),f); 
   f1:=F1[Length(F1)]; 
   d1:=Maximum(d1,Degree(f1));
  od;
  Print(d1,"\n");     
  if d1=1 then 
   K1:=[]; 
   for k in [1..r] do 
    Q:=Eigenvectors(L,b[k]); 
    #Print(k,Q,"\n"); 
    d1:=DiagonalOfMat(Q*b[k]*Q^(-1)); 
    Add(K1,d1);  
   od;
   #Print(K1,"\n");
   J:=DiagonalOfMat(IdentityMat(r));
   Print(J*K1,"\n"); 
   #if not(J*K1=[n,0,0,0,0]) then 
   if not(J*K1=[n,0,0,0]) then 
     G:=SymmetricGroup([2..r]); 
     L3:=Orbit(G,K1[3],Permuted); 
     L4:=Orbit(G,K1[4],Permuted); 
   #  L5:=Orbit(G,K1[5],Permuted); 
     t:=0;
     for m1 in L3 do 
     for m2 in L4 do 
   #  for m3 in L5 do 
   #  Q:=[K1[1],K1[2],m1,m2,m3];
     Q:=[K1[1],K1[2],m1,m2]; 
     J1:=J*Q; 
     Print(J1,"\n");
   #  if J1=[n,0,0,0,0] then
     if J1=[n,0,0,0] then 
      K:=TransposedMat(Q); 
      t:=1;
      break;
     fi;
     od;
   #  if t=1 then 
   #   break;
   #  fi;
   #  od;
     if t=1 then 
      break; 
     fi;
     od; 
   fi;
  fi; 
 fi;
fi; 

return K; 
end;

#####################################
DegreesOfTA:=function(b)
local r,D,i;

r:=Length(b); 
D:=[];
for i in [1..r] do 
  D[i]:=b[i][1][i];
od; 

return D;
end;
#####################################
# Another name for character table of 
# a commutative table algebra is first 
# eigenmatrix. 
#####################################
FirstEigenmatrixOfTA:=function(b)
local P;

P:=CharTableOfTA(b); 

return P;
end;

#####################################
# Note that the second eigenmatrix is 
# computed from the first eigenmatrix.
#####################################
SecondEigenmatrixOfTA:=function(b)
local Q,n,P; 

P:=FirstEigenmatrixOfTA(b);
n:=Sum(P[1]); 
Q:=n*P^(-1);

return Q; 
end; 
#####################################
MultiplicitiesOfTA:=function(b)
local n,P,M,Q;

P:=FirstEigenmatrixOfTA(b);
n:=Sum(P[1]); 
Q:=n*P^(-1); 
M:=Q[1]; 

return M;
end;
#####################################
ClosedSubsetsOfTA:=function(b)
local r,C0,C1,c,C2,C,i,j,k,a,b1,D;

r:=Length(b);
a:=Algebra(Rationals,b);
b1:=Basis(a,b); 
C:=[];
C0:=Combinations([2..r]); 
C1:=[];
for c in C0 do 
  Add(C1,Concatenation([1],c));
od; 
for c in C1 do 
D:=[];
for j in c do 
for k in c do 
  C2:=Coefficients(b1,b1[j]*b1[k]); 
  for i in [1..r] do 
  if not(C2[i]=0) then 
    AddSet(D,i); 
  fi;
  od;
od;
od;

if D=c then
  Add(C,D);
fi; 
od;  

return C;
end;
################################### 
FusionsOfTA:=function(b)
local S,P,D,p,p1,b1,p2,d,a,r;

S:=[];
r:=Length(b);
P:=PartitionsSet([1..r]);
for p in P do 
 if not(p=[]) then 
 D:=[]; 
 for p1 in p do
 b1:=0*b[1]; 
  for p2 in p1 do 
   b1:=b1+b[p2]; 
  od; 
  Add(D,b1); 
 od; 
 a:=Algebra(Rationals,D);
 d:=Dimension(a); 
 if d=Length(p) then 
  Add(S,p); 
 fi; 
 fi;
od;

return S;
end;