#############################################################################
##
#A  matgrp.g                    GAP library                  Martin Schoenert
##
#A  @(#)$Id: matgrp.g,v 3.20.1.6 1996/09/13 11:56:50 sam Exp $
##
#Y  Copyright 1990-1992,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##
##  This file contains  those  functions that mainly deal with matrix groups.
##
##  1995/09/27 sam:
##  For historical reasons, also MeatAxe-like functions can be found in this
##  file.  Now they are, however, independent from 'MatGroupOps'.
##  (They lie in 'MTXOps' now.)
##  The functions use not matrix groups as arguments but module descriptions,
##  always denoted as 'module_descr'.
##  Such a description is a Apply of length 3, the first entry being a Apply of
##  matrices, the second being a field, and the third being the dimension of
##  the matrices.  The interpretation is that of the natural module of the
##  algebra generated by the given matrices over the given field.
##  (The dimension is really needed only in the case of an empty Apply of
##  matrices, which means the algebra generated by the identity matrix.)
##
#H  $Log: matgrp.g,v $
#H  Revision 3.20.1.6  1996/09/13 11:56:50  sam
#H  'MTXOps.EquivalenceTest' needs only the 'vectors' component of the
#H  record returned by 'SpinUpStandard'
#H
#H  Revision 3.20.1.5  1995/11/22  12:12:47  sam
#H  made MeatAxe functions independent of MatGroupOps
#H
#H  Revision 3.20.1.4  1995/03/15  16:42:05  sam
#H  fixed 'InvariantSubspace'
#H
#H  Revision 3.20.1.3  1994/10/10  07:57:19  fceller
#H  added 'MatGroupOps.DerivedSubgroup', 'MatGroupOps.CommutatorSubgroup',
#H  and 'MatGroupOps.CompositionSeries'
#H
#H  Revision 3.20.1.2  1994/08/10  07:23:14  fceller
#H  changed 'MatGroupOps.Subgroup' to avoid changing a parent
#H
#H  Revision 3.20.1.1  1994/08/02  09:30:20  beick
#H  changed 'Peakwords'
#H
#H  Revision 3.20  1994/05/19  13:08:22  sam
#H  added entry 'InvariantForm', fixed minor bug
#H
#H  Revision 3.19  1994/04/21  12:21:34  sam
#H  added 'KroneckerProduct' for matrix groups,
#H  renamed 'Fingerprint' to 'FingerprintEx' (preliminarily)
#H
#H  Revision 3.18  1994/03/04  13:18:23  ahulpke
#H  fixed MatGroupOps.GAP_Intersection for formally different parents
#H
#H  Revision 3.17  1994/03/04  12:41:08  sam
#H  fixed 'CompositionFactors' for case of 0 generators,
#H  bad hack in 'MakePermGroupP'
#H
#H  Revision 3.16  1994/02/16  14:38:40  fceller
#H  added Bettina's MatGroup stuff
#H
#H  Revision 3.15  1994/01/20  12:36:01  sam
#H  moved some dispatchers to 'dispatch.g',
#H  changed 'P.operationImage' to 'P.operation.image'
#H
#H  Revision 3.14  1993/07/16  07:02:36  sam
#H  bad hack in 'Transposed'
#H
#H  Revision 3.13  1993/02/09  14:27:19  martin
#H  made undefined globals local
#H
#H  Revision 3.12  1992/12/16  19:47:27  martin
#H  replaced quoted record names with escaped ones
#H
#H  Revision 3.11  1992/12/03  10:03:28  fceller
#H  changed 'MatGroupOps.PermGroup' to return a bijection
#H
#H  Revision 3.10  1992/05/08  16:21:25  martin
#H  added 'MatGroupOps.RightCoset'
#H
#H  Revision 3.9  1992/05/04  19:04:28  martin
#H  fixed 'MatGroupOps.GAP_Intersection' to assign '<X>.permDomain'
#H
#H  Revision 3.8  1992/04/04  15:27:07  martin
#H  added many more special functions for matrix groups
#H
#H  Revision 3.7  1992/04/03  16:45:12  martin
#H  fixed 'MatricesOps.Group' to check the arguments
#H
#H  Revision 3.6  1992/02/29  13:25:11  jmnich
#H  general library review, some bug fixes
#H
#H  Revision 3.5  1992/02/14  09:46:19  jmnich
#H  changed call of 'Order'
#H
#H  Revision 3.4  1992/01/29  09:09:38  martin
#H  changed 'Order' to take two arguments, group and element
#H
#H  Revision 3.3  1992/01/09  13:25:48  jmnich
#H  added the meataxe functions
#H
#H  Revision 3.2  1992/01/03  15:44:57  martin
#H  changed 'Matrix' to 'Mat'
#H
#H  Revision 3.1  1991/12/06  16:45:52  martin
#H  changed 'MatricesOps.Group' to default to 'GroupElementsOps.Group'
#H
#H  Revision 3.0  1991/11/08  15:09:30  martin
#H  initial revision under RCS
##


#############################################################################
##
#F  IsMatGroup(<obj>) . . . . . . . . . . test if an object is a matrix group
##
IsMatGroup := function ( obj )
    return IsRec( obj )
       and IsBound( obj.isMatGroup )  and obj.isMatGroup;
end;


#############################################################################
##
#F  MatricesOps.Group(<gens>,<id>)  . . . . . . . . . . create a matrix group
##
MatricesOps.Group := function ( Matrices, gens, id )
    local   G,
            d,
            g;

    # check that the generators are all of the same size and invertable
    d := Length(id);
    for g  in gens  do
        if Length(g) <> d  or Length(g[1]) <> d  or RankMat(g) <> d  then
            Error("<gens> must be a Apply of invertable square matrices");
        fi;
    od;

    # make the group record
    G := GroupElementsOps.Group( Matrices, gens, id );

    # Add_ the matrix group tag
    G.isMatGroup     := true;

    # Add_ the known information
    G.dimension         := Length(id);
    G.field             := Field( Flat( Concatenation( gens, [ id ] ) ) );

    # Add_ the operations record
    G.operations        := MatGroupOps;

    # return the group record
    return G;
end;


#############################################################################
##
#V  MatGroupOps . . . . . . . . .  operation record for matrix group category
##
##  'MatGroupOps' is the operation  record  for  matrix  groups.  It contains
##  the domain  functions,  e.g., 'Size'  and 'GAP_Intersection',   and the group
##  functions, e.g., 'Centralizer' and 'SylowSubgroup'.
##
##  'MatGroupOps' is initially a copy of 'GroupOps', and  thus  inherits  the
##  default group  functions.    Currently  we overlay    very few of   those
##  functions.    We should, however,   handle  matrix groups  over small and
##  medium sized finite vector spaces by  treating them as permutation groups
##  over those vector spaces and using the permutation group functions.
##
MatGroupOps := OperationsRecord( "MatGroupOps", GroupOps );


#############################################################################
##
#F  MatGroupOps.Subgroup(<G>,<gens>)  . . . make a subgroup of a matrix group
##
MatGroupOps.Subgroup := function ( G, gens )
    local   S;          # subgroup, result

    # let the default function do the main work
    S := GroupOps.Subgroup( G, gens );

    # Add_ mat group tag and mat group operations record
    if IsBound(S.parent)  then
        S.isMatGroup := true;
        S.operations := MatGroupOps;
        S.field      := G.field;
    fi;

    # return the subgroup
    return S;
end;


#############################################################################
##
#F  MatGroupOps.IsFinite(<G>) . . . . . . .  test if a matrix group is finite
##
MatGroupOps.IsFinite := function ( G )
    if IsFinite( G.field )  then
        return true;
    else
        return GroupOps.IsFinite( G );
    fi;
end;

#############################################################################
##
#F  MatGroupOps.KroneckerProduct( <D>,<G1>,<G2> ) . . . . . Kronecker product
#F                                                           of matrix groups
##
MatGroupOps.KroneckerProduct := function( D, G1, G2 )
    if not Length( G1.generators ) = Length( G2.generators ) then
      Error( "groups must have same number of generators" );
    fi;
    return Group( Apply( [ 1 .. Length( G1.generators ) ],
              x -> KroneckerProduct( G1.generators[x], G2.generators[x] ) ),
                  KroneckerProduct( G1.identity, G2.identity ) );
    end;

