##############################################
#To run these programs in GAP save this file as "*.g" and have GAP read the file in the usual way.
#It is necessary to load the KBMAG package. Information about KBMAG can be found at http://www-gap.mcs.st-and.ac.uk/Packages/kbmag.html.
##############################################
#Input: Free group, relators, prime p
#Output: Dimension of Tor(H_1(G),F_p)

Tor:=function(Freegroup,Relators,Prime)
local AbelInv, list1, list2, x;
AbelInv:=AbelianInvariants(Freegroup/Relators);
list1:=[];
list2:=[];
for x in AbelInv do
     if x<>0 then Add(list1,x);
     fi;
     od;
for x in list1 do
     if x mod Prime = 0 then Add(list2,1);
     fi;
     od;
return Sum(list2);
end;;

###############################################
#P-Primary Rank
#Output: P-Primary Rank of Fp Group

PrimePrimaryRank:=function(Freegroup,Relators,Prime)
local AbelInv, list1, list2, x;
AbelInv:=AbelianInvariants(Freegroup/Relators);
list1:=[];
list2:=[];
for x in AbelInv do
     if x <> 0 then Add(list1,x);
     fi;
     od;
for x in list1 do
     if x mod Prime = 0 then Add(list2,PadicValuation(x,Prime));
     fi;
	     od;
return Sum(list2);	
end;;

###############################################
#The (special) Word Problem
#Output: Reduced word

Reduce_Word:=function(Freegroup,Relators,TestWord,Sublist,Prime)
local Rel_P, GroupGen, comm, G, RG, OR;

Rel_P:=List(Relators,x->x^Prime);
GroupGen:=GeneratorsOfGroup(Freegroup);
comm:=ListX(GroupGen,Relators,Comm);

G:=Freegroup/Concatenation(comm,Rel_P,Sublist);
RG:=KBMAGRewritingSystem(G);
OR:=OptionsRecordOfKBMAGRewritingSystem(RG);
OR.maxeqns:=500000;
OR.tindyint:=100;
MakeConfluent(RG);

return ReducedWord(RG,TestWord);
end;;

#################################################
#Attempts to reduce a generating set
#Output: List of generators

FindBasis:=function(Freegroup,Relators,Prime,Sublist)
local Gen,TestWord,x;
Gen:=Sublist;
for x in Sublist do
     TestWord:=Reduce_Word(Freegroup,Relators,x,Difference(Gen,[x]),Prime);
     if IsOne(TestWord)=true then Gen:=Difference(Gen,[x]);
     fi;
     od;
return Size(Gen);
end;;

#######################################################
#Gives the estimate for H_2
#Output: Upper bound on dimension of H_2
SecondHomologyCoefficients:=function(Freegroup, Relators, Prime, Sublist)

local a,b,c,d,e,f,ff,RPrime;

f:=GeneratorsOfGroup(Freegroup);
ff:=ListX(f,f,Comm);
RPrime:=List(Relators,x->x^Prime);

a:=Tor(Freegroup,Relators,Prime);
b:=PrimePrimaryRank(Freegroup,Concatenation(Relators,ff),Prime);
c:=PrimePrimaryRank(Freegroup,Concatenation(RPrime,ff),Prime);
e:=FindBasis(Freegroup,Relators,Prime,Sublist);
d:=a+b-c+e;
return d;
end;;