# This file contains a system of GAP programs writen by Staszewski and Voelklein
# to compute the cardinality of a Nielsen class using the following structure
# constant formula:
# File rewritten by Adam James January 2011
#
# Let u_1, ..., u_r represent conjugacy classes C_1, ..., C_r of U.  Let
# c_i = |C_i|.
# Then the number of tuples (g_1, ..., g_r) in U with product 1
# such that g_i is U-conjugate to u_i is given by the following formula:
#  c_1 c_2 ... c_r / |U| \sum_{chi \in Irr(U)} chi(C_1) chi (C_2) ... chi(C_r) /
# chi(1)^(r-2).

#
# SubgroupsContainingCC
# returns a list of subgroups all conjugacy classes in tuple.
#
InstallGlobalFunction(SubgroupsContainingCC, function(group, tuple)
  local i,j,cc,reps, ccs, elems, h, subgrps;
  ccs:=List(tuple, i-> Elements(ConjugacyClass(group,i)));
  elems:=Union(ccs);
  h:=Subgroup(group, elems);
  if Size(h) = Size(group) then
    return [group];
  fi;
  subgrps:=IntermediateSubgroups(group,h).subgroups;
  subgrps:=Concatenation([h], subgrps, [group]);
  subgrps:=Sort(subgrps, function(v,w) return Size(v)<Size(w); end);
  return subgrps;
end);


# SubgroupsGenByNielsenTuples(g,ClassTuple)
#
# Returns a list of all subgroups of g generated by a
# tuple in the Nielsen class, taken up to conjugation in g.
#
InstallGlobalFunction(SubgroupsGenByNielsenTuples, function(g,ClassTuple)

# The subroutine OneMoreGenerator is called iteratively, for k=1,2,...
# It successively constructs SubgroupList, which is the
# list that has the following sublist for each conjugacy class
# representative cc[i] of g:  All subgroups of g generated by a
# tuple with product cc[i] from tt[1]^g \times ...\times tt[k]^g,
# taken up to conjugation under the centralizer of cc[i].

# cc is a list of conjugacy class representatives.
  local SL,OneMoreGenerator,numgens,
  ccl,lcc,cc,SubgroupList,k,CC,CCO,t,gg,x,y,h,r,i,j,hh,flag;

  numgens:= Length(ClassTuple);
  ccl:=ConjugacyClasses(g);
  lcc:= Length(ccl);
  cc:=List(ccl,Representative);
  Sort(cc,function(v,w) return Order(v)<Order(w); end); # sort conjugacy classes
  CC:=List(cc,x->Centralizer(g,x));

# The input is a list SubgroupList. Its i-th entry is a list of
# subgroups of g generated by a tuple with product cc[i].
# The output is a NewSubgroupList with the same property.
# In addition, NewSubgroupList[j] contains groups generated
# by a group in SubgroupList[i] plus a conjugate of NewGen.
# The groups in NewSubgroupList[j] are taken up to conjugation
# under the centralizer of cc[j].

# OneMoreGenerator takes as input a list where at index i we have a list of
# subgroups of g generated by a tuple with product cc[i]. The output has the
# same property contains groups generated by groups in our initial list plus
# a conjugate of some new element specified in the function. We then take
# conjugates under the centralizer cc[j]

  OneMoreGenerator :=function(SubgroupList, NewGen)
    local NewSubgroupList,C,t,gg,x,y,h,r,i,j,hh,flag,c;

    NewSubgroupList:=[];
    for j in [1..lcc] do
      NewSubgroupList[j]:=[];
    od;

    for i in [1..lcc] do
      x:= cc[i];
      C:=Centralizer(g,x);
      for gg in SubgroupList[i] do

        for c in DoubleCosetRepsAndSizes(g, Centralizer(g,NewGen),
        Normalizer(C,gg)) do
          t:=NewGen^(c[1]);
          y:=cc[i]*t;

          j:=1;
          while j < lcc+1 do
            r:=RepresentativeAction(g,y,cc[j]);
            if not r=fail then
              break;
            else
              j:=j+1;
            fi;
          od;

          # in some way we are trying to extend the generators for gg
          h:=Group(Concatenation(GeneratorsOfGroup(gg),[t]))^r;
          if not Product(GeneratorsOfGroup(h))=cc[j] then  # why?? Needs to pass
              Print("x");
          fi;

          flag:=0;
          for hh in NewSubgroupList[j] do
            if IsConjugate(CC[j],hh,h) then
              flag:=1;
              break;
            fi;
          od;
          if flag=0 then
            Add(NewSubgroupList[j],h);
          fi;
        od;
      od;
    od;
    return NewSubgroupList;
  end;

# initializing  SubgroupList  ...
  SubgroupList:=[1..lcc];
  for j in [1..lcc] do
    SubgroupList[j]:=[];
  od;

  k:=1;
  while k < lcc+1 do
    if IsConjugate(g,cc[k],ClassTuple[1]) then
      break;
    else
      k:=k+1;
    fi;
  od;

  SubgroupList[k]:=[ Group(cc[k]) ];


# the final iterative application of the subroutine
  for i in [2.. numgens-1] do
    SubgroupList:=OneMoreGenerator(SubgroupList, ClassTuple[i]);
  od;

  k:=1;
  while k < lcc+1 do
    if IsConjugate(g,cc[k]^-1,ClassTuple[numgens]) then
      break;
    else
      k:=k+1;
    fi;
  od;

#  SubgroupList[k] is the desired list of subgroups of g
#  It remains to make it unique up to conjugation

  SL:= [];
  for hh in SubgroupList[k] do
    flag:=0;
    for h in SL  do
      if IsConjugate(g,h,hh) then
        flag:=1;
        break;
      fi;
    od;
    if flag=0 then
      Add(SL,hh);
    fi;
  od;

  Sort(SL,function(v,w) return Size(v)<Size(w); end);
  return SL;
end);

