Quelle  matrep.gi   Sprache: Shell

#############################################################################
##
#W  matrep.gi                  Polycyclic                       Werner Nickel
##

InstallGlobalFunction( "IsMatrixRepresentation",
function( G, matrices )
    local   coll,  conjugates,  d,  I,  i,  j,  conj,  rhs,  k;

    coll := Collector( G );
    conjugates := coll![PC_CONJUGATES];

    d := NumberOfGenerators( coll );

    I := matrices[1]^0;
    for i in [1..d] do
        for j in [i+1..d] do
            conj := matrices[j]^matrices[i];
            if IsBound( conjugates[j] ) and IsBound( conjugates[j][i] ) then
                rhs := I;
                for k in [1,3..Length(conjugates[j][i])-1] do
                    rhs := rhs *
                           matrices[ conjugates[j][i][k] ] ^
                           conjugates[j][i][k+1];
                od;
            else
                rhs := matrices[j];
            fi;
            if conj <> rhs then
                Error( "relation ", [j,i], " not satisfied" );
            fi;
        od;
    od;
    return true;
end );

InstallMethod( ViewObj,
        "for homomorphisms into matrix groups",
        true,
        [ IsHomomorphismIntoMatrixGroup ],
        0,
function( hom )
    local  mapi,  d;

    mapi := MappingGeneratorsImages( hom );
    View( mapi[1] );
    Print( " -> " );

    d := Length(mapi[2][1]);
    Print( "<", Length(mapi[2]), " ", d, "x", d, "-matrices>" );
    return;
end );

