# The following code is used in
# "Norm one tori and Hasse norm principle"
# by Akinari Hoshi, Kazuki Kanai and Aiichi Yamasaki
# Written by Aiichi Yamasaki (2019)
# GAP package HAP version>=1.11.15 is needed

LoadPackage("HAP");

Norm1TorusJ :=function(d,n)
    local I,M1,M2,M,f,Sn,T;
    I:=IdentityMat(d-1);
    Sn:=SymmetricGroup(d);
    T:=TransitiveGroup(d,n);
    M1:=Concatenation(List([2..d-1],x->I[x]),[-List([1..d-1],One)]);
    if d=2 then
        M:=[M1];
    else
        M2:=Concatenation([I[2],I[1]],List([3..d-1],x->I[x]));
        M:=[M1,M2];
    fi;
    f:=GroupHomomorphismByImages(Sn,Group(M),GeneratorsOfGroup(Sn),M);
    return Image(f,T);
end;

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

AbelianizationGen:= function(G)
    local Gab,pi,inv,A,iso,gen,genrep;
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    pi:=NaturalHomomorphismByNormalSubgroup(G,DerivedSubgroup(G));
    Gab:=Image(pi);
    inv:=AbelianInvariantsSNF(Gab);
    A:=AbelianGroup(inv);
    iso:=IsomorphismGroups(A,Gab);
    gen:=List(GeneratorsOfGroup(A),x->Image(iso,x));
    genrep:=List(gen,x->PreImagesRepresentative(pi,x));
    return rec(Gab:=Gab, gen:=gen, genrep:=genrep, inv:=inv, pi:=pi);
end;

FindGenFiniteAbelian:= function(g)
    local e,a,ga,iso;
    e:=AbelianInvariants(g);
    if Length(e)>1 then
        e:=SmithNormalFormIntegerMat(DiagonalMat(e));
        e:=List([1..Length(e)],x->e[x][x]);
        e:=Filtered(e,x->x>1);
    fi;
    a:=AbelianGroup(e);
    ga:=GeneratorsOfGroup(a);
    iso:=IsomorphismGroups(a,g);
    return List(ga,x->Image(iso,x));
end;

EltFiniteAbelian:= function(arg)
    local g,c,gg,F,gF,hom,cF,e;
    g:=arg[1];
    c:=arg[2];
    if Length(arg)=3 then
        gg:=arg[3];
    else
        gg:=GeneratorsOfGroup(g);
    fi;
    F:=FreeGroup(Length(gg));
    gF:=GeneratorsOfGroup(F);
    hom:=GroupHomomorphismByImages(F,g,gF,gg);
    cF:=PreImagesRepresentative(hom,c);
    e:=List(gF,x->ExponentSumWord(cF,x));
    return e;
end;

FirstObstructionN:= function(arg)
    local G,H,Gab,Hab,K,Kinv,mat,v,Habbase,ker1;
    G:=arg[1];
    if Length(arg)=1 then
        H:=Stabilizer(G,1);
    else
        H:=arg[2];
    fi;
    Gab:=AbelianizationGen(G);
    Hab:=AbelianizationGen(H);
    Hab.Hab:=Hab.Gab;
    Unbind(Hab.Gab);
    if DerivedSubgroup(H)=H then
        return rec(ker:=[[],[[],[]]], Hab:=Hab, Gab:=Gab, psi:=[]);
    fi;
    if DerivedSubgroup(G)=G then
        return rec(ker:=[Hab.inv,[Hab.inv,IdentityMat(Length(Hab.inv))]],
            Hab:=Hab, Gab:=Gab, psi:=List(Hab.inv,x->[]));
    fi;
    K:=Image(Hab.pi,Intersection(H,DerivedSubgroup(G)));
    Kinv:=AbelianInvariantsSNF(K);
    mat:=[];
    for v in Hab.genrep do
        Add(mat,EltFiniteAbelian(Gab.Gab,Image(Gab.pi,v),Gab.gen));
    od;
    Habbase:=DiagonalMat(Hab.inv);
    ker1:=List(GeneratorsOfGroup(K),x->EltFiniteAbelian(Hab.Hab,x,Hab.gen));
    ker1:=LatticeBasis(Concatenation(Habbase,ker1));
    ker1:=LatticeBasis(Difference(ker1,Habbase));
    return rec(ker:=[Kinv,[Hab.inv,ker1]],
        Hab:=Hab, Gab:=Gab, psi:=mat);;
