# The following algorithms need
# GAP version>=4.8.7 and GAP package HAP version>=1.11.15.

LoadPackage("HAP");

ConjugacyClassesSubgroups2:= function(G)
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    return ConjugacyClassesSubgroups(G);
end;

BlockSum:= function(l,n)
  local k,ans,i;
  k:=Length(l)/n;
  ans:=[];
  for i in [1..k] do
    Add(ans,Sum([(i-1)*n+1..i*n],x->l[x]));
  od;
  return ans;
end;

Cores := function(arg)
  local v,CRH,CRG,RH,RG,n,RGH,map,G,H,u,k;
  v:=arg[1];
  if Length(arg)=3 then
    CRH:=arg[2];
    CRG:=arg[3];
    k:=Length(CRH.cocyclesBasis[1])/Length(CRG.cocyclesBasis[1]);
    return CRG.cocycleToClass(BlockSum(CRH.classToCocycle(v),k));
  else
    RH:=arg[2];
    RG:=arg[3];
    n:=arg[4];
    G:=RG!.group;
    H:=RH!.group;
    if Length(arg)>5 then
      RGH:=arg[5];
      map:=arg[6];
    else
      if Length(arg)=5 then
        RGH:=arg[5];
      else
        RGH:=ResolutionFiniteSubgroup(RG,H);
      fi;
      map:=HomToIntegers(EquivariantChainMap(RGH,RH,IdentityMapping(H)));
    fi;
  fi;
  u:=Map(map)(v,n);
  k:=Order(G)/Order(H);
  return BlockSum(u,k);
end;

AbelianInvariantsSNF := function(G)
  local n,m,s,l;
  n:=AbelianInvariants(G);
  m:=DiagonalMat(n);
  s:=SmithNormalFormIntegerMat(m);
  return Filtered(DiagonalOfMat(s),x -> x>1);
end;

chooseH := function(G)
  local hh,hhc,hhd,hhdtf,n,ic,i,jc,j;
  hh:=[];
  hhc:=List(ConjugacyClassesSubgroups2(G),Elements);
  hhd:=List(hhc,x->List(x,y->[y,DerivedSubgroup(y)]));
  hhdtf:=List(hhd,x->x[1][1]<>x[1][2]);
  n:=Length(hhd);
  for ic in [n,n-1..2] do
    if hhdtf[ic]=true then
      i:=hhd[ic][1];
      for jc in [1..ic-1] do
        for j in hhd[jc] do
          if i[2]=j[2] and IsSubgroup(i[1],j[1]) then
            hhdtf[jc]:=false;
            break;
          fi;
        od;
      od;
    fi;
  od;
  hh:=List(Filtered([1..n],x->hhdtf[x]),y->hhd[y][1][1]);
  return hh;
end;

chooseHH1trivial := function(G)
  local hh,hhc,kd,hhctf,n,i,j,H,K,HGK,gk,orb,ker,o,im,k,hom;
  hh:=[];
  hhc:=List(ConjugacyClassesSubgroups2(G),Representative);
  kd:=List(hhc,x->[x,DerivedSubgroup(x)]);
  hhctf:=List(kd,x -> x[1]<>x[2]);
  n:=Length(hhc);
  for i in [n,n-1..2] do
    if hhctf[i]=true then
      H:=hhc[i];
      for j in [2..i-1] do
        if hhctf[j]=true then
          K:=hhc[j];
          HGK:=DoubleCosetsNC(G,H,K);
          ker:=kd[j][1];
          for orb in HGK do
            o:=RepresentativesContainedRightCosets(orb);
            o:=List(o,x->RightCoset(H,x));
            im:=[];
            gk:=GeneratorsOfGroup(ker);
            for k in gk do
              Add(im,Product(OrbitsDomain(Group(k),o,OnRight),
                x->(k^Length(x))^(Representative(x[1])^-1)));
            od;
            Apply(im,x->
              Image(NaturalHomomorphismByNormalSubgroup(H,kd[i][2]),x));
            hom:=GroupHomomorphismByImages(ker,H/kd[i][2],gk,im);
            ker:=Kernel(hom);
          od;
          kd[j][1]:=ker;
          if kd[j][1]=kd[j][2] then
            hhctf[j]:=false;
          fi;
        fi;
      od;
    fi;
  od;
  hh:=List(Filtered([1..n],x->hhctf[x]),y->hhc[y]);
  return hh;
end;

HnSubgroupBase := function(tor,im)
  local torbase,im1,v,base1,base2,ans;
  torbase:=DiagonalMat(tor);
  im1:=LatticeBasis(Concatenation(torbase,im));
  im1:=Filtered(im1,x->not x in torbase);
  ans:=[];
  base1:=torbase;
  for v in im1 do
    base2:=LatticeBasis(Concatenation(base1,[v]));
    if base2<>base1 then
      Add(ans,v);
      base1:=base2;
    fi;
  od;
  return ans;
