text:
############################################################################
##
##  recognize.gi                      IRREDSOL                         Burkhard Höfling
##
##  Copyright © 2003–2016 Burkhard Höfling
##

############################################################################
##
#F  IsAvailableIdIrreducibleSolubleMatrixGroup(<G>)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(IsAvailableIdIrreducibleSolubleMatrixGroup,
    function(G)
        return IsMatrixGroup(G) and IsFinite(FieldOfMatrixGroup(G))
            and IsSolvableGroup(G)  and IsIrreducibleMatrixGroup(G)
            and IsAvailableIrreducibleSolubleGroupData(
                DegreeOfMatrixGroup(G), Size(TraceField(G)));
    end);

############################################################################
##
#F  IsAvailableIdAbsolutelyIrreducibleSolubleMatrixGroup(<G>)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(IsAvailableIdAbsolutelyIrreducibleSolubleMatrixGroup,
    function(G)
        return IsMatrixGroup(G) and IsFinite(FieldOfMatrixGroup(G))
            and IsSolvableGroup(G) and IsAbsolutelyIrreducibleMatrixGroup(G)
            and IsAvailableAbsolutelyIrreducibleSolubleGroupData(
                DegreeOfMatrixGroup(G), Size(TraceField(G)));
    end);

############################################################################
##
#F  FingerprintDerivedSeries(<n>)
##
InstallMethod(FingerprintDerivedSeries, "generic method for finite group", true,
    [IsGroup and IsFinite], 0,
    function(G)
        local EncodedPrimePower, EncodedAbInv, primes, max, der;

        EncodedPrimePower := function(n, primes)

            local res, i, p;

            res := -1; # count p^1 as zero

            i := 1;
            while n mod primes[i][1] <> 0 do
                 res := res + primes[i][2];
                 i := i + 1;
            od;
            p := primes[i][1];
            repeat
                    n := n / p;
                    res := res + 1;
            until n mod p <> 0;

            if n <> 1 then
                Error("cannot encode prime power n");
            fi;
            return res;
        end;

        EncodedAbInv := function (inv, primes, max)

            local i, res;
            res := EncodedPrimePower(inv[1], primes);
            for i in [2..Length(inv)] do
                res := res * max + EncodedPrimePower(inv[i], primes);
            od;
            return res;
        end;

        primes := Collected(FactorsInt(Size(G)));
        max := Sum(primes, p -> p[2]);
        der := DerivedSeries(G);
        return List(der{[1..Length(der)-1]}, N -> EncodedAbInv(AbelianInvariants(N), primes, max));
    end);

RedispatchOnCondition(FingerprintDerivedSeries, true,
    [IsGroup], [IsFinite], 0);

############################################################################
##
#M  FingerprintMatrixGroup(<G>)
##
InstallMethod(FingerprintMatrixGroup, "for irreducible FFE matrix group", true,
    [IsMatrixGroup and CategoryCollections(IsFFECollColl)], 0,
    function(G)

    local EncodedPrimeDecomposition, primes, max, ids, len, ccl, cl, id, pos, rep, i, g, q,
        hsize, hstep, counts, fppoly;

    EncodedPrimeDecomposition := function(n, primes)

        local res, i, p;

        res := 0;

        for i in [1..Length(primes)] do
            res := res * (primes[i][2]+1);
            p := primes[i][1];
            while n mod p = 0 do
                n := n / p;
                res := res + 1;
            od;
        od;
        if n <> 1 then
            Error("cannot encode integer n");
        fi;
        return res;
    end;

    primes := Collected(FactorsInt(Size(G)));
    max := Product(primes, p -> p[2]+1);
    q := Size(TraceField(G));
    rep := RepresentationIsomorphism(G);
    ccl := AttributeValueNotSet(ConjugacyClasses, Source(rep));
    hsize := NextPrimeInt(2*Length(ccl));
    hstep := NextPrimeInt(RootInt(hsize));
    ids := EmptyPlist(hsize);
    counts := EmptyPlist(hsize);
    len := 0;
    for cl in ccl do
        g := Representative(cl);
        id := EncodedPrimeDecomposition(Size(cl), primes)
                + max * (EncodedPrimeDecomposition(Order(g), primes)
                        + max * CodeCharacteristicPolynomial(ImageElm(rep, g), q));
        pos := id mod hsize + 1;
        while IsBound(ids[pos]) and ids[pos] <> id do
            pos := pos + hstep;
            if pos > hsize then
                pos := pos - hsize;
            fi;
        od;
        if IsBound(ids[pos]) then
            counts[pos] := counts[pos] + 1;
        else
            ids[pos] := id;
            counts[pos] := 1;
        fi;
    od;
    fppoly := [];
    for id in AsSet(ids) do
        pos := id mod hsize + 1;
        while IsBound(ids[pos]) and ids[pos] <> id do
            pos := pos + hstep;
            if pos > hsize then
                pos := pos - hsize;
            fi;
        od;
        Add(fppoly, [id, counts[pos]]);
    od;
    return fppoly;
end);

############################################################################
##
#F  ConjugatingMatIrreducibleRepOrFail(repG, repH, q, linonly, maxcost, limit_CMIROF)
##
##  computes a record x such that G^x.iso = H and (G^homG)^x.mat = H^homH.
##  If no such matrix exists, the function returns fail.
##
##  maxcost and limit_CMIROF specify when the function will abort the search, returning
##  the string "too expensive".
##
##  This works like the GAP function Morphium but only looks for linear
##  isomorphisms.
##
InstallGlobalFunction(ConjugatingMatIrreducibleRepOrFail,
    function(repG, repH, q, maxcost, limit_CMIROF)

    local EncodedPrimeDecomposition, G, H, ccl, cl, hsize, hstep, ids, id, primes,
        fac, gens, genspos, imglist, g, i, j, k, len,
        cost, weight, mingens, mincost, mingenspos, totalcost, pos, d, ind, imgs, hom, count,
        sizes, D, proj, homG, homH, group, cent, C, centsize, imgreps, matrep, orb, S, occ, cpH,
        moduleG, moduleH, dirr, mat;

    EncodedPrimeDecomposition := function(n, primes)

        local res, i, p;

        res := 0;

        for i in [1..Length(primes)] do
            res := res * (primes[i][2]+1);
            p := primes[i][1];
            while n mod p = 0 do
                n := n / p;
                res := res + 1;
            od;
        od;
        if n <> 1 then
            Error("cannot encode integer n");
        fi;
        return res;
    end;

    G := repG.group;
    H := repH.group;
    primes := Collected(FactorsInt(Size(G)));
    fac := Sum(primes, p -> p[2]);
    if not IsBound(repG.classes) then
        repG.classes := AttributeValueNotSet(ConjugacyClasses, G);
    fi;
    ccl := repG.classes;
    hsize := NextPrimeInt(2*Length(ccl));
    hstep := NextPrimeInt(RootInt(hsize));
    ids := EmptyPlist(hsize);
    occ := EmptyPlist(hsize);
    cost := EmptyPlist(hsize);
    len := 0;
    for i in [1..Length(ccl)] do
        cl := ccl[i];
        id := [Size(cl), Order(Representative(cl)),
            CodeCharacteristicPolynomial(ImageElm(repG.rep, Representative(cl)), q)];
        pos :=((id[3] * fac + EncodedPrimeDecomposition(id[1], primes)) * fac
                + EncodedPrimeDecomposition(id[2], primes)) mod hsize + 1;
        while IsBound(ids[pos]) and ids[pos]<>id do
            pos := pos + hstep;
            if pos > hsize then
                pos := pos - hsize;
            fi;
        od;
        if IsBound(ids[pos]) then
            Add(occ[pos], i);
        else
            ids[pos] := id;
            occ[pos] := [i];
        fi;
    od;

    for i in [1..Length(ids)] do
        if IsBound(ids[i]) then
            cost[i] := ids[i][1]*Length(occ[i]);
        fi;
    od;

    gens := ShallowCopy(MinimalGeneratingSet(G));
    Info(InfoMorph, 1, "ConjugatingMat: MinimalGeneratingSet has ", Length(gens), " generators");

    mincost := 1;
    mingenspos := [];

    for i in [1..Length(gens)] do
        g := gens[i];
        id := [Size(ConjugacyClass(G, g)), Order(g),
            CodeCharacteristicPolynomial(ImageElm(repG.rep, g), q)];
        pos :=((id[3] * fac + EncodedPrimeDecomposition(id[1], primes)) * fac
                + EncodedPrimeDecomposition(id[2], primes)) mod hsize + 1;
        while IsBound(ids[pos]) and ids[pos]<>id do
            pos := pos + hstep;
            if pos > hsize then
                pos := pos - hsize;
            fi;
        od;
        mincost := mincost + cost[pos];
        Add(mingenspos, pos);
    od;
    mingens := gens;

    Info(InfoMorph, 1, "ConjugatingMat: cost of minimal generating system is ", mincost);

    # try to find better generators

    if mincost > maxcost or mincost > 100 then
        weight := EmptyPlist(hsize);
        weight[1] := 0;
        for i in [2..Length(ids)] do
            if IsBound(cost[i]) then
                weight[i] := weight[i-1] + QuoInt(Size(G), cost[i]);
            else
                weight[i] := weight[i-1];
            fi;
        od;

        for i in [1..LogInt(Size(G),2)*Length(MinimalGeneratingSet(G))] do
            gens := [];
            S := TrivialSubgroup(G);
            totalcost := 1;
            genspos := [];
            repeat
                if Random([1..10]) = 1 then
                    g := Random(G);
                    id := [Size(ConjugacyClass(G, g)), Order(g),
                        CodeCharacteristicPolynomial(ImageElm(repG.rep, g), q)];
                    pos :=((id[3] * fac + EncodedPrimeDecomposition(id[1], primes)) * fac
                            + EncodedPrimeDecomposition(id[2], primes)) mod hsize + 1;
                    while IsBound(ids[pos]) and ids[pos]<>id do
                        pos := pos + hstep;
                        if pos > hsize then
                            pos := pos - hsize;
                        fi;
                    od;
                else
                    k := Random(1, weight[Length(weight)]);
                    pos := PositionSorted(weight, k);
                    g := Random(ccl[Random(occ[pos])]);
                fi;
                if not g in S then
                    Add(gens, g);
                    Add(genspos, pos);
                    totalcost := totalcost + cost[pos];
                    S := ClosureGroup(S, g);
                    #remove for redundant generators - we need not check g
                    j := 1;
                    repeat
                        while j < Length(gens) and Size(Subgroup(S, gens{Difference([1..Length(gens)], [j])})) < Size(S) do
                            j := j + 1;
                        od;
                        if j < Length(gens) then
                            totalcost := totalcost-cost[genspos[j]];
                            Remove(gens, j);
                            Remove(genspos, j);
                        fi;
                    until j >= Length(gens);
                    if totalcost > mincost then
                        break;
                    fi;
                fi;
            until Size(S) = Size(G);

            if totalcost < mincost then
                mingens := gens;
                mincost := totalcost;
                mingenspos := genspos;
                Info(InfoMorph, 2, "ConjugatingMat: found better generating system of cost ", mincost, " after ", i, " attempts");
            else
                Info(InfoMorph, 3, "ConjugatingMat: found generating system of cost ", totalcost);
            fi;
        od;
    fi;

    if mincost > maxcost then
        Info(InfoMorph, 1, "ConjugatingMat: too expoensive");
        return "too expensive";
    fi;
    Info(InfoMorph, 1, "ConjugatingMat: using generating system of cost ", mincost);

    if not IsBound(repH.classes) then
        repH.classes := AttributeValueNotSet(ConjugacyClasses, H);
    fi;
    ccl := repH.classes;

    # compute set of possible generator images

    imglist := [];
    cpH := [];

    # sort, most expensive first (one representative will be picked), then ascending cost
    SortParallel(mingenspos, mingens, function(i, j) return cost[i] > cost[j]; end);

    for i in [1..Length(mingens)] do
        id := ids[mingenspos[i]];
        imglist[i] := [];
        count := Length(occ[mingenspos[i]]);
        for j in [1..Length(ccl)] do
            cl := ccl[j];
            if Size(cl) = id[1] and Order(Representative(cl)) = id[2] then
                if not IsBound(cpH[j]) then
                    cpH[j] := CodeCharacteristicPolynomial(ImageElm(repH.rep, Representative(cl)), q);
                fi;
                if cpH[j] = id[3] then
                    if i = 1 then
                        Add(imglist[i], Representative(cl));
                    else
                        Append(imglist[i], AttributeValueNotSet(AsList, cl));
                    fi;
                    count := count - 1;

                    # if H contains more such conjugacy classes, the groups cannot be
                    # isomorphic so it doesn't matter if we miss some possible images
                    if count = 0 then
                        break;
                    fi;
                fi;
            fi;
        od;
        if IsEmpty(imglist[i]) then
            Info(InfoMorph, 1, "ConjugatingMat: no suitable image ccl found, reps cannot be conjugate");
            return fail;
        fi;
    od;

    gens := List(mingens, g -> ImageElm(repG.rep, g));

    dirr := Length(mingens) + 1;

    moduleG := List([1..Length(mingens)], d -> GModuleByMats(gens{[1..d]}, GF(q)));
    repeat
        dirr := dirr - 1;
    until dirr = 0 or not MTX.IsIrreducible(moduleG[dirr]);

    dirr := dirr + 1;

    if dirr > Length(mingens) then
        Error("repG must be irreducible");
    fi;

    sizes := List([1..Length(mingens)], d -> Size(Group(mingens{[1..d]})));

    D := DirectProduct(G, H);
    homG := Embedding(D, 1);
    homH := Embedding(D, 2);
    proj := Projection(D, 2);

    gens := List(mingens, g -> ImageElm(homG, g));

    C := G;
    centsize := [Size(G)];
    for g in mingens do
        C := Centralizer(C, g);
        Add (centsize, Size(C));
    od;

    # do a backtrack search through the possible images
    d := 1;
    ind := [0];
    imgs := [];
    group := [];
    cent := [H];
    matrep := [];
    imgreps := [imglist[1]];
    count := 0;
    repeat
        count := count + 1;
        ind[d] := ind[d] + 1;
        if count > limit_CMIROF then
            Info(InfoMorph, 1, "ConjugatingMat: exceeded limit  ", limit_CMIROF, " attempts");
            return "too expensive";
        elif count mod 1000 = 0 then
            Info(InfoMorph, 2, "count: ", count, " trying indices: ", ind{[1..d]},
                ", lengths: ", List([1..d], t -> Length(imgreps[t])));
        else
            Info(InfoMorph, 3, "step: ", ind{[1..d]});
        fi;

        imgs[d] := imgreps[d][ind[d]];
        matrep[d] := fail;

        if d = 1 then
            group[d] := Group(gens[d]*ImageElm(homH, imgs[d]));
        else
            group[d] := ClosureGroup(group[d-1], gens[d]*ImageElm(homH, imgs[d]));
        fi;
        if d = 1 or Size(group[d]) = sizes[d] then

            # map from gens to imgs is an isomorphism

            if d >= dirr then
                i := d;
                while i > 0 and matrep[i] = fail do
                    matrep[i] := ImageElm(repH.rep, imgs[i]);
                    i := i - 1;
                od;
                moduleH := GModuleByMats(matrep{[1..d]}, GF(q));
                mat :=  MTX.IsomorphismIrred(moduleG[d], moduleH);
            else
                mat := true;
            fi;
            if mat <> fail then
                if d < Length(gens) then
                    d := d + 1;
                    ind[d] := 0;
                    cent[d] := Centralizer(cent[d-1], imgs[d-1]);
                    if Size(cent[d]) = centsize[d] then
                        orb := OrbitsDomain(cent[d], imglist[d], OnPoints);
                        Info(InfoMorph, 3, "ConjugatingMat: ", Length(orb), " orbits found");
                        orb := Filtered (orb, o -> Length(o)=centsize[d]/centsize[d+1]);
                        imgreps[d] := List (orb, o -> o[1]);
                        Info(InfoMorph, 3, "ConjugatingMat: depth ", d, ": ",
                            Length(orb), " orbits of expected length found");
                    else
                        Info(InfoMorph, 3, "ConjugatingMat: centraliser sizes don't match");
                        imgreps[d] := [];
                    fi;
                else
                    hom := GroupGeneralMappingByImagesNC(G, H,
                        List(gens, g -> PreImagesRepresentative(homG, g)), imgs);
                    if not IsGroupHomomorphism(hom) or not IsBijective(hom) then
                        Error("hom is not an isomophism");
                    fi;
                    Info(InfoMorph, 1, "ConjugatingMat: found after ", count, " attempts");
                    Info(InfoMorph, 2, "ConjugatingMat: found indices ", ind);
                    return rec(iso :=hom, mat := mat);
                fi;
            else
                Info(InfoMorph, 3, "step: ", ind{[1..d]}, " linearity fails");
            fi;
        fi;
        while d > 0 and ind[d] = Length(imgreps[d]) do
            d := d - 1;
        od;
    until d = 0;
    Info(InfoMorph, 1, "ConjugatingMat: reps cannot be conjugate,  ", count, " attempts");

    return fail;
end);

