# 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-2020)
# 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
        i:=TransitiveIdentification(Source(epi));
        id:=[n,i];
        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
            i:=TransitiveIdentification(Source(epi));
            id:=[n,i];
            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 ValueOption("info")=true or ValueOption("Info")=true then
        Print(n,"\n");
    fi;
    if n=1 then
        id:=[1,1];
    elif n<=30 then
        i:=TransitiveIdentification(HGgr);
        id:=[n,i];
    else
        id:=fail;
    fi;
    if ValueOption("info")=true or ValueOption("Info")=true then
        Print(id,"\n");
    fi;
    return id;
end;

AlwaysHNPholds:= function(Tid)
    local n,i,tbl,tbl4,tbl6,tbl8,tbl9,tbl10,tbl12,tbl14,tbl15,tbl16;
    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];
    tbl12:=[2,3,4,7,9,10,16,18,20,31,32,33,34,37,40,43,47,52,54,55,56,57,59,
      61,64,65,66,70,71,74,75,77,88,92,93,96,97,100,102,117,122,130,132,133,
      144,168,171,172,174,176,179,188,194,210,214,230,232,234,242,246,255,
      261,271,280];
    tbl14:=[30];
    tbl15:=[9,14];
    tbl16:=[2,3,4,5,7,8,9,10,11,13,15,16,17,18,19,20,21,23,25,27,29,30,32,
      33,34,35,36,37,38,39,40,42,43,44,45,46,47,48,50,51,52,53,54,57,58,61,
      62,63,64,67,68,69,70,71,72,73,74,75,76,78,80,81,82,83,86,87,88,89,90,
      91,92,93,94,96,98,99,100,101,102,105,106,107,109,110,111,112,113,114,
      115,116,117,118,119,120,122,123,127,128,129,130,131,132,133,134,137,
      138,140,141,142,143,148,149,150,151,152,153,154,156,159,160,162,164,
      165,166,167,169,170,171,172,173,174,175,176,177,178,179,180,182,183,
      184,185,192,193,194,197,198,199,200,201,202,203,204,205,206,207,208,
      210,212,214,215,217,218,221,223,225,229,230,231,232,237,239,241,244,
      245,246,248,249,250,251,254,255,256,259,260,261,262,264,265,266,267,
      268,269,270,271,272,274,275,276,277,278,279,280,281,282,283,284,287,
      288,291,292,293,294,295,296,297,298,301,303,304,305,307,308,309,310,
      311,312,313,315,316,317,318,319,320,321,323,324,326,328,329,331,334,
      335,336,340,341,342,343,344,345,347,350,351,352,353,357,358,359,360,
      361,362,363,364,365,368,371,372,373,375,376,377,379,380,381,383,384,
      385,386,387,389,393,394,396,397,398,400,401,402,404,406,407,411,412,
      414,416,417,419,420,423,424,426,427,429,430,433,435,437,440,441,442,
      445,446,448,449,450,452,453,454,455,456,457,459,460,462,463,464,465,
      466,469,470,473,475,477,482,487,489,494,496,497,498,499,501,504,506,
      507,508,509,510,512,513,514,523,528,531,532,535,540,543,544,545,546,
      548,549,555,556,557,558,559,562,566,567,570,571,574,575,576,579,582,
      583,585,586,587,588,589,592,594,595,596,597,598,601,603,604,605,606,
      607,608,609,613,616,618,619,620,622,624,627,628,632,633,635,638,639,
      640,641,642,643,644,645,647,649,651,652,653,654,655,656,657,658,659,
      660,661,665,666,667,668,669,676,677,682,683,684,686,689,690,693,695,
      700,703,706,708,709,713,714,716,717,718,720,721,722,723,724,731,733,
      734,737,738,750,751,752,755,758,766,767,769,771,772,776,779,782,783,
      784,787,793,794,796,797,799,801,802,803,804,805,806,807,809,810,812,
      813,815,820,821,822,823,824,825,829,831,833,834,835,842,845,846,848,
      850,853,854,855,856,857,860,862,869,872,874,875,878,881,884,886,887,
      893,894,895,896,903,905,913,915,916,918,920,921,922,926,927,928,930,
      931,935,948,953,963,966,967,970,976,977,978,981,982,984,985,988,991,
      996,1000,1002,1005,1008,1009,1010,1013,1014,1015,1020,1021,1023,1024,
      1028,1029,1032,1037,1039,1042,1043,1050,1058,1059,1063,1068,1069,1070,
      1071,1072,1080,1082,1083,1084,1088,1089,1092,1095,1097,1098,1099,1101,
      1102,1103,1104,1110,1111,1114,1115,1117,1118,1121,1122,1123,1125,1130,
      1133,1136,1139,1141,1142,1143,1144,1145,1146,1148,1151,1152,1154,1157,
      1159,1163,1167,1168,1171,1172,1175,1182,1183,1187,1188,1190,1196,1201,
      1202,1203,1206,1207,1208,1209,1211,1212,1215,1219,1220,1221,1228,1229,
      1230,1231,1232,1234,1238,1241,1242,1245,1246,1255,1268,1273,1279,1282,
      1283,1284,1285,1288,1290,1295,1301,1302,1303,1305,1306,1308,1309,1327,
      1329,1330,1331,1333,1334,1335,1336,1339,1341,1349,1352,1353,1355,1356,
      1358,1360,1362,1367,1369,1370,1372,1375,1377,1378,1379,1381,1383,1384,
      1386,1389,1393,1398,1403,1404,1406,1407,1412,1414,1415,1418,1419,1424,
      1431,1435,1436,1437,1438,1439,1443,1449,1453,1457,1461,1463,1468,1470,
      1472,1479,1481,1485,1492,1494,1506,1509,1510,1511,1512,1513,1515,1516,
      1517,1518,1522,1525,1526,1543,1544,1545,1546,1547,1548,1549,1550,1551,
      1555,1557,1560,1562,1564,1566,1568,1577,1580,1586,1587,1588,1590,1598,
      1599,1600,1619,1620,1625,1626,1639,1645,1650,1654,1656,1657,1658,1668,
      1669,1672,1676,1695,1696,1697,1703,1705,1707,1716,1717,1747,1749,1752,
      1753,1767,1769,1785,1786,1790,1792,1794,1803,1809,1810,1811,1819,1822,
      1831,1844,1845,1854,1859,1874,1876,1877,1891,1892,1895,1896,1898,1904,
      1908,1911,1917,1919,1921,1927,1936,1937,1941];
    tbl:=[[],[],[],tbl4,[],tbl6,[],tbl8,tbl9,tbl10,[],tbl12,[],tbl14,tbl15,
      tbl16];
    if Tid=fail then
        return fail;
    fi;
    n:=Tid[1];
    i:=Tid[2];
    if IsPrime(n) or n=1 then
        return true;
    elif n>16 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;

