##################################################################
##################################################################
##################################################################
####                                                          ####
####  GAP-Routine                                             ####
####                                                          ####
####                 `FundamentalGroup'                       ####
####                                                          ####
####  Version Nov/2003                                        ####
####  by Frank H. Lutz, TU Berlin, Germany                    ####
####                                                          ####
####  available at                                            ####
####                                                          ####
####     http://www.math.TU-Berlin.de/diskregeom/stellar/     ####
####                                                          ####
##################################################################
##################################################################
##################################################################


##################################################################
##################################################################
##                                                              ##
## The GAP routine  FundamentalGroup  computes a `small'        ## 
## presentation of the fundamental group of a given             ##
## simplicial complex.                                          ##
##                                                              ##
##################################################################
##################################################################


##################################################################
#                                                                #
# How to use the program?                                        #
#                                                                #
# (A) Install GAP on your computer.                              #
#                                                                #
# (B) Write the facets of a simplicial complex                   #
#     in the format                                              #
#                                                                #
#        facets:=[[1,2],[1,4],[2,3],[3,4]];                      #
#                                                                #
#     to a separate                                              #
#                                                                #
###############################################                  #
                                              #                  #
file:=String("FundamentalGroup.testobject");  #                  #
                                              #                  #
###############################################                  #
#                                                                #
# (C) Start GAP and begin computation with                       #
#                                                                #
#        Read("FundamentalGroup");                               #
#                                                                #
# (D) The resulting presentation is printed to the               #
#     output file                                                #
#                                                                #
#        `FundamentalGroup.out'                                  #
#                                                                #
#     Please set:                                                #
#                                                                #
##############                                                   #
             #                                                   #
verbose:=1;  # <--- 0 = output to file only                      #
             #      1 = output to screen and file                #
##############                                                   #
#                                                                #
# (E) gap> quit;                                                 #
#                                                                #
##################################################################



Read("FundamentalGroup.testobject");



### comute defining generators and relators ###

faces:=[];

for k in [1..3] do
    faces[k]:=[];
    for element in facets do
        UniteSet(faces[k],Combinations(element,k));
    od;
od;

tree:=[];
tree_vertices:=ShallowCopy(faces[1][1]);

while Length(tree_vertices) < Length(faces[1]) do
      for edge in faces[2] do
          if edge[1] in tree_vertices and not edge[2] in tree_vertices then
             AddSet(tree,edge);
             AddSet(tree_vertices,edge[2]);
          elif not edge[1] in tree_vertices and edge[2] in tree_vertices then
             AddSet(tree,edge);
             AddSet(tree_vertices,edge[1]);
          fi;
      od;
od;

SubtractSet(faces[2],tree);

PrintTo("FundamentalGroup.out","G:=FreeGroup(");

if Length(faces[2]) > 0 then 
   for k in [1..Length(faces[2])] do
       if k < Length(faces[2]) then
          AppendTo("FundamentalGroup.out","\"a_",faces[2][k][1],"_",faces[2][k][2],"\",");
       else
          AppendTo("FundamentalGroup.out","\"a_",faces[2][k][1],"_",faces[2][k][2],"\");\n");
       fi;
   od;
else
   AppendTo("FundamentalGroup.out","0);\n");
fi;

AppendTo("FundamentalGroup.out","H:=G/[");
for k in [1..Length(faces[3])] do
    flag_comma:=1;
    triangle:=ShallowCopy(faces[3][k]);
    triangle_edges:=Combinations(triangle,2);
    if triangle_edges[1] in tree then
       if triangle_edges[2] in tree then
          if not triangle_edges[3] in tree then
             AppendTo("FundamentalGroup.out","G.",Position(faces[2],triangle_edges[3]));
          else
             flag_comma:=0;
          fi;
       elif triangle_edges[3] in tree then
          AppendTo("FundamentalGroup.out","G.",Position(faces[2],triangle_edges[2]));
       else
          AppendTo("FundamentalGroup.out","G.",Position(faces[2],triangle_edges[3]),"/G.",Position(faces[2],triangle_edges[2]));
       fi;
    elif triangle_edges[2] in tree then
       if triangle_edges[3] in tree then
          AppendTo("FundamentalGroup.out","G.",Position(faces[2],triangle_edges[1]));
       else
          AppendTo("FundamentalGroup.out","G.",Position(faces[2],triangle_edges[1]),"*G.",Position(faces[2],triangle_edges[3]));
       fi;
    elif triangle_edges[3] in tree then
       AppendTo("FundamentalGroup.out","G.",Position(faces[2],triangle_edges[1]),"/G.",Position(faces[2],triangle_edges[2]));
    else
       AppendTo("FundamentalGroup.out","G.",Position(faces[2],triangle_edges[1]),"*G.",Position(faces[2],triangle_edges[3]),"/G.",Position(faces[2],triangle_edges[2]));
    fi;
    if k < Length(faces[3]) and flag_comma = 1 then
       AppendTo("FundamentalGroup.out",",");
    fi;
od;
AppendTo("FundamentalGroup.out","];");

Read("FundamentalGroup.out");