############################################################################
##
#F  ConjugatingMatIrreducibleOrFail(G, H, F)
##
##  G and H must be irreducible matrix groups over the finite field F
##
##  computes a record x such that G^x.mat = H and x.iso is an isomorphism
##  from the source of RepresentationIsomorphism(G) to the source of
##  RepresentationIsomorphism(H). If no such matrix exists, the
##  function returns fail.
##  this works like the GAP function Morphium but only looks for linear
##  isomorphisms.
##
InstallGlobalFunction(ConjugatingMatIrreducibleOrFail, function(G, H, F)

    local q, repG, repH;

    q := Size(F);
    repG := rec(rep := RepresentationIsomorphism(G));
    repG.group := Source(repG.rep);
    repG.classes := ConjugacyClasses(repG.group);

    repH := rec(rep := RepresentationIsomorphism(H));
    repH.group := Source(repH.rep);
    repH.classes := ConjugacyClasses(repH.group);

    return ConjugatingMatIrreducibleRepOrFail(repG, repH, q, infinity, infinity);
end);

############################################################################
##
#F  ConjugatingMatImprimitiveOrFail(G, H, d, F)
##
##  G and H must be irreducible matrix groups over the finite field F
##  H must be block monomial with block dimension d
##
##  computes record x with a single component x.mat such that G^x.mat = H
##  or returns fail if no such x.mat exists
##
##  The function works best if d is small. Irreducibility is only requried
##  if ConjugatingMatIrreducibleOrFail is used
##
InstallGlobalFunction(ConjugatingMatImprimitiveOrFail, function(G, H, d, F)

# H must have minimal block dimension d and must be a block monomial group for that block size

    local n, systemsG, W, hom, r, basis, orb, posbasis, act, i, blocks,
        permW, permH, sys, permG, gensG, C, Cinv, mat;

    n := DegreeOfMatrixGroup(H);
    systemsG := ImprimitivitySystems(G);

    if Minimum(List(systemsG, sys -> Length(sys.bases[1]))) <> d then
        Info(InfoIrredsol, 1, "different minimal block dimension - groups are not conjugate");
        return fail;
    fi;

    if d = n then
        if Size(G) < 100000 then
            return ConjugatingMatIrreducibleOrFail(G, H, F);
        fi;
        hom := NiceMonomorphism(GL(n, Size(F)));
        r := RepresentativeAction(ImagesSource(hom),
                ImagesSet(hom, G), ImagesSet(hom, H));
        if r <> fail then
            Info(InfoIrredsol, 1, " conjugating matrix found");
            return rec( mat := PreImagesRepresentative(hom, r));
        else
            Info(InfoIrredsol, 1, "groups are not conjugate");
            return fail;
        fi;
    else
        W := WreathProductOfMatrixGroup(GL(d, Size(F)), SymmetricGroup(n/d));
        basis := IdentityMat(n, F);
        orb := Orbit(W, basis[1], OnRight);
        posbasis := List(basis, b -> Position(orb, b));
        Assert(1, not fail in posbasis);
        permW := Group(List(GeneratorsOfGroup(W), g -> Permutation(g, orb, OnRight)));

        blocks := [];
        for i in [0,d..n-d] do
            Add(blocks, basis{[i+1..i+d]});
        od;
        act := TransitiveIdentification(Action(H, blocks, OnSubspacesByCanonicalBasis));

        permH := Group(List(GeneratorsOfGroup(H), g -> Permutation(g, orb, OnRight)));
        Assert(1, IsSubgroup(permW, permH));

        for sys in systemsG do
            if Length(sys.bases[1]) = d
                    and TransitiveIdentification(
                        Action(G, sys.bases, OnSubspacesByCanonicalBasis)) = act then
                C := Concatenation(sys.bases);
                Cinv := C^-1;
                gensG := List(GeneratorsOfGroup(G), B -> C*B*Cinv);
                permG := Group(List(gensG, g->Permutation(g, orb, OnRight)));
                Assert(1, IsSubgroup(permW, permG));
                r := RepresentativeAction(permW, permG, permH);
                if r <> fail then
                    # compute the corresponding matrix
                    mat := List(posbasis, k -> orb[k^r]);
                    r := Cinv * mat;
                    Info(InfoIrredsol, 1, " conjugating matrix found");
                    return rec(mat := r);
                fi;
            fi;
        od;
        Info(InfoIrredsol, 1, "groups are not conjugate");
        return fail;
    fi;
end);

