############
# coxeter.g
############
#
# created 31 Jan 03
# by Jon McCammond and Woonjung Choi

# streamlined version of the earlier routineswhich
# uses the reflections (stored in reflections.g) 
# as input but otherwise avoids using chevie at all.

# The routines can be split into three parts
# 1. generic rountines for arbitrary rank
# 2. testing CAT(0) in the std metric for rank 4
# 3. constructing the 3-complex for rank 4
 
Print("...loading the coxeter reflections\n");
Read("reflections.g");
Print("\n To begin type something like createCoxeterInfo(\"A\",4);\n");
Print(" where (\"A\",4) can be replaced with any irreducible Coxeter\n");
Print(" type up to rank ",Length(CGTypeList),"\n\n");

# global variables
if not IsBound(type) then type:=[];fi;
if not IsBound(cw) then cw:=();fi;
if not IsBound(w) then w:=[];fi;
if not IsBound(R) then R:=[];fi;
if not IsBound(L) then L:=[];fi;
if not IsBound(len) then len:=[];fi;
if not IsBound(clen) then clen:=[];fi;
if not IsBound(V) then V:=[];fi;
if not IsBound(P) then P:=[];fi;
if not IsBound(PL) then PL:=[];fi;
if not IsBound(c) then c:=[];fi;
if not IsBound(wc) then wc:=();fi;
if not IsBound(paths) then paths:=[];fi;
if not IsBound(col) then col:=[];fi;
if not IsBound(columns) then columns:=[];fi;
if not IsBound(colPatterns) then colPatterns:=[];fi;
if not IsBound(badTri) then badTri:=[];fi;
if not IsBound(badPatterns) then badPatterns:=[];fi;
if not IsBound(verbose) then verbose:=true;fi;
if not IsBound(nPatterns) then nPatterns:=0;fi;
if not IsBound(Diagonals) then Diagonals:=[];fi;
if not IsBound(tetrahedra) then tetrahedra:=[];fi;


###################################
# Part I: Creating the initial Data
###################################
multiply:=function (L1,L2) # L1,L2 are lists of permutations
    local ans,i;
    ans:=[];
    for i in L1 do Append(ans,i*L2);od;
    return(Set(ans));
end;

kball:=function (L,k) # L is a list of permutations, k is radius
   local ans,i,current,previous,new;
   ans:=[L];
   previous:=[()];
   current:=L;
   for i in [2..k] do
       new:=multiply(current,L);
       SubtractSet(new,Union(current,previous));
       Add(ans,new);
       previous:=current;
       current:=new;
   od;
   return ans;
end;

ksphere:=function (L,k) return kball(L,k)[k];end;

invert:=function(L)
    local i;
    return(List([1..Length(L)],i->L[i]^-1));
end;

lengths:=function(L)
    local i;
    return(List([1..Length(L)],i->Length(L[i])));
end;

levelSets:=function(R,w,n)
    local L,B,l,u,i;          
    u:=(n+(n mod 2))/2;
    l:=(n-(n mod 2))/2;
    L:=[];
    B:=kball(R,u);
    for i in [1..u] do L[i]:=Set(B[i]);od;
    IntersectSet(L[u],invert(L[l])*w);
    for i in [1..u-1] do IntersectSet(L[u-i],multiply(R,L[u-i+1]));od;
    for i in [1..l-1] do L[u+i]:=Set(invert(L[l-i])*w);od;
    L[n]:=[w];
    len:=[1];
    Append(len,lengths(L));
    clen:=[1];
    for i in [2..type[2]] do Add(clen,clen[Length(clen)]+len[i]);od;
    V:=[()];
    Append(V,Flat(L));
    return L; 
end;

coxeterWord:=function(t,n)
    local i,u,l,r,s,ans;
    u:=(n+(n mod 2))/2;
    l:=(n-(n mod 2))/2;
    r:=[];s:=[];
    for i in [1..u] do r[i]:=2*i-1;od;
    for i in [1..l] do s[i]:=2*i;od;
    ans:=[];
    if t="D" then 
        ans:=[1];
        Append(ans,s);
        Append(ans,List([2..Length(r)],i->r[i]));
    elif t="E" then
        ans:=[1];
        Append(ans,List([2..Length(s)],i->s[i]));
        Append(ans,[2]);
        Append(ans,List([2..Length(r)],i->r[i]));
    else
        ans:=r;
        Append(ans,s);
    fi;
    return(ans);
end;

coxeterElement:=function(L,l)
    local ans,i;
    ans:=();
    for i in l do ans:=ans*L[i];od;
    return(ans);
end;

conjugates:=function(R)
    local i,j;
    return(List([1..Length(R)],i->List([1..Length(R)],j->Position(R,R[j]^R[i]))));
end;

wconjugates:=function(R,w)
    local i;
    return(PermList(List([1..Length(R)],i->Position(R,R[i]^w))));
end;

##############temp1
wcon:=function(K,w)
    local i;
    return(PermList(List([1..Length(K)],i->Position(V,w*K[i]*w^-1)-1)));
end;



wcon2:=function(K,w)
    local i;
    return(PermList(List([1..Length(K)],i->(Position(V,K[i]^w)-13)-1)));
end;

##############

wconj2:=function(R,w)
    local i;
    return(PermList(List([1..(Length(V)-2-2*Length(R))],i->(Position(V,V[i+Length(R)+1]^w)-(Length(R)+1)))));
end;




createPoset:=function(L,R)
    local i,po,p,q,r;
    po:=[];
    for i in [1..Length(L)] do
        for p in L[i] do
            for r in R do
                q:=r*p;
                if i=1  then 
                    if q=() then 
                        Add(po,[1,Position(V,p),Position(R,r)]);
                    fi;
                else
                    if (q in L[i-1]) then
                        Add(po,[Position(V,q),Position(V,p),Position(R,r)]);
                    fi;
                fi;
            od;
        od;
    od;
    return(Set(po));
end;

stripLabels:=function(P)
    local po,p;
    po:=[];
    for p in P do Add(po,[p[1],p[2]]);od;
    return(po);
end;