if verbose = 1 then
   Print("\n","Defining generators: ",Length(GeneratorsOfGroup(H)),"\n");
   Print("Defining relators: ",Length(RelatorsOfFpGroup(H)),"\n\n");
fi;



### compute small presentation ###

pi:=SimplifiedFpGroup(H);
#Print("Size = ",Size(pi),"\n");

PrintTo("FundamentalGroup.out","pi_free:=FreeGroup(");

if Length(GeneratorsOfGroup(pi)) > 0 then
   for k in [1..(Length(GeneratorsOfGroup(pi))-1)] do
       AppendTo("FundamentalGroup.out","\"pi.",k,"\",");
   od;
   AppendTo("FundamentalGroup.out","\"pi.",Length(GeneratorsOfGroup(pi)),"\");\n");
else
   AppendTo("FundamentalGroup.out","0);\n");
fi;

AppendTo("FundamentalGroup.out","pi:=pi_free/[\n");
for j in [1..Length(RelatorsOfFpGroup(pi))] do
    string:=String(RelatorsOfFpGroup(pi)[j]);
    split_string:=ShallowCopy(SplitString(string,"*"));
    for k in [1..Length(split_string)] do
        split_element:=ShallowCopy(SplitString(split_string[k],"^"));
        for l in [1..Length(GeneratorsOfGroup(pi))] do
            if split_element[1] = String(GeneratorsOfGroup(pi)[l]) then
               AppendTo("FundamentalGroup.out","pi_free.",l);
               if Length(split_element) = 2 then
                  AppendTo("FundamentalGroup.out","^",split_element[2]);
               fi;
               if k < Length(split_string) then
                  AppendTo("FundamentalGroup.out","*");
               fi;
            fi;
        od;
    od;
    if j < Length(RelatorsOfFpGroup(pi)) then
       AppendTo("FundamentalGroup.out",",\n");
    fi;
od;
AppendTo("FundamentalGroup.out","];");

Read("FundamentalGroup.out");



### compute homology group H_1 ###

PrintTo("FundamentalGroup.out","pi_free:=FreeGroup(");

if Length(GeneratorsOfGroup(pi)) > 0 then
   for k in [1..(Length(GeneratorsOfGroup(pi))-1)] do
       AppendTo("FundamentalGroup.out","\"pi_ab.",k,"\",");
   od;
   AppendTo("FundamentalGroup.out","\"pi_ab.",Length(GeneratorsOfGroup(pi)),"\");\n");
else
   AppendTo("FundamentalGroup.out","0);\n");
fi;

AppendTo("FundamentalGroup.out","pi_ab:=pi_free/[\n");
for k in [1..(Length(GeneratorsOfGroup(pi))-1)] do
    for l in [(k+1)..(Length(GeneratorsOfGroup(pi)))] do
        AppendTo("FundamentalGroup.out","pi_free.",k,"*pi_free.",l,"*pi_free.",k,"^-1*pi_free.",l,"^-1,\n");
    od;
od;
for j in [1..Length(RelatorsOfFpGroup(pi))] do
    string:=String(RelatorsOfFpGroup(pi)[j]);
    split_string:=ShallowCopy(SplitString(string,"*"));
    for k in [1..Length(split_string)] do
        split_element:=ShallowCopy(SplitString(split_string[k],"^"));
        for l in [1..Length(GeneratorsOfGroup(pi))] do
            if split_element[1] = String(GeneratorsOfGroup(pi)[l]) then
               AppendTo("FundamentalGroup.out","pi_free.",l);
               if Length(split_element) = 2 then
                  AppendTo("FundamentalGroup.out","^",split_element[2]);
               fi;
               if k < Length(split_string) then
                  AppendTo("FundamentalGroup.out","*");
               fi;
            fi;
        od;
    od;
    if j < Length(RelatorsOfFpGroup(pi)) then
       AppendTo("FundamentalGroup.out",",\n");
    fi;
od;
AppendTo("FundamentalGroup.out","];");

Read("FundamentalGroup.out");

H_1:=SimplifiedFpGroup(pi_ab);



### print simplified group to file ###

if verbose = 1 then
   Print("Generators (simplified) = ",GeneratorsOfGroup(pi),"\n");
   Print("Relators (simplified) = ",RelatorsOfFpGroup(pi),"\n\n");
fi;

if Length(GeneratorsOfGroup(H)) > Length(RelatorsOfFpGroup(H)) then
   max:=Length(GeneratorsOfGroup(H));
else
   max:=Length(RelatorsOfFpGroup(H));
fi;
count:=0;
while Int(max) > 0 do
      max:=Int(max/10);
      count:=count+1;
od;
PrintTo("FundamentalGroup.out","########################");
for k in [1..count] do
    AppendTo("FundamentalGroup.out","#");
od;
AppendTo("FundamentalGroup.out","\n#\n",
                               "# Defining generators: ",Length(GeneratorsOfGroup(H)),"\n",
                               "# Defining relators: ",Length(RelatorsOfFpGroup(H)),"\n#\n");
for k in [1..count] do
    AppendTo("FundamentalGroup.out","#");