#
# LastOf
#
InstallGlobalFunction(LastOf, function(List)
  return List[Length(List)];
end);

#
# NumberOfNielsenTuples
#
# Computes the number of tuples (u_1, ..., u_r) in U with product 1
# such that u_i is G-conjugate to t_i= ClassTuple[i]. If all classes
# in <ClassTuple> are equal, the function NumberOfNielsenTuplesForSameClasses
# may be faster.
#
InstallGlobalFunction(NumberOfNielsenTuples , function(U, G, ClassTuple)
# We use the distributive law to avoid summing over all combinations of
# conjugacy classes C_1, ..., C_r of U such that the elements of C_i are
# G-conjugate to t_i.
  local Norm, h, ConjClasses,CharTab, TupFusion, N, product,r, i, j, chi, sizes;

  Norm := Normalizer(G, U);
  h := NaturalHomomorphismByNormalSubgroup(Norm, U);
  CharTab:=CharacterTable(Norm);
  ConjClasses := ConjugacyClasses(CharTab);

# We analyse the fusion in Norm for the elements in <ClassTuple>.
# TupFusion[i] is the set of pairs [k,n], where k is the index of some
# conjugacy class of N_G(U)(as found in the character table) also found in U
# which is a subset of ClassTuple[i], and k is the size of this conjugacy class.

  TupFusion := [];
  r:= Length(ClassTuple);
  sizes:=SizesConjugacyClasses(CharTab);
  for i in [ 1 .. r ] do
    TupFusion[i] := [];
    for j in [1.. Length(ConjClasses)] do
      if Order(Image(h, Representative(ConjClasses[j]))) = 1 and
                  IsConjugate(G, Representative(ConjClasses[j]), ClassTuple[i])
      then
        Add(TupFusion[i], [j,sizes[j]]);
      fi;
    od;
    if Length(TupFusion[i]) = 0 then Print("error");
      return 0;
    fi;
  od;

# We use the above structure constant formula.
# We rearrange the sum such that the outer sum runs over the chi's and the
# inner sum over the  possible combinations of r classes of U fusing into
# <ClassTup>. We further use the distributivity law so that the result is
# calculated as a product of sums.
  N:=0;

  for chi in Irr(CharTab) do
    product := 1;
    for i in [ 1 .. r ] do
      product := product * Sum(List(TupFusion[i], x-> chi[x[1]] * x[2]));
    od;

    product := product / chi[1]^(r-2);
    N := N + product;
  od;

  return N / Size(Norm);
end);

#
# TotalNumberOfNielsenTuples
#
# Computes the number of tuples - including commutators
#
InstallGlobalFunction(TotalNumberOfNielsenTuples, function(subgroup, group,
  tuple, genus)
  local i,j,k,s,rep, classreps, comms, debug;

# debugging code
  debug:=false;
  if ValueOption("debug") = true then
    debug:=true;
  fi;

  classreps:=List(ConjugacyClasses(subgroup), Representative);
  s:=0;
# if genus is zero then ignore commutators
  if genus=0 then
      s:=NumberOfNielsenTuples(subgroup, group, tuple);
  else
    # else run over commutators
    comms:= NumberOfCommProducts(group, genus);
    i:=1;
    for rep in classreps do
      s:=s + comms[i]*NumberOfNielsenTuples( subgroup, group,
          Concatenation([rep], tuple));

      # Code for debugging
      if debug then
        Print(" + ", comms[i], " * ");
        Print(NumberOfNielsenTuples( subgroup, group, Concatenation([rep], tuple)));
      fi;

      i:=i+1;
    od;

# Code for debugging
    if debug then
      Print("\n");
    fi;
  fi;

  return s;
end);