end;

H4pFromResolution := function(arg)
  local G,RG,CRG4,hh,H,RH,CRH2,CRH4,n,RGH,map,idH,idlist,Reslist,
    im1,im2,im,I,i,j,tor,torbase,H4,H4gen;
  RG:=arg[1];
  G:=RG!.group;
  if Length(arg)>1 then
    CRG4:=arg[2];
  else
    CRG4:=CR_CocyclesAndCoboundaries(RG,4,true);
  fi;
  if ValueOption("H1trivial")=true or ValueOption("H1Trivial")=true then
    hh:=chooseHH1trivial(G);
  else
    hh:=chooseH(G);
  fi;
  im1:=[];
  Print(Length(hh));Print(List(hh,
    function(H)
      if ID_AVAILABLE(Order(H))<>fail then
        return IdSmallGroup(H);
      else
        return [Order(H),"?"];
      fi;
    end
  ));
  Print("\n");
  idlist:=[];
  Reslist:=[];
  for H in hh do
Print(Position(hh,H));Print("/");Print(Length(hh));Print("\n");
    if H=G then
      CRH2:=CR_CocyclesAndCoboundaries(RG,2,true);
      CRH4:=CRG4;
      n:=Length(CRH2.torsionCoefficients);
      if n>0 then
        I:=IdentityMat(n);
        im2:=List(UnorderedTuples(I,2),
          v->IntegralCupProduct(RG,v[1],v[2],2,2,CRH2,CRH2,CRH4));
      else
        im2:=[];
      fi;
    else
      if ID_AVAILABLE(Order(H))<>fail then
        idH:=IdSmallGroup(H);
        i:=Position(idlist,idH);
        if i<>fail then
          RH:=Reslist[i][1];
          im2:=Reslist[i][2];
        else
          if IsNilpotentGroup(H) then
            RH:=ResolutionNormalSeries(LowerCentralSeries(SmallGroup(idH[1],idH[2])),5);
          else
            RH:=ResolutionFiniteGroup(SmallGroup(idH[1],idH[2]),5);
          fi;
          CRH2:=CR_CocyclesAndCoboundaries(RH,2,true);
          CRH4:=CR_CocyclesAndCoboundaries(RH,4,true);
          n:=Length(CRH2.torsionCoefficients);
          if n>0 then
            I:=IdentityMat(n);
            im2:=List(UnorderedTuples(I,2),
              v->IntegralCupProduct(RH,v[1],v[2],2,2,CRH2,CRH2,CRH4));
            im2:=HnSubgroupBase(CRH4.torsionCoefficients,im2);
            Apply(im2,CRH4.classToCocycle);
          else
            im2:=[];
          fi;
          Add(idlist,idH);
          Add(Reslist,[RH,im2]);
        fi;
      else
        if IsNilpotentGroup(H) then
          RH:=ResolutionNormalSeries(LowerCentralSeries(H),5);
        else
          RH:=ResolutionFiniteGroup(H,5);
        fi;
        CRH2:=CR_CocyclesAndCoboundaries(RH,2,true);
        CRH4:=CR_CocyclesAndCoboundaries(RH,4,true);
        n:=Length(CRH2.torsionCoefficients);
        if n>0 then
          I:=IdentityMat(n);
          im2:=List(UnorderedTuples(I,2),
            v->IntegralCupProduct(RH,v[1],v[2],2,2,CRH2,CRH2,CRH4));
          im2:=HnSubgroupBase(CRH4.torsionCoefficients,im2);
          Apply(im2,CRH4.classToCocycle);
        else
          im2:=[];
        fi;
      fi;
      RGH:=ResolutionFiniteSubgroup(RG,H);
      map:=HomToIntegers(EquivariantChainMap(RGH,RH,IsomorphismGroups(H,RH!.group)));
      im2:=List(im2,i->CRG4.cocycleToClass(Cores(i,RH,RG,4,RGH,map)));
    fi;
    im1:=Concatenation(im1,im2);
  od;
  tor:=CRG4.torsionCoefficients;
  torbase:=DiagonalMat(tor);
  im1:=LatticeBasis(Concatenation(torbase,im1));
  im1:=LatticeBasis(Difference(im1,torbase));
  H4:=AbelianGroup(CRG4.torsionCoefficients);
  H4gen:=GeneratorsOfGroup(H4);
  im:=Group(List(im1,x->Product([1..Length(H4gen)],y->H4gen[y]^x[y])),Identity(H4));
  return [AbelianInvariantsSNF(im),[tor,im1]];
end;