#### Willem's code ##########################################################
##
BindGlobal( "ExtendRep", function( col, new, mats)


    # Here `col' is a from-the-left collector. Let G be the group defined
    # by `col' and H the subgroup generated by the generators with indices
    # `new+1'...`nogens'. Then we assume that the list `mats' defines a
    # representation of H (`new+1' is represented by `mats[1]' and so on.
    # This function extends the representation of the subgroup H to the
    # subgroup generated by H together with the element with index `new'.

    # Elements of `H' are represented as words. For example [ 0, 0, 1, 2 ]
    # is u_3*u_4^2. The length of these words is equal to the number of
    # polycyclic generators of G. So the first entries of such a word
    # will be zero (to be precise: the entries until and including `new').


    local   MakeWord,  EvaluateFunction,  TriangulizeRows,  nogen,
            nulev,  dim,  commutes,  i,  ev,  M,  exrep,  k,  l,
            asbas,  last_inds,  hdeg,  sp,  deg,  ready,  le,  j,  m,
            cc,  cf,  inds,  mons,  tup,  ff,  vecs,  f,  vec,  vecs1,
            B,  B1,  finished,  cls,  done,  changeocc,  vv,  m1,
            num,  isrep,  newmons;

    MakeWord:= function( p )

        # Represents the element `p' of `H' in the form [ ind, exp, ...]
        # e.g., [ 1,0,0,2] is transformed to [1,1,4,2].

        local p1,i;

         p1:=[];
         for i in [1..Length(p)] do
             if p[i]<>0 then Append(p1,[i,p[i]]); fi;
         od;
         return p1;
    end;


    EvaluateFunction:=function( f, a )

        # Here `f' is an element of the dual space of the group algebra of `H'.
        # So it is represented as [ [ i, j ], k ].
        # `a' is a monomial in the group algebra of `H'. We evaluate
        # `f(a)'.

        local p, wd, j, of, M;

        p:= a;

        if f[2] <> 0 then

           # we calculate new^{-f[2]}*a*new^{f[2]}:

           p:= List( [1..NumberOfGenerators( col )], x -> 0 );
           wd:= [ new, -f[2] ];
           for j in [1..Length(a)] do
               if a[j]<>0 then Append( wd, [ j, a[j] ] ); fi;
           od;
           Append( wd, [ new,f[2] ] );
           CollectWordOrFail( col, p, wd );
        fi;

        # We calculate the matrix corresponding to `p'.

        of:= Length(mats)-NumberOfGenerators( col );
        M:= IdentityMat( Length(mats[1]) );
        for j in [1..Length(p)] do
            if of+j >= 1 then
               if p[j] <> 0 then
                  M:= M*(mats[of+j]^p[j]);
               fi;
            fi;
        od;

        # We return the correct entry of the matrix.

        return M[f[1][2]][f[1][1]];

    end;

    TriangulizeRows:= function( vecs )

       # Here `vecs' is a list of integer vectors. This function
       # returns a list of integer vectors spanning the same space
       # over the integers (i.e., the same lattice), in triangular form.
       # The algorithm is similar to the one for Hermite normal form:
       # Suppose we have already taken care of the rows 1..k-1, and suppose that
       # we are dealing with column `col' (here col >= k, and the inequality
       # may be strict because there can be zero columns that do not contribute).
       # Then we look for the minimal entry in column `col' on and below position
       # `k'. We swap the corresponding row to position `k', and subtract it
       # as many times as possible from the rows below. If this produces
       # zeros everywehere then we are happy, and move on. If not then
       # we do this again: move the minimal entry to position `k' etc.
       #
       # The output of this function also contains a second list, in bijection
       # with the rowvectors in the output. The k-th entry of this list contains
       # the position of the first nonzero entry in the k-th row.

       local col, k, i, pos, fac, cols, v, min, c;

       col:=1;
       k:=1;
       cols:= [ ];
       while k <= Length(vecs) do

          # We look for the minimal nonzero element in column `col', where we
          # run through the rows with index >= k.

          min:= 0;
          i:= k-1;

          # First we get the first nonzero entry...
          while min = 0 and i < Length( vecs ) do
             i:=i+1; min:= vecs[i][col];
          od;

          if min = 0 then
             pos:= fail;
          else

             if min < 0 then min:= -min; fi;

             pos:= i;
             while i < Length(vecs) do
                i:=i+1;
                c:= vecs[i][col];
                if c < 0 then c:= -c; fi;
                if c < min and c <> 0 then min:= c; pos:= i; fi;
             od;
          fi;

          if pos = fail then

             # there is no nonzero entry in this column, that means that it
             # will not contribute to the triangular form, we move one column.

             col:= col+1;

          else

             if pos <> k then

                # We swap rows `k' and `pos', so that the minimal value will be
                # in row `k'.

                v:= vecs[k];
                vecs[k]:= vecs[pos];
                vecs[pos]:= v;

             fi;

             # Subtract row `k' as many times as possible.

             for i in [k+1..Length(vecs)] do
                 fac:= (vecs[i][col]-
                              (vecs[i][col] mod vecs[k][col]))/vecs[k][col];
                 vecs[i]:=vecs[i]-fac*vecs[k];
             od;

             # If all entries in the column `col' below position `k' are zero,
             # then we are done. Otherwise we just go through the process again.
             if ForAll( List( [k+1..Length(vecs)], x-> vecs[x][col] ),
                                                             IsZero ) then
                Add( cols, col );
                col:=col+1; k:=k+1;
             fi;

             # Get rid of zero rows...
             vecs:= Filtered( vecs, x -> x <> 0*x );

          fi;

       od;

       return [vecs,cols];

    end;

    nogen:= NumberOfGenerators( col );
    nulev:= List([1..nogen],x->0);
    dim:= Length( mats[1] );

    # We check whether the generator with index `new' commutes with all
    # elements of `H'. In that case we can easily extend the representation.

    commutes:= true;
    for i in [new+1..nogen] do
      ev:= ShallowCopy( nulev );
      CollectWordOrFail( col, ev, [i,-1,new,-1,i,1,new,1] );
      if ev<>0*ev then commutes:= false; break; fi;
    od;

    if commutes then

        # We represent the generator with index `new' by the matrix
        #
        #      / I 0      \
        #      \ 0 E_{12} /
        #
        # where I is the dim x dim identity matrix, and E_{12}
        # is the 2x2 matrix with ones on the diagonal, and a one on pos. (1,2).
        # A generator with index `new+i' is represented by the matrix
        #
        #     / mats[i] 0   \
        #     \   0     I_2 /
        #
        # where I_2 is the 2x2 identity matrix.

      M:= IdentityMat( dim+2 );
      M[dim+1][dim+2]:=1;
      exrep:= [ M ];
      for i in [1..Length(mats)] do
        M:= NullMat( dim+2, dim+2 );
        for k in [1..dim] do
          for l in [1..dim] do
             M[k][l]:=mats[i][k][l];
          od;
        od;
        M[dim+1][dim+1]:=1; M[dim+2][dim+2]:=1;
        Add( exrep, M );
      od;
      return exrep;
    fi;

    # In the other case we compute the space spanned by C_{\rho}. This is
    # the space spanned by the coefficient-functions on the matrix space
    # spanned by all products
    #
    #          mats[1]^k1...mats[s]^ks
    #
    # where k_i\in \Z. So first we calculate a basis of this space.

    asbas:=[ IdentityMat( dim ) ];  # The basis to be.

    last_inds:= [ 1 ];
    # `last_inds' is an array in bijection with `asbas'.
    # If `last_inds[i]=k' then the last letter in the word that defines
    # `asbas[i]' is `k'. (So we only need to multiply with higher
    # elements in order to (maybe) get new basis elements).
    # We basically loop through `asbas' and multiply each element in there
    # with elements with the same or a higher index. If all such products
    # are in `asbas' then we have found our basis. Of course, we only have
    # to try the elements added in the previous round. So `hdeg' is the
    # index where the first of thoses elements is in `asbas'.
    # `deg' will record the maximum degree of a word corresponding to a
    # basis element (for later use).

    hdeg:= 1;
    sp:= MutableBasis( Rationals, asbas );
    deg:= 0;
    ready:= false;
    while not ready do
      deg:= deg + 1;
      i:= hdeg;
      le:= Length( asbas );
      ready:= true;
      while i <= le do

        for j in [ last_inds[i]..Length( mats )] do
          m:= asbas[i]*mats[j];
          if not IsContainedInSpan( sp, m ) then
            ready:= false;
            Add( asbas, m );
            Add( last_inds, j );
            CloseMutableBasis( sp, m );
          fi;
        od;

        i:= i+1;
      od;
      hdeg:= le+1;
    od;

    deg:= deg - 1;

    # Compute the functions. A coefficient function m -> m[j][i] is represented
    # by [i,j]. `cc' will be the list of all such functions.
    # We note that the elements of `asbas' form a discriminating set for
    # the coefficient space. So we represent the coefficcients as vectors using
    # the set `asbas'.

    cc:=[ ];
    sp:= MutableBasis( Rationals, [ List(asbas,m->0)] );

    for i in [1..dim] do
      for j in [1..dim] do

        cf:= List( asbas, m -> m[j][i] );
        if not IsContainedInSpan( sp, cf ) then
          Add( cc, [i,j] );
          CloseMutableBasis( sp, cf );
        fi;
      od;
    od;


    # 'mons' will be a list of all monomials in the group H up to degree 'deg'.

    inds:=[ new+1 .. nogen ];
    mons:=[ ShallowCopy( nulev ) ];
    for i in [1..deg] do
      tup:= UnorderedTuples( inds, i );
      for j in [1..Length(tup)] do
        ev:= ShallowCopy( nulev );
        for k in tup[j] do
          ev[k]:= ev[k]+1;
        od;
        Add( mons, ev );
      od;
    od;

    # 'ff' will be a basis of the subspace of ZH^* spanned by the functions.
    # A function is either coefficient function or new^k applied to a coefficient
    # function. So we represent a function as a list [ [i,j], k].
    # 'vecs' will contain the vectorial representation of the elements of 'ff'
    # relative to the monomials in 'mons'.

    ff:=[]; vecs:=[];
    for i in [1..Length(cc)] do
      f:= [ cc[i], 0 ];
      vec:= List( mons, a -> EvaluateFunction( f, a ) );
      Add( ff, f ); Add( vecs, ShallowCopy( vec ) );
    od;


    while true do

        # We determine the module generated by C_{\rho} (as a subspace
        # of the dual of the vector space spanned by the monomials in `mons').
        # This module is generated by all new^k.f for f\in ff.
        #
        # `vecs1' will be a set of vectors spanning the module over the
        # integers, wheras `vecs' spans the module over the rationals,
        # and `vecs[i]' corresponds exactly to the function `ff[i]' (so we
        # need to keep this information as we
        # do not allow for linear combinations of functions in `ff').

        vecs1:= ShallowCopy( List( vecs, ShallowCopy ) );
        sp:= VectorSpace( Rationals, vecs );
        B:= Basis( sp, vecs );
        k:= 1;
        le:= Length( ff );
        B1:= Basis( sp, vecs1 );

        while k <= le do
            f:= List( ff[k], ShallowCopy );
            finished:= false;
            while not finished do
                f[2]:= f[2]+1; # we let `new' act by increasing
                               # `f[2]' by 1.


                if not f in ff then
                    vec:= List( mons, a -> EvaluateFunction( f, a ) );
                    cf:= Coefficients( B1, vec );
                    if cf <> fail then

                        if not ForAll( cf, IsInt ) then
                            # `vec' lies in the space `sp'
                            # but not in the space over the
                            # integers spanned by `vecs1'.
                            # So we add it, triangularize
                            # and then we get a new basis
                            # `vecs1' that spans the whole
                            # space over the integers.

                            Add( vecs1, vec );
                            vecs1:=TriangulizeRows(vecs1)[1];
                            B1:= Basis( sp, vecs1 );

                        fi;
                        finished:= true;
                        # we are finished with letting `new' act on `f'.
                    else  # we add the new vector, function etc...

                        Add( ff, List( f, ShallowCopy ) );
                        Add( vecs, ShallowCopy( vec ) );
                        Add( vecs1, ShallowCopy( vec ) );
                        sp:= VectorSpace( Rationals, vecs );
                        B1:= Basis( sp, vecs1 );

                    fi;
                else
                    finished:= true;
                fi;

            od;
            k:= k+1;
        od;

        # Now we determine a list of monomials sufficient to "distinguish" the
        # functions. It is of length equal to the dimension of the module.
        # It consists of the monomials corresponding to the columns that
        # come from a call to TriangularizeRows.

        cls:= TriangulizeRows( vecs1 )[2];
        mons:= mons{cls};
        vecs:= List( vecs, u -> u{cls} );
        vecs1:= List( vecs1, u -> u{cls} );

        # We calculate the action of the generators of the group, starting with
        # the new element. `exrep' will contain the matrices of the action.

        # It is possible that the Z-module spanned by `vecs1' is not closed
        # under the action of `new', *over Z*, i.e., that `new.f' is not
        # a Z-linear combination of the elements of `vecs1'. In that case we
        # add the vector we get, triangularize and start again. For that
        # we need the rather complicated loop `while not done do..' etc.

        sp:= VectorSpace( Rationals, vecs );
        B:= Basis( sp, vecs );
        done:= false;

        while not done do

            changeocc:= false;
            B1:= Basis( sp, vecs1 );
            exrep:= [ ];
            M:= [ ];
            for j in [1..Length(vecs1)] do
                vv:= [ ];
                for m in mons do

                    # `ev' will be the element new^{-1}.m.new

                    m1:= [new,-1];
                    Append( m1, MakeWord(m) ); Append( m1, [new,1] );
                    ev:= ShallowCopy( nulev );
                    CollectWordOrFail( col, ev, m1 );

                # Now we calculate the vector corresponding to the function
                # new.f, where f is the function corresponding to the vector
                # vecs1[j]. Now this vector is a linear combination of vectors
                # in `vecs1', i.e., the function is a linear combination of
                # elementary functions (i.e., fcts of the form new^k.c_{ij}).
                # So when evaluating we have to loop over the elements of this
                # linear combination.

                    cf:= Coefficients( B, vecs1[j] );
                    num:= 0;
                    for i in [1..Length(cf)] do
                        if cf[i]<>0 then
                            num:= num +
                                  cf[i]*EvaluateFunction( ff[i], ev );
                        fi;
                    od;
                    Add(vv,num);
                od;

                cf:= Coefficients( B1, vv );
                if not ForAll( cf, x -> IsInt( x ) ) then
                    Add( vecs1, vv );
                    vecs1:= TriangulizeRows(vecs1)[1];
                    changeocc:= true;
                    break;
                else
                    Add( M, Coefficients( B1, vv ) );
                fi;

            od;

            if not changeocc then
                Add( exrep, TransposedMat( M ) );

        # We calculate the action of the "old" generators. Basically works the
        # same as the code for the "new" generator.

                for i in [new+1..nogen] do
                    if changeocc then break; fi;
                    M:= [ ];
                    for j in [1..Length(vecs1)] do
                        vv:= [ ];
                        cf:= Coefficients( B, vecs1[j] );
                        for m in mons do
                            m1:= MakeWord( m ); Append( m1, [i,1] );
                            ev:= ShallowCopy( nulev );
                            CollectWordOrFail( col, ev, m1 );
                            num:= 0;
                            for l in [1..Length(cf)] do
                                if cf[l]<>0 then
                                    num:= num + cf[l]*
                                          EvaluateFunction( ff[l], ev );
                                fi;
                            od;
                            Add(vv,num);
                        od;

                        cf:= Coefficients( B1, vv );
                        if not ForAll( cf, x -> IsInt( x ) ) then
                            Add( vecs1, vv );
                            vecs1:= TriangulizeRows(vecs1)[1];
                            changeocc:= true;
                            break;
                        else
                            Add( M, Coefficients( B1, vv ) );
                        fi;

                    od;

                    Add( exrep, TransposedMat( M ) );

                od;
            fi;

            if not changeocc then done:= true; fi;

        od;


        # If the representation we get is a group representation, then we are
        # happy, if not then we increase the degree.

        isrep:= true;

        for i in [new..nogen] do
            if not isrep then break; fi;
            for j in [i+1..nogen] do

                # We calculate `ev' such that
                # u_j*u_i = u_i*u_j*ev, and we check whether the matrices
                # satisfy this relation.
                ev:= List( [1..nogen], x -> 0 );
                CollectWordOrFail( col, ev, [j,-1,i,-1,j,1,i,1] );
                M:= exrep[1]^0;
                for k in [new..Length(ev)] do
                    if ev[k] <> 0 then
                       M:= M*( exrep[k-new+1]^ev[k] );
                    fi;
                od;
                M:= exrep[j-new+1]*M; M:= exrep[i-new+1]*M;
                if M <> exrep[j-new+1]*exrep[i-new+1] then
                   isrep:= false; break;
                fi;
            od;
        od;

        if not isrep then

           # We increase the degree and compute our new guess for a
           # discriminating set.
           deg:=deg+1;
           newmons:= [];
           tup:= UnorderedTuples( inds, deg );
           for j in [1..Length(tup)] do
               ev:= ShallowCopy( nulev );
               for k in [1..Length(tup[j])] do
                   ev[tup[j][k]]:= ev[tup[j][k]]+1;
               od;
               Add( newmons, ShallowCopy(ev) );
           od;

           for i in [1..Length(vecs)] do
               Append( vecs[i], List( newmons, w
                            -> EvaluateFunction( ff[i], w ) ) );
           od;

           Append( mons, newmons );

        else
           return exrep;
        fi;

    od;  # end of big loop `while true ..etc'

end );