ChooseGiIterator:= function(bG,bH)
    local bGs,Gicandidates,FindGi;
    bGs:=ConjugacyClassesSubgroups(bG);
    Gicandidates:=Filtered(bGs,x->not IsMetacyclic(Representative(x)));
    FindGi:= function(cGi,bH)
        local bG,Gi,His,Hi,flag;
        bG:=ActingDomain(cGi);
        for Gi in Elements(cGi) do
            His:=Reversed(List(ConjugacyClassesSubgroups(Intersection(Gi,bH)),
             Representative));
            for Hi in His do
                if AlwaysHNPholds(CosetRepresentationTid(Gi,Hi))=true then
                    return Gi;
                fi;
            od;
        od;
        return fail;
    end;
    return IteratorByFunctions( rec(
        IsDoneIterator:= function(iter)
            local i;
            for i in [iter!.pos..Length(iter!.Gicandidates)] do
                iter!.found:=iter!.FindGi(iter!.Gicandidates[i],iter!.bH);
                if IsGroup(iter!.found) then
                    iter!.pos:=i;
                    return false;
                fi;
            od;
            return true;
        end,
        NextIterator:= function(iter)
            local i,Gi;
            if IsGroup(iter!.found) then
                Gi:=iter!.found;
                iter!.pos:=iter!.pos+1;
                iter!.found:=false;
                return Gi;
            else
                for i in [iter!.pos..Length(iter!.Gicandidates)] do
                    iter!.found:=iter!.FindGi(iter!.Gicandidates[i],iter!.bH);
                    if IsGroup(iter!.found) then
                        Gi:=iter!.found;
                        iter!.pos:=i+1;
                        iter!.found:=false;
                        return Gi;
                    fi;
                od;
                return fail;
            fi;
        end,
        ShallowCopy:= function(iter)
            return rec(
                bG:=iter!.bG,
                bH:=iter!.bH,
                bGs:=iter!.bGs,
                Gicandidates:=iter!.Gicandidates,
                FindGi:=iter!.FindGi,
                pos:=iter!.pos,
                found:=iter!.found);
        end,
        bG:=bG,
        bH:=bH,
        bGs:=bGs,
        Gicandidates:=Gicandidates,
        FindGi:=FindGi,
        pos:=1,
        found:=false));
end;

MaximalNormalSeries:= function(G)
    local M,ans,Ns,N;
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    M:=G;
    ans:=[G];
    Ns:=NormalSubgroups(G);
    while Order(M)>1 do
        Ns:=Filtered(Ns,x->x<>M and IsSubgroup(M,x));
        M:=Ns[1];
        for N in Ns do
            if Order(N)>Order(M) then
                M:=N;
            fi;
        od;
        Add(ans,M);
    od;
    return ans;
end;

