#############################################################################
#
# This GAP code is a modified version of "tgroupisom.gi" by Bettina Eick and
# Ann-Kristin Engel. It has been used by Alexander Cant and Bettina Eick for
# their research article "Torsion-free nilpotent groups of small Hirsch
# length with isomorphic finite quotients".
#
# This file was last modified on May 5, 2023.
#
#############################################################################


# binomials
B2 := function(x) return x*(x-1)/2; end;
B3 := function(x) return x*(x-1)*(x-2)/6; end;


# Hall multiplication
HallMult := function( t, a, b )
    local f;
    f := a+b;
    f[3] := f[3] + t[1][2][3]*a[2]*b[1];
    f[4] := f[4]
       + t[1][2][4]*a[2]*b[1]
       + t[1][3][4]*a[3]*b[1]
       + t[2][3][4]*a[3]*b[2]
       + t[1][2][3]*t[1][3][4]*a[2]*B2(b[1])
       + t[1][2][3]*t[2][3][4]*B2(a[2])*b[1]
       + t[1][2][3]*t[2][3][4]*a[2]*b[1]*b[2];
    f[5] := f[5]
       + t[3][4][5]*a[4]*b[3]
       + t[2][4][5]*a[4]*b[2]
       + t[2][3][5]*a[3]*b[2]
       + t[1][4][5]*a[4]*b[1]
       + t[1][3][5]*a[3]*b[1]
       + t[1][2][5]*a[2]*b[1]
       + t[2][3][4]*t[3][4][5]*B2(a[3])*b[2]
       + t[2][3][4]*t[2][4][5]*a[3]*B2(b[2])
       + t[1][3][4]*t[3][4][5]*B2(a[3])*b[1]
       + t[1][3][4]*t[1][4][5]*a[3]*B2(b[1])
       + t[2][3][4]*t[3][4][5]*a[3]*b[2]*b[3]
       + t[1][3][4]*t[3][4][5]*a[3]*b[1]*b[3]
       + t[1][3][4]*t[2][4][5]*a[3]*b[1]*b[2]
       + t[1][2][4]*t[3][4][5]*a[2]*b[1]*b[3]
       + t[1][2][4]*t[3][4][5]*a[2]*a[3]*b[1]
       + (B2(t[1][2][3])*t[2][3][4]*t[3][4][5]+t[1][2][3]*t[2][3][5]
            +t[1][2][4]*t[2][4][5])*a[2]*b[1]*b[2]
       + (B2(t[1][2][3])*t[2][3][4]*t[3][4][5]-t[1][2][3]*t[2][3][4]*t[3][4][5]
            +t[1][2][3]*t[1][2][4]*t[3][4][5]+t[1][2][3]*t[2][3][5]
            +t[1][2][4]*t[2][4][5])*B2(a[2])*b[1]
       + (B2(t[1][2][3])*t[1][3][4]*t[3][4][5]+t[1][2][3]*t[1][2][4]*t[3][4][5]
            +t[1][2][3]*t[1][3][5]+t[1][2][4]*t[1][4][5])*a[2]*B2(b[1])
       + t[1][2][3]*t[2][3][4]*t[3][4][5]*B2(a[2])*b[1]*b[3]
       + t[1][2][3]*t[2][3][4]*t[3][4][5]*B2(a[2])*a[3]*b[1]
       + t[1][2][3]*t[2][3][4]*t[2][4][5]*a[2]*b[1]*B2(b[2])
       + t[1][2][3]*t[2][3][4]*t[3][4][5]*a[2]*b[1]*b[2]*b[3]
       + t[1][2][3]*t[2][3][4]*t[3][4][5]*a[2]*a[3]*b[1]*b[2]
       + t[1][2][3]*t[1][3][4]*t[3][4][5]*a[2]*B2(b[1])*b[3]
       + t[1][2][3]*t[1][3][4]*t[3][4][5]*a[2]*a[3]*B2(b[1])
       + (2*t[1][2][3]^2*t[1][3][4]*t[3][4][5]
            +3*t[1][2][3]^2*t[2][3][4]*t[3][4][5]
            +2*t[1][2][3]*t[1][2][4]*t[3][4][5]
            +t[1][2][3]*t[1][3][4]*t[2][4][5])*B2(a[2])*B2(b[1])
       + (2*t[1][2][3]^2*t[2][3][4]*t[3][4][5]
            +t[1][2][3]*t[2][3][4]*t[2][4][5])*B3(a[2])*b[1]
       + (t[1][2][3]^2*t[1][3][4]*t[3][4][5]
            +t[1][2][3]*t[1][3][4]*t[1][4][5])*a[2]*B3(b[1])
       + (t[1][2][3]^2*t[2][3][4]*t[3][4][5]
            +t[1][2][3]*t[2][3][4]*t[2][4][5])*B2(a[2])*b[1]*b[2]
       + (t[1][2][3]^2*t[2][3][4]*t[3][4][5]
            +t[1][2][3]*t[1][3][4]*t[2][4][5])*a[2]*B2(b[1])*b[2]
       + 4*t[1][2][3]^2*t[2][3][4]*t[3][4][5]*B3(a[2])*B2(b[1])
       + 3*t[1][2][3]^2*t[1][3][4]*t[3][4][5]*B2(a[2])*B3(b[1])
       + 2*t[1][2][3]^2*t[2][3][4]*t[3][4][5]*B2(a[2])*B2(b[1])*b[2];
    return f;