################################################
# This is main initial function.  It creates the
# posets and all their associated information.
################################################
createCoxeterInfo:=function(l,n)
    type:=[l,n];
    Print("\ntype = Coxeter group (\"",l,"\",",n,")\n");
          
          R:=CGReflections(l,n);
    Print("R    = ",Length(R)," reflections of W\n");
    
    cw:=coxeterWord(l,n);
    w:=coxeterElement(R,cw);
    Print("w    = coxeter element (",cw,")\n");
    Print("\n");
    
    c:=conjugates(R);
    wc:=wconjugates(R,w);
    Print("c  = conjugacy table for R\n");
    Print("wc = conjugates classes for R rel w\n");
    Print("\n");
    
    L:=levelSets(R,w,n);
    Print("L  = permutations below w (",len,")\n");
    Print("V  = ",Length(V)," permutations listed in L\n\n");
    
    PL:=createPoset(L,R);
    P:=stripLabels(PL);
    Print("P  = poset structure for L\n");
    Print("PL = poset structure for L with labels\n\n");
    
    Print("Other information is stored in len, clen, and cw --but \n");
    Print("these are mostly for internal use.\n\n");
    
    Print("  To calculate the column structure, type createColumns();\n\n");
    
    if n=4 then 
        Print("  Since this group is rank 4, try testGroup(); to see whether\n");
        Print("  the Brady-Krammer complex for the Artin group of this type\n");
        Print("  is CAT(0) using the standard coxeter shape metric.\n");
        Print("\n");
    fi;
end;


###################################################################
# Part II: Testing for CAT(0) in rank 4 (using the standard metric)
###################################################################
caps:=function(po,L)
    local i,j,ans,add;
    ans:=[];
    for i in po do
        if i[1]=L[1] then 
            add:=true;
            for j in [2..Length(L)] do
                if not [L[j],i[2]] in po then add:=false;fi;
            od;
            if add then Add(ans,i[2]);fi;
        fi;
    od;
    return(ans);
end;
                
coverPairs:=function(po,n)
    local i,j,l,ca,c,cp,s;
    cp:=[];s:=[];
    for i in po do if i[1]=n then Add(s,i[2]);fi;od;
    s:=Set(s);
    l:=Length(s);
    for i in [1..l-1] do
        for j in [i+1..l] do
            ca:=caps(po,[s[i],s[j]]);
            if not ca=[] then 
                for c in ca do Add(cp,[s[i],s[j],c]);od;
            fi;
        od;
    od;
    return(Set(cp));
end;

findTris:=function(L)
    local i,j,k,goodtri,badtri,l,e,f;
    goodtri:=[];
    badtri:=[];
    l:=Length(L);
    for i in [1..l-1] do
        for j in [i+1..l] do
            if L[i][2]=L[j][1] then 
                for k in L do
                    if (k[1]=L[i][1]) and (k[2]=L[j][2]) then 
                        e:=Set([L[i][1],L[i][2],L[j][2]]);
                        f:=Set([L[i][3],L[j][3],k[3]]);
                        if Length(f)=3 then Add(badtri,[e,f]);fi;
                    fi;
                od;
            fi;
        od;
    od;
    return(Set(badtri));
end;

testTri:=function(Tr,po)
    local i,bad,new;
    bad:=[];
    for i in Tr do
        if caps(po,i[2])=[] then 
            new:=[Position(R,V[i[1][1]]),Position(R,V[i[1][2]]),Position(R,V[i[1][3]])];
            Add(bad,Set(new));
        fi;
    od;
    return(Set(bad));
end;

reflectionType:=function(n) return(Set(Cycle(wc,n))[1]);end;

pattern:=function(L)
    local i;
    return(List([1..Length(L)],i->reflectionType(L[i])));
end;

#####
# minCycle is to find the minimal cycle of the given one up to numeric order.
# minCycle2 has the position of the 1st coordinate in the original cycle.
# ex) gap> minCycle2([3,2,7,1,3]);
#     [ 4, [ 1, 3, 3, 2, 7 ] ]
#####
minCycle:=function(L)
    local i,j,l,LL,m;
    l:=Length(L);
    LL:=List([1..l],i->L[i]);
    for i in [1..l] do LL[i+l]:=L[i];od;
    m:=0;
    for i in [1..l] do
        j:=1;
        while LL[j+m]=LL[j+i] do
            j:=j+1;
            if j>l then return(List([1..l],i->LL[m+i]));fi;
        od;
        if LL[j+m]>LL[j+i] then m:=i;fi;
    od;
    return(List([1..l],i->LL[m+i])); 
end;

minCycle2:=function(L)
    local i,j,l,LL,m;
    l:=Length(L);
    LL:=List([1..l],i->L[i]);
    for i in [1..l] do LL[i+l]:=L[i];od;
    m:=0;
    for i in [1..l] do
        j:=1;
        while LL[j+m]=LL[j+i] do
            j:=j+1;
            if j>l then return(List([1..l],i->LL[m+i]));fi;
        od;
        if LL[j+m]>LL[j+i] then m:=i;fi;
    od;
    return([m+1,List([1..l],i->LL[m+i])]);
end;



##########################################
# This is the main function to test CAT(0) 
# for 4 generator Artin groups
##########################################
testGroup:=function()
    local i;
    if not type[2]=4 then 
        Print("\nThis is not rank 4\n\n");
        return;
    fi;
    badTri:=testTri(findTris(coverPairs(P,1)),P);
    badPatterns:=[];
    for i in badTri do Add(badPatterns,minCycle(pattern(i)));od;
    badPatterns:=Set(badPatterns);
    if badTri=[] then 
        Print("\nThe type (",type[1],",4) Brady-Krammer complex is CAT(0)\n");
        Print("using the standard metric.\n\n");
        return;
    else 
        Print("\n  The type (",type[1],",4) Brady-Krammer complex is not CAT(0)\n");
        Print("  using the standard metric because some triples of reflections\n");
        Print("  give rise to short geodesic cycles in the unique vertex link\n\n");
        Print("  There are ",Length(badTri)," bad triples of reflections\n");
        Print("  There are ",Length(badPatterns)," bad patterns of reflections\n\n");
        Print("  To see the bad triples type badTri;\n");
        Print("  To see the bad patterns type badPatterns;\n\n");
    fi;
end;

####################################
# Part III: Create Paths and Columns 
####################################


#####
# createVertexPaths(Length(V)) generates all paths with vertex names such as 
#       [[1,8,15,38,50],...]
# createPathLabels(Length(V)) generates all path with reflection names such as
#       [[1,2,3,7],[1,2,4,3],...]
# createVertexPaths(Length(V))[i] & createPathLabels(Length(V))[i] DO NOT match.
#####

createVertexPaths:=function(n)
    local lev,ans,i;
    if n<=clen[2] then return([[1,n]]);fi;
    lev:=type[2];
    while n<=clen[lev] do lev:=lev-1;od;
    ans:=[];
    for i in [clen[lev-1]+1..clen[lev]] do
        if [i,n] in P then Append(ans,createVertexPaths(i));fi;
    od;
    for i in [1..Length(ans)] do Add(ans[i],n);od;
    return ans;
end;