############################################################################
##
#F  RecognitionAISMatrixGroup(G, inds, wantmat, wantgroup, wantiso)
##
##  version of RecognitionIrreducibleSolubleMatrixGroupNC which
##  only works for absolutely irreducible groups G. This version
##  allows to prescribe a set of absolutely irreducible subgroups
##  to which G is compared. This set is described as a subset <inds> of
##  IndicesAbsolutelyIrreducibleSolubleMatrixGroups(n, q), where n is the
##  degree of G and q is the order of the trace field of G. if inds is fail,
##  all groups in the IRREDSOL library are considered.
##
##  WARNING: The result may be wrong if G is not among the groups
##  described by <inds>.
##
InstallGlobalFunction(RecognitionAISMatrixGroup,
    function(G, inds, wantmat, wantgroup, wantiso)

        local allinds, newinds, n, q, order, info, fppos, fpinfo, elminds, pos, t, tinv,
            F, H, systems, rep, i, x;

        F := TraceField(G);
        n := DegreeOfMatrixGroup(G);
        q := Size(F);
        order := Size(G);

        info := rec(); # set up the answer

        allinds := inds = fail;

        inds := SelectionIrreducibleSolubleMatrixGroups(
            n, q, 1, inds, [order], fail, fail);

        Info(InfoIrredsol, 1, Length(inds),
                    " groups of order ", order, " to compare");

        if Length(inds) = 0 then
            Error("panic: no group of order ", order, " in the IRREDSOL library");
        elif Length(inds) > 1 then
            # cut down candidate grops by looking at fingerprints
            if not TryLoadAbsolutelyIrreducibleSolubleGroupFingerprintIndex(n, q) then
                fppos := fail;
            else
                fppos := PositionSet(IRREDSOL_DATA.FP_INDEX[n][q][1], order);

                if fppos <> fail then # fingerprint info available
                    LoadAbsolutelyIrreducibleSolubleGroupFingerprintData(
                        n, q, IRREDSOL_DATA.FP_INDEX[n][q][2][fppos]);

                    fpinfo := IRREDSOL_DATA.FP[n][q][fppos];        # fingerprint info

                    if Length(fpinfo)=2 then # new format using FingerprintDerivedSubgroup
                        if IsBound(fpinfo[1]) then
                            pos := PositionSorted(fpinfo[1],
                                FingerprintDerivedSeries(Source(RepresentationIsomorphism(G))));
                            if pos = fail then
                                Error("panic: FingerprintDerivedSeries not in database");
                            fi;
                        else
                            pos := 1;
                        fi;
                        fpinfo := fpinfo[2][pos];
                        if IsInt(fpinfo) then
                            fpinfo := [[], [[]], [fpinfo]];
                        fi;
                    fi;

                    if Length(fpinfo) <> 3 then
                        Error("panic: error in fingerprint database");
                    fi;

                    Info(InfoIrredsol, 1, "computing fingerprint of group");

                    # which distingushing elements are in G?

                    if IsEmpty(fpinfo[1]) then
                        Info(InfoIrredsol, 1, "all groups have the same fingerprint");
                        # this is enough, same code handles this special case
                        elminds := [];
                    else
                        Info(InfoIrredsol, 1, "computing fingerprint of group");
                        elminds := Filtered([1..Length(fpinfo[1])], i ->
                            fpinfo[1][i] in FingerprintMatrixGroup(G));
                    fi;

                    # find relevant info in database

                    pos := PositionSet(fpinfo[2], elminds);

                    if pos = fail then
                        Error("panic: FingerprintMatrixGroup not found in database");
                    fi;

                    newinds := fpinfo[3][pos];
                    if IsInt(newinds) then
                        newinds := [newinds];
                    fi;
                    if allinds then
                        if not IsSubset(inds, newinds) then
                            Error("panic: fingerprint indices do not match the IRREDSOL library indices");
                        fi;
                        inds := newinds;
                    else
                        inds := Intersection(inds, newinds);
                    fi;

                    Info(InfoIrredsol, 1, Length(inds),
                        " groups in the library have the same fingerprint");
                fi;
            fi;
        fi;

        if Length(inds) > 1 or wantmat or wantiso then
            # rewrite G over smaller field if necc.

            Info(InfoIrredsol, 1, "rewriting group over its trace field");

            t := ConjugatingMatTraceField(G);
            if t <> One(G) then
                tinv := t^-1;
                H := GroupWithGenerators(List(GeneratorsOfGroup(G), g -> tinv * g * t));
                Assert(1, FieldOfMatrixGroup(H) = F);
                SetFieldOfMatrixGroup(H, F);
                SetTraceField(H, F);
                rep := RepresentationIsomorphism(G);
                SetRepresentationIsomorphism(H,
                        GroupHomomorphismByFunction(Source(rep), H,
                            g -> tinv * ImageElm(rep, g)*t, h -> PreImagesRepresentative(rep, t*h*tinv)));
                G := H;

                # we need the imprimitivity systems of G later anyway, so we can rule out groups
                # with different minimal block dimensions as well

                if Length(inds) > 1 then
                    Info(InfoIrredsol, 1, "computing imprimitivity systems of group");
                    systems := ImprimitivitySystems(G);
                    inds := SelectionIrreducibleSolubleMatrixGroups(
                        DegreeOfMatrixGroup(G), Size(F), 1, inds, [Order(G)],
                        [MinimalBlockDimensionOfMatrixGroup(G)], fail);
                fi;
                Info(InfoIrredsol, 1, Length(inds),
                    " groups also have the same minimal block dimension");
            fi;

            # find possible groups in the library

            for i in [1..Length(inds)-1] do
                H := IrreducibleSolubleMatrixGroup(DegreeOfMatrixGroup(G), Size(F), 1, inds[i]);

                if fppos = fail then # compare fingerprints now
                    Info(InfoIrredsol, 1, "comparing fingerprints");
                    if FingerprintMatrixGroup(G) <> FingerprintMatrixGroup(H) then
                        Info(InfoIrredsol, 1, "fingerprint different from id ",
                            IdIrreducibleSolubleMatrixGroup(H));
                        continue;
                    fi;
                fi;

                if wantiso then
                     x := ConjugatingMatIrreducibleOrFail(G, H, F);
                else
                     x := ConjugatingMatImprimitiveOrFail(G, H,
                          MinimalBlockDimensionOfMatrixGroup(G), F);
                fi;
                if x = fail then
                    Info(InfoIrredsol, 1, "group does not have id ", IdIrreducibleSolubleMatrixGroup(H));
                else
                    Assert(1, G^(x.mat) = H);
                    info.id := [DegreeOfMatrixGroup(G), Size(F), 1, inds[i]];
                    if wantgroup then
                        info.group := H;
                    fi;
                    if wantmat then
                        info.mat := t*x.mat;
                    fi;
                    if wantiso then
                        info.iso := x.iso;
                    fi;
                  Info(InfoIrredsol, 1, "group id is ", info.id);
                    return info;
                fi;
            od;
        fi;

        # we are down to the last group - this must be it

        i := inds[Length(inds)];
        info.id := [DegreeOfMatrixGroup(G), Size(F), 1, i];
        Info(InfoIrredsol, 1, "group id is ", info.id);
        if not wantgroup and not wantmat then
            return info;
        fi;
        H := IrreducibleSolubleMatrixGroup(DegreeOfMatrixGroup(G), Size(F), 1, i);
        if wantgroup then
            info.group := H;
        fi;
        if wantmat or wantiso then
            if wantiso then
                x := ConjugatingMatIrreducibleOrFail(G, H, F);
            else
                x := ConjugatingMatImprimitiveOrFail(G, H,
                          MinimalBlockDimensionOfMatrixGroup(G), F);
            fi;
            if x = fail then
                Error("panic: group not found in database");
            fi;
            if wantmat then
                info.mat := t*x.mat;
            fi;
            if wantiso then
                info.iso := x.iso;
            fi;
        fi;
        return info;

    end);

############################################################################
##
#F  RecognitionIrreducibleSolubleMatrixGroup(G, wantmat, wantgroup)
##
##  Let G be an irreducible soluble matrix group over a finite field.
##  This function identifies a conjugate H of G group in the library.
##
##  It returns a record which has the following entries:
##  id:                     contains the id of H (and thus of G),
##                            cf. IdIrreducibleSolubleMatrixGroup
##  mat: (optional)     a matrix x such that G^x = H
##  group: (optional)  the group H
##
##  The entries mat and group are only present if the booleans wantmat and/or
##  wantgroup are true, respectively.
##
##  Currently, wantmat may only be true if G is absolutely irreducible.
##
##  Note that in most cases, the function will be much slower if wantmat
##  is set to true.
##
InstallGlobalFunction(RecognitionIrreducibleSolubleMatrixGroup,
    function(arg)

        local G, wantmat, wantgroup, wantiso, info;

        if Length(arg) < 2 or Length(arg) > 4 then
          Error("RecognitionIrreducibleSolubleMatrixGroup(G [, wantmat [, wantgroup [,wantiso]]])");
        fi;

        # test if G is soluble

        G := arg[1];
        if not IsMatrixGroup(G) or not IsFinite (FieldOfMatrixGroup(G))
                or not IsSolvableGroup(G) then
            Error("G must be a soluble matrix group over a finite field");
        fi;

        wantmat := arg[2];
        if Length(arg) = 2 then
            wantgroup := false;
        else
            wantgroup := arg[3];
            if Length(arg) = 3 then
                wantiso := false;
            else
                wantiso := arg[4];
            fi;
        fi;
        if not IsBool (wantmat) or not IsBool (wantgroup) or not IsBool (wantiso) then
            Error("wantmat, wantgroup and wantiso must be booleans");
        fi;
        info := RecognitionIrreducibleSolubleMatrixGroupNC(G,
            wantmat, wantgroup, wantiso);
        if info = fail then
            Error("This group is not within the scope of the IRREDSOL library");
        fi;
        return info;
    end);