#############################################################################
##
#F  MatGroupOps.MakePermGroupP(<G>) . . . make a isomorphic permutation group
##
##  The difference between this function and the usual  'PermGroup'  is  that
##  the permutation group constructed for <G>  will  be  a  subgroup  of  the
##  corresponding permutation group of <G>\'s parent.
##
#T  This is impossible AT LEAST if the parent is infinite.
#T  Bad hack for dealing with the crystallographic groups library:
#T  In the case of an infinite parent just consider <G> itself.
##
MatGroupOps.MakePermGroupP := function ( G )
    local   P;

    # compute the isomorphic permutation group for the parent group of <G>
    # except if the parent knows that it is infinite.
    P := Parent( G );

    if IsBound( P.isFinite ) and not P.isFinite then

      if not IsBound( G.permGroupP )  then
          G.permDomain := GAP_Union( Orbits( G, G.identity ) );
          G.permGroupP := Operation( G, G.permDomain );
      fi;

    else

      if not IsBound( P.permDomain )  then
          P.permDomain := GAP_Union( Orbits( P, P.identity ) );
          P.permGroupP := Operation( P, P.permDomain );
      fi;

      # compute the isomorphic permutation group for <G>
      if not IsBound( G.permGroupP )  then
          G.permDomain := P.permDomain;
          G.permGroupP := Subgroup( P.permGroupP,
                                Operation( G, P.permDomain ).generators );
      fi;

    fi;
end;


#############################################################################
##
#F  MatGroupOps.Size(<G>) . . . . . . . . . . . . . .  size of a matrix group
##
MatGroupOps.Size := function ( G )

    # compute the isomorphic permutation group for <G> and its parent
    G.operations.MakePermGroupP( G );

    # return the size of the permutation group
    return GAP_Size( G.permGroupP );
end;


#############################################################################
##
#F  MatGroupOps.\in( <obj>, <G> )  . . . . membership test for matrix groups
##
MatGroupOps.\in := function ( obj, G )
    local   l,          # <obj> as a permutation represented by a Apply
            v,          # vector from the operation domain of <P>
            p;          # position of '<v> \^\ <obj>' in the operation domain

    # first a quick test
    if     not IsMat( obj )
        or Length(obj) <> Length(G.identity)
        or Length(obj[1]) <> Length(G.identity[1])
        or RankMat(obj) <> Length(G.identity)
        or not IsSubset( G.field, Field( Flat(obj) ) )
    then
        return false;
    fi;

    # compute the isomorphic permutation group for <G> and its parent
    G.operations.MakePermGroupP( G );

    # try to transform <obj> to a permutation
    l := [];
    for v  in G.permDomain  do
        p := Position( G.permDomain, v ^ obj );
        if p = false  then
            return false;
        fi;
        Add_( l, p );
    od;

    # test if the permutation is in the permutation group
    return PermList(l) in G.permGroupP;
end;


#############################################################################
##
#F  MatGroupOps.GAP_Intersection(<G>,<H>) . . . . . intersection of matrix groups
##
MatGroupOps.GAP_Intersection := function ( G, H )
    local   I,          # intersection of <G> and <H>, result
            P,          # permutation representation of <I>
	    U,          # GAP_Intersection a subgroup
	    i;          # loop variable

    # handle the intersection of two matrix groups with the same parent
    if IsMatGroup(G)  and IsMatGroup(H) then
      if Parent(G) = Parent(H)  then

        # compute the isomorphic permutation groups for <G> and <H>
        G.operations.MakePermGroupP( G );
        H.operations.MakePermGroupP( H );

        # intersect the permutation groups and translate back
        P := GAP_Intersection( G.permGroupP, H.permGroupP );
        I := Subgroup( Parent( G ), Apply( P.generators, gen ->
                    Apply( G.identity,
                        v -> G.permDomain[Position(G.permDomain,v)^gen] ) ));
        if not IsBound( I.permGroupP )  then
            I.permDomain := G.permDomain;
            I.permGroupP := P;
        fi;

      else
	I:=GroupOps.GAP_Intersection(G,H);
	# both groups are subgroups of the same GL(n,K)
	if G.identity=H.identity then
	  U:=TrivialSubgroup(G);
	  for i in I do
	    U:=Closure(U,i);
	  od;
	  I:=Group(U.generators,G.identity);
	  I.elements:=Elements(U);
	fi;
      fi;

    # delegate other cases
    else
        I := GroupOps.GAP_Intersection( G, H );
    fi;

    # return the intersection
    return I;
end;


#############################################################################
##
#F  MatGroupOps.Random(<G>) . . . . . . . .  random element in a matrix group
##
MatGroupOps.Random := function ( G )
    local   rnd;        # random element of '<G>.permGroupP'

    # compute the isomorphic permutation group for <G> and its parent
    G.operations.MakePermGroupP( G );

    # take a random permutation and translate it back
    rnd := Random( G.permGroupP );
    return Apply( G.identity,
                 v -> G.permDomain[Position(G.permDomain,v)^rnd] );
end;


#############################################################################
##
#F  MatGroupOps.Centralizer(<G>,<U>)  . . . . . centralizer in a matrix group
##
MatGroupOps.Centralizer := function ( G, U )
    local    C,         # centralizer of <U> in <G>, result
             P;         # permutation group isomorphic to <C> or <U>

    # compute the isomorphic permutation group for <G> and its parent
    G.operations.MakePermGroupP( G );

    # compute the isomorphic permutation or permutation group for <U>
    if IsMat(U)  then
        P := Permutation( U, G.permDomain );
    else
        U.operations.MakePermGroupP( U );
        P := U.permGroupP;
    fi;

    # compute the centralizer in the permutation group and translate back
    P := Centralizer( G.permGroupP, P );
    C := Subgroup( Parent( G ),
            Apply( P.generators, gen ->
                Apply( G.identity,
                    v -> G.permDomain[Position(G.permDomain,v)^gen] ) ) );
    if not IsBound( C.permGroupP )  then
        C.permDomain := G.permDomain;
        C.permGroupP := P;
    fi;

    # return the centralizer
    return C;
end;


#############################################################################
##
#F  MatGroupOps.Normalizer(<G>,<U>) . . . . . .  normalizer in a matrix group
##
MatGroupOps.Normalizer := function ( G, U )
    local    N,         # normalizer of <U> in <G>, result
             P;         # permutation group isomorphic to <N> or <U>

    # compute the isomorphic permutation group for <G> and its parent
    G.operations.MakePermGroupP( G );

    # compute the isomorphic permutation or permutation group for <U>
    if IsMat(U)  then
        P := Permutation( U, G.permDomain );
    else
        U.operations.MakePermGroupP( U );
        P := U.permGroupP;
    fi;

    # compute the normalizer in the permutation group and translate back
    P := Normalizer( G.permGroupP, P );
    N := Subgroup( Parent( G ),
            Apply( P.generators, gen ->
                Apply( G.identity,
                    v -> G.permDomain[Position(G.permDomain,v)^gen] ) ) );
    if not IsBound( N.permGroupP )  then
        N.permDomain := G.permDomain;
        N.permGroupP := P;
    fi;

    # return the normalizer
    return N;
end;


#############################################################################
##
#F  MatGroupOps.SylowSubgroup(<G>,<p>)  . . . Sylowsubgroup of a matrix group
##
MatGroupOps.SylowSubgroup := function ( G, p )
    local   S,          # <p>-Sylow subgroup of <G>, result
            P;          # permutation group isomorphic to <S>

    # compute the isomorphic permutation group for <G> and its parent
    G.operations.MakePermGroupP( G );

    # compute the Sylow subgroup in the permutation group and translate back
    P := SylowSubgroup( G.permGroupP, p );
    S := Subgroup( Parent( G ),
            Apply( P.generators, gen ->
                Apply( G.identity,
                    v -> G.permDomain[Position(G.permDomain,v)^gen] ) ) );
    if not IsBound( S.permGroupP )  then
        S.permDomain := G.permDomain;
        S.permGroupP := P;
    fi;

    # return the Sylow subgroup
    return S;
end;


#############################################################################
##
#F  MatGroupOps.DerivedSubgroup( <G> )  .  derived subgroup of a matrix group
##
MatGroupOps.DerivedSubgroup := function ( G )
    local   P,  D;

    # compute the isomorphic permutation group for <G> and its parent
    G.operations.MakePermGroupP( G );

    # compute the derived subgroup in the permutation group
    P := DerivedSubgroup( G.permGroupP );
    D := MatGroupOps.Subgroup( Parent( G ),
            Apply( P.generators, gen ->
                Apply( G.identity,
                    v -> G.permDomain[Position(G.permDomain,v)^gen] ) ) );

    if not IsBound( D.permGroupP )  then
        D.permDomain := G.permDomain;
        D.permGroupP := P;
    fi;

    # return the derived subgroup
    return D;
end;


#############################################################################
##
#F  MatGroupOps.CommutatorSubgroup( <G>, <H> )  . . . . . commutator subgroup
##
MatGroupOps.CommutatorSubgroup := function ( G, H )
    local   P,  C;

    # compute the isomorphic permutation group for <G> and <H>
    Parent( G, H );
    G.operations.MakePermGroupP( G );
    H.operations.MakePermGroupP( H );

    # compute the derived subgroup in the permutation group
    P := CommutatorSubgroup( G.permGroupP, H.permGroupP );
    C := MatGroupOps.Subgroup( Parent( G ),
            Apply( P.generators, gen ->
                Apply( G.identity,
                    v -> G.permDomain[Position(G.permDomain,v)^gen] ) ) );

    if not IsBound( C.permGroupP )  then
        C.permDomain := G.permDomain;
        C.permGroupP := P;
    fi;

    # return the derived subgroup
    return C;
