/* ********************************************************************* */
EnumeratedGroup := function(G, p);

    if IsPrime(p) then
        eG, GhomeG := pQuotient(G, p, 63);
        eGhomG := hom< eG -> G | [ (eG.i)@@GhomeG :
                                   i in [1..NumberOfPCGenerators(eG)] ] >;
    else
        CSpc := CosetSpace(G, sub<G|>: CosetLimit := 2000000);
        GhomeG, eG := CosetAction(CSpc);
        eGhomG := hom< eG -> G | [ (eG.i)@@GhomeG :
                                   i in [1..#Generators(eG)] ] >;
    end if;
    return eG, GhomeG, eGhomG;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
RelatorMatrix := function(F,rels);

    Z:=Integers(); M:=[];
    VW:=Hom(RSpace(Z,#rels),RSpace(Z,#Generators(F)));
    for r in rels do     for x in Generators(F) do
                         Append(~M,ExponentSum(r,x));
    end for;             end for;
    M:= VW!M; B,T:=LLL(M);
    return B,T;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
MapPreimage := function(eAmapeB, b);

    eA:=Domain(eAmapeB);
    for a in eA do
    if eAmapeB(a) eq b then return a; break; end if;
    end for;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
NQuotient := function(eA,n);

   eQ, eAhomeQ:=quo<eA| [x^n : x in Generators(eA)]>;
   NQ:=[ MapPreimage(eAhomeQ,a) : a in eQ];
   NQmapeQ:=map<NQ->eQ|x:->eAhomeQ>;
   return NQmapeQ;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
NCoQuotient := function(eA,n);

   NQ :=  [x : x in eA | Order(x) eq GCD(Order(x),n)];
   return NQ;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
AbelianGroupGenerators := function(eA);

    ord:=1;
    for x in Generators(eA) do
        ord:=ord*Order(x);
    end for;
    if ord eq Order(eA) then
        return [x : x in Generators(eA) | x ne Id(eA)];
    else
        Gens:={};
        for i in [1..#Generators(eA)] do
        for x in Transversal(eA,Gens) do
              if Order(x) eq Exponent(quo<eA|Gens>) then
                      Gens:=Gens join {x}; break;
              end if;
        end for;
        end for;
    return SetToSequence(Gens diff {Id(eA)});
    end if;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
Extension := function(v,Cohomology,eAhomA,GhomF,FhomG,RGens,SGens,MGens,p);

    eA:=Domain(eAhomA); A:=Codomain(eAhomA);
    G:=Domain(GhomF); F:=Codomain(GhomF);
    n:=#Generators(F); m:=#Generators(A);
    FA:=FreeGroup(n+m);
    FhomFA:=hom<F->FA | [FA.i : i in [1..n]] >;
    AhomFA:=hom<A->FA | [FA.i : i in [n+1..n+m]] >;
    Erels:= {FhomFA(x) : x in SGens} join
          { AhomFA(LHS(x)*RHS(x)^-1) : x in Relations(A) } join
          { (FA.i,FA.j) : i in [1..n] , j in [n+1..n+m] };
    for j in [1..#RGens] do
        w:=Cohomology[j][v[j]];
        Erels:=Erels join {FhomFA(RGens[j])*AhomFA(eAhomA(w))^-1};
    end for;
    for j in [1+#RGens..#Cohomology] do
        w:=Cohomology[j][v[j]];
        Erels:=Erels join {FhomFA(MGens[j-#RGens])*AhomFA(eAhomA(w))^-1};
    end for;
    E, FAhomE:=quo<FA|Erels>;
    temp:=[F.i : i in [1..n]];
    for i in [1..m] do
        Append(~temp,Id(F));
    end for;
    EhomF:=hom<E->F | temp>;
    EhomG:=hom<E->G | x:-> FhomG(EhomF(x)) >;
    return EhomG, E,  G, eA, FhomFA, AhomFA, FAhomE;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
Cocycle := function(v,Cohomology,eAhomA,GhomF,FhomG,RGens,SGens,MGens,p);

    EhomG, E, G, eA, FhomFA, AhomFA, FAhomE :=
                                            Extension(v,Cohomology,
                                            eAhomA,GhomF,FhomG,RGens,
                                            SGens,MGens,p);
    eE, EhomeE, eEhomE:= EnumeratedGroup(E,p);
    eG, GhomeG, eGhomG:= EnumeratedGroup(G,p);
    eAhomeE:=hom<eA->eE | x:->EhomeE(FAhomE(AhomFA(eAhomA(x)))) >;
    GmapeE:=map<G->eE | x:->EhomeE(FAhomE(FhomFA(GhomF(eGhomG(GhomeG(x)))))) >;
    C:=car<G,G>;
    CmapeA:=map<C->eA | x:->  MapPreimage(eAhomeE,
                 GmapeE(x[1]*x[2])^-1*GmapeE(x[1])*GmapeE(x[2])) >;
    CmapA:=map<C->A | x:->eAhomA(CmapeA(x))>;
    return  CmapA;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
ConstructV:=function(Cohomology);

    Max:=Maximum([#Cohomology[i] : i in [1..#Cohomology]]);
    Z:=quo<Integers()|Max>;
    V:=RSpace(Z,#Cohomology);
    IV:=RSpace(Integers(),#Cohomology);
    V:=[IV!v : v in V];
    for j in [1..#Cohomology] do
        temp:=[v : v in V | v[j] ge #Cohomology[j]];
        for v in temp do
              Exclude(~V,v);
        end for;
    end for;
    V:=[v+IV![1 : i in [1..#Cohomology]] : v in V];
    return V;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
SequenceElementToInteger:=function(S,x);

    for i in [1..#S] do
        if S[i] eq x then return i; break; end if;
    end for;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
InducedPhiOnA:=function(Phi,i,V,Cohomology,AhomeA,eAhomA,UsefulOrders);

    A:=Domain(Phi); eA:=Domain(eAhomA);
    eAPhieA:=hom<eA->eA | x:->AhomeA(Phi(eAhomA(x)))>;
    eAphieA:=[];
    uo:=#UsefulOrders;
    for j in [1..uo] do
         NQ,NQuotient:=quo<eA|[x^UsefulOrders[j] : x in Generators(eA)]>;
         compo:=map<Cohomology[j]->NQ|x:->NQuotient(x)>;
         eAphieA[j]:=map<Cohomology[j]->Cohomology[j] | x:->
              MapPreimage(compo,NQuotient(eAPhieA(x)))>;
    end for;
    for j in [1+uo..#Cohomology] do
         eAphieA[j]:=eAPhieA;
    end for;
    elt1:=[eAphieA[j](Cohomology[j][V[i][j]]) : j in [1..#Cohomology]];
    elt2:=V![SequenceElementToInteger(Cohomology[j],elt1[j]) :
                                                  j in [1..#Cohomology]];
     return SequenceElementToInteger(V,elt2);

end function;
/* ******************************************************************** */



/* ******************************************************************** */
IthComponentOfX:=function(i,x,SchurMultiplier,Mgens);

    S:=[Mgens[j] : j in [1..#Mgens] | j ne i];
    eQ,SchurMultiplierhomeQ:=quo<SchurMultiplier|S>;
    for e in [0..Exponent(SchurMultiplier)-1] do
        if SchurMultiplierhomeQ(Mgens[i]^e*x^-1) eq Id(eQ) then return e;
        break; end if;
    end for;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
SolveYeqXB:=function(Y,B);

    for j in [1..#Y] do
    Append(~B,Y[j]);
    end for;
    B:= Hom(RSpace(Integers(),#Y+1),RSpace(Integers(),#Y))!B;
    B1,T:=LLL(B);
    return RSpace(Integers(),#Y)![-T[#Y+1][j] : j in [1..#Y]];

end function;
/* ******************************************************************** */



/* ******************************************************************** */
InducedPhiOnG:=function(Phi,i,V,Cohomology,UsefulOrders,FhomG,GhomF,FhomDG,
            DGhomF, DGhomeDG,eDGhomDG, SchurMultiplier,RGens,Mgens,B,eA);

    G:=Domain(Phi); F:=Domain(FhomDG); eDG:=Domain(eDGhomDG);
    GhomeDG:=hom<G->eDG|x:-> DGhomeDG(FhomDG(GhomF(x)))>;
    eDGhomG:=hom<eDG->G|x:-> FhomG(DGhomF(eDGhomDG(x)))>;
    FPhiF:=hom<F->F|x:-> GhomF(Phi(FhomG(x)))>;
    eDGPhieDG:=hom<eDG->eDG|x:-> GhomeDG(Phi(eDGhomG(x)))>;
    B:=[B[i][j] : i, j in [1..#RGens]];

    M:=[]; N:=[];
    for k in [1..#RGens] do
    Y:=[];
         for x in Generators(F) do
         Append(~Y,ExponentSum(FPhiF(RGens[k]),x));
         end for;
         X:=SolveYeqXB(Y,B);
         Append(~M,X);
    end for;
    M:=Hom(RSpace(Integers(),#RGens),RSpace(Integers(),#RGens))!M;
    for k in [1..#Mgens] do
         for j in [1..#Mgens] do
         Append(~N,IthComponentOfX(j,eDGPhieDG(Mgens[k]),SchurMultiplier,
                                                  Mgens));
         end for;
    end for;
    N:=Hom(RSpace(Integers(),#Mgens),RSpace(Integers(),#Mgens))!N;

    elt:=[Cohomology[j][V[i][j]] : j in [1..#Cohomology]];
    elt1:=[];
    for k in [1..#RGens] do
         w:=Id(eA);
         for j in [1..#RGens] do
         w:=w*elt[j]^M[k][j];
         end for;
         NQ,NQuotient:=quo<eA|[x^UsefulOrders[k] : x in Generators(eA)]>;
         compo:=map<Cohomology[k]->NQ | x:->NQuotient(x)>;
         w:=MapPreimage(compo,NQuotient(w));
    Append(~elt1,w);
    end for;
    for k in [1..#Mgens] do
          w:=Id(eA);
           for j in [1..#Mgens] do
            w:=w*elt[j+#RGens]^N[k][j];
           end for;
    Append(~elt1,w);
    end for;

    elt2:=V![SequenceElementToInteger(Cohomology[j],elt1[j]) :
                                    j in [1..#Cohomology]];
    return SequenceElementToInteger(V,elt2);

end function;
/* ******************************************************************** */



/* ******************************************************************** */
InducedPhi:=function(Phi,V,Cohomology,AhomeA,eAhomA,UsefulOrders,
                     FhomG,GhomF,FhomDG, DGhomF, DGhomeDG, eDGhomDG,
                     SchurMultiplier,RGens,Mgens,B,eA);

if Domain(Phi).1 in Domain(AhomeA) then
   return map<[1..#V]->[1..#V] | x:->
           InducedPhiOnA(Phi,x,V,Cohomology,AhomeA,eAhomA,UsefulOrders)>;
else
   return map<[1..#V]->[1..#V] | x:->
           InducedPhiOnG(Phi,x,V,Cohomology,UsefulOrders,FhomG,GhomF,FhomDG,
                       DGhomF, DGhomeDG, eDGhomDG, SchurMultiplier,RGens,
                                                               Mgens,B,eA)>;
end if;

end function;
/* ******************************************************************** */




/* ******************************************************************** */
/*                  THE MAIN FUNCTION                                   */
/* ******************************************************************** */
SecondCohomology := function(G,A,Phi,p,t);

    n:=#Generators(G);
    F:=FreeGroup(n);
    GhomF:=hom<G->F | [F.i : i in [1..n]] >;
    FhomG:=hom<F->G | [G.i : i in [1..n]] >;
    eA, AhomeA, eAhomA := EnumeratedGroup(A,p);
    expA:=Exponent(eA);

    Grels:={ GhomF(LHS(x)*RHS(x)^-1) : x in Relations(G) };
    Grels:=SetToSequence(Grels);
    B,T:=RelatorMatrix(F,Grels);
    RGens:=[]; SGens:=[]; DGrels:={x^expA : x in Grels};
    for i in [1..n] do
        w:=Identity(F);
        for j in [1..#Grels] do
            w:=w*Grels[j]^T[i][j];
        end for;
        Append(~RGens,w);
        DGrels:=DGrels join {w};
    end for;
    for x in Generators(F) do          for r in Grels do
            Append(~SGens,(r,x));
            DGrels:=DGrels join {(r,x)};
    end for;                            end for;

    DGrels:=DGrels diff {Id(F)};
    DG, FhomDG := quo<F | DGrels>; delete DGrels;
    DGhomF:=hom<DG->F | [F.i : i in [1..n]]>;
    eDG, DGhomeDG, eDGhomDG := EnumeratedGroup(DG,p);
    QuotientedSchurMultiplier:=sub<eDG|DGhomeDG(FhomDG(Grels))>;
    QuotientedSchurMultiplier;
    GAbelianRels:= {(F.i,F.j) : i,j in [1..n] | i lt j };
    for i in [1 .. n] do
        w:=Identity(F);
        for j in [1..n] do
             w:=w*F.j^B[j][i];
        end for;
       GAbelianRels:= GAbelianRels join {w};
    end for;
    GAbelian:= quo<F|GAbelianRels>;
    eGAbelian, GAbelianhomeGAbelian, eGAbelianhomGAbelian :=
            EnumeratedGroup(GAbelian,p);
    FhomGAbelian:=hom<F->GAbelian | [GAbelian.i : i in [1..n]] >;
    FhomeGAbelian:=hom<F->eGAbelian | x:->
                                     GAbelianhomeGAbelian(FhomGAbelian(x)) >;
    eA, AhomeA, eAhomA := EnumeratedGroup(A,p);
    Cohomology:=[];
    COHOMOLOGY:=[];
    UsefulOrders:=[];
    for i in [1..n] do
        temphom:=NQuotient(eA,Order(FhomeGAbelian(F.i)));
        Append(~Cohomology, Domain(temphom));
        Append(~COHOMOLOGY, Codomain(temphom));
        Append(~UsefulOrders,Order(FhomeGAbelian(F.i)));
    end for;

    Mgens:=[x : x in AbelianGroupGenerators(QuotientedSchurMultiplier)];
    MGens:=[DGhomF(eDGhomDG(x)) : x in
                      AbelianGroupGenerators(QuotientedSchurMultiplier)];
    for x in Mgens do
        Append(~Cohomology, NCoQuotient(eA,Order(x)));
        Append(~COHOMOLOGY, sub<eA|NCoQuotient(eA,Order(x))> );

    end for;
    EXTENSIONS:=[]; COCYCLES:=[];
    V:=ConstructV(Cohomology);
    for v in V do
       EhomG := Extension(v,Cohomology,eAhomA,GhomF,
                    FhomG,RGens,SGens,MGens,p);
       Append(~EXTENSIONS,EhomG);
       if t eq 1 then
           CmapA := Cocycle(v,Cohomology,eAhomA,GhomF,
                    FhomG,RGens,SGens,MGens,p);
           Append(~COCYCLES,CmapA);
       end if;
    end for;
    INVARIANTS:=AbelianInvariants(DirectProduct(COHOMOLOGY));

    INDUCEDPHI:=InducedPhi(Phi,V,Cohomology,AhomeA,eAhomA,UsefulOrders,
                           FhomG,GhomF,FhomDG,DGhomF,DGhomeDG,eDGhomDG,
                           QuotientedSchurMultiplier, RGens,Mgens,B,eA);

    return COHOMOLOGY, COCYCLES, EXTENSIONS, INDUCEDPHI, INVARIANTS, V;

end function;
/* ******************************************************************** */