############################################################################
##
#F  RecognitionIrreducibleSolubleMatrixGroupNC(G, wantmat, wantgroup, wantiso)
##
##  version of RecognitionIrreducibleSolubleMatrixGroup which does not check
##  its arguments and returns fail if G is not within the scope of the
##  IRREDSOL library
##
InstallGlobalFunction(RecognitionIrreducibleSolubleMatrixGroupNC,
    function(G, wantmat, wantgroup, wantiso)

        local moduleG, module, info, newinfo, conjinfo, gens, perm_pow, basis,
            gensG, repG, repH, H, d, e, q, gal;

        # reduce to the absolutely irreducible case

        repG := RepresentationIsomorphism (G);
        gensG := List(GeneratorsOfGroup(Source(repG)), x -> ImageElm(repG, x));
        moduleG := GModuleByMats (gensG, FieldOfMatrixGroup(G));

        if not MTX.IsIrreducible (moduleG) then
            Error("G must be irreducible over FieldOfMatrixGroup(G)");
        elif MTX.IsAbsolutelyIrreducible (moduleG) then
            Info(InfoIrredsol, 1, "group is absolutely irreducible");
            info := RecognitionAISMatrixGroup(G, fail, wantmat, wantgroup, wantiso);
            if info = fail then
                return fail;
            else
                newinfo := rec(id := [DegreeOfMatrixGroup(G), Size(TraceField(G)), 1, info.id[4]]);
                if wantgroup then
                    newinfo.group := info.group;
                fi;
                if wantmat then
                    newinfo.mat := info.mat;
                fi;
                if wantiso then
                    newinfo.iso := info.iso;
                fi;
             fi;
        else
            Info(InfoIrredsol, 1, "reducing to the absolutely irreducible case");
            module := GModuleByMats (gensG, GF(CharacteristicOfField (G)^MTX.DegreeSplittingField (moduleG)));
            repeat
                basis := MTX.ProperSubmoduleBasis (module);
                module := MTX.InducedActionSubmodule (module, basis);
            until MTX.IsIrreducible (module);
            Assert(1, MTX.IsAbsolutelyIrreducible (module));

            # construct absolutely irreducible group

            H := Group(MTX.Generators (module));
            SetIsAbsolutelyIrreducibleMatrixGroup(H, true);
            SetIsIrreducibleMatrixGroup(H, true);
            SetIsSolvableGroup(H, true);
            SetSize(H, Size(G));
            e := LogInt(Size(TraceField (H)), CharacteristicOfField (H));
            d := DegreeOfMatrixGroup(G)/MTX.Dimension (module);
            Assert(1, e mod d = 0);

            # construct representation isomorphism for H

            repH := GroupGeneralMappingByImagesNC (Source (repG), H,
                    GeneratorsOfGroup(Source(repG)),
                    MTX.Generators (module));
            SetIsGroupHomomorphism (repH, true);
            SetIsBijective (repH, true);
            SetRepresentationIsomorphism (H, repH);

            # recognize absolutely irreducible component

            info := RecognitionAISMatrixGroup(H, fail, wantmat, false, wantiso);
            if info = fail then
                return fail;
            fi;

            # translate results back

            q := CharacteristicOfField(H)^(e/d);
            SetTraceField (G, GF(q));
            Assert(1, q^d = info.id[2]);
            perm_pow := PermCanonicalIndexIrreducibleSolubleMatrixGroup(
                    DegreeOfMatrixGroup(G), q, d, info.id[4]);

            newinfo := rec(
                id := [DegreeOfMatrixGroup(G), q, d, info.id[4]^(perm_pow.perm^perm_pow.pow)]);
            if wantgroup then

                newinfo.group := CallFuncList(IrreducibleSolubleMatrixGroup, newinfo.id);
            fi;

            if wantmat then
                Info(InfoIrredsol, 1, "determining conjugating matrix");
                # raising to gal-th power is the field automorphism which maps H to a conjugate of
                # the absolutely irreducible group which is used to construct the irreducible group
                # we want to find
                gal := q^perm_pow.pow;
                if gal = 1 then
                     gens := List(MTX.Generators(module), mat -> mat^info.mat);
                     Info(InfoIrredsol, 1, "already got right Galois conjugate");
                     if wantiso then
                        if Source (RepresentationIsomorphism (newinfo.group)) <> Range (info.iso) then
                            Error("sources of isomorphisms don't match");
                        fi;
                        newinfo.iso := info.iso;
                    fi;
               else
                    # construct Galois conjugate of module

                    Info(InfoIrredsol, 1, "computing and recognizing Galois conjugate");

                    module := GModuleByMats (
                        List(MTX.Generators (module), mat -> List(mat, row -> List(row, x -> x^gal))),
                        MTX.Dimension (module), MTX.Field(module));

                    H := Group(MTX.Generators (module));
                    SetIsAbsolutelyIrreducibleMatrixGroup(H, true);
                    SetIsIrreducibleMatrixGroup(H, true);
                    SetIsSolvableGroup(H, true);
                    SetSize(H, Size(G));

                    # construct representation isomorphism for H

                    repH := GroupGeneralMappingByImagesNC (Source (repG), H,
                            GeneratorsOfGroup(Source(repG)),
                            MTX.Generators (module));
                    SetIsGroupHomomorphism (repH, true);
                    SetIsBijective (repH, true);
                    SetRepresentationIsomorphism (H, repH);

                    # recognize absolutely irreducible component

                    conjinfo := RecognitionAISMatrixGroup(H, [perm_pow.min], true, false, wantiso);
                    if conjinfo = fail then
                        Error("panic: internal error, Galois conjugate isn't in the library");
                    elif conjinfo.id <> [info.id[1], info.id[2], 1, perm_pow.min] then
                        Error("panic: internal error, didn't find correct Galois conjugate");
                    fi;
                    gens := List(MTX.Generators (module), mat -> mat^conjinfo.mat);
                    if wantiso then
                        if Source (RepresentationIsomorphism (newinfo.group)) <> Range (conjinfo.iso) then
                            Error("sources of isomorphisms don't match");
                        fi;
                        newinfo.iso := conjinfo.iso;
                    fi;
                fi;

                basis := CanonicalBasis (AsVectorSpace (GF(q), GF(q^d)));
                      # it is important to use CanonicalBasis here, in order to be sure
                      # that the result is the same as during the construction of
                      # the corresponding irreducible group

                # now construct a module over GF(q)
                 module := GModuleByMats (
                     List(gens, mat -> BlownUpMat (basis, mat)),
                     MTX.Dimension (moduleG),
                     MTX.Field (moduleG));

                 newinfo.mat := MTX.Isomorphism (moduleG, module);
                    if newinfo.mat = fail then
                        Error("panic: no conjugating mat found");
                    fi;
            fi;

        fi;
        Info(InfoIrredsol, 1, "irreducible group id is", newinfo.id);
        return newinfo;
    end);

############################################################################
##
#M  IdIrreducibleSolubleMatrixGroup(<G>)
##
##  see the IRREDSOL manual
##
InstallMethod(IdIrreducibleSolubleMatrixGroup,
    "for irreducible soluble matrix group over finite field", true,
    [IsFFEMatrixGroup
        and IsIrreducibleMatrixGroup and IsSolvableGroup], 0,

    function(G)

        local info;

        info := RecognitionIrreducibleSolubleMatrixGroupNC(G, false, false, false);
        if info = fail then
            Error("group is beyond the scope of the IRREDSOL library");
        else
            return info.id;
        fi;
    end);

RedispatchOnCondition(IdIrreducibleSolubleMatrixGroup, true,
    [IsFFEMatrixGroup],
    [IsIrreducibleMatrixGroup and IsSolvableGroup], 0);

############################################################################
##
#V  IRREDSOL_DATA.MS_GROUP_INDEX
##
##  translation table for Mark Short's irredsol library
##
IRREDSOL_DATA.MS_GROUP_INDEX :=
[ ,
  [ , [ [ 2, 1 ], [ 1, 1 ] ], [ [ 2, 1 ], [ 1, 1 ], [ 1, 5 ], [ 2, 2 ], [ 1,
                  4 ], [ 1, 3 ], [ 1, 2 ] ],,
        [ [ 2, 1 ], [ 1, 11 ], [ 2, 2 ], [ 1, 5 ], [ 1, 4 ], [ 2, 3 ],
             [ 1, 8 ], [ 1, 9 ], [ 2, 4 ], [ 1, 3 ], [ 1, 2 ], [ 1, 10 ],
             [ 1, 7 ], [ 2, 5 ], [ 1, 14 ], [ 1, 1 ], [ 1, 6 ], [ 1, 13 ],
             [ 1, 12 ] ],,
        [ [ 2, 1 ], [ 1, 7 ], [ 1, 10 ], [ 1, 18 ], [ 2, 2 ], [ 1, 8 ],
             [ 1, 6 ], [ 2, 3 ], [ 1, 14 ], [ 2, 4 ], [ 1, 16 ], [ 1, 4 ],
             [ 1, 9 ], [ 1, 5 ], [ 2, 5 ], [ 1, 17 ], [ 1, 21 ], [ 1, 22 ],
             [ 1, 13 ], [ 1, 3 ], [ 1, 2 ], [ 1, 15 ], [ 2, 6 ], [ 1, 12 ],
             [ 1, 23 ], [ 1, 1 ], [ 1, 20 ], [ 1, 11 ], [ 1, 19 ] ],,,,
        [ [ 2, 1 ], [ 2, 2 ], [ 2, 3 ], [ 1, 21 ], [ 1, 10 ], [ 2, 4 ],
             [ 1, 26 ], [ 1, 6 ], [ 2, 5 ], [ 1, 20 ], [ 1, 17 ], [ 2, 6 ],
             [ 1, 25 ], [ 1, 5 ], [ 1, 7 ], [ 2, 7 ], [ 1, 14 ], [ 2, 8 ],
             [ 1, 16 ], [ 1, 29 ], [ 2, 9 ], [ 1, 22 ], [ 1, 9 ], [ 1, 8 ],
             [ 1, 24 ], [ 2, 10 ], [ 1, 13 ], [ 1, 30 ], [ 1, 4 ], [ 1, 18 ],
             [ 1, 19 ], [ 2, 11 ], [ 1, 23 ], [ 1, 2 ], [ 1, 3 ], [ 1, 12 ],
             [ 2, 12 ], [ 1, 15 ], [ 1, 28 ], [ 1, 1 ], [ 1, 11 ], [ 1, 27 ] ],,
        [ [ 1, 26 ], [ 2, 1 ], [ 1, 21 ], [ 1, 23 ], [ 2, 2 ], [ 1, 25 ],
             [ 1, 27 ], [ 2, 3 ], [ 1, 43 ], [ 1, 18 ], [ 1, 29 ], [ 1, 11 ],
             [ 2, 4 ], [ 1, 14 ], [ 1, 20 ], [ 1, 22 ], [ 1, 13 ], [ 1, 19 ],
             [ 1, 30 ], [ 2, 5 ], [ 1, 51 ], [ 1, 50 ], [ 2, 6 ], [ 1, 40 ],
             [ 1, 38 ], [ 1, 32 ], [ 1, 8 ], [ 1, 9 ], [ 2, 7 ], [ 1, 44 ],
             [ 1, 17 ], [ 1, 28 ], [ 1, 15 ], [ 1, 12 ], [ 1, 49 ], [ 1, 48 ],
             [ 2, 8 ], [ 1, 36 ], [ 1, 42 ], [ 1, 5 ], [ 1, 10 ], [ 1, 7 ],
             [ 1, 6 ], [ 1, 47 ], [ 1, 39 ], [ 1, 35 ], [ 2, 9 ], [ 1, 31 ],
             [ 1, 16 ], [ 1, 52 ], [ 1, 37 ], [ 1, 3 ], [ 1, 2 ], [ 1, 46 ],
             [ 2, 10 ], [ 1, 34 ], [ 1, 41 ], [ 1, 1 ], [ 1, 45 ], [ 1, 33 ],
             [ 1, 24 ], [ 1, 4 ] ] ],
  [ , [ [ 3, 1 ], [ 1, 1 ] ], [ [ 1, 5 ], [ 3, 1 ], [ 1, 4 ], [ 1, 3 ],
             [ 1, 2 ], [ 3, 2 ], [ 1, 7 ], [ 1, 1 ], [ 1, 6 ] ],,
        [ [ 1, 12 ], [ 1, 15 ], [ 1, 16 ], [ 1, 6 ], [ 3, 1 ], [ 1, 3 ],
             [ 1, 10 ], [ 1, 11 ], [ 1, 9 ], [ 3, 2 ], [ 1, 19 ], [ 1, 4 ],
             [ 1, 5 ], [ 1, 14 ], [ 1, 13 ], [ 3, 3 ], [ 1, 18 ], [ 1, 8 ],
             [ 1, 2 ], [ 1, 7 ], [ 1, 17 ], [ 1, 1 ] ] ],
  [ , [ [ 4, 1 ], [ 2, 3 ], [ 4, 2 ], [ 2, 1 ], [ 1, 5 ], [ 2, 2 ], [ 1, 2 ],
             [ 1, 3 ], [ 1, 4 ], [ 1, 1 ] ],
        [ [ 4, 1 ], [ 2, 20 ], [ 4, 2 ], [ 2, 11 ], [ 2, 9 ], [ 2, 8 ],
             [ 4, 3 ], [ 2, 12 ], [ 2, 16 ], [ 2, 17 ], [ 4, 4 ], [ 1, 73 ],
             [ 1, 11 ], [ 1, 10 ], [ 2, 4 ], [ 1, 9 ], [ 1, 12 ], [ 1, 34 ],
             [ 1, 32 ], [ 2, 7 ], [ 2, 6 ], [ 1, 33 ], [ 2, 5 ], [ 2, 10 ],
             [ 4, 5 ], [ 1, 71 ], [ 2, 18 ], [ 2, 15 ], [ 1, 72 ], [ 2, 24 ],
             [ 2, 23 ], [ 1, 13 ], [ 1, 14 ], [ 1, 7 ], [ 1, 6 ], [ 2, 3 ],
             [ 1, 36 ], [ 1, 41 ], [ 2, 2 ], [ 1, 27 ], [ 1, 42 ], [ 1, 26 ],
             [ 1, 35 ], [ 2, 19 ], [ 2, 14 ], [ 1, 69 ], [ 4, 6 ], [ 1, 70 ],
             [ 1, 5 ], [ 1, 49 ], [ 1, 38 ], [ 1, 37 ], [ 2, 22 ], [ 2, 25 ],
             [ 1, 57 ], [ 1, 65 ], [ 2, 26 ], [ 1, 8 ], [ 1, 28 ], [ 1, 39 ],
             [ 1, 30 ], [ 1, 40 ], [ 1, 29 ], [ 1, 44 ], [ 1, 43 ], [ 1, 21 ],
             [ 2, 1 ], [ 1, 67 ], [ 1, 68 ], [ 2, 13 ], [ 1, 76 ], [ 1, 3 ],
             [ 1, 4 ], [ 1, 2 ], [ 1, 31 ], [ 1, 59 ], [ 1, 63 ], [ 1, 64 ],
             [ 1, 62 ], [ 2, 21 ], [ 1, 61 ], [ 1, 19 ], [ 1, 45 ], [ 1, 24 ],
             [ 1, 48 ], [ 1, 46 ], [ 1, 47 ], [ 1, 58 ], [ 1, 66 ], [ 1, 75 ],
             [ 1, 1 ], [ 1, 22 ], [ 1, 23 ], [ 1, 60 ], [ 1, 25 ], [ 1, 54 ],
             [ 1, 55 ], [ 1, 56 ], [ 1, 74 ], [ 1, 20 ], [ 1, 18 ], [ 1, 52 ],
             [ 1, 51 ], [ 1, 53 ], [ 1, 17 ], [ 1, 16 ], [ 1, 50 ], [ 1, 15 ] ] ]
     ,
  [ , [ [ 5, 1 ], [ 1, 1 ] ], [ [ 5, 1 ], [ 5, 2 ], [ 1, 11 ], [ 1, 4 ], [ 1,
                  12 ], [ 5, 3 ], [ 1, 3 ], [ 1, 5 ], [ 1, 6 ], [ 5, 4 ],
             [ 1, 7 ], [ 1, 2 ], [ 1, 8 ], [ 1, 10 ], [ 1, 1 ], [ 1, 9 ] ] ],
  [ , [ [ 6, 1 ], [ 3, 2 ], [ 3, 4 ], [ 6, 2 ], [ 2, 4 ], [ 2, 3 ],
             [ 1, 15 ], [ 3, 5 ], [ 2, 5 ], [ 1, 9 ], [ 1, 10 ], [ 6, 3 ],
             [ 2, 7 ], [ 2, 8 ], [ 2, 2 ], [ 3, 1 ], [ 1, 11 ], [ 2, 11 ],
             [ 1, 16 ], [ 3, 3 ], [ 2, 1 ], [ 1, 6 ], [ 1, 7 ], [ 2, 6 ],
             [ 2, 12 ], [ 1, 21 ], [ 1, 20 ], [ 1, 14 ], [ 1, 13 ], [ 1, 8 ],
             [ 1, 5 ], [ 1, 17 ], [ 1, 19 ], [ 1, 3 ], [ 1, 4 ], [ 1, 2 ],
             [ 2, 10 ], [ 1, 12 ], [ 1, 1 ], [ 1, 18 ] ] ],
  [ , [ [ 7, 1 ], [ 1, 1 ] ] ] ];