end;


#############################################################################
##
#F  MatGroupOps.CompositionSeries( <G> )  . . . . . . . .  composition series
##
MatGroupOps.CompositionSeries := function( G )
    local   P,  C,  U,  i;

    # compute the isomorphic permutation group for <G> and its parent
    G.operations.MakePermGroupP(G);

    # compute the derived subgroup in the permutation group
    P := CompositionSeries( G.permGroupP );
    C := Apply( P, U -> MatGroupOps.Subgroup( Parent(G),
            Apply( U.generators, gen ->
                Apply( G.identity,
                    v -> G.permDomain[Position(G.permDomain,v)^gen] ) ) ) );

    for i  in [ 1 .. Length(C) ]  do
        U := C[i];
        if not IsBound(U.permGroupP)  then
            U.permDomain := G.permDomain;
            U.permGroupP := P[i];
        fi;
    od;

    # return the composition series
    return C;
end;


#############################################################################
##
#F  MatGroupOps.ConjugacyClasses(<G>) . . conjugacy classes of a matrix group
##
MatGroupOps.ConjugacyClasses := function ( G )
    local   classes,    # conjugacy classes of <G>, result
            pclasses,   # conjugacy classes of '<G>.permGroupP'
            class,      # one conjugacy class in <pclasses>
            rep;

    # compute the isomorphic permutation group for <G> and its parent
    G.operations.MakePermGroupP( G );

    # compute the conjugacy classes in the permutation group
    pclasses := ConjugacyClasses( G.permGroupP );

    # translate every conjugacy class back
    classes := [];
    for class in pclasses  do
        rep := Representative( class );
        rep := Apply( G.identity,
                     v -> G.permDomain[Position(G.permDomain,v)^rep] );
        Add_( classes, ConjugacyClass( G, rep ) );
    od;

    # return the classes
    return classes;
end;


#############################################################################
##
#F  MatGroupOps.PermGroup(<G>)  . . . . convert a matrix group to a permgroup
##
MatGroupOps.PermGroup := function ( G )
    local   P;

    # construct the permutation group
    P := Operation( G, GAP_Union( Orbits( G, G.identity ) ) );

    # construct the bijection
    P.bijection := GroupHomomorphismByImages( P, G,
                                              P.operation.genimages,
                                              G.generators );
    P.bijection.isMapping           := true;
    P.bijection.isGroupHomomorphism := true;
    P.bijection.isInjective         := true;
    P.bijection.isMonomorphism      := true;
    P.bijection.isSurjective        := true;
    P.bijection.isEpimorphism       := true;
    P.bijection.isBijection         := true;
    P.bijection.isIsomorphism       := true;

    # return the permutation group
    return P;
end;


#############################################################################
##
#F  MatGroupOps.Stabilizer(<G>,<d>,<opr>) . . .  stabilizer in a matrix group
##
MatGroupOps.Stabilizer := function ( G, d, opr )
    local   S,          # stabilizer of <d> in <G>, result
            P;          # permutation group isomorphic to <S>

    # special case to find the stabilizer of a vector
    if IsVector(d)  and opr = OnPoints  then

        # compute the isomorphic permutation group for <G> and its parent
        G.operations.MakePermGroupP( G );

        # test whether we can handle this case
        if not d in G.permDomain  then
            return GroupOps.Stabilizer( G, d, opr );
        fi;

        # translate the vector to a point
        d := Position( G.permDomain, d );

        # find the stabilizer in the permutation group and translate back
        P := Stabilizer( G.permGroupP, d );
        S := Subgroup( Parent( G ),
            Apply( P.generators, gen ->
                Apply( G.identity,
                    v -> G.permDomain[Position(G.permDomain,v)^gen] ) ) );
        if not IsBound( S.permGroupP )  then
            S.permDomain := G.permDomain;
            S.permGroupP := P;
        fi;

    # delegate other cases
    else
        S := GroupOps.Stabilizer( G, d, opr );
    fi;

    # return the stabilizer
    return S;
end;

                      
#############################################################################
##
#F  MatGroupOps.RepresentativeOperation(<G>,<d>,<e>,<opr>)  .  representative
#F                                               of a point in a matrix group
##
MatGroupOps.RepresentativeOperation := function ( G, d, e, opr )
    local   rep;        # representative taking <d> to <e>, result

    # special case to find a conjugating element
    if d in G  and e in G  and opr = OnPoints  then

        # compute the isomorphic permutation group for <G> and its parent
        G.operations.MakePermGroupP( G );

        # translage the two matrices
        d := Permutation( d, G.permDomain );
        e := Permutation( e, G.permDomain );

        # find a conjugating permutation and translate back
        rep := RepresentativeOperation( G.permGroupP, d, e );
        if rep <> false  then
            rep := Apply( G.identity,
                         v -> G.permDomain[Position(G.permDomain,v)^rep] );
        fi;

    # special case to find a matrix taking one vector to another
    elif IsVector(d)  and IsVector(e)  and opr = OnPoints  then

        # compute the isomorphic permutation group for <G> and its parent
        G.operations.MakePermGroupP( G );

        # make sure we can handle this case
        if not d in G.permDomain  or not e in G.permDomain  then
            return GroupOps.RepresentativeOperation( G, d, e, opr );
        fi;

        # translate the two vector to points
        d := Position( G.permDomain, d );
        e := Position( G.permDomain, e );

        # find a representative and translate back
        rep := RepresentativeOperation( G.permGroupP, d, e );
        if rep <> false  then
            rep := Apply( G.identity,
                         v -> G.permDomain[Position(G.permDomain,v)^rep] );
        fi;

    # delegate other cases
    else
        rep := GroupOps.RepresentativeOperation( G, d, e, opr );
    fi;

    # return the representative
    return rep;
end;


#############################################################################
##
#F  MatGroupOps.Order(<G>,<g>)  . . . . . . . . . . . . . . order of a matrix
##
MatGroupOps.Order := MatricesOps.Order;


#############################################################################
##
#F  MatGroupOps.InvariantForm( <G> )
##
MatGroupOps.InvariantForm := function( G )
    return MatricesOps.InvariantForm( G.generators );
    end;

#############################################################################
##
#F  MatGroupOps.RightCoset(<U>,<g>) . . . . . . right coset in a matrix group
#V  RightCosetMatGroupOps operations record of right cosets in a matrix group
##
##  'MatGroupOps.RightCoset'  is the  function to create a  right coset in  a
##  matrix   group.   It  computes  a  special  element  (namely  the  matrix
##  corresponding to the  smallest element  of the corresponding coset in the
##  permutation group) of the coset and stores <U> together with this special
##  element as  representative in a record, and  enters the operations record
##  'RightCosetMatGroupOps'.
##
##  'RightCosetMatGroupOps'  is the operations  record of  right cosets in  a
##  matrix    group.     It    inherits   the    default    functions    from
##  'RightCosetGroupOps', and  overlays  the comparison function,  using  the
##  fact  that  matrix  group   cosets  have  a  special  unique  element  as
##  representative.
##
MatGroupOps.RightCoset := function ( U, g )
    local   C,          # right coset of <U> and <g>, result
            p;          # permutation corresponding to <g>

    # compute the isomorphic permutation group for <G> and its parent
    U.operations.MakePermGroupP( U );

    # compute the isomorphic permutation for <g>
    p := Permutation( g, U.permDomain );

    # compute the right coset in the permutation group
    C := U.permGroupP.operations.RightCoset( U.permGroupP, p );

    # take its representative and translate it back
    p := C.representative;
    g := Apply( U.identity, v -> U.permDomain[Position(U.permDomain,v)^p] );

    # make the domain
    C := rec( );
    C.isDomain          := true;
    C.isRightCoset      := true;

    # enter the identifying information
    C.group             := U;
    C.representative    := g;
    C.special           := g;

    # enter knowledge
    if IsBound( U.isFinite )  then
        C.isFinite      := U.isFinite;
    fi;
    if IsBound( U.size )  then
        C.size          := U.size;
    fi;

    # enter the operations record
    C.operations        := RightCosetMatGroupOps;

    # return the coset
    return C;
end;

RightCosetMatGroupOps := Copy( RightCosetGroupOps );