createPathLabels:=function(n)
    local i,j,ans,p;
    paths:=createVertexPaths(n);
    ans:=[];
    for i in [1..Length(paths)] do
        ans[i]:=[];
        for j in [2..Length(paths[i])] do
            Add(ans[i],Position(R,V[paths[i][j-1]]^-1*V[paths[i][j]]));
        od;
    od;
    return Set(ans);
end;

# This creates the vertex paths with the reflections
createPathLabels2:=function(n)
    local i,j,ans,p;
    paths:=createVertexPaths(n);
    ans:=[];
    for i in [1..Length(paths)] do
        ans[i]:=[];
        ans[i][1]:=paths[i];
        ans[i][2]:=[];
        for j in [2..Length(paths[i])] do
            Add(ans[i][2],Position(R,V[paths[i][j-1]]^-1*V[paths[i][j]]));
        od;
    od;
    return Set(ans);
end;

shiftCol:=function(s)
    local i,j,ans;
    ans:=List([2..Length(s)],j->s[j]);
    Add(ans,s[1]^wc);
    return(ans);
end;

columnNumbers:=function(L)
    local i;
    return(List([1..Length(L)],i->paths[L[i]][1]));
end;


######################################
# This is the function which finds the
# columns in the Brady-Krammer complex
######################################
createColumns:=function() 
    local shift,ans,cy,i;
    paths:=createPathLabels(Length(V));
    shift:=List([1..Length(paths)],i->Position(paths,shiftCol(paths[i])));
    col:=PermList(shift);
    cy:=CycleStructurePerm(col);
    columns:=[];
    for i in [1..Length(paths)] do Add(columns,minCycle(columnNumbers(Cycle(col,i))));od;
    columns:=Set(columns);
    colPatterns:=[];
    for i in columns do Add(colPatterns,minCycle(pattern(i)));od;
    colPatterns:=Set(colPatterns);
    if verbose then
    Print("\n");
    Print(" There are ",Length(paths)," simplices\n");
    for i in [1..Length(cy)] do
        if IsBound(cy[i]) then
            if cy[i]>1 then
                Print(" There are ",cy[i]," columns each containing ",
                      i+1," simplices\n");
            else
                Print(" There is 1 column containing ",i+1," simplices\n");
            fi;
        fi;
    od;
    Print(" There are ",Length(colPatterns)," column patterns\n");
    Print("\nType paths; to see the reflections creating the simplices\n");
    Print("Type columns; to see the sequences creating the columns\n");
    Print("Type colPatterns; to see the list of column patterns\n");
    Print("\n");
    
    Print("If this type has diagram symmetries, the column patterns can\n");
    Print("be simplified by typing simplifyPatterns(s); where s is a\n");
    Print("list of length ",type[2]," that the i-th generater should\n");
    Print("be replaced with.  The list of patterns returned will be\n");
    Print("cyclically minimal and without repetitions.\n\n");
    fi;
end;

simplify:=function(L,s)
    local i,j,ans;
    ans:=List([1..Length(L)],i->[]);
    for i in [1..Length(L)] do 
        for j in L[i] do 
            Add(ans[i],s[j]);
        od;
    od;
    return ans;
end;

simplifyPatterns:=function(s)
    local i,ans;
    ans:=simplify(colPatterns,s);
    for i in [1..Length(ans)] do
        ans[i]:=minCycle(ans[i]);
    od;
    return(Set(ans));
end;



###########################
# Part VI create tetrahedra
###########################

### z is the diagram symmetries.
z:=function(t)
    local ans;
    if t="A" then
        ans:=[1,2,2,1];
    elif t="B" then
        ans:=[1,2,3,4];
    elif t="D" then 
        ans:=[1,1,2,1];
    elif t="F" then
        ans:=[1,2,2,1];
    elif t="H" then
        ans:=[1,2,3,4];        
    fi;
    return(ans);
end;


## simpPatterns generates minimal cycles of the given cycles of repeating 4 digits
#  with the original position of the first coordinate of the output.
#  ex : gap> simpPatterns([[2,3,1,4,2,3,1,4],[2,2,1,4,2,2,1,4]]);
##     " [ [ 3, [ 1, 4, 2, 3, 1, 4, 2, 3 ] ], [ 3, [ 1, 4, 2, 2, 1, 4, 2, 2 ] ]" 
simpPatterns:=function(b)
    local i,ans;
    ans:=List([1..Length(b)],i->[]);
    for i in [1..Length(b)] do
        ans[i]:=[minCycle2([b[i][1],b[i][2],b[i][3],b[i][4]])[1],minCycle(b[i])];
    od;
    return(ans);
end;

### CreateDiagonalSets gives list of sets for diagonals of each columns
# this is a subloop for createCapset.
createDiagonalSets:=function()
    local i,j,k,a,b,c,d,l;  
    verbose:=false;
    createColumns();
    verbose:=true;
    a:=List([1..Length(columns)],i->pattern(columns[i]));
    b:=simplify(a,z(type[1]));
    c:=simpPatterns(b);
    d:=List([1..2*(Length(columns))],i->[]); # couples for the diagonal in coulumn
    for i in [1..Length(columns)] do
        if (Length(columns[i]) mod 2)=0 then
            l:=Length(columns[i])/2;
            if c[i][1] in [1,3] then        
                for j in [1..l] do
                    Add(d[2*i-1],[columns[i][2*j-1],columns[i][2*j]]);   
                od;
                for k in [1..l-1] do
                    Add(d[2*i],[columns[i][2*k],columns[i][2*k+1]]);
                od;     
                Add(d[2*i],[columns[i][Length(columns[i])],columns[i][1]]); 
            fi;   
            if c[i][1] in [2,4] then        
                for j in [1..l] do
                    Add(d[2*i],[columns[i][2*j-1],columns[i][2*j]]);   
                od;
                for k in [1..l-1] do
                    Add(d[2*i-1],[columns[i][2*k],columns[i][2*k+1]]);
                od;   
                Add(d[2*i-1],[columns[i][Length(columns[i])],columns[i][1]]);
            fi; 
        elif (Length(columns[i]) mod 2)=1 then
            l:=Length(columns[i]);
            for j in [1..l-1] do
                Add(d[2*i-1],[columns[i][j],columns[i][j+1]]);
            od;
            Add(d[2*i-1],[columns[i][l],columns[i][1]]);
            for k in d[2*i-1] do
                Add(d[2*i],k);
            od;
        fi;
    od;    
    return([a,b,c,d]);
end;


### this creates capset which is the list whose ith coordinate is every numbers j such that 
# d[i] cap d[j] is nonempty. 
# this is a subloop for createDiagonals.
createCapset:=function()
     local d,i,j,capset;
    d:=createDiagonalSets()[4];
    capset:=List([1..2*(Length(columns))],i->[]);
    for i in [1..2*(Length(columns))] do
         for j in [1..2*(Length(columns))] do
             if Length(Intersection(d[i],d[j]))>0 then
                Add(capset[i],j);
             fi;
         od;
    od;
    return(capset);