end;

FirstObstructionDnr:= function(arg)
    local G,H,Gab,Hab,HG,HGrep,Dnrgen,h,x,Dnr,Dnrinv,Habbase,Dnrmat;
    G:=arg[1];
    if Length(arg)=1 then
        H:=Stabilizer(G,1);
    else
        H:=arg[2];
    fi;
    Gab:=AbelianizationGen(G);
    Hab:=AbelianizationGen(H);
    Hab.Hab:=Hab.Gab;
    Unbind(Hab.Gab);
    if DerivedSubgroup(H)=H then
        return rec(Dnr:=[[],[[],[]]], Hab:=Hab, Gab:=Gab);
    fi;
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    HG:=RightCosets(G,H);
    HGrep:=List(HG,Representative);
    Dnrgen:=[];
    for x in HGrep do
        for h in GeneratorsOfGroup(Intersection(H,H^x)) do
            Add(Dnrgen,Image(Hab.pi,Comm(h,x^-1)));
        od;
    od;
    Dnr:=Group(Dnrgen,Identity(Hab.Hab));
    Dnrinv:=AbelianInvariantsSNF(Dnr);
    Habbase:=DiagonalMat(Hab.inv);
    Dnrmat:=List(Dnrgen,x->EltFiniteAbelian(Hab.Hab,x,Hab.gen));
    Dnrmat:=LatticeBasis(Concatenation(Habbase,Dnrmat));
    Dnrmat:=LatticeBasis(Difference(Dnrmat,Habbase));
    return rec(Dnr:=[Dnrinv,[Hab.inv,Dnrmat]],
        Hab:=Hab, Gab:=Gab);
end;

FirstObstructionDr:= function(arg)
    local G,Gv,H,Gab,Hab,HGGv,HGGvrep,Hwi,Hwiab,Gvab,psi2i,i,psi2iimage,Hw,
          psi2,ker,phi1i,phi1iimage,phi1,Dr,Drinv,Habbase,Drmat;
    G:=arg[1];
    Gv:=arg[2];
    if Length(arg)=2 then
        H:=Stabilizer(G,1);
    else
        H:=arg[3];
    fi;
    Gab:=AbelianizationGen(G);
    Hab:=AbelianizationGen(H);
    Hab.Hab:=Hab.Gab;
    Unbind(Hab.Gab);
    if DerivedSubgroup(H)=H then
        return rec(Dr:=[[],[[],[]]], Hab:=Hab, Gab:=Gab);
    fi;
    HGGv:=DoubleCosets(G,H,Gv);
    HGGvrep:=List(HGGv,Representative);
    Hwi:=List(HGGvrep,x->Intersection(Gv^(x^(-1)),H));
    Hwiab:=List(Hwi,AbelianizationGen);
    Gvab:=AbelianizationGen(Gv);
    psi2i:=[];
    for i in [1..Length(HGGv)] do
        psi2iimage:=List(Hwiab[i].genrep,x->x^HGGvrep[i]);
        psi2iimage:=List(psi2iimage,x->Image(Gvab.pi,x));
        Add(psi2i,GroupHomomorphismByImages(Hwiab[i].Gab,Gvab.Gab,Hwiab[i].gen,
            psi2iimage));
    od;
    Hw:=DirectProduct(List(Hwiab,x->x.Gab));
    psi2:=GroupHomomorphismByFunction(Hw,Gvab.Gab,x->
     Product([1..Length(HGGv)],i->Image(psi2i[i],Image(Projection(Hw,i),x))));
    ker:=Kernel(psi2);
    phi1i:=[];
    for i in [1..Length(HGGv)] do
        phi1iimage:=List(Hwiab[i].genrep,x->Image(Hab.pi,x));
        Add(phi1i,GroupHomomorphismByImages(Hwiab[i].Gab,Hab.Hab,Hwiab[i].gen,
            phi1iimage));
    od;
    phi1:=GroupHomomorphismByFunction(Hw,Hab.Hab,x->
     Product([1..Length(HGGv)],i->Image(phi1i[i],Image(Projection(Hw,i),x))));
    Dr:=Image(phi1,ker);
    Drinv:=AbelianInvariantsSNF(Dr);
    Habbase:=DiagonalMat(Hab.inv);
    Drmat:=List(GeneratorsOfGroup(Dr),x->EltFiniteAbelian(Hab.Hab,x,Hab.gen));
    Drmat:=LatticeBasis(Concatenation(Habbase,Drmat));
    Drmat:=LatticeBasis(Difference(Drmat,Habbase));
    return rec(Dr:=[Drinv,[Hab.inv,Drmat]],
        Hab:=Hab, Gab:=Gab);