MakeImmutable(IRREDSOL_DATA.MS_GROUP_INDEX);

############################################################################
##
#F  IdIrreducibleSolubleMatrixGroupIndexMS(<n>, <p>, <k>)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(IdIrreducibleSolubleMatrixGroupIndexMS,
    function(n, p, k)
        if IsInt(n) and n > 1 and IsPosInt(p) and IsPrimeInt(p)
                and IsPosInt(k) then
            if IsBound(IRREDSOL_DATA.MS_GROUP_INDEX[n]) then
                if IsBound(IRREDSOL_DATA.MS_GROUP_INDEX[n][p]) then
                    if IsBound(IRREDSOL_DATA.MS_GROUP_INDEX[n][p][k]) then
                        return Concatenation([n, p], IRREDSOL_DATA.MS_GROUP_INDEX[n][p][k]);
                    else
                        Error("k is out of range");
                    fi;
                else
                    Error("p is out of range");
                fi;
            else
                Error("n is out of range");
            fi;
        else
            Error("n, p, k must be integers, n > 1, and p must be a prime");
        fi;
    end);

############################################################################
##
#F  IndexMSIdIrreducibleSolubleMatrixGroup(<n>, <q>, <d>, <k>)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(IndexMSIdIrreducibleSolubleMatrixGroup,
    function(n, p, d, k)
        local pos;

        if IsInt(n) and n > 1 and IsPosInt(p) and IsPrimeInt(p)
                and IsPosInt(k) and IsPosInt(d) and n mod d = 0 then
            if IsBound(IRREDSOL_DATA.MS_GROUP_INDEX[n]) then
                if IsBound(IRREDSOL_DATA.MS_GROUP_INDEX[n][p]) then
                    pos := Position (IRREDSOL_DATA.MS_GROUP_INDEX[n][p], [d, k]);
                    if pos <> fail then
                        return [n, p, pos];
                    else
                        Error("inadmissible value for k");
                    fi;
                else
                    Error("p is out of range");
                fi;
            else
                Error("n is out of range");
            fi;
        else
            Error("n, p, d, k must be integers, n > 1, p must be a prime, and d must divide n");
        fi;
    end);

############################################################################
##
#E
##

max_stars_repo_name: fingolfin/irredsol

max_issues_repo_name: fingolfin/irredsol

max_forks_repo_name: fingolfin/irredsol


text:
#############################################################################
##
#W  tree.gi                 automgrp package                   Yevgen Muntyan
#W                                                             Dmytro Savchuk
##
#Y  Copyright (C) 2003 - 2018 Yevgen Muntyan, Dmytro Savchuk
##

#############################################################################
##
#M  TopDegreeOfTree (<obj>)
##
InstallMethod(TopDegreeOfTree, "for [IsActingOnTree]", [IsActingOnTree],
function(obj)
  local si;
  si := SphericalIndex(obj);
  if not IsEmpty(si.start) then
    return si.start[1];
  else
    return si.period[1];
  fi;
end);

#############################################################################
##
#M  DegreeOfLevel (<obj>)
##
InstallMethod(DegreeOfLevel, "for [IsActingOnTree, IsPosInt]",
              [IsActingOnTree, IsPosInt],
function(obj, k)
  return AG_DegreeOfLevelInSphericalIndex(SphericalIndex(obj), k);
end);

#############################################################################
##
#M  IsActingOnRegularTree(<obj>)
##
InstallMethod(IsActingOnRegularTree, "for [IsActingOnTree]",
              [IsActingOnTree],
function(obj)
  local si;
  si := SphericalIndex(obj);
  return IsEmpty(si.start) and Length(si.period) = 1;
end);

#############################################################################
##
#M  DegreeOfTree(<obj>)
##
InstallMethod(DegreeOfTree, "for [IsActingOnTree]",
              [IsActingOnTree],
function(obj)
  if not IsActingOnRegularTree(obj) then Error("The object ",obj," is not acting on a regular tree"); fi;
  return SphericalIndex(obj).period[1];
end);

#############################################################################
##
#M  IsActingOnBinaryTree(<obj>)
##
InstallMethod(IsActingOnBinaryTree, "for [IsActingOnTree]",
              [IsActingOnTree],
function(obj)
  return IsActingOnRegularTree(obj) and DegreeOfTree(obj) = 2;
end);

#############################################################################
##
#M  FixesLevel (<obj>, <k>)
##
##  The most stupid method - just for case.
##
InstallMethod(FixesLevel, "for [IsActingOnTree, IsPosInt]",
              [IsActingOnTree, IsPosInt],
function(obj, k)
  local lev, v;

  if HasIsSphericallyTransitive(obj) then
    if IsSphericallyTransitive(obj) then
      return false;
    fi;
  fi;

  lev := Tuples([1..DegreeOfTree(obj)], k);
  for v in lev do
    if not FixesVertex(obj, v) then
      return false;
    else
      Info(InfoAutomGrp, 3, "IsSphericallyTransitive(obj): false");
      Info(InfoAutomGrp, 3, "  obj fixes vertex", v);
      Info(InfoAutomGrp, 3, "  obj = ", obj);
      SetIsSphericallyTransitive(obj, false);
    fi;
  od;

  return true;
end);

#E

ext: gi

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

ImportAll(paradigms.smp); # for AParSMP
ImportAll(paradigms.vector); # for VRCLR
Import(approx); # for CeilingRat.

Load(spiral.paradigms.dram.hacks); # for TTensor3 etc.
Import(hacks);

#_swrap := (nt) -> let(t:=nt.firstTag(), ScratchWrap(t.size,t.nsgmts,t.linesize));

#_nodeFitsTensorDoesnt := (nt) -> let (
#    k := nt.getTag(1).size,
#    m := nt.params[1].dims()[2],
#    n := nt.params[2],
#
#    2*m <= k and m*n >= k
#);

#_isALStoreTag := (nt) -> nt.isTag(1,ALStore) or nt.isTag(1,ALStoreCx);

_isADramTag := (nt) -> nt.isTag(1,ADram);
_isATileTag := (nt) -> nt.isTag(1,ATile);
_isATileRdTag := (nt) -> nt.isTag(1,ATileRd);
_isATileWrTag := (nt) -> nt.isTag(1,ATileWr);
_isACubeTag := (nt) -> nt.isTag(1,ACube);
_isACubeRdTag := (nt) -> nt.isTag(1,ACubeRd);
_isACubeWrTag := (nt) -> nt.isTag(1,ACubeWr);
_isTL_I := (nt) -> nt.params[1]=1 and nt.params[2]=1 and nt.params[3]>1 and nt.params[4]=1;
_isTL_L := (nt) -> nt.params[1]>1 and nt.params[2]>1 and nt.params[3]=1 and nt.params[4]=1;
_isTL_I_L := (nt) -> nt.params[1]>1 and nt.params[2]>1 and nt.params[3]>1 and nt.params[4]=1;
_isReduced1D := (nt) -> let(p := Length(nt.params[1].params), When(p=3, (Length(nt.params[1].params[1]) = 1), true));

_DetermineTag := function(t)
	local ll,tag,k;
	ll := When(Length(t.params)=3,Length(t.params[1]),1); #if MDDFT,x,1

	if(ll = 2 or ll = 1) then
		k := Int(Log2Int(t.getTag(1).rb)/2);
		k := 2^k;
		if(not (k^2 = t.getTag(1).rb)) then
			PrintLine("WARNING: Tile <-> Row-Buffer cannot be mathced perfectly!");
		fi;
		tag := ATile(k,t.getTag(1).m);
	else if(ll = 3) then
		k := Int(Log2Int(t.getTag(1).rb)/3);
		k := 2^k;
		if(not (k^3 = t.getTag(1).rb)) then
			PrintLine("WARNING: Cube <-> Row-Buffer cannot be mathced perfectly!");
		fi;
		tag := ACube(k,t.getTag(1).m);
	else
		PrintLine("I cannot handle an algorithm with ",ll," stages for now, sorry..");
		Error("");
	fi; fi;

	return tag;