od;
AppendTo("FundamentalGroup.out","########################\n\n\n\n#####\n#\n",
                               "# Generators (simplified) = ",GeneratorsOfGroup(pi),"\n");
if Length(RelatorsOfFpGroup(pi)) <= 1 then
   AppendTo("FundamentalGroup.out","# Relators (simplified) = ",RelatorsOfFpGroup(pi),"\n#\n");
else
   AppendTo("FundamentalGroup.out","# Relators (simplified) = [ ",RelatorsOfFpGroup(pi)[1],",\n");
   for k in [2..(Length(RelatorsOfFpGroup(pi))-1)] do
       AppendTo("FundamentalGroup.out","#                           ",RelatorsOfFpGroup(pi)[k],",\n");
   od;
   AppendTo("FundamentalGroup.out","#                           ",RelatorsOfFpGroup(pi)[Length(RelatorsOfFpGroup(pi))]," ]\n#\n");
fi;
AppendTo("FundamentalGroup.out","#####\n\n");

AppendTo("FundamentalGroup.out","pi_free:=FreeGroup(");

if Length(GeneratorsOfGroup(pi)) > 0 then
   for k in [1..(Length(GeneratorsOfGroup(pi))-1)] do
       AppendTo("FundamentalGroup.out","\"pi.",k,"\",");
   od;
   AppendTo("FundamentalGroup.out","\"pi.",Length(GeneratorsOfGroup(pi)),"\");\n");
else
   AppendTo("FundamentalGroup.out","0);\n");
fi;

AppendTo("FundamentalGroup.out","pi:=pi_free/[");
for j in [1..Length(RelatorsOfFpGroup(pi))] do
    string:=String(RelatorsOfFpGroup(pi)[j]);
    split_string:=ShallowCopy(SplitString(string,"*"));
    for k in [1..Length(split_string)] do
        split_element:=ShallowCopy(SplitString(split_string[k],"^"));
        for l in [1..Length(GeneratorsOfGroup(pi))] do
            if split_element[1] = String(GeneratorsOfGroup(pi)[l]) then
               AppendTo("FundamentalGroup.out","pi_free.",l);
               if Length(split_element) = 2 then
                  AppendTo("FundamentalGroup.out","^",split_element[2]);
               fi;
               if k < Length(split_string) then
                  AppendTo("FundamentalGroup.out","*");
               fi;
            fi;
        od;
    od;
    if j < Length(RelatorsOfFpGroup(pi)) then
       AppendTo("FundamentalGroup.out",",\n");
    fi;
od;
AppendTo("FundamentalGroup.out","];");


if verbose = 1 then
   Print("Generators (abelian) = ",GeneratorsOfGroup(H_1),"\n");
   Print("Relators (abelian) = ",RelatorsOfFpGroup(H_1),"\n\n");
fi;

AppendTo("FundamentalGroup.out","\n\n\n\n#####\n#\n",
                               "# Generators (abelian) = ",GeneratorsOfGroup(H_1),"\n");
if Length(RelatorsOfFpGroup(H_1)) <= 1 then
   AppendTo("FundamentalGroup.out","# Relators (abelian) = ",RelatorsOfFpGroup(H_1),"\n#\n");
else
   AppendTo("FundamentalGroup.out","# Relators (abelian) = [ ",RelatorsOfFpGroup(H_1)[1],",\n");
   for k in [2..(Length(RelatorsOfFpGroup(H_1))-1)] do
       AppendTo("FundamentalGroup.out","#                        ",RelatorsOfFpGroup(H_1)[k],",\n");
   od;
   AppendTo("FundamentalGroup.out","#                        ",RelatorsOfFpGroup(H_1)[Length(RelatorsOfFpGroup(H_1))]," ]\n#\n");
fi;
AppendTo("FundamentalGroup.out","#####\n\n");

AppendTo("FundamentalGroup.out","pi_ab:=pi_free/[");
for k in [1..(Length(GeneratorsOfGroup(pi))-1)] do
    for l in [(k+1)..(Length(GeneratorsOfGroup(pi)))] do
        AppendTo("FundamentalGroup.out","pi_free.",k,"*pi_free.",l,"*pi_free.",k,"^-1*pi_free.",l,"^-1,\n");
    od;
od;
for j in [1..Length(RelatorsOfFpGroup(pi))] do
    string:=String(RelatorsOfFpGroup(pi)[j]);
    split_string:=ShallowCopy(SplitString(string,"*"));
    for k in [1..Length(split_string)] do
        split_element:=ShallowCopy(SplitString(split_string[k],"^"));
        for l in [1..Length(GeneratorsOfGroup(pi))] do
            if split_element[1] = String(GeneratorsOfGroup(pi)[l]) then
               AppendTo("FundamentalGroup.out","pi_free.",l);
               if Length(split_element) = 2 then
                  AppendTo("FundamentalGroup.out","^",split_element[2]);
               fi;
               if k < Length(split_string) then
                  AppendTo("FundamentalGroup.out","*");
               fi;
            fi;
        od;
    od;
    if j < Length(RelatorsOfFpGroup(pi)) then
       AppendTo("FundamentalGroup.out",",\n");
    fi;
od;
AppendTo("FundamentalGroup.out","];");