KerResH3Z:= function(G,H)
    local MNS,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
        MNS:=MaximalNormalSeries(G);
        if Length(MNS)>2 then
            RG:=ResolutionNormalSeries(MNS,4);
        else
            RG:=ResolutionFiniteGroup(G,4);
        fi;
    fi;
    CRG3:=CR_CocyclesAndCoboundaries(RG,3,true);
    H3Z:=CRG3.torsionCoefficients;
    if H3Z=[] then
        return [[],[[],[]]];
    fi;
    torbase:=DiagonalMat(H3Z);
    kerbase:=IdentityMat(Length(H3Z));
    if ValueOption("iterator")=true or ValueOption("Iterator")=true then
        Gis:=ChooseGiIterator(G,H);
    else
        Gis:=ChooseGi(G,H);
    fi;
    for Gi in Gis do
        if IsNilpotent(Gi) then
            RGi:=ResolutionNormalSeries(LowerCentralSeries(Gi),4);
        elif IsSolvable(Gi) then
            RGi:=ResolutionNormalSeries(DerivedSeries(Gi),4);
        else
            MNS:=MaximalNormalSeries(Gi);
            if Length(MNS)>2 then
                RGi:=ResolutionNormalSeries(MNS,4);
            else
                RGi:=ResolutionFiniteGroup(G,4);
            fi;
        fi;
        ker:=ResHnZ(RG,RGi,3,CRG3).Ker;
        kerbase:=LatticeIntersection(kerbase,Union(ker[2][2],torbase));
        kerbase:=LatticeBasis(kerbase);
        if kerbase=torbase then
            break;
        fi;
    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;

ConjugacyClassSubgroupsNGHOrbitRep:= function(arg)
    local G,H,Gv,NGH,NGGv,NGHNGGv,reps;
    G:=arg[1];
    Gv:=arg[2];
    if Length(arg)=2 then
        H:=Stabilizer(G,1);
    else
        H:=arg[3];
    fi;
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    NGH:=Normalizer(G,H);
    NGGv:=Normalizer(G,Gv);
    reps:=List(DoubleCosets(G,NGGv,NGH),Representative);
    return List(reps,x->Gv^x);
end;

ConjugacyClassesSubgroupsNGHOrbitRep:= function(arg)
    local cs,G,H;
    cs:=arg[1];
    G:=ActingDomain(cs[1]);
    if Length(arg)=1 then
        H:=Stabilizer(G,1);
    else
        H:=arg[2];
    fi;
    return List(List(cs,Representative),x->ConjugacyClassSubgroupsNGHOrbitRep(G,x,H));
end;

MinConjugacyClassesSubgroups:= function(cslist)
    local G,cs,K,MK,M,flag,ans;
    if cslist=[] then
        return [];
    fi;
    G:=ActingDomain(cslist[1]);
    ans:=[];
    for cs in cslist do
        K:=Representative(cs);
        MK:=MaximalSubgroups(K);
        flag:=true;
        for M in MK do
            if ConjugacyClassSubgroups(G,M) in cslist then
                flag:=false;
                break;
            fi;
        od;
        if flag=true then
            Add(ans,cs);
        fi;
    od;
    return ans;
end;

IsInvariantUnderAutG:= function(cslist)
    local G,AutG,InAutG,OutGRep,todolist,todolistold,cs,Imcs;
    if cslist=[] then
        return true;
    fi;
    G:=ActingDomain(cslist[1]);
    AutG:=AutomorphismGroup(G);
    InAutG:=InnerAutomorphismsAutomorphismGroup(AutG);
    OutGRep:=List(RightCosets(AutG,InAutG),Representative);
    todolist:=cslist;
    while todolist<>[] do
        todolistold:=todolist;
        for cs in todolistold do
            if not cs in todolist then
                continue;
            fi;
            Imcs:=Set(OutGRep,x->ConjugacyClassSubgroups(G,Image(x,Representative(cs))));
            if not IsSubset(Set(cslist),Imcs) then
                return false;
            else
                todolist:=Difference(todolist,Imcs);
            fi;
        od;
    od;
    return true;
end;

AbelianInvariantsGoverH := function(lg,lh)
    local A,Ag,e,Gg,G,Hg,H;
    if(lg[1]<>lh[1]) then
        return false;
    fi;
    A:=AbelianGroup(lg[1]);
    Ag:=GeneratorsOfGroup(A);
    e:=Identity(A);
    Gg:=List(lg[2],x->Product([1..Length(Ag)],y->Ag[y]^x[y]));
    G:=Group(Gg,e);
    Hg:=List(lh[2],x->Product([1..Length(Ag)],y->Ag[y]^x[y]));
    H:=Group(Hg,e);
    if IsSubgroup(G,H) then
        return AbelianInvariantsSNF(G/H);
    else
        return fail;
    fi;
end;

AllSubgroups2 := function(G)
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    return AllSubgroups(G);
end;

StemExtensions:= function(G)
    local ScG,ScGg,K,Ks,ans,m,pi,cG,cGg,iso,GG,GGg,Gg,epi,n,i,id;
    ScG:=SchurCoverG(G);
    ScGg:=GeneratorsOfGroup(ScG.SchurCover);
    K:=Kernel(ScG.epi);
    Ks:=AllSubgroups2(K);
    ans:=[];
    for m in Ks 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
            i:=TransitiveIdentification(Source(epi));
            id:=[n,i];
            Add(ans,rec(StemExtension:=Source(epi), epi:=epi, Tid:=id));
        else
            Add(ans,rec(StemExtension:=Source(epi), epi:=epi));
        fi;
    od;
    return ans;
end;