mast := function(RG,RH,RI,RJ)
  local I,J,Ig,Jg,iso,HI,RHI,m,eqchmap;
  I:=RI!.group;
  J:=RJ!.group;
  Ig:=First(I,x->Order(x)=Order(I));
  Jg:=First(J,x->Order(x)=Order(J));
  iso:=GroupHomomorphismByImages(I,J,[Ig],[Jg]);
  RHI:=ResolutionFiniteDirectProduct(RH,RI);
  HI:=RHI!.group;
  m:=GroupHomomorphismByFunction(HI,RG!.group,
    x->Image(Projection(HI,1),x)*Image(iso,Image(Projection(HI,2),x)));
  eqchmap:=EquivariantChainMap(RHI,RG,m);
  return eqchmap;
end;

i2ast := function(RHI,RI)
  local i2,eqchmap;
  i2:=Embedding(RHI!.group,2);
  eqchmap:=EquivariantChainMap(RI,RHI,i2);
  return eqchmap;
end;

pr2ast := function(RI,RHI)
  local pr2,eqchmap;
  pr2:=Projection(RHI!.group,2);
  eqchmap:=EquivariantChainMap(RHI,RI,pr2);
  return eqchmap;
end;

sHI := function(xi,RHI,RI,n)
  local i2star,pr2star;
  i2star:=HomToIntegers(i2ast(RHI,RI));
  pr2star:=HomToIntegers(pr2ast(RI,RHI));
  return xi-Map(pr2star)(Map(i2star)(xi,n),n);
end;

HnZtoHnminus1QoverZfromResolution := function(arg)
  local v,R,n,p,M,Mp,l,i,j,w,basis,cols,u,c,ansp,ans;
  v:=arg[1];
  R:=arg[2];
  n:=arg[3];
  M:=[];
  for i in [1..R!.dimension(n)] do
    l:=List([1..R!.dimension(n-1)],x->0);
    w:=R!.boundary(n,i);
    for j in w do
      l[AbsInt(j[1])]:=l[AbsInt(j[1])]+SignInt(j[1]);
    od;
    Add(M,l);
  od;
  M:=TransposedMat(M);
  if Length(arg)>3 then
    p:=arg[4];
    Mp:=List(p,x->M[x]);
  else
    p:=[1..R!.dimension(n-1)];
    Mp:=M;
  fi;
  basis:=NormalFormIntMat(Mp,6);
  cols:=List([1..basis.rank],i->First([1..R!.dimension(n)],
    j->basis.normal[i][j]<>0));
  ansp:=List([1..Length(p)],x->0);
  u:=v;
  for i in [1..basis.rank] do
    c:=u[cols[i]]/basis.normal[i][cols[i]];
    ansp:=ansp+c*basis.rowtrans[i];
    u:=u-c*basis.normal[i];
  od;
  if Set(u)<>[0] then
    return fail;
  fi;
  ans:=List([1..R!.dimension(n-1)],x->0);
  for i in [1..Length(p)] do
    ans[p[i]]:=ansp[i];
  od;
  return ans;
end;

Hnminus1QoverZtoHnZfromResolution := function(v,R,n)
  local M,l,i,j,w;
  M:=[];
  for i in [1..R!.dimension(n)] do
    l:=List([1..R!.dimension(n-1)],x->0);
    w:=R!.boundary(n,i);
    for j in w do
      l[AbsInt(j[1])]:=l[AbsInt(j[1])]+SignInt(j[1]);
    od;
    Add(M,l);
  od;
  M:=TransposedMat(M);
  return v*M;
end;

HS := function(f,RHI,n)
  local fQoverZ,d1,d2;
  if n < 2 then
    Print(
     "ERROR: n must be at least 2. \n" );
    return fail;
  fi;
  d1:=RHI!.dimension(n-1)-RHI!.dimension(n-2);
  d2:=RHI!.dimension(n-1)-RHI!.dimension(n-3);
  fQoverZ:=HnZtoHnminus1QoverZfromResolution(f,RHI,n,[1..d2]);
  return fQoverZ{[d1+1..d2]};
end;

GeneratorOfCyclicGroup := function(C)
  return First(C,x->Order(x)=Order(C));
end;

chooseHI := function(G)
  local I,HI,HItf,n,i,j,ic,jc;
  I:=List(Filtered(ConjugacyClassesSubgroups2(G),
    x->IsCyclic(Representative(x))),Elements);
  HI:=List(I,x->List(x,y->[Centralizer(G,y),y]));
  HItf:=List(HI,x->true);
  n:=Length(HI);
  for ic in [n,n-1..2] do
    if HItf[ic]=true then
      i:=HI[ic][1];
      if GroupCohomology(i[1],3)=[] then
        HItf[ic]:=false;
      fi;
      for jc in [1..ic-1] do
        for j in HI[jc] do
          if i[1]=j[1] and IsSubgroup(i[2],j[2]) then
            HItf[jc]:=false;
            break;
          fi;
        od;
      od;
    fi;
  od;
  HI:=List(Filtered([1..n],x->HItf[x]),y->HI[y][1]);
  Apply(HI,x->[x[1],Group(GeneratorOfCyclicGroup(x[2]))]);
  HI:=Filtered(HI,x->Order(x[2])>1);
  return HI;