RightCosetMatGroupOps.\= := function ( C, D )
    local   isEql;

    # compare a right coset with minimal representative
    if IsRightCoset( C )  and IsBound( C.special )  then

        # with another right coset with minimal representative
        if IsRightCoset( D )  and IsBound( D.special )  then
            if C.group = D.group  then
                isEql := C.special = D.special;
            else
                isEql := RightCosetGroupOps.\=( C, D );
            fi;

        # with a subgroup, which is a special right coset
        elif IsGroup( D )  then
            if C.group = D  then
                isEql := C.special = D.identity;
            else
                isEql := RightCosetGroupOps.\=( C, D );
            fi;

        # with something else
        else
            isEql := RightCosetGroupOps.\=( C, D );
        fi;

    # compare a subgroup, which is a special right coset
    elif IsGroup( C )  then

        # with a right coset with minimal representative
        if IsRightCoset( D )  and IsBound( D.special )  then
            if C = D.group  then
                isEql := C.identity = D.special;
            else
                isEql := RightCosetGroupOps.\=( C, D );
            fi;

        # with something else
        else
            isEql := RightCosetGroupOps.\=( C, D );
        fi;

    # compare something else
    else
        isEql := RightCosetGroupOps.\=( C, D );
    fi;

    # return the result
    return isEql;
end;


#############################################################################
##
#F  MatGroup( <gens>, <field>[, <identity>] ) . . . . . create a matrix group
##
MatGroup := function( arg )
    local   gens, m, idmat;

    if Length( arg ) = 2 then
        if arg[1] = [] then
            Error( "sorry, need at least one element" );
        fi;
        idmat := arg[1][1] ^ 0;
    elif Length( arg ) = 3 then
        idmat := arg[3];
    else
        Error( "usage: MatGroup( <generators>, <field>[, <identity>] )" );
    fi;

    gens := [];
    for m in arg[1] do
        if m <> idmat then  Add_( gens, m );  fi;
    od;

    return rec(
        generators := gens,
        field      := arg[2],
        identity   := idmat,
        dimension  := Length( idmat ),
        isMatGroup := true,
        isGroup    := true,
        isDomain   := true,
        operations := MatGroupOps
    );
end;


#############################################################################
##
#F  RandomMatGroup( <dim>, <field>, <numgens> ) .  create random matrix group
##
RandomMatGroup := function( dim, field, k )
    local   group, mats, id, i;

    id := IdentityMat( dim, field );
    mats := [1..k];
    for i in [1..k] do
        mats[i] := RandomInvertableMat( dim, field );
    od;
    return MatGroup( mats, field, id );
end;


#############################################################################
##
#F  MatGroupOps.Transposed( <matgroup> )  . . . . . . . . . . . . . . . . . .
##
MatGroupOps.Transposed := function( group )
    local   mats, m;

    mats := [];
    for m in group.generators do
        Add_( mats, TransposedMat( m ) );
    od;
    return MatGroup( mats, group.field, group.identity );
end;

#############################################################################
##
##  Here the MeatAxe functions start
##
MTXOps := rec();

#############################################################################
##
#F  MTXOps.SpinUp( <vectors>, <matrices> )
##
##  returns a record that describes a semi-echelonized basis for the module
##  spanned by the vectors in the Apply <vectors> under the action of the
##  matrices in the Apply <matrices>.
##
MTXOps.SpinUp := function( vectors, matrices )

    local ech,       # corresponding base in echelon form
          dim,       # length of vectors
          heads,     # info about leading columns
          vector,    # loop over seed vectors
          j,         # loop over ...	
          localseed, # Apply of vectors to be processed for actual seed vector
          v,         # actual vector (reduced)
          gen;       # loop over 'matrices'

    ech   := [];
    dim   := Length( vectors[1] );
    heads := 0 * [ 1 .. dim ];

    # Loop over the seed vectors.
    for vector in vectors do

      # Spin up the space generated by the seed vector.
      localseed:= [ vector ];

      for v in localseed do

        # Reduce the vector with the basis vectors.
        v:= ShallowCopy( v );
        for j in [ 1 .. dim ] do
          if heads[j] <> 0 then
            v:= v - v[j] * ech[ heads[j] ];
          fi;
        od;

        j:= DepthVector( v );
        if j <= dim then

          # We found a new basis vector.
          Add_( ech, v / v[j] );
          heads[j]:= Length( ech );

          # Store the images under the generators in 'localseed'.
          for gen in matrices do
            Add_( localseed, v * gen );
          od;

        fi;
      od;

    od;

    # Return the basis record.
    return rec( 
                vectors            := ech,
                heads              := heads
               );
    end;

#############################################################################
##
#F  MTXOps.InvariantSubspace( <module_descr> )
#F  MTXOps.InvariantSubspace( <module_descr>, <matrix> )
#F  . . . . . . . . . . . . . . . . . . . . . . .  find an invariant subspace
##
##  This function tries to find an invariant subspace of the module as
##  suggested  by  Norton's criterion for irreducibility.  However some basic
##  calculations  have  to  be  made  if  reducibility  is  concluded  in the
##  transposed case.
##
MTXOps.InvariantSubspace := function( arg )
    local module_descr,
          m, subm, sdim, ns, tsubm, tsdim, tns,
          base, tbase, enum, line, wgts, perm, i, j,
          pivots, nonpivots, zero, F, dim, matlist;

    if Length( arg ) = 1 then
        module_descr := arg[1];
        F:= module_descr[2];
        m     := SmallCorankMatrixRecord( module_descr );
    elif Length( arg ) = 2 then
        module_descr := arg[1];
        F:= module_descr[2];
        m     := arg[2];
        ns    := NullspaceMat( m );
        if ns = [] then
            Error( "sorry, <m> has to be singular" );
        fi;
        m := rec(
            matrix    := m,
            corank    := Length( ns ),
            nullspace := RowSpace( ns, F ),
            corankGcd := Length( ns )
        );
    else
      Error( "usage: MTXOps.InvariantSubspace(<module_descr>[,<mat>])" );
    fi;

    enum   := LineEnumeration( m.nullspace );
    line   := 1;
    matlist:= module_descr[1];
    dim:= module_descr[3];
    while line <= enum.numberLines do
        subm := MTXOps.SpinUp( [ enum.line( line ) ], matlist );
        sdim := Length( subm.vectors );
        line := line + 1;
        if 0 < sdim and sdim < dim then
          return RowSpace( subm.vectors, F, 0 * subm.vectors[1] );
        fi;
    od;

    dim:= module_descr[3];

    # try to find a proper submodule in the transposed representation

    matlist:= Apply( matlist, TransposedMat );
    tns     := NullspaceMat( TransposedMat( m.matrix ) );
    tsubm   := MTXOps.SpinUp( [ tns[1] ], matlist ); 
    tsdim   := Length( tsubm.vectors );

      # Check whether one nontrivial vector in the right nullspace of the
      # matrix spans a proper submodule (under action from the left).

      if tsdim < dim then

        # We found an invariant subspace for the action of the transposed
        # matrices, and now compute from this the submodule for the original
        # matrices.

        tbase := Copy( tsubm.vectors );
        TriangulizeMat( tbase );
        pivots:= Apply( tbase, DepthVector );
        nonpivots:= GAP_Difference( [ 1 .. dim ], pivots );
  
        # determine the submodule base for the normal case
  
        zero := 0 * tbase[1];
        base := [];
  
        for i in [ 1 .. dim - tsdim ] do
          base[i]:= Copy( zero );
          base[i][ nonpivots[i] ]:= One( F );
          for j in [ 1 .. tsdim ] do
            base[i][ pivots[j] ]:= - tbase[j][ nonpivots[i] ];
          od;
        od;

        return RowSpace( base, F, zero );

    else
        return m.corankGcd;
    fi;
end;


# #############################################################################
# ##
# #F  MTXOps.IsInvariantSubspace( <matgroup>, <rowspace> ) . . . . . . . .
# #F  . . . . . . . . . . . . . . . . . . . . . . test if is invariant subspace
# ##
# MTXOps.IsInvariantSubspace := function( group, vs )
#     return Submodule( RowModule( group ), vs ).abelianGroup = vs;
# end;
# 
# 
# #############################################################################
# ##
# #F  IsInvariantSubspace( <domain>, <object> ) . . . . . . . . . . . . . . . .
# ##
# IsInvariantSubspace := function( D, obj )
#     local   isinv;
# 
#     if IsDomain( D ) and IsBound( D.operations.IsInvariantSubspace ) then
#         isinv := D.operations.IsInvariantSubspace( D, obj );
#     else
#         Error( "sorry, can't test on invariant subspace for <domain>" );
#     fi;
#     return isinv;
# end;


#############################################################################
##
#F  MTXOps.IrreducibilityTest( <module_descr>[, <matrix>] )  . . . . . .
#F  . . . . . . . . . . . . . . . . . .  test whether a module is irreducible
##
MTXOps.IrreducibilityTest := function( arg )
    local   subs;

    if   Length( arg ) = 1 then
        subs := MTXOps.InvariantSubspace( arg[1] );
    elif Length( arg ) = 2 then
        subs := MTXOps.InvariantSubspace( arg[1], arg[2] );
    else
        Error( "usage: IrreducibilityTest( <module_descr>[, <mat>] )" );
    fi;

    if IsInt( subs ) then
        return rec(
            isIrreducible := true,
            isReducible   := false,
            degree        := subs
        );
    else
        return rec(
            isIrreducible     := false,
            isReducible       := true,
            invariantSubspace := subs
        );
    fi;