end;
#  NewRulesFor(TCompose, rec(
#  # (AB) -> (A) fence (B)
#  	AB_tile := rec(
#  	)
#  ));

NewRulesFor(MDDFT, rec(
    MDDFT_tSPL_RowCol_break_2D := rec(
        info := "tSPL MDDFT_n -> MDDFT_n/d, MDDFT_d",

        applicable := (self, t) >> Length(t.params[1]) > 1 and not _isADramTag(t) and _isATileTag(t),
        freedoms := t -> [ [1..Length(t.params[1])-1] ],

        child := (t, fr) -> let(
            newdims := SplitAt(t.params[1], fr[1]),
            rot := t.params[2],
            [ TTensor(
                MDDFT(newdims[1], rot),
                MDDFT(newdims[2], rot)
            ).withTags(t.getTags())]
        ),

        apply := (t, C, Nonterms) -> C[1],
        switch := false
    ),

	MDDFT_tSPL_RowCol_break_3D := rec(
        info := "tSPL MDDFT_n -> MDDFT_n/d, MDDFT_d",

        applicable := (self, t) >> Length(t.params[1]) > 1 and not _isADramTag(t) and _isACubeTag(t),
        freedoms := t -> [ [1..Length(t.params[1])-1] ],

        children := (t) -> let(
            newdims := t.params[1],

            [[ TTensor3(
                MDDFT([newdims[1]]),
                MDDFT([newdims[2]]),
                MDDFT([newdims[3]])
            ).withTags(t.getTags())]]
        ),

        apply := (t, C, Nonterms) -> C[1],
        switch := false
    ),

	MDDFT_tSPL_RowCol_push := rec(
        info := "tSPL MDDFT_n -> MDDFT_n/d, MDDFT_d",

        applicable := (self, t) >> Length(t.params[1]) > 1 and _isADramTag(t),
        freedoms := t -> [ [1..Length(t.params[1])-1] ],

        child := (t, fr) -> let(
			tag := _DetermineTag(t),
			newTag := Concat([tag], Drop(t.getTags(), 1)),
            [ t.withoutFirstTag().withTags(newTag) ]
        ),

        apply := (t, C, Nonterms) -> C[1],
        switch := false
    )

));

NewRulesFor(TTwiddle, rec(
    TTwiddle_dram := rec(
        minSize := false,
        maxSize := false,
        applicable := (self, nt) >> ((IsBool(self.minSize) and not self.minSize) or nt.params[1] >= self.minSize) and
            ((IsBool(self.maxSize) and not self.maxSize) or nt.params[1] <= self.maxSize),
        forTransposition := false,
        apply := (t, C, Nonterms) -> TwiddleROM(t.params[1], t.params[2], t.params[3])
    )
));

#######################################################################################
#   tSPL DFT rules
NewRulesFor(DFT, rec(
    #F DFT_CT: 1965
    #F   General Cooley-Tukey Rule
    #F   DFT_n = (DFT_n/d tensor I_d) * diag * (I_n/d tensor F_d) * perm
    #F
    #F Cooley/Tukey:
    #F   An Algorithm for the Machine Calculation of Complex Fourier Series.
    #F   Mathematics of Computation, Vol. 19, 1965, pp. 297--301.
    #F
    DFT_tSPL_CT_tiled := rec(
    info          := "tSPL DFT(mn,k) -> DFT(m, k%m), DFT(n, k%n)",

    maxSize       := false,

    applicable    := (self, nt) >> nt.params[1] > 2
        and (self.maxSize = false or nt.params[1] <= self.maxSize)
		and _isATileTag(nt)
        and not IsPrime(nt.params[1])
		and IsPosInt(Sqrt(nt.params[1]))
        and nt.hasTags(),

    children      := nt -> let(n := Sqrt(nt.params[1]), m := n,

	[[
        TCompose([
            TGrp(TCompose([
                TTensorI(DFT(m, nt.params[2] mod m), n, AVec, AVec),
                TTwiddle(m*n, n, nt.params[2])
            ])),
            TGrp(TTensorI(DFT(n, nt.params[2] mod n), m, APar, AVec))
        ]).withTags(nt.getTags())
    ]]),

    apply := (nt, c, cnt) -> c[1],

    switch := false
    ),

	DFT_tSPL_push_tiled := rec(
        info := "tSPL DFT(mn,k) -> DFT(m, k%m), DFT(n, k%n)",

		maxSize       := false,

		applicable := (self, t) >> t.params[1] > 2 and _isADramTag(t),

        children := (t) -> let(
			tag := _DetermineTag(t),
			newTag := Concat([tag], Drop(t.getTags(), 1)),
            [[ t.withoutFirstTag().withTags(newTag) ]]
        ),

        apply := (t, C, Nonterms) -> C[1],
        switch := false
    )
));

############################## 3D Rules ##############################

#   (A x B) rules
NewRulesFor(TTensor3, rec(
#   (A x B x C) -> (A x B x I)(I x I x C)
    AxBxI_IxIxC := rec(
        info := "(A x B x C) -> (A x B x I)(I x I x C)",
        forTransposition := false,
        applicable := nt -> true,
        #inplace := false,
        children := (self, nt) >> let(#inp := When(self.inplace, Inplace, x->x),
            [[ TCompose([
                TTensorI3(nt.params[1], nt.params[2], I(nt.params[3].dims()[2])),
                TTensorI3(I(nt.params[1].dims()[2]), I(nt.params[2].dims()[2]), nt.params[3])
            ]).withTags(nt.getTags()) ]]),
        apply := (nt, c, cnt) -> c[1],
#D        isApplicable := P -> true,
#D        allChildren := P -> [[TCompose([TTensorI(P[1], P[2].dims()[1], AVec, AVec), TTensorI(P[2], P[1].dims()[2], APar, APar)], P[3])]],
#D        rule := (P, C) -> C[1]
    ),
	#   (A x B x C) -> (A x I x I)(I x B x C)
    AxIxI_IxBxC := rec(
        info := "(A x B x C) -> (A x I x I)(I x B x C)",
        forTransposition := false,
        applicable := nt -> true,
        #inplace := false,
        children := (self, nt) >> let(#inp := When(self.inplace, Inplace, x->x),
            [[ TCompose([
				TTensorI3(nt.params[1], I(nt.params[2].dims()[2]), I(nt.params[3].dims()[2])),
				TTensorI3(I(nt.params[1].dims()[2]), nt.params[2], nt.params[3])
            ]).withTags(nt.getTags()) ]]),
        apply := (nt, c, cnt) -> c[1],
	),

	# (A x B x C) -> (A x I x I)(I x B x I)(I x I x C)
	AxIxI_IxBxI_IxIxC := rec(
        info := "(A x B x C) -> (A x I x I)(I x B x I)(I x I x C)",
        forTransposition := false,
        applicable := nt -> true,
        #inplace := false,
        children := (self, nt) >> let(#inp := When(self.inplace, Inplace, x->x),
            [[ TCompose([
				TTensorI3(nt.params[1], I(nt.params[2].dims()[2]), I(nt.params[3].dims()[2])),
				TTensorI3(I(nt.params[1].dims()[2]), nt.params[2], I(nt.params[3].dims()[2])),
				TTensorI3(I(nt.params[1].dims()[2]), I(nt.params[2].dims()[2]), nt.params[3])
            ]).withTags(nt.getTags()) ]]),
        apply := (nt, c, cnt) -> c[1],
	)
));

NewRulesFor(TTensorI3, rec(
	# (In x In x An) -> In3 (In2/k2 x (Ik x Ik x DFTn)) In3
	IxIxA_cube := rec(
		forTransposition := false,
		applicable := nt -> nt.hasTags()
							and _isACubeTag(nt)
							and IsIIA(nt.params)
							and IsPosInt(nt.params[3].dims()[2] / nt.getTag(1).k / nt.getTag(1).k)
							and (Length(nt.params[3].params[1]) = 1) # reduced to 1D-DFT kernel
							and (nt.dims()[2] > nt.getTag(1).m),

		children := nt -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[3].dims()[2],
			RdTag := Concat([ACubeRd(k,m)], Drop(nt.getTags(), 1)),
			WrTag := Concat([ACubeWr(k,m)], Drop(nt.getTags(), 1)),

			[[
				TL(1,1,n*n*n,1).withTags(WrTag),
				TTensorI3(I(k),I(k),nt.params[3]).withTags(nt.getTags()), #.withTags(Drop(nt.getTags(), 1)) ,
				TL(1,1,n*n*n,1).withTags(RdTag)
			]]
		),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[3].dims()[2],

			MemFence(
				c[1] *
				DTensor(
					I(n*n/k/k),
					(c[2])
				) *
				c[3]
			)
		)
	),

	# (An x In x In) -> Ln3_n (In2/k2 x (Ik x Ik x DFTn)) Ln3_n2
	AxIxI_cube := rec(
		forTransposition := false,
		applicable := nt -> nt.hasTags()
							and _isACubeTag(nt)
							and IsAII(nt.params)
							and IsPosInt(nt.params[1].dims()[2] / nt.getTag(1).k / nt.getTag(1).k)
							and (Length(nt.params[1].params[1]) = 1) # reduced to 1D-DFT kernel
							and (nt.dims()[2] > nt.getTag(1).m),

		children := nt -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[1].dims()[2],
			RdTag := Concat([ACubeRd(k,m)], Drop(nt.getTags(), 1)),
			WrTag := Concat([ACubeWr(k,m)], Drop(nt.getTags(), 1)),

			[[
				TL(n*n*n,n,1,1).withTags(WrTag) ,
				TTensorI3(I(k),I(k),nt.params[1]).withTags(nt.getTags()), #.withTags(Drop(nt.getTags(), 1)) ,
				TL(n*n*n,n*n,1,1).withTags(RdTag)
			]]
		),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[1].dims()[2],

			MemFence(
				c[1] *
				DTensor(
					I(n*n/k/k),
					(c[2])
				) *
				c[3]
			)
		)
	),

	# (In x An x In) -> In x Ln2_n (In2/k2 x (Ik x Ik x DFTn)) In x Ln2_n
	IxAxI_cube := rec(
		forTransposition := false,
		applicable := nt -> nt.hasTags()
							and _isACubeTag(nt)
							and IsIAI(nt.params)
							and IsPosInt(nt.params[2].dims()[2] / nt.getTag(1).k / nt.getTag(1).k)
							and (Length(nt.params[2].params[1]) = 1) # reduced to 1D-DFT kernel
							and (nt.dims()[2] > nt.getTag(1).m),

		children := nt -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[2].dims()[2],
			RdTag := Concat([ACubeRd(k,m)], Drop(nt.getTags(), 1)),
			WrTag := Concat([ACubeWr(k,m)], Drop(nt.getTags(), 1)),

			[[
				TL(n*n,n,n,1).withTags(WrTag) ,
				TTensorI3(I(k),I(k),nt.params[2]).withTags(nt.getTags()), #.withTags(Drop(nt.getTags(), 1)) ,
				TL(n*n,n,n,1).withTags(RdTag)
			]]
		),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[2].dims()[2],

			MemFence(
				c[1] *
				DTensor(
					I(n*n/k/k),
					(c[2])
				) *
				c[3]
			)
		)
	),

	IxIxA_cube_base := rec(
		forTransposition := false,
		applicable := nt -> nt.hasTags()
							and _isACubeTag(nt)
							and IsIIA(nt.params)
							and (Length(nt.params[3].params[1]) = 1) # reduced to 1D-DFT kernel
							and (nt.dims()[2] <= nt.getTag(1).m),

		children := nt -> [[ ]],

		apply := (nt, c, cnt) -> let(
			n := nt.params[3].dims()[2],
			k2 := nt.params[2].dims()[2],
			k1 := nt.params[1].dims()[2],

			CompKern(I(k1),I(k2),DFT(n))
		)

	)

));