BindGlobal( "RepresentationForPcpCollector", function( col )

  local n,m,mats,i;

  n:= NumberOfGenerators( col );
  m:=IdentityMat(2);
  m[1][2]:=1;
  mats:=[m];
  for i in [2..n] do
    mats:= ExtendRep( col, n-i+1, mats );
  od;
  return mats;
end );

##
##
#### End of Willem's code ###################################################

InstallMethod( UnitriangularMatrixRepresentation,
        "for torsion free fin. gen. nilpotent pcp-groups",
        true,
        [ IsPcpGroup and IsNilpotentGroup ],
        0,
function( tgroup )
    local  coll,  mats,  mgroup,  phi;

    ## Does the group have power relations?
    if not IsTorsionFree( tgroup ) then
        Error("there are power relations in the collector of the pcp-group");
        ## Here we could compute the upper central series and construct an
        ## isomorphism to a group defined along the upper central series.
    fi;

    coll := Collector( tgroup );
    mats := LowerUnitriangularForm( RepresentationForPcpCollector( coll ) );

    mgroup := Group( mats, mats[1]^0 );
    UseIsomorphismRelation( tgroup, mgroup );

    phi := GroupHomomorphismByImagesNC( tgroup, mgroup,
                   GeneratorsOfGroup(tgroup),
                   GeneratorsOfGroup(mgroup) );
    SetIsBijective( phi, true );
    SetIsHomomorphismIntoMatrixGroup( phi, true );
    # FIXME: IsHomomorphismIntoMatrixGroup should perhaps be
    # a plain filter not a property. Especially since no methods
    # for it are installed.

    return phi;
end );