end;


#############################################################################
##
#F  MTXOps.AbsoluteIrreducibilityTest( <module_descr>[, <matrix>] )  . .
#F  . . . . . . . . . . . . . test whether a module is absolutely irreducible
##
MTXOps.AbsoluteIrreducibilityTest := function( arg )
    local   subs, deg, module_descr, mats, field, m;

    if Length( arg ) = 1 then
        module_descr := arg[1];
        subs  := MTXOps.InvariantSubspace( arg[1] );
    elif Length( arg ) = 1 then
        module_descr := arg[1];
        subs  := MTXOps.InvariantSubspace( arg[1], arg[2] );
    else
      Error( "usage: AbsoluteIrreducibilityTest( <module_descr>[,<mat>] )" );
    fi;

    if IsInt( subs ) then
        deg := subs;
        if deg = 1 then
            return rec(
                isIrreducible         := true,
                isReducible           := false,
                isAbsoluteIrreducible := true,
                degree                := deg
            );
        fi;

        if IsFinite( arg[1] ) then

            # construct the new module description over the bigger field

            field := GAP_GF( GAP_Size( module_descr[2] ) ^ deg );
            module_descr:= [ module_descr[1], field, module_descr[3] ];

            if Length( arg ) = 1 then
                subs := MTXOps.InvariantSubspace( module_descr );
            else
                subs := MTXOps.InvariantSubspace( module_descr, arg[2] );
            fi;

            if IsInt( subs ) then
                return rec(
                    isIrreducible         := true,
                    isReducible           := false,
                    isAbsoluteIrreducible := true,
                    degree                := deg * subs
                );
            else
                return rec(
                    isIrreducible         := true,
                    isReducible           := false,
                    isAbsoluteIrreducible := false,
                    degree                := deg,
                    invariantSubspace     := subs
                    
                );
            fi;

        else
            Error( "sorry, don't know how to extend infinte fields" );
        fi;
    else
        return rec(
            isIrreducible         := false,
            isReducible           := true,
            isAbsoluteIrreducible := false,
            invariantSubspace     := subs
        );
    fi;
end;


#############################################################################
##
#F  MTXOps.MatRepresentation( <gens>, <vs> ) .  action of <gens> on <vs>
##
MTXOps.MatRepresentation := function( gens, vs )
    local    base, new, mat, m, b;

    # compute a base for <vs>
    base := Base(vs);

    # compute the new operation
    new := [];
    for mat  in gens  do
        m := [];
        for b  in base  do
            Add_( m, Coefficients( vs, b*mat ) );
        od;
        Add_( new, m );
    od;
    return rec( range := new, base := base );

end;

#############################################################################
##
#F  MTXOps.CompositionFactors( <module_descr> )  . . composition factors
##
MTXOps.CompositionFactors := function( module_descr )
    local   factors,  bases,  groups,  decompose, F, idmat;

    F:= module_descr[2];
    idmat:= IdentityMat( module_descr[3], F );

    # catch trivial case
    # 1995/09/27 sam
    # The 'groups' component was 'group' before,
    # so it seems the trivial case was never touched seriously.
    if 0 = Length( module_descr[1] )  then
        return rec( 
            bases   := Apply( idmat, x -> [x] ),
            groups  := Apply( [ 1 .. module_descr[3] ],
                             x -> [ [ [[ One( F )]] ], F, 1 ] ),
            factors := Apply( [ 1 .. module_descr[3] ], x -> [] )
        );
    fi;

    decompose := function ( m, b )
        local   descr, subs, vs, rep;

        # construct an invariant subspace
        descr:= [ m, F, Length( m[1] ) ];
        subs := MTXOps.InvariantSubspace( descr );

        # if group operates irreducible,  Add_ generators to <factors>
        if IsInt(subs) then
            Add_( factors, m );
            Add_( groups,  descr );
            Add_( bases,   b );

        # otherwise decompose again
        else

            # compute operation on the factor space
            vs  := VectorSpace( IdentityMat( descr[3], F ), F ) mod subs;
            rep := MTXOps.MatRepresentation( m, vs );
            decompose( rep.range, rep.base * b );

            # compute the operation on the subspace
            rep := MTXOps.MatRepresentation( m, subs );
            decompose( rep.range, rep.base * b );
        fi;
    end;

    factors := [];
    bases   := [];
    groups  := [];
    decompose( module_descr[1], idmat );

    return rec(
        groups  := groups,
        factors := factors,
        bases   := bases
    );
end;


# 1995/09/27 sam
# The following is used nowhere in the library.
# Anyhow, it was no good idea to choose a group as argument.
#
# #############################################################################
# ##
# #F  MTXOps.GeneratorsConsituents( <grp> )  . . . . constituents of <grp>
# ##
# ##
# ##  'MTXOps.GeneratorsConstituents' returns the   Apply of generators  of
# ##  the constituents of <grp>.
# ##
# MTXOps.GeneratorsConstituents := function( arg )
# 
#     local   grp,           # the matrixgroup
#             facs,          # the composition factors
#             consti,        # the constituents
#             equi,          # boolean
#             i, j; 
# 
#     # check arguments
#     if Length(arg) = 1  then
#         grp  := arg[1];
#         facs := CompositionFactors(grp);
#     elif Length(arg) = 2  then
#         grp  := arg[1];
#         facs := arg[2];
#     else
#         Error( "usage: GeneratorsConstituents( <grp>[, <factors>] )" );
#     fi;
# 
#     # calculate constituens 
#     consti := [1];
#     for i  in [ 2 .. Length(facs.factors) ]  do
#         j := 0;
#         equi := false;
#         while not equi and j < Length(consti)  do
#             j := j+1;
#             equi := EquivalenceTest ( facs.groups[consti[j]],
#                                       facs.groups[i] ).areequivalent;
#         od;
#         if not equi  then
#             Add_( consti, i );
#         fi;
#     od;
# 
#     return Apply( consti, x -> facs.factors[x] );
# end;


#############################################################################
##
#F  MTXOps.SpinedBase( <module_descr>, <base>, <dim> ) . . . . spin base
##
##  'SpinedBase'  spinns  each (normed) vector  of the vector space generated
##  by  <base>  as module  of  <group>. It returns a  Apply  of bases  for all
##  spinned vectorspaces, which dimensions are less or equal than <dim>
##
MTXOps.SpinedBase := function( module_descr, base, dim )

    local   elem, elements, result, v, base, space, Apply, dimension, w, gen, 
            z, field, i;

    field := module_descr[2];

    # look for normed element
    elements := NormedVectors( VectorSpace( base, field ) );

    # spinn up
    result := [];
    if Length( module_descr[1] ) = 0 and dim = 1  then
        for v  in elements  do
            Add_( result, [v] );
        od;
        return result;
    fi;
    for v  in elements  do
        base  := [ v ];
        space := VectorSpace( base, field );
        Apply  := [ v ];
        dimension := 1;
        while 0 < Length(Apply) and dimension <= dim  do
            w := Apply[1];
            for gen in module_descr[1] do
                z := w * gen;
                if not z  in space  then
                    Add_( base, z );
                    Add_( Apply, z );
                    dimension := dimension + 1;
                    space := VectorSpace( base, field );
                fi;
            od;
            Apply := Filtered( Apply, x -> x <> w );
        od;
        if dimension = dim  then
            TriangulizeMat( base );
            Add_( result, base );
        fi;
    od;
    return result;
end; 


#############################################################################
##
#F  InfoPeak2( <arg> )  . . . . . . . . . . . . . peak word debug information
##
if not IsBound(InfoPeak2) then InfoPeak2 := Ignore;   fi;


#############################################################################
##
#F  MTXOps.Peakwords( <group>, <compositionfactors>) . . . . . peakwords
##
##  Peakwords returns a record consisting of the entry  'peaks'  which  gives  
##  normalized kernstable peakword to each constituent  of  <group>  and  the
##  entry 'dims' containing the dimensions of the constituents of <group>.
##
if not IsBound(PEAKWORD_LIMIT)  then PEAKWORD_LIMIT := 50;  fi;

MTXOps.NewElement := function( facslist, matlist ) 

    local found, i, j, h, k, new; 

    found := false; 
    i := 1; 
    j := 1; 

    # look for a group element which is not in Apply and Add_ it
    while not found  do
        new := matlist[ i ] * matlist[ j ];
        if  not new  in matlist  then
            found := true;
            Add_ ( matlist, new );

            # calculate compfacs and Add_ them
            for h  in [ 1..Length( facslist ) ]  do
                Add_( facslist[h], facslist[h][i]*facslist[h][j] );
            od;
        else
            if i = Length(matlist)  then
                i := 1;
                j := j + 1;
            else
                i := i + 1;
            fi;
        fi;
    od;