# TODO
#
# NumberOfCommProducts
#
# Computes the number of ways elements in a given group can be written as a
# product of g commutators. Outputs a list which at index i is the number of
# ways an element of ConjugacyClasses[i] can be written a product of g comms
#
InstallGlobalFunction(NumberOfCommProducts, function(group, genus)
  local i, j, k, n, tbl, ccs, iden, centsize, irrs, c, tablecomms, comms;

    tbl:=CharacterTable(group);
    ccs:=ConjugacyClasses(group);
    iden := IdentificationOfConjugacyClasses(tbl);
    ccs:=List(ccs, x -> Representative);
    centsize:= SizesCentralizers(tbl)[1]^(2*genus -1);
    irrs:=Irr(tbl);
    tablecomms:=[];
    comms:=[];
  if genus = 0 then
    comms:=List(ccs, x -> 1);
  else
    for c in [1..Length(ccs)] do
      tablecomms[c] := centsize*Sum(List(irrs, x -> x[c]/(x[1]^(2*genus-1))));
    od;
    for c in [1..Length(ccs)] do
      comms[iden[c]]:= tablecomms[c];
    od;
  fi;
  return comms;
end);



#
# NumberOfGeneratingTuples
#
# computes the number of tuples in the Nielsen class given by <ClassTuple>
# which generate the whole group <G>.
#

InstallGlobalFunction(NumberOfGeneratingTuples , function(G, ClassTuple, genus)
  local ContributingSubgrps, numberOfGeneratingTuples,n,h , j, k, U;

# <ContributingSubgrps> is the list of all subgroups of G
# generated by tuples in the Nielsen class.
# These subgroups are listed in ascending order.
# if the genus is positive then we don't bother worrying about generation since
# this is a given (pretty much).
  ContributingSubgrps := SubgroupsGenByNielsenTuples(G, ClassTuple);

# if the whole group cannot be generated, we return 0.
  if  Length(ContributingSubgrps)  = 0 or not Size(LastOf(ContributingSubgrps))
    = Size(G) then
    Print("There are no generating tuples\n");
    return 0;
  fi;

  numberOfGeneratingTuples := [];
  for k in [ 1 ..  Length(ContributingSubgrps)-1 ] do

    U := ContributingSubgrps[k];
    Print("We handle the ", k, "-th contributing subgroup, which has order ",
    Size(U), "\n");

    numberOfGeneratingTuples[k] := TotalNumberOfNielsenTuples(U, G,
    ClassTuple, genus);

#   make sure that the number of generating tuples per subgroup is positive
#   and calculate using inclusion exclusion the contributing factors for each
#   subgroup chain
    if numberOfGeneratingTuples[k]<0 then
      numberOfGeneratingTuples[k]:=-
      numberOfGeneratingTuples[k];
    else
      for j in  [1 ..  k-1]  do
        for h in IsomorphicSubgroups( U, ContributingSubgrps[j] ) do
          if IsConjugate(G, Image(h),  ContributingSubgrps[j])  then
              numberOfGeneratingTuples[k] := numberOfGeneratingTuples[k] -
              numberOfGeneratingTuples[j] * Index(U,Normalizer(U,Image(h)));
          fi;
        od;
      od;
    fi;
  od;

  n:= TotalNumberOfNielsenTuples(G, G, ClassTuple, genus);
  for k in [ 1 ..  Length(ContributingSubgrps)-1] do
    n:= n - numberOfGeneratingTuples[k] *
    Index(G,Normalizer(G,ContributingSubgrps[k]));
##      Print(numberOfGeneratingTuples[k] *
##      Size(Center(ContributingSubgrps[k])) / Size(ContributingSubgrps[k]), "\n");
  od;
  Print(" The number of GeneratingNielsenTuples mod G-conjugation is ",
  n*Size(Center(G))/Size(G), "\n");
# Experiment so that sizes are not mod G-conjugation
  return n*Size(Center(G))/Size(G);
end);
#End