end;


#############################################################################
#                         Isomorphisms 211
#############################################################################

ZpIsom211 := function( )
    local s, t, h, w, i, j, res, l, r;

    s := List([1..5], x -> List([1..5], y -> List([1..5], z ->
    Indeterminate(Rationals,
        Concatenation("s_{", String(x), String(y), String(z),"}")))));
    s[2][3][4] := 0;
    s[1][2][5] := 0;
    s[1][3][5] := 0;
    s[1][4][5] := 0;
    s[2][3][5] := 0;
    s[2][4][5] := 0;
    s[3][4][5] := 0;

    t := List([1..5], x -> List([1..5], y -> List([1..5], z ->
    Indeterminate(Rationals,
        Concatenation("t_{", String(x), String(y), String(z),"}")))));
    t[2][3][4] := 0;
    t[1][2][5] := 0;
    t[1][3][5] := 0;
    t[1][4][5] := 0;
    t[2][3][5] := 0;
    t[2][4][5] := 0;
    t[3][4][5] := 0;

    # the isomorphism
    h := NullMat(5,5);
    w := [ [2,1], [3,1],[3,2], [4,1],[4,2],[4,3]];
    for i in [1..4] do
        for j in [1..4] do
            if not [i,j] in w then
                h[i][j] := Indeterminate(Rationals,
                   Concatenation("m_{",String(i),String(j),"}"));
            fi;
        od;
    od;

    res := [];

    # evaluate rel: h2 h1 = h1 h2 h3^t123 h4^t124
    l := HallMult(s, h[2], h[1]);
    r := HallMult(s, h[1], h[2]);
    r := HallMult(s, r, t[1][2][3]*h[3]);
    r := HallMult(s, r, t[1][2][4]*h[4]);
    Add( res, l-r );

    # evaluate rel: h3 h1 = h1 h3 h4^t134
    l := HallMult(s, h[3], h[1]);
    r := HallMult(s, h[1], h[3]);
    r := HallMult(s, r, t[1][3][4]*h[4]);
    Add( res, l-r );

    # evaluate rel: h3 h2 = h2 h3
    l := HallMult(s, h[3], h[2]);
    r := HallMult(s, h[2], h[3]);
    Add( res, l-r );

    # summarize
    res := Flat(res);
    res := Filtered(res, x -> x <> 0*x);
    res := Permuted(res, (2,3));

    # adjust
    res[3] := Value( res[3], [t[1][2][3]], [(res[1]+t[1][2][3]*h[3][3])/h[3][3]] );

    return res;
end;


#############################################################################
#                       Isomorphisms 311
#############################################################################