end; 


# this function creates diagonals based on the intersection of 
# their diagonals without regarding column types.
# this is a subloop for createDiagonals.
createDia:=function()
     local i,j,k, 
           dia,
           m,
           capset,
           A,B,C,
           checked;
    capset:=createCapset();
    dia:=List([1..2*(Length(columns))],i->-1);
    A:=[1..Length(capset)];
    checked:=List([1..Length(capset)], i->0);   # not checked;
    m:=5;
    for i in A do
        if checked[i] = 0 then  # not-checked, thus do it.
            B:=capset[i];
            for j in capset[i] do
                C:=UnionSet(B,capset[j]);
                while Length(Difference(C,B))>0 do
                    for k in Difference(C,B) do
                        B:=C;
                        C:=UnionSet(B,capset[k]);
                     od;
                od;
            od;
            for j in B do
                checked[j]:=1;   # checked and assigned
                dia[j]:=m;
            od;
            m:=m+1;            
        fi;
    od;
    return(dia);
end;


# createTypeSets creates list of length L which maps each coordinate of L 
# to the coresponding element of M.
# this is a subloop for createDiagonals. 
createTypeSets:=function(M,L)
   local se,i,j;
    se:=List([1..Length(L)],i->[]);
    for i in [1..Length(columns)] do
        for j in [1..Length(L)] do
            if M[i]=L[j] then
               Add(se[j],i);
            fi;    
        od;
    od;
    return(se);
end;



# this creates diagonals for each columns.
createDiagonals:=function()
    local dia,simple,diagonals,minE,minF,minI,minJ,diaJ,minK,diaK,minII,
          diagonals1,minN,diaset,diaset1,idx,
          a,b,c,d,e,f,i,j,k,n,p,q,r,s,u,x,y,l,m,t,
          A,B,C,D,E,F,H,I,J,K,L,M,N,O;

    dia:=createDia();        
    simple:=simplifyPatterns(z(type[1]));    
    c:=createDiagonalSets()[3];               # columns in terms of generators with shift.
    diaset:=createDiagonalSets()[4];          # couples for each diagonal of each columns.
    diagonals:=StructuralCopy(dia);         
    A:=List([1..Length(columns)],i->0);       # index list for columns as patterns.
    for i in [1..Length(columns)] do
         for j in [1..Length(simple)] do
            if simple[j]=c[i][2] then
               A[i]:=j;
            fi;
         od;
    od;
    B:=[];C:=[];                              # B: columns with more than one generator.
    for k in [1..Length(simple)] do       # C: columns with one generator.
         if Length(Set(simple[k]))>1 then
            Add(B,k);            
         fi;
    od;
    Append(C,Difference([1..Length(simple)],B));
    idx := List([1..Length(columns)], i->0);
    H:=createTypeSets(A,C);                   # columns for each type with one generator.
    for x in [1..Length(H)] do
        t:=1;     
        while t>0 do
            diagonals1:=StructuralCopy(diagonals);   # for comparison.
            I:=[];                                   # I: for a tetra with one diagonal value.
            for y in H[x] do
                if diagonals1[2*y-1]=diagonals1[2*y] then
                   Add(I,diagonals1[2*y]); 
                fi;
    #Print("I =",I);   
            od;
            if Length(I)>0 then           # if I not empty, we get all same value for the type.
                minI:=Minimum(I);
                for l in H[x] do
                    for m in [1..Length(dia)] do
                        if diagonals[m]=diagonals1[2*l-1] then
                           diagonals[m]:=minI;                   
                        fi;
                        if diagonals[m]=diagonals1[2*l] then
                           diagonals[m]:=minI;
                        fi;   
                    od;
                od;
    #           Print(minI);
    #           Print(diagonals);
            fi; 
            if Length(I)=0 then                       
                for a in H[x] do         # finding same values and matching the pairs.    
                J:=[];K:=[];
                    for b in H[x] do
                        if diagonals[2*a-1]=diagonals[2*b-1] then
                           Append(J,[2*a,2*b]);
                        fi;
                        if diagonals[2*a-1]=diagonals[2*b] then
                           Append(J,[2*a,2*b-1]);
                        fi;
    #Print("J =",J);
                        if Length(J)>2 then
                             diaJ:=[];
                             for f in J do
                                  Add(diaJ,diagonals[f]);
                             od;
                             minJ:=Minimum(diaJ);
                             for d in J do
                                 for e in [1..Length(dia)] do
                                      if diagonals[e]=diagonals1[d] then
                                         diagonals[e]:=minJ;
                                      fi;
                                 od;                           
                             od; 
                             J:=[J[1],J[2]];  
                        fi;
                        diagonals1:=StructuralCopy(diagonals);     
                        if diagonals[2*a]=diagonals[2*b-1] then
                           Append(K,[2*a-1,2*b]);
                        fi;
                        if diagonals[2*a]=diagonals[2*b] then
                           Append(K,[2*a-1,2*b-1]);
                        fi;
    #Print("J =",J);
                        if Length(K)>2 then
                             diaK:=[];
                             for f in K do
                                  Add(diaK,diagonals[f]);
                             od;
                             minK:=Minimum(diaK);
                             for d in K do
                                 for e in [1..Length(dia)] do
                                      if diagonals[e]=diagonals1[d] then
                                         diagonals[e]:=minK;
                                      fi;
                                 od;                           
                             od;
                             K:=[K[1],K[2]];
                        fi;
                        diagonals1:=StructuralCopy(diagonals);     
                    od;
                od;
            fi;  
                   
            if diagonals=diagonals1 then     
               t:=0;
            fi;
            for i in [1..Length(H)] do
                L:=[];M:=[];
                for j in H[i] do
                    L:=[diagonals[2*j-1],diagonals[2*j]];
                    Add(M,Length(Set(L)));          
                od;
                if Maximum(M)-Minimum(M)>0 then   # in case one pattern has more 
                    t:=1;                         # than one type set of diagonals.
                fi;
            od;
        od;
        N:=[];                           # in case same type of tetra with diff. order
        for i in H[x] do                 # of diagonals, we can change them.
            Add(N,diagonals[2*i-1]);        
        od;
        if Length(Set(N))>1 then
            minN:=Minimum(N);
            O:=List([1..Length(H[x])],i->0);
            for i in H[x] do 
                if diagonals[2*i-1]<>minN then
                   O[Position(H[x],i)]:=1;
                fi;
            od;
            diaset1:=StructuralCopy(diaset);
            for j in [1..Length(O)] do
                if O[j]=1 then
                    diaset[2*H[x][j]-1]:=diaset1[2*H[x][j]];           
                    diaset[2*H[x][j]]:=diaset1[2*H[x][j]-1];
                    diagonals[2*H[x][j]-1]:=diagonals1[2*H[x][j]];
                    diagonals[2*H[x][j]]:=diagonals1[2*H[x][j]-1];
                    idx[H[x][j]]:=1;
                fi;               
            od;              
        fi;
    od; 
    D:=createTypeSets(A,B);        # columns with more than one generator.
    for n in [1..Length(D)] do
        E:=[];
        F:=[]; 
        for p in D[n] do            
            Add(E,diagonals[2*p-1]);
            Add(F,diagonals[2*p]); 
        od; 
        if Length(Intersection(E,F))>0 then
             Append(E,F);
             minE:=Minimum(E);
             for q in E do
                for r in [1..Length(dia)] do
                   if diagonals[r]=q then
                      diagonals[r]:=minE;
                   fi;
                od;
             od;
        elif Length(Intersection(E,F))=0 then
             minE:=Minimum(E);
             minF:=Minimum(F);
             for q in E do
                 for r in [1..Length(dia)] do
                    if diagonals[r]=q then
                       diagonals[r]:=minE;
                    fi;
                 od;
             od;  
             for s in F do
                 for u in [1..Length(dia)] do
                    if diagonals[u]=s then
                       diagonals[u]:=minF;
                    fi;
                 od;
             od;  
        fi;
    od; 
    for i in [1..Length(columns)] do
        if diagonals[2*i-1]<>diagonals[2*i] and idx[i]=1 then
            c[i][1]:=2;
        fi;
    od;
    return([c,A,B,C,D,H,diagonals,diaset,idx]);