end;

chooseHISubgroup := function(G)
  local I,HI,HItf,n,i,j,ic,jc,idH,Hs3trivial27and125and343,Hs3trivial81,
    Hs3trivial625and2401,Hs3trivial243,Hs3trivial3125;
  Hs3trivial27and125and343:=[1..4];
  Hs3trivial81:=Concatenation([[1..6],[8..10],[13,14]]);
  Hs3trivial625and2401:=Concatenation([[1..6],[9,10],[13,14]]);
  Hs3trivial243:=Concatenation([[1],[3..12],[19..27],[33],[35],[43..47],[49,50],[66]]);
  Hs3trivial3125:=Concatenation([[1],[15..17],[24..29],[42],[44],[51..57],[59,60],[76]]);
  I:=List(Filtered(ConjugacyClassesSubgroups2(G),
    x->IsCyclic(Representative(x))),Elements);
  HI:=List(I,x->List(x,y->[Centralizer(G,y),y]));
  HItf:=List(HI,x->true);
  n:=Length(HI);
  for ic in [n,n-1..2] do
    if HItf[ic]=true then
      i:=HI[ic][1];
      if GroupCohomology(i[1],3)=[] then
        HItf[ic]:=false;
      fi;
      if ID_AVAILABLE(Order(i[1]))<>fail then
        idH:=IdSmallGroup(i[1]);
        if idH[1] in [9,25,49] then
          HItf[ic]:=false;
        elif idH[1] in [27,125,343] and idH[2] in Hs3trivial27and125and343 then
          HItf[ic]:=false;
        elif idH[1]=81 and idH[2] in Hs3trivial81 then
          HItf[ic]:=false;
        elif idH[1] in [625,2401] and idH[2] in Hs3trivial625and2401 then
          HItf[ic]:=false;
        elif idH[1]=243 and idH[2] in Hs3trivial243 then
          HItf[ic]:=false;
        elif idH[1]=3125 and idH[2] in Hs3trivial3125 then
          HItf[ic]:=false;
        fi;
      fi;
      for jc in [1..ic-1] do
        for j in HI[jc] do
          if i[1]=j[1] and IsSubgroup(i[2],j[2]) then
            HItf[jc]:=false;
            break;
          fi;
        od;
      od;
    fi;
  od;
  HI:=List(Filtered([1..n],x->HItf[x]),y->HI[y][1]);
  Apply(HI,x->[x[1],Group(GeneratorOfCyclicGroup(x[2]))]);
  HI:=Filtered(HI,x->Order(x[2])>1);
  return HI;
end;

IsUnramifiedH3 := function(RG,v)
  local G,CRG4,u,hi,k,H,I,RH,RI,RJ,RHI,m,mstar,mu,smu,HSmu,ZHSmu,CRH3,c;
  G:=RG!.group;
  if Length(v)=RG!.dimension(4) then
    u:=v;
  else
    CRG4:=CR_CocyclesAndCoboundaries(RG,4,true);
    u:=CRG4!.classToCocycle(v);
  fi;
  if ValueOption("subgroup")=true or ValueOption("Subgroup")=true then
    hi:=chooseHISubgroup(G);
  else
    hi:=chooseHI(G);
  fi;
  for k in hi do
    H:=k[1];
    I:=k[2];
    if H=G then
      RH:=RG;
    elif IsNilpotentGroup(H) then
      RH:=ResolutionNormalSeries(LowerCentralSeries(H),5);
    else
      RH:=ResolutionFiniteGroup(H,5);
    fi;
    RI:=ResolutionAbelianGroup([Order(I)],5);
    RJ:=ResolutionAbelianGroup(I,5);
    m:=mast(RG,RH,RI,RJ);
    RHI:=m!.source;
    mstar:=HomToIntegers(m);
    mu:=Map(mstar)(u,4);
    smu:=sHI(mu,RHI,RI,4);
    Print(Position(hi,k));Print("/");Print(Length(hi));Print("\n");
    HSmu:=HS(smu,RHI,4);
    ZHSmu:=Hnminus1QoverZtoHnZfromResolution(HSmu,RH,3);
    CRH3:=CR_CocyclesAndCoboundaries(RH,3,true);
    c:=[CRH3!.torsionCoefficients,CRH3!.cocycleToClass(ZHSmu)];
Print(c);Print("\n");
    if c[2]=fail then
      return fail;
    elif LatticeBasis(c)<>[c[1]] then
      return false;
    fi;
  od;
  return true;
end;