##
## GAP -- Groups, Algorithms and Programming Version 4.4.9 [17]
## Implementation to compute G\otimes G for G a nilpotent group
##
## Input:  G a nilpotent group represented as a PcpGroup
## Output: G\otimes G as a PcpGroup
##
## Dependencies: GAP Packages Polycyclic [11] and nq [23] 
##
TensorSquare := 
    function(G)

    local FP,         # Free product G*G^\varphi 
          lg,rg,      # Generators of the left and right base groups of 
                      #   G*G^\varphi
          pairs,      # Generators of the left and right base group paired up
          lcj, lcomm, # Helper functions for left conjugation and commutation
          nurel,      # Helper function to create relation pairs for nu(G)
          mingp,      # Positions of the minimal generating set of G relative
                      #   to the polycyclic generating set (Igs)
          class,      # Nilpotency class of G
          nu,         # Finite presentation of nu(G)
          LG, RG,     # Images of the generators of the left and right base 
                      #   groups in nu(G) as a finitely presented group
          epi,        # Epimorphism from nu(G)-> polycyclic presentation of
                      #   nu(G)
          LSG, RSG;   # left and right base groups as subgroups of the
                      #   polycyclic presentation of nu(G)
        
    ## Check to see if G is nilpotent and a Pcp group
    ##
    if not (IsPcpGroup(G) and IsNilpotent(G)) then
        Error("Error: G must be a nilpotent pcp group \n");
    fi;

    ## Define functions for left conjugation and commutation
    ##
    lcj   := function (a,b) return a*b*a^-1; end;
    lcomm := function(a,b) return lcj(a,b)*b^-1; end;
 
    ## Define a function to create the relation for \nu(G)        
    ##    {^X}[a,b]=[{^x}a,({^x^\phi}b)^\varphi] 
    ## where X=x (if i=1) or X=x^\phi (if i=2).  
    ##
    nurel := function(x,a,b,i)
        return lcj(x[i],lcomm(a,b)) / lcomm(lcj(x[1],a),lcj(x[2],b)); 
    end;

    ## Create the free product G*G^\varphi and record the class of G
    ##
    FP    := FreeProduct(G,G);
    class := NilpotencyClassOfGroup(G);

    ## Record the positions in the polycyclic generating set (Igs) of 
    ## the minimal generators of the group
    ##
    mingp := List(MinimalGeneratingSet(G), i->Position(Igs(G),i)); 
    if fail in mingp then 
        Error("Error: Minimal generating set must be in the Igs");
    fi;

    ## Image of the Igs (polycyclic generating set) of each base group
    ## in the underlying free group of the free product FP. And pair
    ## them up as needed to create the relations of nu(G).
    ##
    lg    := List(Igs(G), g->UnderlyingElement(Image(Embedding(FP,1),g)));
    rg    := List(Igs(G), g->UnderlyingElement(Image(Embedding(FP,2),g)));
    pairs := List([1..Length(lg)],i->[lg[i],rg[i]]);

    ## Create a finite presentation of nu(G)
    ##
    nu    := FreeGroupOfFpGroup(FP)/Concatenation(RelatorsOfFpGroup(FP),
                 ListX(pairs,lg{mingp},rg{mingp},[1,2], nurel) );

    ## Record the images of the generators of the left and right base 
    ## groups in nu(G) as a finitely presented group.
    ##
    LG    := List(lg, g->ElementOfFpGroup(FamilyObj(nu.1),g));
    RG    := List(rg, g->ElementOfFpGroup(FamilyObj(nu.1),g));

    ## Create an isomorphism from the finite presentation of nu(G) to
    ## the polycyclic presentation for nu(G) using the known 
    ## nilpotency bound on nu(G).
    ##
    epi   := NqEpimorphismNilpotentQuotient(nu, class+1);

    ## Obtain the images of the left and right base groups in the
    ## polycyclic presentation of nu(G)
    ##
    LSG   := Subgroup(Image(epi),List(LG,g->Image(epi,g)));
    RSG   := Subgroup(Image(epi),List(RG,g->Image(epi,g)));

    ## Return [G,G^\phi] which is isomorphic to the tensor square 
    ## 
    return CommutatorSubgroup(LSG,RSG);

    end;