end;


# columns with diagonals and tetrahedra with all labels. 
createSpces:=function()
    local c,diag, diagonals,A,i,j,k,l,
          ntetra, nDiag;
    diag := createDiagonals();
    c := diag[1];
    diagonals := diag[7];
#    Print("diagonals =",diagonals, "\n");
    A:=List([1..Length(columns)],i->[]);    
    for i in [1..Length(columns)] do
       for j in [1..4] do
           A[i][j]:=c[i][2][j];
       od;
       for k in [5,6] do
           A[i][k]:=diagonals[(2*i)-6+k];
       od;
   od;
    tetrahedra := [];
    ntetra:=Length(Set(diag[2]));
    nDiag := Length(diag[2]);
    for l in [1..ntetra] do
        # find the first occurence of l in diag[2]
        for k in [1..nDiag] do
            if l = diag[2][k] then
                break;
            fi;
        od;
        Add(tetrahedra, A[k]);
    od;
    return([A,tetrahedra]);
end;

createNodes1:=function()
    local A,B,C,D,E,N1,
          i,j;
    A:=Set(z(type[1]));
    B:=[];
    for i in A do
        Add(B,Position(z(type[1]),i));
    od;
    D:=List([1..Length(B)],i->[]);
    for i in A do
        for j in [1..4] do
            if z(type[1])[j]=i then
                Add(D[i],j);
            fi;
        od;
    od;
    E:=List([1..Length(D)],i->[]);
    for i in [1..Length(D)] do
        for j in [1..Length(D[i])] do
            E[i][j]:=Position(V,R[D[i][j]]);
        od;
    od;
    N1:=[];
    C:=[];
    for j in B do
        Add(C,Position(V,R[j]));
    od;
    return([C,B,D,E]);
end;



basics:=function()
    local a,A,i,j;
    a:=Cycles(wconjugates(R,w),[1..Length(R)]);
    A:=[0,0,0,0];
    for i in [1..4] do
        for j in [1..Length(a)] do
            if i in a[j] then A[i]:=j; fi;
        od;
    od;
    return([a,A]);
end;

genSet:=function()
    local a,A,D,i,j,k,cyc;
    D:=createNodes1()[3];
    a:=basics()[1];
    A:=basics()[2];
    cyc:=List([1..Length(D)],i->[]);
    for i in [1..Length(D)] do
        for j in D[i] do
            for k in a[A[j]] do
                Add(cyc[i],k);
            od;
        od;
    od;
    return(cyc);
end;

genSetNodes:=function()
    local cyc,nodesSet1,nodesSet2,i,j; 
    cyc:=genSet();
    nodesSet1:=StructuralCopy(cyc);
    nodesSet2:=StructuralCopy(cyc);
    for i in [1..Length(cyc)] do
        for j in [1..Length(cyc[i])] do
            nodesSet1[i][j]:=Position(V,R[cyc[i][j]]);
        od;
    od;
    for i in [1..Length(cyc)] do
        for j in [1..Length(cyc[i])] do
            nodesSet2[i][j]:=Position(V,w*R[cyc[i][j]]);
        od;
    od;
    return([nodesSet1,nodesSet2]);
end;


# subloop for nodes2       
mutableCopy:=function(obj)
    local mcp, i;
    
    mcp := [];
    if IsList(obj) then
        for i in [1..Length(obj)] do
            mcp[i] := mutableCopy(obj[i]);
        od;
    else
        mcp := obj;
    fi;
    
    return(mcp);
end;


# subloop for createNodes2
nodes2:=function()
    local a,b,c,i,j,A;
    a:=wconj2(R,w);
    A:=Cycles(a,[1..(Length(V)-2-(2*(Length(R))))]);
    b:=mutableCopy(A);
    for i in [1..Length(A)] do
        for j in [1..Length(A[i])] do
            b[i][j]:=A[i][j]+Length(R)+1;
        od;
    od;
    c:=[];
    for i in [1..Length(b)] do
        Add(c,b[i][1]);
    od;
    return([b,c]);
end;


createNodes2:=function()
    local N,A,B,C,
          i,j,k;
    N:=nodes2()[2];
    C:=[];
    for i in N do
        A:=[]; B:=[];
        for j in [1..(Length(P)-2*Length(R))/2] do
            if i=P[j+Length(R)][2] then
                Add(A,P[j+Length(R)][1]);
            fi;
        od;
        for k in [1..(Length(P)-2*Length(R))/2] do
            if i=P[Length(P)/2+k][1] then
                Add(B,P[Length(P)/2+k][2]);
            fi;
        od;
        C[Position(N,i)]:=[A,B];
    od;
    return([N,C]);
end;