end;

MTXOps.IsPeakWord := function( facs, equiv )

    local   kernels, index, help, length, i, j, k;

    kernels := Apply( facs, x -> Length( NullspaceMat(x) ) );

    i := DepthVector( kernels );
    if i > Length( kernels )  then
        return false;
    fi;
    index := i;

    # find index in equiv
    for i  in [ 1..Length(equiv) ]  do
        if index  in equiv[i]  then
            j := i;
        fi;
    od;

    # test whether all other composition factors have trivial kernel
    help := GAP_Difference( [ 1..Length(facs) ], equiv[j] );
    if Length( help ) > 0  then
        help := Apply( help, x -> kernels[x] );
        if Length( Set( help ) ) > 1 or help[1] <> 0  then
            return false;
        fi;
    fi;

    # test whether all equivalent composition factors have equal kernel
    help := Apply( equiv[j], x -> kernels[x] );
    if Length( Set( help ) ) > 1 then
        return false;
    fi;
    length := help[1];

    # test if word is kernstable
    help := Apply( equiv[j], x -> Length( NullspaceMat( facs[x]^2 ) ) );
    if Length( Set( help ) ) > 1 or length <> help[1] then
        return false;
    fi;

    return rec( index    := j,
                equilist := equiv[j],
                length   := length );
end;

        
MTXOps.IsNormalizedPeak := function( dims, facs, facslist, peak, q )

    local j, matgrp;
   
    # catch trivial case
    if peak.length = 1  then
        return true;
    fi;

    j := peak.index;

    # change dims[j] if necessary
    if IsBool( dims[j] )  then 
        dims[j] := [ peak.length ];
    elif IsList( dims[j] )  then
        dims[j] := [ Gcd( dims[j][1], peak.length ) ];
    fi;
 
    # get dimension if possible and necessary
    if IsList( dims[j] ) and q^dims[j][1] <= 2^16  then
        dims[j] := Length( MTXOps.CompositionFactors(
                          [ facslist[peak.equilist[1]],
                            GAP_GF(q^dims[j][1]),
                            Length( facslist[ peak.equilist[1] ][1] ) ]
                            ).factors );
    fi;

    # test wether it is a normalized peakword
    if IsInt( dims[j] )  then
        if peak.length = dims[j]  then
            return true;
        fi;
    fi;
    return false;
end;
         
MTXOps.Peakwords := function( descr_module, compfacs )

    local field,                 # field of matgrp
          q,                     # char of field 
          matlist, l,            # gens of matgrp and their number         
          facslist,              # Apply of composition factors
          equiv,                 # gives constituents
          const,                 # constituents
          length,                # number of constituents
          peakwords,             # peakwords
          dims,                  # Apply of dimensions of kernels
          dimstest,              # number of unsuccesful searches
          count,                 # number of found peakwords
          v,                     # a vector corresponding to a element of the
                                 # groupalgebra
          facs, fac,             # composition factor to v
          peak,                  # true if it is a peakword to peak.index
          elem,                  # an element of the groupalgebra
          found, i, j, h;

    # trivial case
    if 0 = Length( descr_module[1] )  then
        return [ NullMat( descr_module[3], descr_module[3],
                          Zero( descr_module[2] ) ) ];
    fi;

    field    := descr_module[2];
    q        := field.char;
    matlist  := Copy( descr_module[1] );
    l        := Length( matlist );
    facslist := Copy( compfacs.factors );

    # set up equivalence array    
    const := [];
    equiv := [];
    for i  in [ 1..Length(compfacs.factors) ]  do
        j := 1;
        found := false;
        while j <= Length(const)  and not found  do
            if EquivalenceTest( compfacs.groups[i], const[j] ).areequivalent 
            then
                Add_( equiv[j], i );
                found := true;
            else
                j := j + 1;
            fi;
        od;
        if not found  then
            Add_( const, compfacs.groups[i] );
            Add_( equiv, [i] );
        fi;
    od;

    length    := Length ( equiv );
    peakwords := Apply ( [ 1..length ], x -> false ); 
    dims      := Apply ( [ 1..length ], x -> false ); 
    dimstest  := Apply ( [ 1..length ], x -> 0 );
    count     := 0;                                

    InfoPeak2("#I  start searching for peakwords: ", length," constituents\n");

    i := 0;
    while count < length  do

        # v is the coeffientvektor to an element via matlist 
        v := CoefficientsQadic( i, q );
        MakeVecFFE( v, field.one );

        # generate a new matrix if necessary 
        if Length(v) > Length(matlist)  then
            MTXOps.NewElement( facslist, matlist );
        fi;

        # calculate facs of the element represented by <v>
        facs := [ ];
        for j  in [ 1..Length ( facslist ) ]  do
            fac := field.zero * facslist[ j ][ 1 ];
            for h  in [ 1..Length( v ) ]  do
                fac := fac + facslist[ j ][ h ] * v[ h ];
            od;
            Add_( facs, fac );
        od;

        # look if the element is a peakword for a constituent
        peak := MTXOps.IsPeakWord( facs, equiv );
        if IsRec( peak )  then

            # its a peakword to the jth constituent 
            j := peak.index;
            if IsBool(peakwords[j]) and dimstest[j] < PEAKWORD_LIMIT  then

                # test whether it is normalized
                if MTXOps.IsNormalizedPeak( dims, facs,
                                                       facslist, peak, q )  
                then
                    peakwords[j] := v;
                    count        := count + 1;
                    InfoPeak2("#I  found peakword to ",j,"th constituent \n");
                else
                    dimstest[j]  := dimstest[j]+1; 
                    InfoPeak2("#I  unsuccessful search for peakword to ",j);
                    InfoPeak2("th constituent \n");
                fi;
            fi;
        fi;
        i := i + 1;
    od;

    # calculate matrices from exponent vectors
    for i  in [ 1..length ]  do
        if not IsBool ( peakwords[ i ] )  then
            elem := field.zero * matlist[ 1 ];
            for j  in [ 1..Length( peakwords[ i ] ) ]  do
                elem := elem + peakwords[ i ][ j ] * matlist[ j ];
            od;
            peakwords[ i ] := elem;
        fi;
    od;
    return rec( peaks := peakwords,
                dims  := Apply( const, x -> x[3] ) );
end;    

#############################################################################
##
#F  MTXOps.EquivalenceTest( <module_descr1>, <module_descr2> )
##
MTXOps.EquivalenceTest := function( module_descr1, module_descr2 )
    local   fp1, fp2, module, minpos, ns1, ns2, crk, result, base, line,
            enum, y1, x, x_inv, i;

   if    module_descr1[3] <> module_descr2[3]
      or module_descr1[2] <> module_descr2[2]
      or Length( module_descr1[1] ) <> Length( module_descr2[1] ) then
        return rec( areequivalent := false );
    fi;

    # first check the trivial cases

    if module_descr1[1] = module_descr2[1] then
        return rec(
            areequivalent        := true,
            transformationMatrix := IdentityMat( module_descr1[3],
                                                 module_descr1[2] ) );
    elif module_descr1[3] = 1 then
        return rec( areequivalent := false );
    fi;

    fp1 := FingerprintEx( module_descr1 );
    fp2 := FingerprintEx( module_descr2 );

    if fp1.coranks <> fp2.coranks then
        return rec( areequivalent := false );
    fi;
    if fp1.failed then
        Error( "sorry, fingerprint failed" );
    fi;

    minpos := Position( fp1.coranks,
                        Minimum( Filtered( fp1.coranks, x -> x <> 0 ) ) );

    ns1 := fp1.nullspaces[minpos];
    ns2 := fp2.nullspaces[minpos];
    crk := fp1.coranks[minpos];

    result := rec( areequivalent := false );

    enum := LineEnumeration( ns1 );
    line := 1;
    repeat
        base := SpinUpStandard( [ enum.line( line ) ],
                                module_descr1[1] ).vectors;
        line := line + 1;
    until Length( base ) = module_descr1[3] or line > enum.numberLines;

    if Length( base ) = module_descr1[3] then
        y1   := base;
        enum := LineEnumeration( ns2 );
        line := 1;

        repeat

            base := SpinUpStandard( [ enum.line( line ) ],
                                    module_descr2[1] ).vectors;
            line := line + 1;

            if Length( base ) = module_descr2[3] then
                x     := base^-1 * y1;
              # x_inv := x ^ -1;
                i     := 1;
                result.areequivalent := true;
                while i <= Length( module_descr2[1] )
                      and result.areequivalent do
                    if x * module_descr1[1][i]
                       <> module_descr2[1][i] * x then
                        result.areequivalent := false;
                    fi;
                    i := i + 1;
                od;
                if result.areequivalent then
                    result.transformationMatrix := x;
                fi;
            fi;

        until result.areequivalent or line > enum.numberLines;

    else
        Error( "sorry, failed finding a cyclic generator" );
    fi;
    return result;