NewRulesFor(TL, rec(
# In^2_{read} -> (In/k x L^n_k x Ik)(In/k x L^n_n/k x Ik)
	I_cubeRd := rec(
		applicable := nt -> nt.hasTags()
							and _isACubeRdTag(nt)
							and _isTL_I(nt)
							and IsPosInt(RootInt(nt.params[3],3)) #n3 -> n
							and IsPosInt(RootInt(nt.params[3],3) / nt.getTag(1).k),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := RootInt(nt.params[3],3),
			#mem := MemRdPrm(fTensor(fId(n/k),L(n*n,n*n/k/k),fId(k))) * MemRdPrm(fTensor(fId(n),L(n,k),fId(n))),
			mem := MemRdPrm(fTensor(fId(n/k),L(n*n,n*n/k/k),fId(k))) * MemRdPrm(fTensor(fId(n),L(n,k),fId(n/k),fId(k))),
			loc := DTensor(I(n*n/k/k), LocalRdPrm(fTensor(L(n*k,k*k), fId(k)))),

			loc*mem

		)
	),
# In^2_{write} -> (In/k x L^n_k x Ik)(In/k x L^n_n/k x Ik)
	I_cubeWr := rec(
		applicable := nt -> nt.hasTags()
							and _isACubeWrTag(nt)
							and _isTL_I(nt)
							and IsPosInt(RootInt(nt.params[3],3)) #n3 -> n
							and IsPosInt(RootInt(nt.params[3],3) / nt.getTag(1).k),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := RootInt(nt.params[3],3),
			#mem := MemRdPrm(fTensor(fId(n/k),L(n*n,n*n/k/k),fId(k))) * MemRdPrm(fTensor(fId(n),L(n,k),fId(n))),
			mem := MemRdPrm(fTensor(fId(n/k),L(n*n,n*n/k/k),fId(k))) * MemRdPrm(fTensor(fId(n),L(n,k),fId(n/k),fId(k))),
			loc := DTensor(I(n*n/k/k), LocalRdPrm(fTensor(L(n*k,k*k), fId(k)))),

			TransposedSPL(loc*mem)
		)
	),
# Ln^2_n_{read} -> (In/k x L^nk_k)(L^n^2/k_n/k x Ik)
	Ln3n2_cubeRd := rec(
		applicable := nt -> nt.hasTags()
							and _isACubeRdTag(nt)
							and _isTL_L(nt)
							and IsPosInt(RootInt(nt.params[1],3)) #n3 -> n
							and IsPosInt(RootInt(nt.params[1],3) / nt.getTag(1).k),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := RootInt(nt.params[1],3),
			#mem := 	MemRdPrm(fTensor(L(n*n*n/k,n*n/k/k), fId(k))) * MemRdPrm(fTensor(fId(n),L(n,k),fId(n))),
			mem := 	MemRdPrm(fTensor(L(n*n*n/k,n*n/k/k), fId(k))) * MemRdPrm(fTensor(fId(n),L(n,k),fId(n/k),fId(k))),
			loc := DTensor(I(n*n/k/k), LocalRdPrm(fTensor(L(n*k*k,k*k)))),

			loc*mem
		)
	),
# Ln^2_n_{write} -> (L^n^2/k_n x Ik)(In/k x L^nk_n)
	Ln3n_cubeWr := rec(
		applicable := nt -> nt.hasTags()
							and _isACubeWrTag(nt)
							and _isTL_L(nt)
							and IsPosInt(RootInt(nt.params[1],3)) #n3 -> n
							and IsPosInt(RootInt(nt.params[1],3) / nt.getTag(1).k),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := RootInt(nt.params[1],3),
			#mem := 	MemRdPrm(fTensor(L(n*n*n/k,n*n/k/k), fId(k))) * MemRdPrm(fTensor(fId(n),L(n,k),fId(n))),
			mem := 	MemRdPrm(fTensor(L(n*n*n/k,n*n/k/k), fId(k))) * MemRdPrm(fTensor(fId(n),L(n,k),fId(n/k),fId(k))),
			loc := DTensor(I(n*n/k/k), LocalRdPrm(fTensor(L(n*k*k,k*k)))),

			TransposedSPL(loc*mem)
		)

	),
# Ln^2_n_{read} -> (In/k x L^nk_k)(L^n^2/k_n/k x Ik)
	InLn2n_cubeRd := rec(
		applicable := nt -> nt.hasTags()
							and _isACubeRdTag(nt)
							and _isTL_I_L(nt)
							and IsPosInt(RootInt(nt.params[1],2)) #n2 -> n
							and IsPosInt(RootInt(nt.params[1],2) / nt.getTag(1).k),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := RootInt(nt.params[1],2),
			#mem := MemRdPrm(fTensor(fId(n/k), L(n*n,n/k), fId(k))) * MemRdPrm(fTensor(fId(n/k), L(n,n/k), fId(n*k))),
			mem := MemRdPrm(fTensor(fId(n/k), L(n*n,n/k), fId(k))) * MemRdPrm(fTensor(fId(n/k), L(n,n/k), fId(n),fId(k))),
			loc := DTensor(I(n*n/k/k), LocalRdPrm(L(n*k*k,k*k))) * DTensor(I(n*n/k/k), LocalRdPrm(fTensor(fId(n/k),L(k*k,k),fId(k)))),

			loc*mem
		)
	),

	InLn2n_cubeWr := rec(
		applicable := nt -> nt.hasTags()
							and _isACubeWrTag(nt)
							and _isTL_I_L(nt)
							and IsPosInt(RootInt(nt.params[1],2)) #n2 -> n
							and IsPosInt(RootInt(nt.params[1],2) / nt.getTag(1).k),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := RootInt(nt.params[1],2),
			#mem := MemRdPrm(fTensor(fId(n/k), L(n*n,n/k), fId(k))) * MemRdPrm(fTensor(fId(n/k), L(n,n/k), fId(n*k))),
			mem := MemRdPrm(fTensor(fId(n/k), L(n*n,n/k), fId(k))) * MemRdPrm(fTensor(fId(n/k), L(n,n/k), fId(n),fId(k))),
			loc := DTensor(I(n*n/k/k), LocalRdPrm(L(n*k*k,k*k))) * DTensor(I(n*n/k/k), LocalRdPrm(fTensor(fId(n/k),L(k*k,k),fId(k)))),

			TransposedSPL(loc*mem)
		)
	)

));

############################## 2D Rules ##############################

NewRulesFor(TTensorI, rec(

# ================================================================
# nt.params; => [DFT_matix, Size_of_I, <APar,Avec>, <APar,Avec>]
# So,
# nt.params[1].dims()[2] -> Size of DFT matrix
# nt.params[2] -> Size of I matrix
# IsParPar actually checks nt.params[3]=APar and nt.params[4]=APar
#
# nt.getTag(1) -> Tag object. So access its methods via "."
#
# ================================================================
#  tSPL I x L x I
#  TL(m, n, k, j)
#  I(k), L(m, n), I(j)
# ================================================================

# (In x An) -> In^2 (In/k x Ik x An) In^2
# (In x An)_{Tile(k,m)} -> In^2_{TileWr(k,m)} (In/k x Ik x An) In^2_{TileRd(k,m)} for k|n and m>kn
	IxA_tile := rec(
		forTransposition := false,
		applicable := nt -> nt.hasTags()
							and _isATileTag(nt)
							and IsParPar(nt.params)
							and IsPosInt(nt.params[1].dims()[2] / nt.getTag(1).k)
							and _isReduced1D(nt) # reduced to 1D-DFT kernel
							#and (Length(nt.params[1].params[1]) = 1) # reduced to 1D-DFT kernel
							and (nt.params[1].dims()[2]*nt.params[2] > nt.getTag(1).m),

		children := nt -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[1].dims()[2],
			RdTag := Concat([ATileRd(k,m)], Drop(nt.getTags(), 1)),
			WrTag := Concat([ATileWr(k,m)], Drop(nt.getTags(), 1)),

			[[
				TL(1,1,n*n,1).withTags(WrTag) ,
				TTensorI(nt.params[1],k,APar,APar).withTags(nt.getTags()), #.withTags(Drop(nt.getTags(), 1)) ,
				TL(1,1,n*n,1).withTags(RdTag)
			]]
		),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[1].dims()[2],

			MemFence(
				c[1] *
				DTensor(
					I(n/k),
					(c[2])
				) *
				c[3]
			)
		)
	),

# (An x In) -> L^{n^2}_n (In/k x Ik x An) L^{n^2}_n
	AxI_tile := rec(
		forTransposition := false,
		applicable := nt -> nt.hasTags()
							and _isATileTag(nt)
							and IsVecVec(nt.params)
							and IsPosInt(nt.params[1].dims()[2] / nt.getTag(1).k)
							and _isReduced1D(nt) # reduced to 1D-DFT kernel
							#and (Length(nt.params[1].params[1]) = 1) # reduced to 1D-DFT kernel
							and (nt.params[1].dims()[2]*nt.params[2] > nt.getTag(1).m),

		children := nt -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[1].dims()[2],
			RdTag := Concat([ATileRd(k,m)], Drop(nt.getTags(), 1)),
			WrTag := Concat([ATileWr(k,m)], Drop(nt.getTags(), 1)),

			[[
				TL(n*n,n,1,1).withTags(WrTag) ,
				TTensorI(nt.params[1],k,APar,APar).withTags(nt.getTags()), #.withTags(Drop(nt.getTags(), 1)) ,
				TL(n*n,n,1,1).withTags(RdTag)
			]]
		),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[1].dims()[2],

			MemFence(
				c[1] *
				DTensor(
					I(n/k),
					(c[2])
				) *
				c[3]
			)
		)
	),
# L^{n^2}_n (In x An) -> L^{n^2}_n (In/k x Ik x An) In^2
	L_IxA_tile := rec(
		forTransposition := false,
		applicable := nt -> nt.hasTags()
							and _isATileTag(nt)
							and IsVecPar(nt.params)
							and IsPosInt(nt.params[1].dims()[2] / nt.getTag(1).k)
							and _isReduced1D(nt) # reduced to 1D-DFT kernel
							#and (Length(nt.params[1].params[1]) = 1) # reduced to 1D-DFT kernel
							and (nt.params[1].dims()[2]*nt.params[2] > nt.getTag(1).m),

		children := nt -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[1].dims()[2],
			RdTag := Concat([ATileRd(k,m)], Drop(nt.getTags(), 1)),
			WrTag := Concat([ATileWr(k,m)], Drop(nt.getTags(), 1)),

			[[
				TL(n*n,n,1,1).withTags(WrTag) ,
				TTensorI(nt.params[1],k,APar,APar).withTags(nt.getTags()), #.withTags(Drop(nt.getTags(), 1)) ,
				TL(1,1,n*n,1).withTags(RdTag)
			]]
		),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[1].dims()[2],

			MemFence(
				c[1] *
				DTensor(
					I(n/k),
					(c[2])
				) *
				c[3]
			)
		)
	),