elimNodes :=function()
    local nodes,edges,nodesSet1,nodesSet2,edgetypes,typeSet,typeSet2,typeSet3,Nodes2,
          i,j,k,l,numGen,gptype;
    gptype := type[1];
    numGen := Length(Set(z(gptype)));
    nodes:=createNodes2()[1];
    edges:=createNodes2()[2];
    nodesSet1:=genSetNodes()[1];
    nodesSet2:=genSetNodes()[2];
    edgetypes:=[];
    for i in [1..Length(nodes)] do 
        if gptype in ["A","D","F"] then
            edgetypes[i] := [[0,0],[0,0]];
        elif gptype in ["B","H"] then
            edgetypes[i]:= [[0,0,0,0],[0,0,0,0]];
        fi;
    od;
    for i in [1..Length(nodes)] do
        for j in [1,2] do           
            for l in edges[i][j] do
                for k in [1..numGen] do
                    if l in nodesSet1[k] then
                    edgetypes[i][1][k]:=edgetypes[i][1][k]+1;                
                    elif l in nodesSet2[k] then
                    edgetypes[i][2][k]:=edgetypes[i][2][k]+1;
                    fi;
                od;
            od;
        od;
    od;   
    typeSet:=Set(edgetypes);
    typeSet2:=[];
    for i in [1..Length(typeSet)] do
        Add(typeSet2,Set(typeSet[i]));
    od;
    typeSet3:=Set(typeSet2);
    Nodes2:=[];
    for i in [1..Length(typeSet3)] do
        if Length(typeSet3[i])=2 then
            Add(Nodes2,nodes[Position(edgetypes,typeSet3[i])]);
        elif Length(typeSet3[i])=1 then
            Add(Nodes2,nodes[Position(edgetypes,[typeSet3[i][1],typeSet3[i][1]])]);
        fi;
    od;
    return([edgetypes,typeSet,typeSet2,typeSet3,Nodes2]);
end;

nodesCheck:=function()
    return(Concatenation(createNodes1()[1],elimNodes()[5]));
end;



#
# find all nodes which are connected to 's' in a graph, 'graph',
#    and returns nodes in the first element (= b)
#    and their locations defined in 'graph' in the second element (= idx_graph)
#
findEdges := function(s, graph)
    local n, i, b, idx_graph;
    n := Length(graph[1]);  # edges
    b := [];
    idx_graph := [];
    for i in [1..n] do
        if graph[2][i] = 0 then  # not used
            if graph[1][i][1] = s then
                Add(b, graph[1][i][2]);
                Add(idx_graph, i);
            elif graph[1][i][2] = s then
                Add(b, graph[1][i][1]);
                Add(idx_graph, i);
            fi;              
        fi;
    od;
    return([b, idx_graph]);
end;

#
# as a recursive function, this will find all loops containing 'path'
#    and store them in 'paths'.
#
# 'path' is a collection of nodes traveled and can be just a starting node.
#
findClosedLoop := function(path, paths, graph, nVertices)
    local s_edges, n, i, ret, npath;
    
    # traveled too far, no loop.
    if nVertices+1 = Length(path) then
        return(-1);
    fi;
    
    # traveled so far
    npath := Length(path);
    
    # next nodes available -> edges
    s_edges := findEdges(path[npath], graph);
    n := Length(s_edges[1]);
    for i in [1..n] do
        
        # move forward a step
        Add(path, s_edges[1][i]);
        
        # and see if a loop is found
        if path[1] = s_edges[1][i] then   # loop found
            Append(paths, [path]);
            return(0);
        else
            graph[2][s_edges[2][i]] := 1;  # used
            ret := findClosedLoop(path, paths, graph, nVertices);
            graph[2][s_edges[2][i]] := 0;  # restore
        fi;
        
        # move backward a step (to its original location)
        path := path{[1..npath]};
    od;
    
    return(-1);
end;
findClosedLoop2 := function(path, paths, graph, nVertices)
    local s_edges, n, i, ret, npath;
    
    # traveled too far, no loop.
    if nVertices+1 = Length(path) then
        return(-1);
    fi;
    

    # traveled so far
    npath := Length(path);
    
    # next nodes available -> edges
    s_edges := findEdges(path[npath], graph);
    n := Length(s_edges[1]);
    for i in [1..n] do
        
        # move forward a step
        Add(path, s_edges[1][i]);
        
        # and see if a loop is found
        if path[1] = s_edges[1][i] then   # loop found
            Append(paths, [path]);
            return(0);
        elif npath+1 > Length(Set(path)) then
            return(-1);
        else
#            Add(vtces, s_edges[1][i]);
            graph[2][s_edges[2][i]] := 1;  # used
            ret := findClosedLoop2(path, paths, graph, nVertices);
            graph[2][s_edges[2][i]] := 0;  # restore
        fi;
        
        # move backward a step (to its original location)
        path := path{[1..npath]};
    od;
    
    return(-1);
end;

#
# find All Closed Loops from a given graph
#   need to simplify using minCycle...
#
findAllClosedLoops := function(graph)
    local i, path, paths, nVertices, vertices, graph_flag, grp, nGrp;
    
    grp := graph[1];  # only graph
    nGrp := Length(grp);

    graph_flag := List([1..nGrp], i->0);  # traveld or not    
    
    vertices := Set(Flat(grp));    # all vertices
    nVertices := Length(vertices);
    vertices := Set(Flat(grp{[1..nGrp]}{[1]}));  # check only a side of edge
    
#    vtces:=[];
    paths := [];
    for i in vertices do
        findClosedLoop2([i], paths, [grp, graph_flag], nVertices+1);
    od;
    
    return paths;
end;


# subloop for findGraphForNode.
# pathLabelsToGraph will give a graph 'graphs' for the specific node and 
# find spces for all edge in the graph 'spcs'   
pathLabelsToGraph := function(pathLabels, nodeNum, numGen)
    local i, n, graphs, spcs;
    
    n := Length(pathLabels);
    if (2 <= nodeNum) and (nodeNum <= numGen+1) then  # in rank one
        graphs := [];
        spcs := [];
        for i in [1..n] do
            if pathLabels[i][1][2] = nodeNum then
                Add(graphs, pathLabels[i][1]{[3,4]});
                Add(spcs, pathLabels[i][2]);
            fi;
        od;
    else  # in rank two
        graphs := [];
        spcs := [];
        for i in [1..n] do
            if pathLabels[i][1][3] = nodeNum then
                Add(graphs, pathLabels[i][1]{[2,4]});
                Add(spcs, pathLabels[i][2]);
            fi;
        od;        
    fi;
      
    return( [graphs, spcs] );
end;


#
# find a graph for a node
#
findGraphForNode := function(node)
    return(pathLabelsToGraph(createPathLabels2(Length(V)), node, Length(R)));
end;

#
# find all closed loops for a given node
#    simply find a graph for the node and call findAllClosedLoops();
#
findAllClosedLoopsForNode := function(node)
    return(findAllClosedLoops(findGraphForNode(node)));
end;