end;

#############################################################################
##
#F  EquivalenceTest( <domain>, <domain> ) . . . . . .  equivalence of domains
##
EquivalenceTest := function( D1, D2 )
    local   equiv;

    if IsDomain( D1 ) and IsBound( D1.operations.EquivalenceTest ) then
        equiv := D1.operations.EquivalenceTest( D1, D2 );
    else
        equiv := MTXOps.EquivalenceTest( D1, D2 );
    fi;
    return equiv;
end;


#############################################################################
##
#F  AddNextMatrixFunction( <function> ) . . . . . set the NextMatrix function
##
NEXTMATRIX := false;
AddNextMatrixFunction := function( func )
   NEXTMATRIX:= func;
end;

#############################################################################
##
#F  NextMatrix( <module_descr>, <info> )  . . . . . . .  generate next matrix
##
NextMatrix := function( module_descr, info )
    local   m;

    if NEXTMATRIX = false then
        AddNextMatrixFunction( NextMatrix2 );
    fi;
    m := NEXTMATRIX( module_descr, info );
    return m;
end;

#############################################################################
##
#F  SmallCorankMatrixRecord( <module_descr> ) . . . random matrix of low rank
##
SmallCorankMatrixRecord := G -> RandomMatrixRecord( G, 1 );

#############################################################################
##
#F  RandomMatrixRecord( <module_descr> ) . . . . . . . . choose random matrix
#F  RandomMatrixRecord( <module_descr>, <integer> )
##
##  returns
##
##      rec(
##          matrix    := <matrix>,
##          corank    := <integer>,
##          nullspace := <rowspace>,
##          corankGcd := <integer>
##      )
##
RandomMatrixRecord := function( arg )
    local   module_descr, maxcrk, m, crk, gcd, ns, minmat, mincrk, minns, zero, info;

    if   Length( arg ) = 1 then   module_descr := arg[1];  maxcrk := module_descr[3];
    elif Length( arg ) = 2 then   module_descr := arg[1];  maxcrk := arg[2];
    else Error( "usage: RandomMatrixRecord( <matmodule_descr>[, <maxcorank>] )" );
    fi;

    zero := NullMat( module_descr[3], module_descr[3], module_descr[2] );

    if module_descr[1] = [] or module_descr[3] = 1 then
        return rec(
            matrix    := zero,
            corank    := module_descr[3],
            nullspace := RowSpace( module_descr[3], module_descr[2] ),
            corankGcd := module_descr[3]
        );
    fi;

    info := rec();
    m    := NextMatrix( module_descr, info );

    ns  := NullspaceMat( m );
    crk := Length( ns );
    gcd := crk;
    if 0 < crk and crk <= maxcrk then
        return rec(
            matrix    := m,
            corank    := crk,
            nullspace := RowSpace( ns, module_descr[2], zero[1] ),
            corankGcd := gcd
        );
    elif 0 < crk then
        minmat := m;
        mincrk := crk;
        minns  := ns;
    else
        minmat := 0 * m;
        mincrk := module_descr[3];
        minns  := RowSpace( module_descr[3], module_descr[2] ).generators;
        gcd    := module_descr[3];
    fi;

    m := NextMatrix( module_descr, info );
    while m <> false do
        ns  := NullspaceMat( m );
        crk := Length( ns );
        gcd := Gcd( gcd, crk );
        if 0 < crk and crk <= maxcrk then
            return rec(
                matrix    := m,
                corank    := crk,
                nullspace := VectorSpace( ns, module_descr[2], zero[1] ),
                corankGcd := gcd
            );
        elif 0 < crk and crk < mincrk then
            minmat := m;
            mincrk := crk;
            minns  := ns;
        fi;

        m := NextMatrix( module_descr, info );
    od;

    return rec(
        matrix    := minmat,
        corank    := mincrk,
        nullspace := VectorSpace( minns, module_descr[2], zero[1] ),
        corankGcd := gcd
    );
end;


#############################################################################
##
#F  FingerprintEx( <module_descr> )
#F  FingerprintEx( <module_descr>, <integer> )
##
##  returns
##
##      rec(
##          matrices   := [ <matrix> ],
##          coranks    := [ <integer> ],
##          nullspaces := [ <rowspace> ],
##          corankGcd  := <integer>,
##          failed     := <boolean>
##      )
##
FingerprintEx := function( arg )
    local module_descr,
          maxcrk, corank, nullspace, matrix, m, zero, info, gcd, i;

    if   Length( arg ) = 1 then   module_descr := arg[1];  maxcrk := 0;
    elif Length( arg ) = 2 then   module_descr := arg[1];  maxcrk := arg[2];
    else Error( "usage: FingerprintEx( <module_descr>[, <maxcorank>] )" );
    fi;

    # if the module is easy to handle do it now. Remember that a corank of
    # one (although possibly trivial) means the fingerprint did not fail.

    zero := NullMat( module_descr[3], module_descr[3], module_descr[2] );

    if module_descr[1] = [] or module_descr[3] = 1 then
        return rec(
            coranks    := [ 1 ],
            nullspaces := [ RowSpace( 1, module_descr[2] ) ],
            matrices   := [ zero ],
            failed     := false
        );
    fi;

    corank    := [];
    nullspace := [];
    matrix    := [];

    gcd  := module_descr[3];
    info := rec();
    m    := NextMatrix( module_descr, info );
    i    := 0;

    while m <> false do

        i := i + 1;
        nullspace[i] := RowSpace( NullspaceMat( m ), module_descr[2],
                                  zero[1] );
        corank[i] := Dimension( nullspace[i] );
        matrix[i] := m;
        gcd := Gcd( gcd, corank[i] );

        if 0 < corank[i] and corank[i] <= maxcrk then
            return rec(
                coranks    := corank,
                nullspaces := nullspace,
                matrices   := matrix,
                corankGcd  := gcd,
                failed     := false
            );
        fi;

        m := NextMatrix( module_descr, info );
    od;

    return rec(
        coranks    := corank,
        nullspaces := nullspace,
        matrices   := matrix,
        corankGcd  := gcd,
        failed     := maxcrk <> 0 or Maximum( corank ) = 0
    );
end;


#############################################################################
##
#F  ClassicNextMatrix( <module_descr>, <info> )
##
##  This  function  implements  the  classical  fingerprint  methods  for two
##  matrices  as  proposed  by  R.A.   Parker  in  his  first paper about the
##  meat-axe system.
##
ClassicNextMatrix := function( module_descr, info )
    local   m;

    if info = rec() then
        info.R := [];
        info.number := 0;

        if Length( module_descr[1] ) >= 2 then
            info.R[1] := module_descr[1][1];
            info.R[2] := module_descr[1][2];
        elif Length( module_descr[1] ) = 1 then
            info.R[1] := module_descr[1][1];
            info.R[2] := IdentityMat(module_descr[3],module_descr[2]);
        else
            info.R[1] := IdentityMat(module_descr[3],module_descr[2]);
            info.R[2] := IdentityMat(module_descr[3],module_descr[2]);
        fi;

        info.R[3] := info.R[1] * info.R[2];
        info.R[4] := info.R[1] + info.R[2];
        info.R[5] := info.R[3] + info.R[4];
    fi;

    info.number := info.number + 1;

    if   info.number = 1 then
        m := info.R[5];
    elif info.number = 2 then
        info.R[6] := info.R[3] * info.R[2];
        info.R[7] := info.R[5] + info.R[6];
        m         := info.R[7];
    elif info.number = 3 then
        info.R[8] := info.R[2] * info.R[7];
        info.R[9] := info.R[1] + info.R[8];
        m         := info.R[9];
    elif info.number = 4 then
        info.R[10] := info.R[2] + info.R[9];
        m          := info.R[10];
    elif info.number = 5 then
        info.R[11] := info.R[3] + info.R[10];
        m          := info.R[11];
    elif info.number = 6 then
        info.R[12] := info.R[1] + info.R[11];
        m          := info.R[12];
    else
        m := false;
    fi;

    return m;
end;