ZpIsom311 := function()
    local s, t, h, w, i, j, res, l, r;

    s := List([1..5], x -> List([1..5], y -> List([1..5], z ->
    Indeterminate(Rationals,
        Concatenation("s_{", String(x), String(y), String(z),"}")))));
    s[1][2][3] := 0;
    s[1][3][4] := 0;
    s[2][3][4] := 0;
    s[2][4][5] := 0;
    s[3][4][5] := 0;

    t := List([1..5], x -> List([1..5], y -> List([1..5], z ->
    Indeterminate(Rationals,
        Concatenation("t_{", String(x), String(y), String(z),"}")))));
    t[1][2][3] := 0;
    t[1][3][4] := 0;
    t[2][3][4] := 0;
    t[2][4][5] := 0;
    t[3][4][5] := 0;

    # the isomorphism
    h := NullMat(5,5);
    w := [ [2,1], [3,1],[3,2], [4,1],[4,2],[4,3], [5,1], [5,2], [5,3], [5,4] ];
    for i in [1..5] do
        for j in [1..5] do
            if not [i,j] in w then
                h[i][j] := Indeterminate(Rationals,
                   Concatenation("m_{",String(i),String(j),"}"));
            fi;
        od;
    od;

    res := [];

    # evaluate rel: h2 h1 = h1 h2 h4^t124 h5^t125
    l := HallMult(s, h[2], h[1]);
    r := HallMult(s, h[1], h[2]);
    r := HallMult(s, r, t[1][2][4]*h[4]);
    r := HallMult(s, r, t[1][2][5]*h[5]);
    Add( res, l-r );

    # evaluate rel: h3 h1 = h1 h3 h4^t134 h5^t135
    l := HallMult(s, h[3], h[1]);
    r := HallMult(s, h[1], h[3]);
    r := HallMult(s, r, t[1][3][4]*h[4]);
    r := HallMult(s, r, t[1][3][5]*h[5]);
    Add( res, l-r );

    # evaluate rel: h4 h1 = h1 h4 h5^t145
    l := HallMult(s, h[4], h[1]);
    r := HallMult(s, h[1], h[4]);
    r := HallMult(s, r, t[1][4][5]*h[5]);
    Add( res, l-r );

    # evaluate rel: h3 h2 = h2 h3 h5^t[2][3][5]
    l := HallMult(s, h[3], h[2]);
    r := HallMult(s, h[2], h[3]);
    r := HallMult(s, r, t[2][3][5]*h[5]);
    Add( res, l-r );

    # evaluate rel: h4 h2 = h2 h4
    l := HallMult(s, h[4], h[2]);
    r := HallMult(s, h[2], h[4]);
    Add( res, l-r );

    # evaluate rel: h3 h4 = h4 h3
    l := HallMult(s, h[4], h[2]);
    r := HallMult(s, h[2], h[4]);
    Add( res, l-r );

    # summarize
    res := Flat( res );
    res := Filtered( res, x -> x <> 0*x );
    res := Permuted( res, (2,5,3,4) );

    # adjust
    res[5] := Value( res[5], [t[1][2][4]], [(res[1]+t[1][2][4]*h[4][4])/h[4][4]] );

   return res;
end;


############################################################################
#                    Isomorphisms 2111
############################################################################

ZpIsom2111 := function()
    local t, s, h, w, i, j, res, l, r;

    s := List([1..5], x -> List([1..5], y -> List([1..5], z ->
    Indeterminate(Rationals,
        Concatenation("s_{", String(x), String(y), String(z),"}")))));
    s[2][3][4] := 0;
    s[2][4][5] := 0;
    s[3][4][5] := 0;

    t := List([1..5], x -> List([1..5], y -> List([1..5], z ->
    Indeterminate(Rationals,
        Concatenation("t_{", String(x), String(y), String(z),"}")))));
    t[2][3][4] := 0;
    t[2][4][5] := 0;
    t[3][4][5] := 0;

    # the isomorphism
    h := NullMat(5,5);
    w := [ [2,1], [3,1],[3,2], [4,1],[4,2],[4,3], [5,1],[5,2],[5,3],[5,4] ];
    for i in [1..5] do
        for j in [1..5] do
            if not [i,j] in w then
                h[i][j] := Indeterminate(Rationals,
                   Concatenation("m_{",String(i),String(j),"}"));
            fi;
        od;
    od;

    res := [];

    # evaluate rel: h2 h1 = h1 h2 h3^t123 h4^t124 h5^t125
    l := HallMult(s, h[2], h[1]);
    r := HallMult(s, h[1], h[2]);
    r := HallMult(s, r, t[1][2][3]*h[3]);
    r := HallMult(s, r, t[1][2][4]*h[4]);
    r := HallMult(s, r, t[1][2][5]*h[5]);
    Add( res, l-r );

    # evaluate rel: h3 h1 = h1 h3 h4^t134 h5^t135
    l := HallMult(s, h[3], h[1]);
    r := HallMult(s, h[1], h[3]);
    r := HallMult(s, r, t[1][3][4]*h[4]);
    r := HallMult(s, r, t[1][3][5]*h[5]);
    Add( res, l-r );

    # evaluate rel: h3 h2 = h2 h3 h5^t235
    l := HallMult(s, h[3], h[2]);
    r := HallMult(s, h[2], h[3]);
    r := HallMult(s, r, t[2][3][5]*h[5]);
    Add( res, l-r );

    # evaluate rel: h4 h1 = h1 h4 h5^t145
    l := HallMult(s, h[4], h[1]);
    r := HallMult(s, h[1], h[4]);
    r := HallMult(s, r, t[1][4][5]*h[5]);
    Add( res, l-r );

    # evaluate rel: h4 h2 = h2 h4
    l := HallMult(s, h[4], h[2]);
    r := HallMult(s, h[2], h[4]);
    Add( res, l-r );

    # evaluate rel: h4 h3 = h3 h4
    l := HallMult(s, h[4], h[3]);
    r := HallMult(s, h[3], h[4]);
    Add( res, l-r );

    # summarize
    res := Flat(res);
    res := Filtered(res, x -> x <> 0*x);
    res := Permuted(res, (2,5,6,4)(3,7));

    # adjust
    res[5] := Value( res[5], [t[1][2][3]], [(res[1]+t[1][2][3]*h[3][3])/h[3][3]] );
    res[6] := Value( res[6], [t[1][3][4]], [(res[2]+t[1][3][4]*h[4][4])/h[4][4]] );
    res[7] := Value( res[7], [t[1][2][3]], [(res[1]+t[1][2][3]*h[3][3])/h[3][3]] );
    res[7] := Value( res[7], [t[1][2][4]], [(res[5]+t[1][2][4]*h[4][4])/h[4][4]] );

    return res;