#
# simpify loops by removing repetions and only taking simple loops
#
simplifyLoops := function(loops)
    local nLoops, minLoops, simplifiedLoops, i, k;
    
    nLoops := Length(loops);
    
    minLoops := [];
    for i in [1..nLoops] do
        
        # breaking loops by removing the closing vertice
        k := Length(loops[i])-1;
        
        Add(minLoops, minCycle(loops[i]{[1..k]}));
    od;
    minLoops := Set(minLoops);   # no repeated loops
    
    nLoops := Length(minLoops);
    simplifiedLoops := [];
    for i in [1..nLoops] do
        # find only simple loops
        if Length(minLoops[i]) = Length(Set(minLoops[i])) then
            Add(minLoops[i], minLoops[i][1]);  # closing back the loop
            Add(simplifiedLoops, minLoops[i]);
        fi;
    od;
    
    return simplifiedLoops;
end;

#
#
#
findAllSimpleClosedLoopsForNode := function(node)
    return(simplifyLoops(findAllClosedLoopsForNode(node)));
end;
    

#
# this is to find a match of Simplex responding to 
#    each edge of a link in a column and its location
#
findSubMatch:=function(spx, dblCol)
    local nspx, ndblCol, i, j, matched;    
    
    nspx := Length(spx);
    ndblCol := Length(dblCol);
    
    for i in [1..(ndblCol-nspx+1)] do
        if spx[1] = dblCol[i] then   # first match
            matched := 1;
            for j in [2..nspx] do    #  complete match
                if spx[j] <> dblCol[i+j-1] then  # match broken
                    matched := 0;
                    break;
                fi;
            od;
            if matched = 1 then
                return (i);
            fi;
        fi;
    od;
    return (0);  # not found
end;

#
# find spx from list of columns
#   and return which column and its location in the column
#
findMatch:=function(spx, cols)
    local nCols, i, dblCol, pos;

    nCols := Length(cols);
    for i in [1..nCols] do
        dblCol := Flat([cols[i], cols[i]]);  # makes double column
        pos := findSubMatch(spx, dblCol);
        # pos := PositionSublist(dblCol, spx);
        if pos > 0 then  # found
            return([i, pos]);
        fi;
    od;
end;  

#
# find label (=[type, location]) for each edge
#
findEdgeLabel:=function(edge, graphs, cols)
    local i, col_loc, diag, spx, pos;
    
    # find spx corresponding to edge in graphs
    for i in [1..Length(graphs[1])] do
        if Set(edge) = Set(graphs[1][i]) then
            spx := graphs[2][i];
        fi;
    od;
    
    # findMatch -> (col, loc)
    col_loc := findMatch(spx, cols);  # col_loc[1] = col, col_loc[2] = location
    
    diag := createDiagonals();
    
    pos := (1 + col_loc[2] - diag[1][col_loc[1]][1] - 1) mod 4 + 1;
    
    return([diag[2][col_loc[1]], pos]);
    
end;


#
#
#
findLabelForEdge := function(edge, els)
    local i;
    # find label corresponding to edge
    for i in [1..Length(els)] do
        if Set(edge) = Set(els[i][1]) then
            return(els[i][2]);
        fi;
    od;
end;


#
#
#
findAllEdgeLabels:=function(node)
    local grp, nGrp, graphs, i, eLabel;
    
    verbose := false;
    createColumns();
    verbose := true;
    
    grp := findGraphForNode(node);
    nGrp := Length(grp[1]);
    
    graphs := [];
    Print("nGraph = ", nGrp, "\n");
    for i in [1..nGrp] do  # for each edge
        eLabel := findEdgeLabel(grp[1][i], grp, columns);
        
        # if node is in rank 1
        if node <= Length(R) + 1 then
            Add(graphs, [grp[1][i], eLabel{[1,2]}]);
        else
            if eLabel[2] mod 2 = 0 then  # even
                Add(graphs, [grp[1][i], [eLabel[1], 6]]);
            else                         # odd
                Add(graphs, [grp[1][i], [eLabel[1], 5]]);
            fi;
        fi;
        
        Print(".\c");
    od;
    Print("\n");
    
    return(graphs);
end;
    

simpleLoopToEdgeLabels := function(loop, els)
    local nEdges, lbl, edge, i;
    
    nEdges := Length(loop);
    
    lbl := [];
    for i in [1..(nEdges-1)] do
        edge := loop{[i, i+1]};
        Add(lbl, findLabelForEdge(edge, els));
    od;
    
    return (lbl);
end;


simpleLoopsToEdgeLables := function(loops, els)
    local nLoops, EL, i;

    nLoops := Length(loops);
    
    EL := [];
    for i in [1..nLoops] do
        Add(EL, simpleLoopToEdgeLabels(loops[i], els));
    od;
    
    return (EL);
end;


countEdgeLabels := function(el)
    local nEL, tetraedges, countEL, i;
    
    nEL := Length(el);
   
    # constant
    tetraedges := 6;  
    
    countEL := List([1..nPatterns * tetraedges], i->0);
    for i in [1..nEL] do
        countEL[(el[i][1]-1)*tetraedges + el[i][2]] := countEL[(el[i][1]-1)*tetraedges + el[i][2]]+1;
    od;
    
    return countEL;
end;


countAllEdgeLabels := function(EL)
    local i, nEL, countELs;
    
    nEL := Length(EL);
    countELs := [];
    for i in [1..nEL] do
        Add(countELs, countEdgeLabels(EL[i]));
    od;
    
    return Set(countELs);
end;

findAllIneqes := function()
    local nodes, allCountELs, graph, loops, ELs, countELs, i,Patterns;
    
    nodes := nodesCheck();   
    createColumns();
    createSpces();
    nPatterns:=Length(tetrahedra);

    Print(" There are ", nPatterns," simplices.\n"); 
    Print(" The simplices are\n ",tetrahedra,".\n\n");
    
    Print(" There Are ",Length(nodes)," nodes to check which are ",nodes,".\n\n");

    allCountELs := [];
    for i in [1..Length(nodes)] do
        graph := findGraphForNode(nodes[i]);
        loops := findAllSimpleClosedLoopsForNode(nodes[i]);
        ELs := simpleLoopsToEdgeLables(loops, findAllEdgeLabels(nodes[i]));
        countELs := countAllEdgeLabels(ELs);
        allCountELs := UnionSet(countELs, allCountELs);
    od;

    return allCountELs;
end;

# compareList compares two elements r and s 
# and gives true when r[i] is not more than s[i] for all i. 
# subloop of elimRedundant.
# that is, s is redundant in our case.
compareList := function(r,s)
    local i;
    
    for i in [1..Length(r)] do
        if r[i]>s[i] then
            return false;
        fi;
    od;
    return true;    
end;