end;

MaximalSubgroups2:= function(G)
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    return SortedList(MaximalSubgroups(G));
end;

SchurCoverG:= function(G)
    local epi,iso,ScG,ScGg,GG,GGg,Gg,n,i,id;
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    epi:=EpimorphismSchurCover(G);
    iso:=IsomorphismPermGroup(Source(epi));
    ScG:=Source(epi);
    ScGg:=GeneratorsOfGroup(ScG);
    GG:=Range(iso);
    GGg:=List(ScGg,x->Image(iso,x));
    Gg:=List(ScGg,x->Image(epi,x));
    epi:=GroupHomomorphismByImages(GG,G,GGg,Gg);
    n:=NrMovedPoints(Source(epi));
    if n>=2 and n<=30 and IsTransitive(Source(epi),[1..n]) then
        for i in [1..NrTransitiveGroups(n)] do
            if Order(TransitiveGroup(n,i))=Order(Source(epi)) and
             IsConjugate(SymmetricGroup(n),
              TransitiveGroup(n,i),Source(epi)) then
                id:=[n,i];
                break;
            fi;
        od;
        return rec(SchurCover:=Source(epi), epi:=epi, Tid:=id);
    else
        return rec(SchurCover:=Source(epi), epi:=epi);
    fi;
end;

MinimalStemExtensions:= function(G)
    local ScG,ScGg,K,MK,ans,m,pi,cG,cGg,iso,GG,GGg,Gg,epi,n,i,id;
    ScG:=SchurCoverG(G);
    ScGg:=GeneratorsOfGroup(ScG.SchurCover);
    K:=Kernel(ScG.epi);
    MK:=MaximalSubgroups2(K);
    ans:=[];
    for m in MK do
        pi:=NaturalHomomorphismByNormalSubgroup(ScG.SchurCover,m);
        cG:=Range(pi);
        cGg:=List(ScGg,x->Image(pi,x));
        iso:=IsomorphismPermGroup(Range(pi));
        GG:=Range(iso);
        GGg:=List(cGg,x->Image(iso,x));
        Gg:=List(ScGg,x->Image(ScG.epi,x));
        epi:=GroupHomomorphismByImages(GG,G,GGg,Gg);
        n:=NrMovedPoints(Source(epi));
        if n>=2 and n<=30 and IsTransitive(Source(epi),[1..n]) then
            for i in [1..NrTransitiveGroups(n)] do
                if Order(TransitiveGroup(n,i))=Order(Source(epi)) and
                 IsConjugate(SymmetricGroup(n),
                  TransitiveGroup(n,i),Source(epi)) then
                    id:=[n,i];
                    break;
                fi;
            od;
            Add(ans,rec(MinimalStemExtension:=Source(epi), epi:=epi, Tid:=id));
        else
            Add(ans,rec(MinimalStemExtension:=Source(epi), epi:=epi));
        fi;
    od;
    return ans;
end;