#############################################################################
##
#F  ExtendedClassicNextMatrix( <module_descr>, <info> )
##
##  This  function  implements  the  classical  fingerprint  methods  for two
##  matrices  as  proposed  by  R.A.   Parker  in  his  first paper about the
##  meat-axe  system, extended for the case when more generating matrices are
##  given.
##
ExtendedClassicNextMatrix := function( module_descr, info )
    local   m;

    if info = rec() then

        info.R      := [];
        info.number := 0;
        info.next   := false;
        if Length( module_descr[1] ) >= 2 then
            info.R[1] := module_descr[1][1];   info.1 := 1;
            info.R[2] := module_descr[1][2];   info.2 := 2;
        elif Length( module_descr[1] ) = 1 then
            info.R[1] := IdentityMat(module_descr[3],module_descr[2]);
            info.2 := 0;
            info.R[2] := module_descr[1][1];   info.1 := 1;
        else
            info.R[1] := IdentityMat(module_descr[3],module_descr[2]);
            info.1 := 0;
            info.R[2] := IdentityMat(module_descr[3],module_descr[2]);
            info.2 := 0;
        fi;
        info.R[3] := info.R[1] * info.R[2];
        info.R[4] := info.R[1] + info.R[2];
        info.R[5] := info.R[3] + info.R[4];

    elif info.next then

        info.2 := info.2 + 1;
        if info.2 > Length( module_descr[1] ) then
            info.1 := info.1 + 1;
            info.2 := info.1 + 1;
        fi;
        if info.2 > Length( module_descr[1] ) then
            return false;
        fi;

        info.R      := [];
        info.number := 0;
        info.next   := false;
        info.R[1] := module_descr[1][info.1];
        info.R[2] := module_descr[1][info.2];
        info.R[3] := info.R[1] * info.R[2];
        info.R[4] := info.R[1] + info.R[2];
        info.R[5] := info.R[3] + info.R[4];

    fi;

    info.number := info.number + 1;

    if   info.number = 1 then
        m := info.R[5];
    elif info.number = 2 then
        info.R[6] := info.R[3] * info.R[2];
        info.R[7] := info.R[5] + info.R[6];
        m         := info.R[7];
    elif info.number = 3 then
        info.R[8] := info.R[2] * info.R[7];
        info.R[9] := info.R[1] + info.R[8];
        m         := info.R[9];
    elif info.number = 4 then
        info.R[10] := info.R[2] + info.R[9];
        m          := info.R[10];
    elif info.number = 5 then
        info.R[11] := info.R[3] + info.R[10];
        m          := info.R[11];
    elif info.number = 6 then
        info.R[12] := info.R[1] + info.R[11];
        m          := info.R[12];
    else
        info.next := true;
        m := NextMatrix( module_descr, info );
    fi;

    return m;
end;


#############################################################################
##
#F  NextMatrix1( <module_descr>, <info> )
##
##
NextMatrix1 := function( module_descr, info )
    local   m, g, i;

    if info = rec() then
        if module_descr[1] = [] then
            m := IdentityMat(module_descr[3],module_descr[2]);
        else
            m := module_descr[1][1];
            for i in [2..Length( module_descr[1] )] do
                m := m * module_descr[1][i];
            od;
            info.product := m;
            for i in [1..Length( module_descr[1] )] do
                m := m + module_descr[1][i];
            od;
        fi;

        info.matrix    := m;
        info.generator := 0;
        info.nextgen   := true;

        return m;
    fi;

    if info.nextgen then
        info.generator := info.generator + 1;
        if info.generator > Length( module_descr[1] ) then
            return false;
        fi;

        m := info.matrix;
        g := module_descr[1][info.generator];

        m := m * g;
        for i in [1..info.generator] do
            m := m + module_descr[1][i];
        od;
        info.matrix  := m;
        info.nextgen := false;
    else
        m := info.matrix;
        m := m + info.product;
        info.matrix  := m;
        info.nextgen := true;
    fi;

    return m;
end;


#############################################################################
##
#F  NextMatrix2( <module_descr>, <info> )
##
##
NextMatrix2 := function( module_descr, info )
    local   m, g, ew, ew2, i;

    if info = rec() then
        if module_descr[1] = [] then
            m := IdentityMat(module_descr[3],module_descr[2]);
        else
            m := module_descr[1][1];
            for i in [2..Length( module_descr[1] )] do
                m := m * module_descr[1][i];
            od;
            info.product := m;
            for i in [1..Length( module_descr[1] )] do
                m := m + module_descr[1][i];
            od;
        fi;

        if IsFinite( module_descr[2] ) then
            ew     := module_descr[2].one;
            ew2    := module_descr[2].zero;
            info.eigenvalues := [ ew2 ];
            for i in [2..Order( module_descr[2], module_descr[2].root )+1] do
                ew := ew * module_descr[2].root;
                info.eigenvalues[i] := ew - ew2;
                ew2 := ew;
            od;
        else
            info.eigenvalues := [ 0, -3, 1, 1, 2, 1, 1 ];
        fi;

        info.matrix    := m;
        info.generator := 0;
        info.nextgen   := true;
        info.nextev    := 1;

        return m;
    fi;

    if info.nextgen then
        if info.nextev = 1 then

            info.generator := info.generator + 1;
            if info.generator > Length( module_descr[1] ) then
                return false;
            fi;

            m := info.matrix;
            g := module_descr[1][info.generator];

            m := m * g;
            for i in [1..info.generator] do
                m := m + module_descr[1][i];
            od;
            info.matrix := m;
            info.nextev := info.nextev + 1;
            if info.nextev > Length( info.eigenvalues ) then
                info.nextgen := false;
                info.nextev  := 1;
            fi;
        else
            if info.nextev = 2 then
                info.matrix2 := info.matrix;
            fi;

            m := Copy( info.matrix2 );
            for i in [1..module_descr[3]] do
                m[i][i] := m[i][i] + info.eigenvalues[info.nextev];
            od;
            info.matrix2 := m;

            info.nextev := info.nextev + 1;
            if info.nextev > Length( info.eigenvalues ) then
                info.nextgen := false;
                info.nextev  := 1;
            fi;
        fi;
    else
        if info.nextev = 1 then
            m := info.matrix;
            m := m + info.product;
            info.matrix := m;
            info.nextev := info.nextev + 1;
            if info.nextev > Length( info.eigenvalues ) then
                info.nextgen := true;
                info.nextev  := 1;
            fi;
        else
            if info.nextev = 2 then
                info.matrix2 := info.matrix;
            fi;

            m := Copy( info.matrix2 );
            for i in [1..module_descr[3]] do
                m[i][i] := m[i][i] + info.eigenvalues[info.nextev];
            od;
            info.matrix2 := m;

            info.nextev := info.nextev + 1;
            if info.nextev > Length( info.eigenvalues ) then
                info.nextgen := true;
                info.nextev  := 1;
            fi;
        fi;
    fi;

    return m;
end;


#############################################################################
##
#F  NextMatrix3( <module_descr>, <info> )
##
##
NextMatrix3 := function( module_descr, info )
    local   m;

    if info = rec() then
        info.R      := [];
        info.number := 0;
        info.next   := false;
        if Length( module_descr[1] ) >= 2 then
            info.R[1] := module_descr[1][1];   info.1 := 1;
            info.R[2] := module_descr[1][2];   info.2 := 2;
        elif Length( module_descr[1] ) = 1 then
            info.R[1] := IdentityMat(module_descr[3],module_descr[2]);
            info.2 := 0;
            info.R[2] := module_descr[1][1];   info.1 := 1;
        else
            info.R[1] := IdentityMat(module_descr[3],module_descr[2]);
            info.1 := 0;
            info.R[2] := IdentityMat(module_descr[3],module_descr[2]);
            info.2 := 0;
        fi;

    elif info.next then

        info.2 := info.2 + 1;
        if info.2 > Length( module_descr[1] ) then
            info.1 := info.1 + 1;
            info.2 := info.1 + 1;
        fi;
        if info.2 > Length( module_descr[1] ) then
            return false;
        fi;

        info.R      := [];
        info.number := 0;
        info.next   := false;
    fi;

    info.number := info.number + 1;

    if info.number = 1 then
        info.R[3] := info.R[1] * info.R[2];
        info.R[4] := info.R[1] + info.R[2];
        info.R[5] := info.R[3] + info.R[4];
        m := info.R[5];
    elif info.number = 2 then
        info.R[6] := info.R[3] * info.R[2];
        info.R[7] := info.R[5] + info.R[6];
        m := info.R[7];
    elif info.number = 3 then
        info.R[8] := info.R[2] * info.R[7];
        info.R[9] := info.R[1] + info.R[8];
        m := info.R[9];
    elif info.number = 4 then
        info.R[10] := info.R[2] + info.R[9];
        m := info.R[10];
    elif info.number = 5 then
        info.R[11] := info.R[3] + info.R[10];
        m := info.R[11];
    elif info.number = 6 then
        info.R[12] := info.R[1] + info.R[11];
        m := info.R[12];
    elif 7 <= info.number and info.number <= 18 then
        m := info.R[info.number-6]
             - IdentityMat(module_descr[3],module_descr[2]);
    elif 19 <= info.number and info.number <= 30 then
        m := info.R[info.number-18]
             + IdentityMat(module_descr[3],module_descr[2]);
    else
        info.next := true;
        m := NextMatrix( module_descr, info );
    fi;

    return m;
end;