end;


############################################################################
#                    Isomorphisms 212
############################################################################

ZpIsom212 := function()
    local t, s, h, w, i, j, res, l, r, rev, l1, l2;

    t := List([1..5], x -> List([1..5], y -> List([1..5], z ->
    Indeterminate(Rationals,
        Concatenation("t_{", String(x), String(y), String(z),"}")))));
    t[1][3][5] := 0;
    t[2][3][4] := 0;
    t[1][4][5] := 0;
    t[2][4][5] := 0;
    t[3][4][5] := 0;

    s := List([1..5], x -> List([1..5], y -> List([1..5], z ->
    Indeterminate(Rationals,
        Concatenation("s_{", String(x), String(y), String(z),"}")))));
    s[1][3][5] := 0;
    s[2][3][4] := 0;
    s[1][4][5] := 0;
    s[2][4][5] := 0;
    s[3][4][5] := 0;

    # the isomorphism
    h := NullMat(5,5);
    w := [ [3,1],[3,2], [4,1],[4,2],[4,3], [5,1],[5,2],[5,3] ];
    for i in [1..5] do
        for j in [1..5] do
            if not [i,j] in w then
                h[i][j] := Indeterminate(Rationals,
                   Concatenation("m_{",String(i),String(j),"}"));
            fi;
        od;
    od;

    res := [];

    # evaluate rel: h2 h1 = h1 h2 h3^t123 h4^t124 h5^t125
    l := HallMult(s, h[2], h[1]);
    r := HallMult(s, h[1], h[2]);
    r := HallMult(s, r, t[1][2][3]*h[3]);
    r := HallMult(s, r, t[1][2][4]*h[4]);
    r := HallMult(s, r, t[1][2][5]*h[5]);
    Add( res, l-r );

    # evaluate rel: h3 h1 = h1 h3 h4^t134
    l := HallMult(s, h[3], h[1]);
    r := HallMult(s, h[1], h[3]);
    r := HallMult(s, r, t[1][3][4]*h[4]);
    Add( res, l-r );

    # evaluate rel: h3 h2 = h2 h3 h5^t235
    l := HallMult(s, h[3], h[2]);
    r := HallMult(s, h[2], h[3]);
    r := HallMult(s, r, t[2][3][5]*h[5]);
    Add( res, l-r );

    # summarize
    res := Flat(res);
    res := Filtered(res, x -> x <> 0*x);
    res := Permuted(res, (2,6,4)(3,7,5) );

    # adjust
    res[6] := Value( res[6], [t[1][2][3]], [(res[1]+t[1][2][3]*h[3][3])/h[3][3]] );
    res[7] := Value( res[7], [t[1][2][3]], [(res[1]+t[1][2][3]*h[3][3])/h[3][3]] );

    return res;
end;