ResHnZ:= function(arg)
    local RG,RH,n,G,H,inj,map,mapZ,CRGn,CRHn,HnG,HnH,m,res,null,ker,Hng,Hnggen,
          Hnh,Hnhgen,resHnggen,torbase,im,coker,hom,cokergen,cokergen1;
    RG:=arg[1];
    RH:=arg[2];
    n:=arg[3];
    G:=RG!.group;
    H:=RH!.group;
    inj:=GroupHomomorphismByFunction(H,G,x->x);
    map:=EquivariantChainMap(RH,RG,inj);
    mapZ:=HomToIntegers(map);
    if Length(arg)>=4 then
        CRGn:=arg[4];
    else
        CRGn:=CR_CocyclesAndCoboundaries(RG,n,true);
    fi;
    if Length(arg)=5 then
        CRHn:=arg[5];
    else
        CRHn:=CR_CocyclesAndCoboundaries(RH,n,true);
    fi;
    HnG:=CRGn.torsionCoefficients;
    HnH:=CRHn.torsionCoefficients;
    if HnG=[] then
        if HnH=[] then
            return rec(HnGZ:=[],HnHZ:=HnH,Res:=[],Ker:=[[],[[],[]]],
                       Coker:=[[],[[],[]]]);
        else
            return rec(HnGZ:=[],HnHZ:=HnH,Res:=[],Ker:=[[],[[],[]]],
                       Coker:=[HnH,[HnH,IdentityMat(Length(HnH))]]);
        fi;
    fi;
    if HnH=[] then
        return rec(HnGZ:=HnG,HnHZ:=[],
            Res:=List(HnG,x->[]),Ker:=[HnG,[HnG,IdentityMat(Length(HnG))]],
                      Coker:=[[],[[],[]]]);
    fi;
    m:=List(IdentityMat(Length(HnG)),x->
            CRHn.cocycleToClass(mapZ!.mapping(CRGn.classToCocycle(x),n)));
    null:=NullspaceIntMat(m);
    Hng:=AbelianGroup(HnG);
    Hnggen:=GeneratorsOfGroup(Hng);
    Hnh:=AbelianGroup(HnH);
    Hnhgen:=GeneratorsOfGroup(Hnh);
    resHnggen:=List(m,x->Product([1..Length(Hnhgen)],y->Hnhgen[y]^x[y]));
    res:=GroupHomomorphismByImages(Hng,Hnh,Hnggen,resHnggen);
    ker:=Kernel(res);
    im:=Image(res);
    null:=List(GeneratorsOfGroup(ker),x->EltFiniteAbelian(Hng,x,Hnggen));
    torbase:=DiagonalMat(HnG);
    null:=LatticeBasis(Concatenation(torbase,null));
    null:=LatticeBasis(Difference(null,torbase));
    hom:=NaturalHomomorphismByNormalSubgroup(Hnh,im);
    coker:=Image(hom);
    if Order(coker)=1 then
        return rec(HnGZ:=HnG,HnHZ:=HnH,Res:=m,
            Ker:=[AbelianInvariantsSNF(ker),[HnG,null]],Coker:=[[],[HnH,[]]]);
    fi;
    cokergen:=FindGenFiniteAbelian(coker);
    cokergen1:=List(cokergen,x->Representative(PreImages(hom,x)));
    cokergen1:=List(cokergen1,x->EltFiniteAbelian(Hnh,x,Hnhgen));
    return rec(HnGZ:=HnG,HnHZ:=HnH,Res:=m,
        Ker:=[AbelianInvariantsSNF(ker),[HnG,null]],
              Coker:=[AbelianInvariants(coker),[HnH,cokergen1]]);
end;

CosetRepresentationTid:= function(G,H)
    local Gg,HG,HGg,HGgr,n,i,id;
    Gg:=GeneratorsOfGroup(G);
    HG:=RightCosets(G,H);
    HGg:=List(Gg,x->Permutation(x,HG,OnRight));
    HGgr:=Group(HGg,());
    n:=Index(G,H);
    if n=1 then
        id:=[1,1];
    elif n<=30 then
        for i in [1..NrTransitiveGroups(n)] do
            if Order(TransitiveGroup(n,i))=Order(HGgr) and
              IsConjugate(SymmetricGroup(n),TransitiveGroup(n,i),HGgr) then
                id:=[n,i];
                break;
            fi;
        od;
    else
        id:=fail;
    fi;
    return id;