# elimRedundant gets rid of redundant elements from L using compareList.
# subloop of simplifiedIneqes.
elimRedundant := function(L) 
    local simplified, nL,flag, i, j;
    
    simplified :=[];
    nL := Length(L);
    flag := List([1..nL],i->0);
    
    for i in [1..nL] do
        for j in [1..nL] do 
            if  (not i=j) and compareList(L[j],L[i]) then
                if not L[i]=L[j] then                
                    flag[i] :=j;
                    break;
                elif L[i]=L[j] then
                    if i>j then
                        flag[i] := j;
                        break;
                    fi;
                fi;
            fi;
        od;
    od;
    for i in [1..nL] do 
        if flag[i] =0 then
            Add(simplified,L[i]);
        fi;
    od;
    return ([flag,simplified]);
end;

# simplifiedIneqes returns inequations after eliminate redundants.
simplifiedIneqes := function()
    local ineq;
    ineq := elimRedundant(findAllIneqes());
    Print("We have ",Length(ineq[2])," inequalities.\n\n");
    return ineq;
end;

# numEntries produces a sequence in which each entry indicates 
# the number of nonzero entries of corresponding entry of given list.
# subloop of sortList
numEntries := function(L)
    local i,j,ans;
    
    ans := List([1..Length(L)],i->0);
    for i in [1..Length(L)] do
        for j in L[i] do 
            if not j=0 then
                ans[i] := ans[i]+1;
            fi;
        od;
    od;
    return ans;
end;

# sortList produces 1)the number of elements with each number of nonzero 
#entries,  2) those elements sorted by the number of nonzero entries. 
# subloop of sortIneqes.
sortList := function(L,numEntries)
    local i,j, sortLength,sortNum,sortL;
    
    sortLength := Maximum(numEntries);
    sortNum := List([1..sortLength],i->[0,0]);
    sortL := List([1..sortLength],i->[]);
    
    for i in [1..sortLength] do
        sortNum[i][1] :=i;
    od;
    for i in [1..Length(L)] do
        for j in [1..sortLength] do
            if numEntries[i]=j then
                sortNum[j][2] := sortNum[j][2]+1;
                Add(sortL[j],L[i]);
            fi;
        od;
    od;
    
    return([sortNum,sortL]);
end;

# sortIneqes sorts inequations up to the number of nonzero entries.
sortIneqes := function()
    local ineqes;
    ineqes := simplifiedIneqes()[2];
    return(sortList(ineqes,numEntries(ineqes)));
end;



# dihedralAngles generates a sequence of two six-tuples which indicates 
# edges for the angle. subloop for simplifiedDihAng.
dihedralAngles := function(tetra)
    local t, dihAng, i ;
    
    t := tetra;
    dihAng := List([1..6],i->[]);
    
    Add(dihAng[1],[t[1],t[2],t[3],t[4],t[5],t[6]]); # this is the same way each tetrahedron counted. not same to Jon's way.
    Add(dihAng[2],[t[2],t[3],t[4],t[1],t[6],t[5]]);
    Add(dihAng[3],[t[3],t[4],t[1],t[2],t[5],t[6]]);
    Add(dihAng[4],[t[4],t[1],t[2],t[3],t[6],t[5]]);
    Add(dihAng[5],[t[5],t[3],t[6],t[1],t[4],t[2]]);
    Add(dihAng[6],[t[6],t[1],t[5],t[3],t[4],t[2]]);
    
    for i in [1..6] do  # edges for the reverse direction  which produces same angle.
        Add(dihAng[i],[dihAng[i][1][1],dihAng[i][1][4],dihAng[i][1][3],dihAng[i][1][2],dihAng[i][1][6],dihAng[i][1][5]]);
    od;
    
    return dihAng;
end;

# simplifiedDihAng is to find same dihedral angles by comparing 
# the surounding edges obtained from dihedralAngles. 
simplifiedDihAng := function(dihAng)
    local modifiedDihAng, simplifiedDihAng,
          player,index, rest,flag,
          i;
    
    modifiedDihAng:=[];
    
    for i in [1..6] do
        Add(modifiedDihAng,Set(dihAng[i]));
    od; 
    
    simplifiedDihAng := List([1..6],i->0);
    simplifiedDihAng[1]:=1;
    player := 1; # the standard one for matching
    index := 1;  
    rest := [2,3,4,5,6]; # not matched yet
    
    while  Length(rest)>0 do
        flag :=[]; 
        for i in rest do
            if modifiedDihAng[i]=modifiedDihAng[player] then
                simplifiedDihAng[i] := index;
                Add(flag,i);
            fi;
        od;
        rest := Difference(rest,flag);
        index := index+1;
        if Length(rest)>0 then
            player := Minimum(rest);
        fi;
    od;
            
    return simplifiedDihAng;
end;

# dihAngForGp produces a sequence of dihedral angles of tetrahedra 
# for the group using tetrahedra(tetrahedra is generated in createSpces)
dihAngForGp := function()
    local tetra, n, dihAngForGp, ang, max,
          i,j;
    
    tetra := tetrahedra;
    n := Length(tetra);
    dihAngForGp := List([1..6*n],i->0);
    
    for i in [1..n] do
        max := Maximum(dihAngForGp);
        for j in [1..6] do
            ang := simplifiedDihAng(dihedralAngles(tetra[i]));
            dihAngForGp[6*(i-1)+j] := ang[j]+max;
        od;
    od;
    
    Print("The dihedral angles for the group is  ",dihAngForGp," .\n\n");
    return dihAngForGp;
end;

reducedIneq := function(ineq, dihAngForGp)
    local lengIneq, numAngles, reducedIneq,
          i, j;
    
    lengIneq := Length(ineq);
    numAngles := Maximum(dihAngForGp);
    reducedIneq := List([1..numAngles],i->0);
    
    for i in [1..numAngles] do
        for j in [1..lengIneq] do
            if dihAngForGp[j] = i then
                reducedIneq[i] := reducedIneq[i]+ineq[j];
            fi;
        od;
    od;
    
    return reducedIneq;
end;

        
reducedIneqes := function()
    local ineqes, numIneq, reducedIneqes, dihAngForTheGp, i;
    
    ineqes := simplifiedIneqes()[2];
    numIneq := Length(ineqes);
    reducedIneqes := [];
    dihAngForTheGp := dihAngForGp();
    
    for i in [1..numIneq] do
        Add(reducedIneqes, reducedIneq(ineqes[i],dihAngForTheGp));
    od;
    
    return reducedIneqes;
end;


simplifiedReducedIneqes := function()
    local ineq;
    ineq := elimRedundant(reducedIneqes());
#    ineq := elimRedundant(red);   
    Print("We have ",Length(ineq[2])," inequalities.\n\n");
    return ineq;
end;