# (In x An) L^{n^2}_n -> In^2 (In/k x Ik x An) L^{n^2}_n
	IxA_L_tile := rec(
		forTransposition := false,
		applicable := nt -> nt.hasTags()
							and _isATileTag(nt)
							and IsParVec(nt.params)
							and IsPosInt(nt.params[1].dims()[2] / nt.getTag(1).k)
							and _isReduced1D(nt) # reduced to 1D-DFT kernel
							#and (Length(nt.params[1].params[1]) = 1) # reduced to 1D-DFT kernel
							and (nt.params[1].dims()[2]*nt.params[2] > nt.getTag(1).m),

		children := nt -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[1].dims()[2],
			RdTag := Concat([ATileRd(k,m)], Drop(nt.getTags(), 1)),
			WrTag := Concat([ATileWr(k,m)], Drop(nt.getTags(), 1)),

			[[
				TL(1,1,n*n,1).withTags(WrTag) ,
				TTensorI(nt.params[1],k,APar,APar).withTags(nt.getTags()), #.withTags(Drop(nt.getTags(), 1)) ,
				TL(n*n,n,1,1).withTags(RdTag)
			]]
		),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := nt.params[1].dims()[2],

			MemFence(
				c[1] *
				DTensor(
					I(n/k),
					(c[2])
				) *
				c[3]
			)
		)
	)
));

# Termination Rules, i.e. computation fits in local memory
NewRulesFor(TTensorI, rec(
	IxA_tile_base := rec(
		applicable := nt -> nt.hasTags()
							and _isATileTag(nt)
							and IsParPar(nt.params)
							and (nt.params[1].dims()[2] * nt.params[2] <= nt.getTag(1).m),

		children := nt -> [[ ]],

		apply := (nt, c ,cnt) -> let(
			k := nt.params[2],
			n := nt.params[1].dims()[2],

			CompKern(I(k), DFT(n))
		)
	),

	AxI_tile_base := rec(
		applicable := nt -> nt.hasTags()
							and _isATileTag(nt)
							and IsVecVec(nt.params)
							and (nt.params[1].dims()[2] * nt.params[2] <= nt.getTag(1).m),

		children := nt -> [[ ]],

		apply := (nt, c ,cnt) -> let(
			k := nt.params[2],
			n := nt.params[1].dims()[2],

			# NOTE: Check the L index!
			LocalWrPrm(fTensor(L(k*n,n))) * CompKern(I(k), DFT(n)) * LocalRdPrm(fTensor(L(n*k,k)))
			# TPrm and TL for splhdl
			#LocalWrPrm(Tensor(L(k*n,n))) * CompKern(I(k), DFT(n)) * LocalRdPrm(Tensor(L(n*k,k)))
		)
	),

	L_IxA_tile_base := rec(
		applicable := nt -> nt.hasTags()
							and _isATileTag(nt)
							and IsVecPar(nt.params)
							and (nt.params[1].dims()[2] * nt.params[2] <= nt.getTag(1).m),

		children := nt -> [[ ]],

		apply := (nt, c ,cnt) -> let(
			k := nt.params[2],
			n := nt.params[1].dims()[2],

			# NOTE: Check the L index!
			# TTensorI(DFTn,k,AVec,APar) -> L(n*k,???) * (Ik x DFTn)
			#
			#
			#
			LocalWrPrm(fTensor(L(n*k,k))) * CompKern(I(k),DFT(n))
			# TPrm and TL for splhdl
			#LocalWrPrm(Tensor(L(n*k,k))) * CompKern(I(k),DFT(n))
		)
	),

	IxA_L_tile_base := rec(
		applicable := nt -> nt.hasTags()
							and _isATileTag(nt)
							and IsParVec(nt.params)
							and (nt.params[1].dims()[2] * nt.params[2] <= nt.getTag(1).m),

		children := nt -> [[ ]],

		apply := (nt, c ,cnt) -> let(
			k := nt.params[2],
			n := nt.params[1].dims()[2],

			# NOTE: Check the L index!
			CompKern(I(k),DFT(n)) * LocalRdPrm(fTensor(L(k*n,n)))
			# TPrm and TL for splhdl
			#CompKern(I(k),DFT(n)) * LocalRdPrm(Tensor(L(k*n,n)))
		)
	)
));

NewRulesFor(TL, rec(
# In^2_{read} -> (In/k x L^n_k x Ik)(In/k x L^n_n/k x Ik)
	I_tileRd := rec(
		applicable := nt -> nt.hasTags()
							and _isATileRdTag(nt)
							and _isTL_I(nt)
							and IsPosInt(Sqrt(nt.params[3]))
							and IsPosInt(Sqrt(nt.params[3]) / nt.getTag(1).k),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := Sqrt(nt.params[3]),

			DTensor(I(n/k), LocalRdPrm(fTensor(L(n,k),fId(k)))) * MemRdPrm(fTensor(fId(n/k), L(n,n/k), fId(k)))
			# TPrm and TL for splhdl
			#DTensor(I(n/k), LocalRdPrm(Tensor(L(n,k),I(k)))) * MemRdPrm(fTensor(fId(n/k), L(n,n/k), fId(k)))

		)
	),
# In^2_{write} -> (In/k x L^n_k x Ik)(In/k x L^n_n/k x Ik)
	I_tileWr := rec(
		applicable := nt -> nt.hasTags()
							and _isATileWrTag(nt)
							and _isTL_I(nt)
							and IsPosInt(Sqrt(nt.params[3]))
							and IsPosInt(Sqrt(nt.params[3]) / nt.getTag(1).k),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := Sqrt(nt.params[3]),

			MemWrPrm(fTensor(fId(n/k), L(n,k), fId(k))) * DTensor(I(n/k), LocalWrPrm(fTensor(L(n,n/k),fId(k))))
			# TPrm and TL for splhdl
			#MemWrPrm(fTensor(fId(n/k), L(n,k), fId(k))) * DTensor(I(n/k), LocalWrPrm(Tensor(L(n,n/k),I(k))))
		)
	),
# Ln^2_n_{read} -> (In/k x L^nk_k)(L^n^2/k_n/k x Ik)
	L_tileRd := rec(
		applicable := nt -> nt.hasTags()
							and _isATileRdTag(nt)
							and _isTL_L(nt)
							and IsPosInt(Sqrt(nt.params[1]))
							and IsPosInt(Sqrt(nt.params[1]) / nt.getTag(1).k),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := Sqrt(nt.params[1]),

			DTensor(I(n/k), LocalRdPrm(fTensor(L(n*k,k)))) * MemRdPrm(fTensor(L(n*n/k,n/k), fId(k)))
			# TPrm and TL for splhdl
			#DTensor(I(n/k), LocalRdPrm(Tensor(L(n*k,k)))) * MemRdPrm(fTensor(L(n*n/k,n/k), fId(k)))

		)
	),
# Ln^2_n_{write} -> (L^n^2/k_n x Ik)(In/k x L^nk_n)
	L_tileWr := rec(
		applicable := nt -> nt.hasTags()
							and _isATileWrTag(nt)
							and _isTL_L(nt)
							and IsPosInt(Sqrt(nt.params[1]))
							and IsPosInt(Sqrt(nt.params[1]) / nt.getTag(1).k),

		apply := (nt, c, cnt) -> let(
			k := nt.getTag(1).k,
			m := nt.getTag(1).m,
			n := Sqrt(nt.params[1]),

			MemWrPrm(fTensor(L(n*n/k,n), fId(k))) * DTensor(I(n/k), LocalWrPrm(fTensor(L(n*k,n))))
			# TPrm and TL for splhdl
			#MemWrPrm(fTensor(L(n*n/k,n), fId(k))) * DTensor(I(n/k), LocalWrPrm(Tensor(L(n*k,n))))
		)

	)
));

#==============================================================
# First version of TL rules
# Generates children for permutations that is to be further optimized
#==============================================================
#	NewRulesFor(TL, rec(
#	# In^2_{read} -> (In/k x L^n_k x Ik)(In/k x L^n_n/k x Ik)
#		I_tileRd := rec(
#			applicable := nt -> nt.hasTags()
#								and _isATileRdTag(nt)
#								and _isTL_I(nt)
#								and IsPosInt(Sqrt(nt.params[3]))
#								and IsPosInt(Sqrt(nt.params[3]) / nt.getTag(1).k),
#
#			children := nt -> let(
#				k := nt.getTag(1).k,
#				m := nt.getTag(1).m,
#				n := Sqrt(nt.params[3]),
#
#				[[
#					TL(n, k, n/k, k).withTags(Drop(nt.getTags(), 1)) ,
#					TL(n, n/k, n/k, k).withTags(Drop(nt.getTags(), 1))
#				]]
#			),
#
#			apply := (nt, c, cnt) -> c[1] * c[2]
#		),
#	# In^2_{write} -> (In/k x L^n_k x Ik)(In/k x L^n_n/k x Ik)
#		I_tileWr := rec(
#			applicable := nt -> nt.hasTags()
#								and _isATileWrTag(nt)
#								and _isTL_I(nt)
#								and IsPosInt(Sqrt(nt.params[3]))
#								and IsPosInt(Sqrt(nt.params[3]) / nt.getTag(1).k),
#
#			children := nt -> let(
#				k := nt.getTag(1).k,
#				m := nt.getTag(1).m,
#				n := Sqrt(nt.params[3]),
#
#				[[
#					TL(n, k, n/k, k).withTags(Drop(nt.getTags(), 1)) ,
#					TL(n, n/k, n/k, k).withTags(Drop(nt.getTags(), 1))
#				]]
#			),
#
#			apply := (nt, c, cnt) -> c[1] * c[2]
#		),
#	# Ln^2_n_{read} -> (In/k x L^nk_k)(L^n^2/k_n/k x Ik)
#		L_tileRd := rec(
#			applicable := nt -> nt.hasTags()
#								and _isATileRdTag(nt)
#								and _isTL_L(nt)
#								and IsPosInt(Sqrt(nt.params[1]))
#								and IsPosInt(Sqrt(nt.params[1]) / nt.getTag(1).k),
#
#			children := nt -> let(
#				k := nt.getTag(1).k,
#				m := nt.getTag(1).m,
#				n := Sqrt(nt.params[1]),
#
#				[[
#					TL(n*k, k, n/k, 1).withTags(Drop(nt.getTags(), 1)) ,
#					TL(n*n/k, n/k, 1, k).withTags(Drop(nt.getTags(), 1))
#				]]
#			),
#
#			apply := (nt, c, cnt) -> c[1] * c[2]
#		),
#	# Ln^2_n_{write} -> (L^n^2/k_n x Ik)(In/k x L^nk_n)
#		L_tileWr := rec(
#			applicable := nt -> nt.hasTags()
#								and _isATileWrTag(nt)
#								and _isTL_L(nt)
#								and IsPosInt(Sqrt(nt.params[1]))
#								and IsPosInt(Sqrt(nt.params[1]) / nt.getTag(1).k),
#
#			children := nt -> let(
#				k := nt.getTag(1).k,
#				m := nt.getTag(1).m,
#				n := Sqrt(nt.params[1]),
#
#				[[
#					TL(n*n/k, n, 1, k).withTags(Drop(nt.getTags(), 1)),
#					TL(n*k, n, n/k, 1).withTags(Drop(nt.getTags(), 1))
#				]]
#			),
#
#			apply := (nt, c, cnt) -> c[1] * c[2]
#
#		)
#	));
#