end;

AlwaysHNPholds:= function(Tid)
    local n,i,tbl,tbl4,tbl6,tbl8,tbl9,tbl10,tbl14,tbl15;
    tbl4:=[2,4];
    tbl6:=[4,12];
    tbl8:=[2,3,4,9,11,13,14,15,19,21,22,31,32,37,38];
    tbl9:=[2,5,7,9,11,14,23];
    tbl10:=[7,26,32];
    tbl14:=[30];
    tbl15:=[9,14];
    tbl:=[[],[],[],tbl4,[],tbl6,[],tbl8,tbl9,tbl10,[],[],[],tbl14,tbl15];
    if Tid=fail then
        return fail;
    fi;
    n:=Tid[1];
    i:=Tid[2];
    if IsPrime(n) or n=1 then
        return true;
    elif n=12 or n>15 then
        return fail;
    elif i in tbl[n] then
        return false;
    else
        return true;
    fi;
end;

IsMetacyclic:= function(G)
    local p;
    if Order(G)=1 then
        return true;
    fi;
    for p in Set(Factors(Order(G))) do
        if not IsCyclic(SylowSubgroup(G,p)) then
            return false;
        fi;
    od;
    return true;
end;

ChooseGi:= function(bG,bH)
    local bGs,Gicandidates,Gis,cGi,Gi,His,Hi,flag;
    bGs:=ConjugacyClassesSubgroups(bG);
    Gicandidates:=Filtered(bGs,x->not IsMetacyclic(Representative(x)));
    Gis:=[];
    for cGi in Gicandidates do
        for Gi in Elements(cGi) do
            His:=Reversed(List(ConjugacyClassesSubgroups(Intersection(Gi,bH)),
             Representative));
            flag:=false;
            for Hi in His do
                if AlwaysHNPholds(CosetRepresentationTid(Gi,Hi))=true then
                    Add(Gis,Gi);
                    flag:=true;
                    break;
                fi;
            od;
            if flag=true then
                break;
            fi;
        od;
    od;
    return Gis;
end;

KerResH3Z:= function(G,H)
    local RG,CRG3,H3Z,torbase,kerbase,Gis,Gi,RGi,ker,H3,H3g,K;
    if IsNilpotent(G) then
        RG:=ResolutionNormalSeries(LowerCentralSeries(G),4);
    elif IsSolvable(G) then
        RG:=ResolutionNormalSeries(DerivedSeries(G),4);
    else
        RG:=ResolutionFiniteGroup(G,4);
    fi;
    CRG3:=CR_CocyclesAndCoboundaries(RG,3,true);
    H3Z:=CRG3.torsionCoefficients;
    if H3Z=[] then
        return [[],[[],[]]];
    fi;
    torbase:=DiagonalMat(H3Z);
    kerbase:=IdentityMat(Length(H3Z));
    Gis:=ChooseGi(G,H);
    for Gi in Gis do
        if IsNilpotent(Gi) then
            RGi:=ResolutionNormalSeries(LowerCentralSeries(Gi),4);
        elif IsSolvable(Gi) then
            RGi:=ResolutionNormalSeries(DerivedSeries(Gi),4);
        else
            RGi:=ResolutionFiniteGroup(Gi,4);
        fi;
        ker:=ResHnZ(RG,RGi,3,CRG3).Ker;
        kerbase:=LatticeIntersection(kerbase,Union(ker[2][2],torbase));
        kerbase:=LatticeBasis(kerbase);
    od;
    kerbase:=LatticeBasis(Difference(kerbase,torbase));
    H3:=AbelianGroup(H3Z);
    H3g:=GeneratorsOfGroup(H3);
    K:=Group(List(kerbase,x->Product([1..Length(x)],y->H3g[y]^x[y])),Identity(H3));
    return [AbelianInvariantsSNF(K),[H3Z,kerbase]];
end;