Quelle  meataxe.gi   Sprache: unbekannt

#############################################################################
##
##  This file is part of GAP, a system for computational discrete algebra.
##  This file's authors include Derek Holt, Sarah Rees, Alexander Hulpke.
##
##  Copyright of GAP belongs to its developers, whose names are too numerous
##  to list here. Please refer to the COPYRIGHT file for details.
##
##  SPDX-License-Identifier: GPL-2.0-or-later
##
##  This file contains the 'Smash'-MeatAxe modified for GAP4 and using the
##  standard MeatAxe interface.  It defines the MeatAxe SMTX.
##

InstallGlobalFunction(GModuleByMats,function(arg)
local l,f,dim,m;
  if Length(arg)<>2 and Length(arg)<>3 then
    Error("Usage: GModuleByMats(<mats>,[<dim>,]<field>)");
  fi;
  l:=arg[1];
  f:=arg[Length(arg)];
  if Length(l)>0 and Characteristic(l[1])<>Characteristic(f) then
      Error("matrices and field do not fit together");
  fi;
  l:=List(l,i->ImmutableMatrix(f,i));

  if ForAny(l,i->NrRows(i)<>NrCols(i)) or
    Length(Set(l,NrRows))>1 then
    Error("<l> must be a list of square matrices of the same dimension");
  fi;
  m:=rec(field:=f,
         isMTXModule:=true);
  if Length(l)>0 then
    dim:=NrRows(l[1]);
    if Length(arg) = 3 then
      if dim <> arg[2] then
        Error("matrices and dim do not fit together");
      fi;
    fi;
  elif Length(arg)=2 then
    Error("if no generators are given the dimension must be given explicitly");
  else
    dim:=arg[2];
    l:=[ ImmutableMatrix(f, IdentityMat(dim,f) ) ];
    m.smashMeataxe:=rec(isZeroGens:=true);
  fi;
  m.dimension:=dim;
  m.generators:=MakeImmutable(l);
  m.IsOverFiniteField:= Size(f)<>infinity and IsFFECollCollColl(l);
  return m;
end);

# variant of Value: if we evaluate the polynomial `f` at a matrix `x`, then it
# is usually beneficial to first factor `f` and evaluate at the factors
BindGlobal("SMTX_Value",function(f,x,one)
local fa;
  fa:=Factors(f);
  if Length(fa)>1 then
    return Product(List(Collected(fa),y->Value(y[1],x,one)^y[2]),one);
  else
    return Value(f,x,one);
  fi;
end);

#############################################################################
##
#F  TrivialGModule( g, F ) . . . trivial G-module
##
##  g is a finite group, F a field, trivial smash G-module computed.
InstallGlobalFunction(TrivialGModule,function(g, F)
local mats;
  mats:=List(GeneratorsOfGroup(g),i->[[One(F)]]);
  return GModuleByMats(mats,F);
end);

#############################################################################
##
#F  InducedGModule( g, h, m ) . . . calculate an induced G-module
##
## h should be a subgroup of a finite group g, and m a smash
## GModule for h.
## The induced module for g is calculated.
InstallGlobalFunction(InducedGModule,function(g, h, m)
   local  gensh, mats, ghom, gdim, hdim, F, index, gen, genim,
         gensim, r, i, j, k, l, elt, im;

   if IsGroup(g) = false then
      return Error("First argument is not a group.");
   fi;
   if SMTX.IsMTXModule(m) = false then
      return Error("Second argument is not a meataxe module.");
   fi;

   gensh:=GeneratorsOfGroup(h);
   mats:=SMTX.Generators(m);
   if Length(gensh) <> Length(mats) then
      Error("m does not have same number of generators as h = G1");
   fi;

   hdim:=SMTX.Dimension(m);
   F:=SMTX.Field(m);
   if Characteristic(F)=0 then
       ghom:=GroupHomomorphismByImagesNC(h,Group(mats),gensh,mats);
   else
       ghom:=GroupHomomorphismByImages(h,GL(hdim,F),gensh,mats);
   fi;

   # set up transveral
   r:=RightTransversal(g, h);
   index:=Length(r);

   gdim:=index*hdim;

   # Now calculate images of generators.
   gensim:=[];
   for gen in GeneratorsOfGroup(g) do
      genim:=NullMat(gdim, gdim, F);
      for i in [1..index] do
         j:=PositionCanonical(r, r[i]*gen);
         elt:=r[i]*gen/r[j];
         im:=Image(ghom, elt);
         # Now insert hdim x hdim matrix im in the correct place in the genim.
         for k in [1..hdim] do
            for l in [1..hdim] do
               genim[ (i-1)*hdim+k][ (j-1)*hdim+l]:=im[k][l];
            od;
         od;
      od;
      Add(gensim, genim);
   od;

   return GModuleByMats(gensim, F);

end);

#############################################################################
##
#F  NaturalGModule ( g[, F] )
##
## g is a matrix group, F a field.
## The corresponding natural module is output.
InstallGlobalFunction(NaturalGModule,function(group, field...)
  if not IsMatrixGroup(group) then
    Error("<group> must be a matrix group");
  fi;
  if Length(field) = 0 then
    field := DefaultFieldOfMatrixGroup(group);
  elif Length(field) = 1 then
    field := field[1];
  else
    Error("too many arguments");
  fi;
  return GModuleByMats(GeneratorsOfGroup(group), DimensionOfMatrixGroup(group), field);
end);

#############################################################################
##
#F PermutationGModule( g, F) . permutation module
##
## g is a permutation group, F a field.
## The corresponding permutation module is output.
InstallGlobalFunction(PermutationGModule,function(g, F)
   local gens, deg;
   gens:=GeneratorsOfGroup(g);
   deg:=LargestMovedPoint(gens);
   return GModuleByMats(List(gens,g->PermutationMat(g,deg,F)),F);
end);

###############################################################################
##
#F  TensorProductGModule( m1, m2 )  . . tensor product of two G-modules
##
## TensorProductGModule calculates the tensor product of smash
## modules m1 and m2.
## They are assumed to be modules over the same algebra so, in particular,
## they  should have the same number of generators.
##
InstallGlobalFunction(TensorProductGModule,function( m1, m2)
   local mat1, mat2, F1, F2,  gens, i, l;

   mat1:=SMTX.Generators(m1); mat2:=SMTX.Generators(m2);
   F1:=SMTX.Field(m1); F2:=SMTX.Field(m2);
   if F1 <> F2 then
      Error("GModules are defined over different fields.\n");
   fi;
   l:=Length(mat1);
   if l <> Length(mat2) then
      Error("GModules have different numbers of generators.");
   fi;

   gens:=[];
   for i in [1..l] do
      gens[i]:=KroneckerProduct(mat1[i], mat2[i]);
   od;

   return GModuleByMats(gens, F1);
end);

###############################################################################
##
#F  WedgeGModule( module ) . . . . . wedge product of a G-module
##
## WedgeGModule calculates the wedge product of a G-module.
## That is the action on antisymmetrix tensors.
##
InstallGlobalFunction(WedgeGModule,function( module)
   local mats, mat, newmat, row, F, gens, dim, nmats, i, j, k, m, n, x;

   mats:=SMTX.Generators(module);
   F:=SMTX.Field(module);
   nmats:=Length(mats);
   dim:=SMTX.Dimension(module);

   gens:=[];
   for i in [1..nmats] do
      mat:=mats[i];
      newmat:=[];
      for j in [1..dim] do
         for k in [1..j - 1] do
            row:=[];
            for m in [1..dim] do
               for n in [1..m - 1] do
                  x:=mat[j,m] * mat[k,n] - mat[j,n] * mat[k,m];
                  Add(row, x);
               od;
            od;
            Add(newmat, row);
         od;
      od;
      Add(gens, newmat);
   od;

   return GModuleByMats(gens, F);
end);

SMTX.Setter:=function(string)
  MakeImmutable(string);
  return function(module,obj)
    if not IsBound(module.smashMeataxe) then
      module.smashMeataxe:=rec();
    fi;
    module.smashMeataxe.(string):=obj;
  end;
end;

SMTX.IsMTXModule:=function(module)
  return IsBound(module.isMTXModule) and
         IsBound(module.field) and
         IsBound(module.generators) and
         IsBound(module.dimension);
end;

SMTX.IsZeroGens:=function(module)
  return IsBound(module.smashMeataxe)
     and IsBound(module.smashMeataxe.isZeroGens)
     and module.smashMeataxe.isZeroGens=true;
end;

SMTX.Dimension:=function(module)
  return module.dimension;
end;

SMTX.Field:=function(module)
  return module.field;
end;

SMTX.Generators:=function(module)
  if SMTX.IsZeroGens(module) then
    return [];
  else
    return module.generators;
  fi;
end;

SMTX.SetIsIrreducible:=function(module,b)
  module.IsIrreducible:=b;
end;

SMTX.HasIsIrreducible:=function(module)
  return IsBound(module.IsIrreducible);
end;


SMTX.IsAbsolutelyIrreducible:=function(module)
  if not IsBound(module.IsAbsolutelyIrreducible) then
    if not SMTX.IsIrreducible(module) then
      return false;
    fi;
    module.IsAbsolutelyIrreducible:=SMTX.AbsoluteIrreducibilityTest(module);
  fi;
  return module.IsAbsolutelyIrreducible;
end;

SMTX.SetIsAbsolutelyIrreducible:=function(module,b)
  module.IsAbsolutelyIrreducible:=b;
end;

SMTX.HasIsAbsolutelyIrreducible:=function(module)
  return IsBound(module.IsAbsolutelyIrreducible);
end;

SMTX.SetSmashRecord:=SMTX.Setter("dummy");
SMTX.Subbasis:=SMTX.Getter("subbasis");
SMTX.SetSubbasis:=SMTX.Setter("subbasis");
SMTX.AlgEl:=SMTX.Getter("algebraElement");
SMTX.SetAlgEl:=SMTX.Setter("algebraElement");
SMTX.AlgElMat:=SMTX.Getter("algebraElementMatrix");
SMTX.SetAlgElMat:=SMTX.Setter("algebraElementMatrix");
SMTX.AlgElCharPol:=SMTX.Getter("characteristicPolynomial");
SMTX.SetAlgElCharPol:=SMTX.Setter("characteristicPolynomial");
SMTX.AlgElCharPolFac:=SMTX.Getter("charpolFactors");
SMTX.SetAlgElCharPolFac:=SMTX.Setter("charpolFactors");
SMTX.AlgElNullspaceVec:=SMTX.Getter("nullspaceVector");
SMTX.SetAlgElNullspaceVec:=SMTX.Setter("nullspaceVector");
SMTX.AlgElNullspaceDimension:=SMTX.Getter("ndimFlag");
SMTX.SetAlgElNullspaceDimension:=SMTX.Setter("ndimFlag");

SMTX.CentMat:=SMTX.Getter("centMat");
SMTX.SetCentMat:=SMTX.Setter("centMat");
SMTX.CentMatMinPoly:=SMTX.Getter("centMatMinPoly");
SMTX.SetCentMatMinPoly:=SMTX.Setter("centMatMinPoly");

SMTX.FGCentMat:=SMTX.Getter("fieldGenCentMat");
SMTX.SetFGCentMat:=SMTX.Setter("fieldGenCentMat");
SMTX.FGCentMatMinPoly:=SMTX.Getter("fieldGenCentMatMinPoly");
SMTX.SetFGCentMatMinPoly:=SMTX.Setter("fieldGenCentMatMinPoly");

SMTX.SetDegreeFieldExt:=SMTX.Setter("degreeFieldExt");


#############################################################################
##
#F  SMTX.OrthogonalVector( subbasis ) single vector othogonal to a submodule,
##  N.B. subbasis is assumed to consist of normed vectors,
##  submodule is assumed proper.
##
SMTX.OrthogonalVector:=function( subbasis )
   local zero, one, v, i, j, k, x, dim, len;
   subbasis:=ShallowCopy(subbasis);
   Sort(subbasis);
   subbasis:=Reversed(subbasis);
   # Now subbasis is in order so that the vector whose leading coefficient
   # comes furthest to the left comes first.
   len:=Length(subbasis);
   dim:=Length(subbasis[1]);
   i:= 1;
   v:=[];
   one:=One(subbasis[1][1]);
   zero:=Zero(one);
   for i in [1..dim] do
      v[i]:=zero;
   od;
   i:=1;
   while i <= len and subbasis[i][i] = one do
      i:= i + 1;
   od;
   v[i]:=one;
   for j in Reversed([1..i-1]) do
      x:=zero;
      for k in [j + 1..i] do
         x:=x + v[k] * subbasis[j][k];
      od;
      v[j]:=-x;
   od;

   return v;
end;

BindGlobal( "SubGModLeadPos", function(sub)
local leadpos;
   ## As in SpinnedBasis, leadpos[i] gives the position of the first nonzero
   ## entry (which will always be 1) of sub[i].
   leadpos:=List(sub, PositionNonZero);
   if not IsDuplicateFree(leadpos) then
      Error("Subbasis isn't normed.");
   fi;
  return leadpos;
end );

#############################################################################
##
#F  SpinnedBasis( v, matrices, F, [ngens] ) . . . .
##
## The first argument v  can either be a vector over the module on
## which matrices act or a subspace.
##
## SpinnedBasis computes a basis for the submodule defined by the action of the
## matrix group generated by the list matrices on v.
## F is the field over which we act.
## It is returned as a list of normed vectors.
## If the optional fourth argument is present, then only the first ngens
## matrices in the list are used.
SMTX.SpinnedBasis:=function( v, matrices, F, ngens... )
   local   zero, ldim, step, ans, dim, subdim, leadpos, w, i, j, k, l, m;

   if Length(ngens) > 1 then
      Error("Usage:  SpinnedBasis( v, matrices, F, [ngens] )");
   fi;
   if Length(ngens) = 1 then
      ngens:=ngens[1];
      if ngens <= 0 or ngens > Length(matrices) then
         ngens:=Length(matrices);
      fi;
   else
      ngens:=Length(matrices);
   fi;
   ans:=[];
   zero:=ZeroOfBaseDomain(matrices[1]);
   if IsList(v) and Length(v)=0 then
     return [];
   elif IsMatrix(v) then
     v:= TriangulizedMat(v);
     ans:=Filtered(v,x->not IsZero(x));
   elif IsList(v) and IsVectorObj(v[1]) then
     v:=TriangulizedMat(Matrix(F,v));
     ans:=Filtered(List(v),x->not IsZero(x));
   else
     # single vector (as vector or list)
     ans:=[v];
   fi;


   #ans:=ShallowCopy(Basis(VectorSpace(F,v)));

   ans:=Filtered(ans,x->not IsZero(x));
   if Length(ans)=0 then return ans; fi;
   ans:=List(ans, v -> ImmutableVector(F,v));
   dim:=Length(ans[1]);
   subdim:=Length(ans);
   ldim:=subdim;
   step:=10;
   leadpos:=SubGModLeadPos(ans);
   for i in [1..Length(ans)] do
     w:=ans[i];
     j:=w[leadpos[i]];
     if not IsOne(j) then
       ans[i]:=j^-1*ans[i];
     fi;
   od;

   i:=1;
   while i <= subdim do
      for l in [1..ngens] do
         m:=matrices[l];
         # apply generator m to submodule generator i
         w:=ShallowCopy(ans[i] * m);
         # try to express w in terms of existing submodule generators
         for  j in [1..subdim] do
            k:=w[leadpos[j]];
            if k <> zero then
               AddRowVector(w,ans[j],-k);
            fi;
         od;

         j := PositionNonZero(w);
         if j <= dim then
            # we have found a new generator of the submodule
            subdim:=subdim + 1;
            leadpos[subdim]:=j;
            MultVector(w,w[j]^-1);
            Add( ans, w );
            if subdim = dim then
               ans:=ImmutableMatrix(F,ans);
               return ans;
            fi;
            if subdim-ldim>step then
              Info(InfoMeatAxe,4,"subdimension ",subdim);
              ldim:=subdim;
              if ldim>10*step then step:=step*3;fi;
            fi;
         fi;
      od;
      i:=i + 1;
   od;

   Sort(ans);
   ans:=Reversed(ans); # To bring it into semi-echelonised form.
   ans:=ImmutableMatrix(F,ans);
   return ans;
end;

SMTX.SubGModule:=function(module, subspace)
## The submodule of module generated by <subspace>.
  return SMTX.SpinnedBasis(subspace, SMTX.Generators(module),
                                    SMTX.Field(module));
end;

SMTX.SubmoduleGModule:=SMTX.SubGModule;

# internal helper to create a new induced/quotient/sub module over
# a given module, taking MTX.IsZeroGens into account
SMTX._NewGModule:=function(module, newgens)
   local F;
   F:=SMTX.Field(module);
   if SMTX.IsZeroGens(module) then
     return GModuleByMats([],Length(newgens[1]),F);
   else
     return GModuleByMats(newgens,F);
   fi;
end;

#############################################################################
##
#F  SMTX.SubQuotActions(module,sub,typ) . . .
##  generators of sub- and quotient-module and original module wrt new basis
##
##  IT IS ASSUMED THAT THE GENERATORS OF SUB ARE NORMED.
##
##  this function is used to compute all submodule/quotient stuff, as
##  indicated by  typ: 1=Sub, 2=Quotient, 4=Common
##  The function returns a record with components 'smatrices', 'qmatrices',
##  'nmatrices' and 'nbasis' if applicable.
##
##  See the description for 'SMTX.InducedAction' for
##  description of the matrices
##
SMTX.SubQuotActions:=function(arg)
local module, sub, typ, matrices, dim, subdim, F,
      s, c, q, leadpos, zero, zerov, smatrices, newg, im, newim, k, subi,
      qmats, smats, nmats, sr, qr, g, h, erg, i, j;

  if Length(arg) = 3 then
    # module,sub,typ
    module:=arg[1];
    sub:=arg[2];
    typ:=arg[3];
    dim:=SMTX.Dimension(module);
    F:=SMTX.Field(module);
    matrices:=module.generators;
  elif Length(arg) = 6 then
    # matrices,sub,dim,subdim,F,typ
    #
    # old variant, supported for backwards compatibility with
    # the polycyclic package and possibly private code that directly
    # calls this function
    matrices:=arg[1];
    sub:=arg[2];
    dim:=arg[3];
    subdim:=arg[4];
    F:=arg[5];
    typ:=arg[6];
    Assert(0, Length(sub) = subdim);
  else
    Error("wrong number of arguments");
  fi;

  Assert(0, ForAll(sub, x -> dim = Length(x)));

  s:=(typ mod 2)=1; # subspace indicator
  typ:=QuoInt(typ,2);
  q:=(typ mod 2)=1; # quotient indicator
  c:=typ>1; # common indicator

  zero:=Zero(F);
  leadpos:=SubGModLeadPos(sub);
  subdim:=Length(sub);

  if subdim*2<dim and not (q or c) then
    # the subspace dimension is small and we only want the subspace action:
    # performing a base change is too expensive

    zerov:=ListWithIdenticalEntries(subdim,zero);
    ConvertToVectorRep(zerov,F);

    smatrices:=[];
    for g in matrices do
      newg:=[];
      for i in [1..subdim] do
        im:=ShallowCopy(sub[i] * g);
        newim:=ShallowCopy(zerov);
        for j in [1..subdim] do
          k:=im[leadpos[j]];
          if k<> zero then
            newim[j]:=k;
            AddRowVector(im,sub[j],-k);
          fi;
        od;

        # Check that the vector is now zero - if not, then sub was
        # not the basis of a submodule
        if not IsZero(im) then return fail; fi;
        Add(newg, newim);
      od;
      Add(smatrices,ImmutableMatrix(F,newg));
    od;
    return rec(smatrices:=smatrices);
  else
    # we want the quotient or all or the subspace dimension is big enough to
    # merit a basechange

    # first extend the basis
    sub:=List(sub);
    Append(sub,List(One(matrices[1])){Difference([1..dim],leadpos)});
    sub:=ImmutableMatrix(F,sub);
    subi:=sub^-1;
    qmats:=[];
    smats:=[];
    nmats:=[];
    sr:=[1..subdim];qr:=[subdim+1..dim];
    for g in matrices do
      g:=sub*g*subi;
      if s then
        h:=ExtractSubMatrix(g,sr,sr);
        h:=ImmutableMatrix(F,h);
        Add(smats,h);
      fi;
      if q then
        h:=ExtractSubMatrix(g,qr,qr);
        h:=ImmutableMatrix(F,h);
        Add(qmats,h);
      fi;
      if c then Add(nmats,g);fi;
    od;
    erg:=rec();
    if s then
      erg.smatrices:=smats;
    fi;
    if q then
      erg.qmatrices:=qmats;
    fi;
    if c then
      erg.nmatrices:=nmats;
    fi;
    if q or c then
      erg.nbasis:=sub;
    fi;
    return erg;
  fi;
end;



#############################################################################
##
##  SMTX.NormedBasisAndBaseChange(sub)
##
##  returns a list [bas,change] where bas is a normed basis for <sub> and
##  change is the base change from bas to sub (the basis vectors of bas
##  expressed in coefficients for sub)
SMTX.NormedBasisAndBaseChange:=function(sub)
local l,m,d;
  l:=Length(sub);
  d:=Length(sub[1]);
  m:= IdentityMat(d,One(sub[1][1]));
  sub:=List([1..l],i->Concatenation(ShallowCopy(sub[i]),m[i]));
  TriangulizeMat(sub);
  m:=d+l;
  return [sub{[1..l]}{[1..d]},sub{[1..l]}{[d+1..m]}];
end;

#############################################################################
##
#F  SMTX.InducedActionSubmoduleNB( module, sub ) . . . . construct submodule
##
## module is a module record, and sub is a list of generators of a submodule.
## IT IS ASSUMED THAT THE GENERATORS OF SUB ARE NORMED.
## (i.e. each has leading coefficient 1 in a unique place).
## SMTX.InducedActionSubmoduleNB( module, sub ) computes the submodule of
## module for which sub is the basis.
## If sub does not generate a submodule then fail is returned.
SMTX.InducedActionSubmoduleNB:=function( module, sub )
   local   ans;

   if Length(sub) = 0 then
      return List(module.generators,i->[[]]);
   fi;

   ans:=SMTX.SubQuotActions(module,sub,1);

   if ans=fail then
     return fail;
   fi;

   return SMTX._NewGModule(module, ans.smatrices);
end;

# Ditto, but allowing also unnormed modules
SMTX.InducedActionSubmodule:=function(module,sub)
local nb,ans;

   nb:=SMTX.NormedBasisAndBaseChange(sub);
   sub:=nb[1];
   nb:=nb[2];

   if Length(sub) = 0 then
      return List(module.generators,i->[[]]);
   fi;

   ans:=SMTX.SubQuotActions(module,sub,1);

   if ans=fail then
     return fail;
   fi;

   # conjugate the matrices to correspond to given sub
   return SMTX._NewGModule(module, List(ans.smatrices,i->i^nb));;
end;

SMTX.ProperSubmoduleBasis:=function(module)
  if SMTX.IsIrreducible(module) then
    return fail;
  fi;
  return SMTX.Subbasis(module);
end;


#############################################################################
##
#F  SMTX.InducedActionFactorModule( module, sub [,compl] )
##
## module is a module record, and sub is a list of generators of a submodule.
## (i.e. each has leading coefficient 1 in a unique place).
## Qmodule is returned, where qmodule
## is the quotient module.
##
SMTX.InducedActionFactorModule:=function(module, sub, compl...)
local ans;

   sub:=TriangulizedMat(sub);

   if Length(sub) = SMTX.Dimension(module) then
      return List(module.generators,i->[[]]);
   fi;

   ans:=SMTX.SubQuotActions(module,sub,2);

   if ans=fail then
     return fail;
   fi;

   if Length(compl)=1 then
     # compute basechange
     sub:=Concatenation(sub,compl[1]);
     sub:=sub*Inverse(ans.nbasis);
     ans.qmatrices:=List(ans.qmatrices,i->i^sub);
   fi;

   return SMTX._NewGModule(module, ans.qmatrices);
end;

#############################################################################
##
#F  SMTX.InducedActionFactorModuleWithBasis( module, sub )
##
# FIXME: this function is never used and undocumented. Keep it or remove it?
SMTX.InducedActionFactorModuleWithBasis:=function(module,sub)
local ans, qmodule;

   sub:=TriangulizedMat(sub);

   if Length(sub) = SMTX.Dimension(module) then
      return List(module.generators,i->[[]]);
   fi;

   ans:=SMTX.SubQuotActions(module,sub,2);

   if ans=fail then
     return fail;
   fi;

   # fetch new basis
   sub:=ans.nbasis{[Length(sub)+1..module.dimension]};
   qmodule:=SMTX._NewGModule(module, ans.qmatrices);
   return [qmodule,sub];

end;

#############################################################################
##
#F  SMTX.InducedAction( module, sub, typ )
##  generators of sub- and quotient-module and original module wrt new basis
##  and new basis
##
## module is a module record, and sub is a list of generators of a submodule.
## IT IS ASSUMED THAT THE GENERATORS OF SUB ARE NORMED.
## (i.e. each has leading coefficient 1 in a unique place).
## SMTX.InducedAction computes the submodule and quotient
## and the original module with its matrices written wrt to the basis used
## to compute smodule and qmodule.
## [smodule, qmodule, nmodule] is returned,
## where smodule is the submodule and qmodule the quotient module.
## The matrices of nmodule have the form  A  0  where  A  and  B  are the
##                                        C  B
## corresponding matrices of smodule and qmodule resepctively.
## If sub is not the basis of a submodule then fail is returned.
SMTX.InducedAction:=function(module, sub, typ... )
local ans,erg;

   if Length(typ)>0 then
     typ:=typ[1];
   else
     typ:=7;
   fi;

   erg:=SMTX.SubQuotActions(module,sub,typ);

   if erg=fail then
     return fail;
   fi;

   ans:=[];

   if IsBound(erg.smatrices) then
     Add(ans,SMTX._NewGModule(module, erg.smatrices));
   fi;
   if IsBound(erg.qmatrices) then
     Add(ans,SMTX._NewGModule(module, erg.qmatrices));
   fi;
   if IsBound(erg.nmatrices) then
     Add(ans,SMTX._NewGModule(module, erg.nmatrices));
   fi;
   if IsBound(erg.nbasis) then
     Add(ans,erg.nbasis);
   fi;

   return ans;

end;

#############################################################################
##
#F  SMTX.InducedActionSubMatrixNB( mat, sub ) . . . . construct submodule
##
##  as InducedActionSubmoduleNB but for a matrix.
# FIXME: this function is never used and undocumented. Keep it or remove it?
SMTX.InducedActionSubMatrixNB:=function( mat, sub )
local module, ans;

   if Length(sub) = 0 then
      return [];
   fi;

   module:=GModuleByMats([mat], DefaultFieldOfMatrix(mat));
   ans:=SMTX.SubQuotActions(module,sub,1);

   if ans=fail then
     return fail;
   else
     return ans.smatrices[1];
   fi;

end;

# Ditto, but allowing also unnormed modules
# FIXME: this function is never used and undocumented. Keep it or remove it?
SMTX.InducedActionSubMatrix:=function(mat,sub)
local nb, module, ans;

   nb:=SMTX.NormedBasisAndBaseChange(sub);
   sub:=nb[1];
   nb:=nb[2];

   if Length(sub) = 0 then
      return [];
   fi;

   module:=GModuleByMats([mat], DefaultFieldOfMatrix(mat));
   ans:=SMTX.SubQuotActions(module,sub,1);

   if ans=fail then
     return fail;
   else
    # conjugate the matrices to correspond to given sub
     return ans.smatrices[1]^nb;
   fi;

end;

#############################################################################
##
#F  SMTX.InducedActionFactorMatrix( mat, sub [,compl] )
##
##  as InducedActionFactor, but for a matrix.
##
SMTX.InducedActionFactorMatrix:=function(mat, sub, compl...)
local module, ans;

   module:=GModuleByMats([mat], DefaultFieldOfMatrix(mat));
   sub:=TriangulizedMat(sub);

   if Length(sub) = SMTX.Dimension(module) then
      return [];
   fi;

   ans:=SMTX.SubQuotActions(module,sub,2);

   if ans=fail then
     return fail;
   fi;

   if Length(compl)=1 then
     # compute basechange
     sub:=Concatenation(sub,compl[1]);
     sub:=sub*Inverse(ans.nbasis);
     ans.qmatrices:=List(ans.qmatrices,i->i^sub);
   fi;

   return ans.qmatrices[1];
end;

SMTX.SMCoRaEl:=function(matrices,ngens,newgenlist,dim,F)
local g1,g2,coefflist,M,pol;
  g1:=Random(1, ngens);
  g2:=g1;
  while g2=g1 and ngens>1 do
     g2:=Random(1, ngens);
  od;
  ngens:=ngens + 1;
  matrices[ngens]:=matrices[g1] * matrices[g2];
  Add(newgenlist, [g1, g2]);
  # Take a random linear sum of the existing generators as new generator.
  # Record the sum in coefflist
  coefflist:=[Random(F)];
  M:=coefflist[1]*matrices[1];
  for g1 in [2..ngens] do
     g2:=Random(F);
     if IsOne(g2) then
       M:=M + matrices[g1];
     elif not IsZero(g2) then
       M:=M + g2 * matrices[g1];
     fi;
     Add(coefflist, g2);
  od;
  Info(InfoMeatAxe,3,"Evaluated random element in algebra.");
  pol:=CharacteristicPolynomialMatrixNC(F,M,1);
  return [M,coefflist,pol];
end;

# how many random elements should we try before (temporarily ) giving up?
# This number is set relatively high to minimize the chance of an unlucky
# random run in functions such as composition series computation.
SMTX.RAND_ELM_LIMIT:=5000;

#############################################################################
##
#F  SMTX.IrreduciblityTest( module ) try to reduce a module over a finite
##                                      field
##
## 27/12/2000.
## New version incorporating Ivanyos/Lux method of handling one difficult case
## for proving reducibility.
## (See G.Ivanyos and K. Lux, `Treating the exceptional cases of the meataxe',
##  Experimental Mathematics 9, 2000, 373-381.
##
## module is a module record
## IsIrreducible( ) attempts to decide whether module is irreducible.
## When it succeeds it returns true or false.
## We choose at random elements of the group algebra of the group.
## If el is such an element, we define M, p, fac, N, e and v as follows:-
## M is the matrix corresponding to el, p is its characteristic polynomial,
## fac an irreducible factor of p, N the nullspace of the matrix fac(M),
## ndim the dimension of N, and v a vector in N.
## If we can find the above such that ndim = deg(fac) then we can test
## conclusively for irreducibility. Then, in the case where irreducibility is
## proved, we store the information as fields for the module, since it may be
## useful later (e.g. to test for absolute irreducibility, equivalence with
## another module).
## These  fields are accessed by the functions
## AlgEl()(el), AlgElMat(M), AlgElCharPol(p),
## AlgElCharPolFac(fac), AlgElNullspaceDimension(ndim), and
## AlgElNullspaceVec(v).
##
## If we cannot find such a set with ndim = deg(fac) we may nonetheless prove
## reducibility  by finding a submodule. However we can never prove
## irreducibility without such a set (and hence the algorithm could run
## forever, but hopefully this will never happen!)
## Where reducibility is proved, we set the field .subbasis
## (a basis for the submodule, normed in the sense that the first non-zero
## component of each basis vector is 1, and is in a different position from
## the first non-zero component of every other basis vector).
## The test for irreducibility is based on the meataxe method (but in the
## meataxe, ndim is always very small, usually 1. The modification here is put
## in to enable the method to work over modules with large centralizing fields).
## We simply spin v. If we do not get  the whole space, we have a submodule,
## on the other hand, if we do get the whole space, we calculate the
## nullspace NT of the transpose of fac(M), spin that under the group
## generated by the transposes of the generating matrices, and thus either
## find the transpose of a submodule or conclusively prove irreducibility.
##
## This function can also be used to get a random submodule. Therefore it
## is not an end-user function but only called internally
SMTX.IrreducibilityTest:=function( module )
   local matrices, tmatrices, ngens, ans,  M, mat, g1, g2, maxdeg,
         newgenlist, coefflist, zero,
         N, NT, v, subbasis, fac, sfac, pol, orig_pol, q, dim, ndim, i,
         l, trying, deg, facno, bestfacno, F, count, R, rt0,idmat,
         pfac1, pfac2, quotRem, pfr, idemp, M2, mat2, mat3;

   rt0:=Runtime();
   #Info(InfoMeatAxe,1,"Calling MeatAxe. All times will be in milliseconds");
   if not SMTX.IsMTXModule(module) then
      Error("Argument of IsIrreducible is not a module.");
   fi;

   if not module.IsOverFiniteField then
      Error("Argument of IsIrreducible is not over a finite field.");
   fi;
   matrices:=ShallowCopy(module.generators);
   dim:=SMTX.Dimension(module);
   ngens:=Length(module.generators);
   F:=SMTX.Field(module);
   zero:=Zero(F);
   R:=PolynomialRing(F);

   # Now compute random elements M of the group algebra, calculate their
   # characteristic polynomials, factorize, and apply the irreducible factors
   # to M to get matrices with nontrivial nullspaces.
   # tmatrices will be a list of the transposed generators if required.

   tmatrices:=[];
   trying:=true;
   # trying will become false when we have an answer
   maxdeg:=1;
   newgenlist:=[];
   # Do a small amount of preprocessing to increase the generator set.
   for i in [1..1] do
      g1:=Random(1, ngens);
      g2:=g1;
      while g2=g1 and Length(matrices) > 1 do
         g2:=Random(1, ngens);
      od;
      ngens:=ngens + 1;
      matrices[ngens]:=matrices[g1] * matrices[g2];
      Add(newgenlist, [g1, g2]);
   od;
   Info(InfoMeatAxe,4,"Done preprocessing. Time = ",Runtime()-rt0,".");
   count:=0;

   # Main loop starts - choose a random element of group algebra on each pass
   while trying  do
      count:=count + 1;
      if count mod SMTX.RAND_ELM_LIMIT = 0 then
         Error("Have generated ",SMTX.RAND_ELM_LIMIT,
                "random elements and failed to prove\n",
                "or disprove irreducibility. Type return to keep trying.");
      fi;
      maxdeg:=Minimum(maxdeg * 2,dim);
      # On this pass, we only consider irreducible factors up to degree maxdeg.
      # Using higher degree factors is very time consuming, so we prefer to try
      # another element.
      # To choose random element, first add on a new generator as a product of
      # two randomly chosen unequal existing generators
      # Record the product in newgenlist.
      Info(InfoMeatAxe,3,"Choosing random element number ",count);

      M:=SMTX.SMCoRaEl(matrices,ngens,newgenlist,dim,F);

      ngens:=Length(matrices);

      coefflist:=M[2];
      pol:=M[3];
      M:=ImmutableMatrix(F,M[1]);
      idmat:=M^0;

      orig_pol:=pol;
      Info(InfoMeatAxe,3,"Evaluated characteristic polynomial. Time = ",
           Runtime()-rt0,".");
      # Now we extract the irreducible factors of pol starting with those
      # of low degree
      deg:=0;
      fac:=[];
      # The next loop is through the degrees of irreducible factors
      while DegreeOfLaurentPolynomial(pol) > 0 and deg < maxdeg and trying do
         repeat
            deg:=deg + 1;
            if deg > Int(DegreeOfLaurentPolynomial(pol) / 2) then
               fac:=[pol];
            else
               fac:=Factors(R, pol: factoroptions:=rec(onlydegs:=[deg]));
               fac:=Filtered(fac,i->DegreeOfLaurentPolynomial(i)=deg);
               Info(InfoMeatAxe,3,Length(fac)," factors of degree ",deg,
                    ", Time = ",Runtime()-rt0,".");
            fi;
         until fac <> [] or deg = maxdeg;

         if fac <> [] then
            if DegreeOfLaurentPolynomial(fac[1]) = dim then
               # In this case the char poly is irreducible, so the
               # module is irreducible.
               ans:=true;
               trying:=false;
               bestfacno:=1;
               v:=ListWithIdenticalEntries(dim,zero);
               v[1]:=One(F);
               ndim:=dim;
            fi;
            # Otherwise, first see if there is a non-repeating factor.
            # If so it will be decisive, so delete the rest of the list
            l:=Length(fac);
            facno:=1;
            while facno <= l and trying do
               if facno = l  or  fac[facno] <> fac[facno + 1] then
                  fac:=[fac[facno]]; l:=1;
               else
                  while facno < l and fac[facno] = fac[facno + 1] do
                     facno:=facno + 1;
                  od;
               fi;
               facno:=facno + 1;
            od;
            # Now we can delete repetitions from the list fac
            sfac:=Set(fac);

            if DegreeOfLaurentPolynomial(fac[1]) <> dim then
               # Now go through the factors and attempt to find a submodule
               facno:=1; l:=Length(sfac);
               while facno <= l and trying do
                  mat:=SMTX_Value(sfac[facno], M,idmat);
                  Info(InfoMeatAxe,5,"Evaluated matrix on factor. Time = ",
                       Runtime()-rt0,".");
                  N:=NullspaceMat(mat);
                  v:=N[1];
                  ndim:=Length(N);

                  Info(InfoMeatAxe,5,"Evaluated nullspace. Dimension = ",
                       ndim,". Time = ",Runtime()-rt0,".");
                  subbasis:=SMTX.SpinnedBasis(v, module.generators, F);
                  Info(InfoMeatAxe,5,"Spun up vector. Dimension = ",
                       Length(subbasis),". Time = ",Runtime()-rt0,".");
                  if Length(subbasis) < dim then
                     # Proper submodule found
                     trying:=false;
                     ans:=false;
                     SMTX.SetSubbasis(module, subbasis);
                  elif ndim = deg then
                     trying:=false;
                     # if we transpose and find no proper submodule, then the
                     # module is definitely irreducible.
                     mat:=TransposedMat(mat);
                     if Length(tmatrices)=0 then
                        tmatrices := List(module.generators, TransposedMat);
                     fi;
                     Info(InfoMeatAxe,5,"Transposed matrices. Time = ",
                          Runtime()-rt0,".");
                     NT:=NullspaceMat(mat);
                     Info(InfoMeatAxe,5,"Evaluated nullspace. Dimension = ",
                          Length(NT),". Time = ",Runtime()-rt0, ".");
                     subbasis:=SMTX.SpinnedBasis(NT[1],tmatrices,F);
                     Info(InfoMeatAxe,5,"Spun up vector. Dimension = ",
                          Length(subbasis),". Time = ",Runtime()-rt0, ".");
                     if Length(subbasis) < dim then
                        # subbasis is a basis for a submodule of the transposed
                        # module, and the orthogonal complement of this is a
                        # submodule of the original module. So we find a vector
                        # v in that, and then spin it. Of course we won't
                        # necessarily get the full orthogonal complement
                        # that way, but we'll certainly get a proper submodule.
                        v:=SMTX.OrthogonalVector(subbasis);
                        SMTX.SetSubbasis(module,
                          SMTX.SpinnedBasis(v,module.generators,F));
                        ans:=false;
                     else
                        ans:=true;
                        bestfacno:=facno;
                     fi;
                  fi;
                  if trying and deg>1 and count>2 then
                     Info(InfoMeatAxe,3,"Trying Ivanyos/Lux Method");
                     # first find the appropriate idempotent
                     pfac1:=sfac[facno];
                     pfac2:=orig_pol;
                     while true do
                       quotRem := QuotRemLaurpols(pfac2, sfac[facno], 3);
                       if not IsZero(quotRem[2]) then
                         break;
                       fi;
                       pfac2 := quotRem[1];
                       pfac1:=pfac1*sfac[facno];
                     od;
                     pfr:=GcdRepresentation(pfac1, pfac2);
                     idemp:=QuotRemLaurpols(pfr[2]*pfac2, orig_pol, 2);
                     # Now another random element in the group algebra.
                     # and a random vector in the module
                     g2:=Random(F);
                     if IsOne(g2) then
                       M2:=matrices[1];
                     else
                       M2:=g2 * matrices[1];
                     fi;
                     for g1 in [2..ngens] do
                        g2:=Random(F);
                        if IsOne(g2) then
                          M2:=M2 + matrices[g1];
                        elif not IsZero(g2) then
                          M2:=M2 + g2 * matrices[g1];
                        fi;
                     od;
                     Info(InfoMeatAxe,5,
                         "Evaluated second random element in algebra.");
                     v:=Random(FullRowSpace(F,dim));
                     mat2:=SMTX_Value(idemp, M,idmat);
                     mat3:=mat2*M2*mat2;
                     v:=v*(M*mat3 - mat3*M);
                     # This vector might lie in a proper subspace!
                     subbasis:=SMTX.SpinnedBasis(v, module.generators, F);
                     Info(InfoMeatAxe,5,"Spun up vector. Dimension = ",
                       Length(subbasis),". Time = ",Runtime()-rt0,".");
                     if Length(subbasis) < dim and Length(subbasis) <> 0  then
                       # Proper submodule found
                       trying:=false;
                       ans:=false;
                       SMTX.SetSubbasis(module, subbasis);
                    fi;
                  fi;
                  facno:=facno + 1;
               od; # going through irreducible factors of fixed degree.
               # If trying is false at this stage, then we don't have
               # an answer yet, so we have to go onto factors of the next degree.
               # Now divide p by the factors used if necessary
               if trying and deg < maxdeg then
                  for q in fac do
                     pol:=Quotient(R, pol, q);
                  od;
               fi;
            fi;           #DegreeOfLaurentPolynomial(fac[1]) <> dim
         fi;             #fac <> []
      od; # loop through degrees of irreducible factors

      # if we have not found a submodule and trying is false, then the module
      # must be irreducible.
      if trying = false and ans = true then
         SMTX.SetAlgEl(module, [newgenlist, coefflist]);
         SMTX.SetAlgElMat(module, M);
         SMTX.SetAlgElCharPol(module, orig_pol);
         SMTX.SetAlgElCharPolFac(module, sfac[bestfacno]);
         SMTX.SetAlgElNullspaceVec(module, v);
         SMTX.SetAlgElNullspaceDimension(module, ndim);
      fi;

   od;  # main loop

   Info(InfoMeatAxe,4,"Total time = ",Runtime()-rt0," milliseconds.");
   return ans;

end;


SMTX.IsIrreducible:=function(module)
  if not IsBound(module.IsIrreducible) then
    module.IsIrreducible:=SMTX.IrreducibilityTest(module);
  fi;
  return module.IsIrreducible;
end;

#############################################################################
##
#F SMTX.RandomIrreducibleSubGModule( module ) . . .
## find a basis for a random irreducible
## submodule of module, and return that basis and the submodule, with all
## the irreducibility flags set.
## Returns false if module is irreducible.
SMTX.RandomIrreducibleSubGModule:=function( module )
   local  ranSub, subbasis, submodule, subbasis2, submodule2,
   F, el, M, fac, N, i, matrices, genpair;

   if not SMTX.IsMTXModule(module) then
      return Error("Argument of RandomIrreducibleSubGModule is not a module.");
   elif SMTX.HasIsIrreducible(module) and SMTX.IsIrreducible(module) then
      return false;
   fi;

   # now call an irreducibility test that will compute a new subbasis

   i:=SMTX.IrreducibilityTest(module);

   if i then
     # we just found out it is irreducible
     SMTX.SetIsIrreducible(module,true);
     return false;
   elif not SMTX.HasIsIrreducible(module) then
     # or store reducibility
     SMTX.SetIsIrreducible(module,false);
   fi;

   subbasis:=SMTX.Subbasis(module);
   submodule:=SMTX.InducedActionSubmoduleNB(module, subbasis);
   ranSub:=SMTX.RandomIrreducibleSubGModule(submodule);
   if ranSub = false then
      # submodule has been proved irreducible in a call to this function,
      # so the flags have been set.
      return [ subbasis, submodule];
   else
      # ranSub[1] is given in terms of the basis for the submodule,
      # but we want it in terms of the basis of the original module.
      # So we multiply it by subbasis.
      # Then we need our basis to be normed.
      # this is done by triangulization
      F:=SMTX.Field(module);
      subbasis2:=ranSub[1] * subbasis;
      subbasis2:=TriangulizedMat(subbasis2);

      # But now since we've normed the basis subbasis2,
      # the matrices of the submodule ranSub[2] are given with respect to
      # the wrong basis.  So we have to recompute the submodule.
      submodule2:=SMTX.InducedActionSubmoduleNB(module, subbasis2);
      # Unfortunately, although it's clear that this submodule is
      # irreducible, we'll have to reset the flags that IsIrreducible sets.
      # AH Why can't we keep irreducibility?

      # Some will be the same # as in ranSub[2], but some are affected by
      # the base change, or at least part of it, since the flags gets
      # screwed up by the base change.
      # We need to set the following flags:
      el:=SMTX.AlgEl(ranSub[2]);
      SMTX.SetAlgEl(submodule2,el);
      SMTX.SetAlgElCharPol(submodule2,SMTX.AlgElCharPol(ranSub[2]));
      fac:=SMTX.AlgElCharPolFac(ranSub[2]);
      SMTX.SetAlgElCharPolFac(submodule2,fac);
      SMTX.SetAlgElNullspaceDimension(submodule2,
             SMTX.AlgElNullspaceDimension(ranSub[2]));

      # Only the actual algebra element and its nullspace have to be recomputed
      # This code is essentially from IsomorphismGModule
      matrices:=ShallowCopy(submodule2.generators);
      for genpair in el[1] do
         Add(matrices, matrices[genpair[1]] * matrices[genpair[2]]);
      od;
      M:=ImmutableMatrix(F,Sum([1..Length(matrices)], i-> el[2][i] * matrices[i]));
      SMTX.SetAlgElMat(submodule2,M);
      N:=NullspaceMat(SMTX_Value(fac,M,M^0));
      SMTX.SetAlgElNullspaceVec(submodule2,N[1]);
      return [subbasis2, submodule2];
   fi;

end;

#############################################################################
##
#F  SMTX.GoodElementGModule( module ) . .  find good group algebra element
##                                       in an irreducible module
##
## module is a module that is already known to be irreducible.
## GoodElementGModule finds a group algebra element with nullspace of
## minimal possible dimension. This dimension is 1 if the module is absolutely
## irreducible, and the degree of the relevant field extension otherwise.
## This is needed for testing for equivalence of modules.
SMTX.GoodElementGModule:=function( module )
local matrices, M, mat,  N, newgenlist, coefflist,
      fac, pol, oldpol,  q, deg, i, l,
      trying, dim, mindim, F, R, count, rt0, idmat;

   rt0:=Runtime();
   if not SMTX.IsMTXModule(module) or not SMTX.IsIrreducible(module) then
     ErrorNoReturn("Argument is not an irreducible module.");
   fi;
   if not SMTX.HasIsAbsolutelyIrreducible(module) then
      SMTX.IsAbsolutelyIrreducible(module);
   fi;
   if  SMTX.IsAbsolutelyIrreducible(module) then
     mindim:=1;
   else
     mindim:=SMTX.DegreeFieldExt(module);
   fi;

   if SMTX.AlgElNullspaceDimension(module) = mindim then return; fi;
   # This is the condition that we want. If it holds already, then there is
   # nothing else to do.

   dim:=SMTX.Dimension(module);
   matrices:=ShallowCopy(module.generators);
   F:=SMTX.Field(module);
   R:=PolynomialRing(F);

   # Now compute random elements el of the group algebra, calculate their
   # characteristic polynomials, factorize, and apply the irreducible factors
   # to el to get matrices with nontrivial nullspaces.

   trying:=true;
   count:=0;
   newgenlist:=[];
   while trying do
      count:=count + 1;
      if count mod SMTX.RAND_ELM_LIMIT = 0 then
         Error("Have generated ",SMTX.RAND_ELM_LIMIT,
                " random elements and failed ",
                "to find a good one. Type return to keep trying.");
      fi;
      Info(InfoMeatAxe,5,"Choosing random element number ",count,".");

      M:=SMTX.SMCoRaEl(matrices,Length(matrices),newgenlist,dim,F);

      coefflist:=M[2];
      pol:=M[3];
      M:=M[1];
      idmat:=M^0;

      Info(InfoMeatAxe,4,"Evaluated characteristic polynomial. Time = ",
           Runtime()-rt0,".");
      # That is necessary in case p is defined over a smaller field that F.
      oldpol:=pol;
      # Now we extract the irreducible factors of pol starting with those
      # of low degree
      deg:=0;
      fac:=[];
      while deg  <= mindim and trying do
         repeat
            deg:=deg + 1;
            if deg > mindim then
               fac:=[pol];
            else
               fac:=Factors(R, pol: factoroptions:=rec(onlydegs:=[deg]));
               fac:=Filtered(fac,i->DegreeOfLaurentPolynomial(i)<=deg);
               Info(InfoMeatAxe,4,Length(fac)," factors of degree ",deg,
                    ", Time = ",Runtime()-rt0,".");
            fi;
         until fac <> [];
         l:=Length(fac);
         if trying and deg <= mindim then
            i:=1;
            while i <= l and trying do
               mat:=SMTX_Value(fac[i], M,idmat);
               Info(InfoMeatAxe,5,"Evaluated matrix on factor. Time = ",
                    Runtime()-rt0,".");
               N:=NullspaceMat(mat);
               Info(InfoMeatAxe,5,"Evaluated nullspace. Dimension = ",
                    Length(N),". Time = ",Runtime()-rt0,".");
               if Length(N) = mindim then
                  trying:=false;
                  SMTX.SetAlgEl(module, [newgenlist, coefflist]);
                  SMTX.SetAlgElMat(module, M);
                  SMTX.SetAlgElCharPol(module, oldpol);
                  SMTX.SetAlgElCharPolFac(module, fac[i]);
                  SMTX.SetAlgElNullspaceVec(module, N[1]);
                  SMTX.SetAlgElNullspaceDimension(module, Length(N));
               fi;
               i:=i + 1;
            od;
         fi;

         if trying then
            for q in fac do
               pol:=Quotient(R, pol, q);
            od;
         fi;
      od;
   od;
   Info(InfoMeatAxe,5,"Total time = ",Runtime()-rt0," milliseconds.");

end;

#############################################################################
##
#F  EnlargedIrreducibleGModule(module, mat) . .add a generator to a module that
#
# 2bdef!


#############################################################################
##
#F  SMTX.FrobeniusAction(F, A, v [, basis]) . . action of matrix A on
##                                          . . Frobenius block of vector v
##
## FrobeniusAction(A, v) computes the Frobenius block of the dxd matrix A
## generated by the length - d vector v, and returns it.
## It is based on code of MinPolCoeffsMat.
## The optional fourth argument is for returning the basis for this block.
##
SMTX.FrobeniusAction:=function( fld, A, v, basis... )
local   L, d, p, M, one, zero, R, h, w, i, j, nd, ans;

   if Length(basis) = 0  then
      basis:=0;
   elif Length(basis) = 1  then
      basis:=basis[1];
   else
      return Error("usage: SMTX.FrobeniusAction(<F>, <A>, <v>, [, <basis>] )");
   fi;
   one :=OneOfBaseDomain(A);
   zero:=ZeroOfBaseDomain(A);
   d:=NrRows( A );
   M:=ListWithIdenticalEntries(NrCols(A) + 1,zero);
   M:=ImmutableVector(fld,M);
   v:=ImmutableVector(fld,v);
   A:=ImmutableMatrix(fld,A);

   # L[i] (length d) will contain a vector with head entry 1 at position i,
   # which is in the current block.
   # R[i] (length d + 1 but (d + 1) - entry always 0) is vector expressing
   # L[i] in terms of the basis of the block.
   L:=[];
   R:=[];

   # <j> - 1 gives the power of <A> we are looking at
   j:=1;

   # spin vector around and construct polynomial
   repeat

      # compute the head of <v>
      h:=PositionNonZero(v);

      # start with appropriate polynomial x^(<j> - 1)
      p:=ShallowCopy( M );
      p[j]:=one;

      # divide by known left sides
      w:=ShallowCopy( v );
      while h <= d and IsBound( L[h] ) do
         AddVector(p, R[h], -w[h]);
         AddVector(w, L[h], -w[h]);
         h:=PositionNonZero(w, h);
      od;

      # if <v> is not the zero vector try next power
      if h <= d  then
         #AH replaced Copy by ShallowCopy as only vector is used
         if basis <> 0 then basis[j]:=ShallowCopy(v); fi;

         R[h]:=p * w[h]^-1;
         L[h]:=w * w[h]^-1;
         j:=j + 1;
         v:=v * A;
      fi;
   until h > d;

   nd:=Length(p);
   while 0 < nd  and p[nd] = zero  do
      nd:=nd - 1;
   od;
   nd:=nd - 1;
   ans:=[];
   for i in [1..nd - 1] do
      ans[i]:=ListWithIdenticalEntries(nd,zero);
      ans[i][i + 1]:=one;
   od;
   ans[nd]:=-p{[1..nd]};

   return ans;
end;

#############################################################################
##
#F SMTX.CompleteBasis(matrices,basis) . complete a basis under a group action
##
##  CompleteBasis( matrices, basis ) takes the partial basis 'basis' of the
##  underlying space of the (irreducible) module defined by matrices, and
##  attempts to extend it to a complete basis which is a direct sum of
##  translates of the original subspace under group elements. It returns
##  true or false according to whether it succeeds.
##  It is called by IsAbsolutelyIrreducible()
##
SMTX.CompleteBasis:=function( matrices, basis )
local  L, d, subd, subd0, h, v, w, i, bno, gno, vno, newb, ngens;

   subd:=Length(basis);
   subd0:=subd;
   d:=Length( basis[1] );
   if d = subd then
      return true;
   fi;
   # L is list of normalized generators of the subspace spanned by basis.
   L:=[];
   ngens:=Length(matrices);

   # First find normalized generators for subspace itself.
   for i in [1..subd] do
      v:=basis[i];
      h:=PositionNonZero(v);
      w:=ShallowCopy(v);
      while h <= d and IsBound( L[h] )  do
         AddVector(w, L[h], -w[h]);
         h:=PositionNonZero(w, h);
      od;
      if h <= d then
         L[h]:=w * w[h]^-1;
      else
         return Error("Initial vectors are not linearly independent.");
      fi;
   od;

   # Now start translating
   bno:=1; gno:=1; vno:=1;
   while subd < d do
      # translate vector vno of block bno by generator gno
      v:= basis[ (bno - 1) * subd0 + vno] * matrices[gno];
      h:=PositionNonZero(v);
      w:=ShallowCopy(v);
      while h <= d and IsBound( L[h] )  do
         AddVector(w, L[h], -w[h]);
         h:=PositionNonZero(w, h);
      od;
      if h <= d then
         # new generator (and block)
         if vno = 1 then
            newb:=true;
         elif newb = false then
            return false;
         fi;
         L[h]:=w * w[h]^-1;
         subd:=subd + 1;
         basis[subd]:=v;
      else
         # in existing subspace
         if vno = 1 then
            newb:=false;
         elif newb = true then
            return false;
         fi;
      fi;
      vno:=vno + 1;
      if vno > subd0 then
         vno:=1;
         gno:=gno + 1;
         if gno > ngens then
            gno:=1;
            bno:=bno + 1;
         fi;
      fi;
   od;

   return true;
end;

#############################################################################
##
#F SMTX.AbsoluteIrreducibilityTest( module ) . . decide if an irreducible
##                    module over a  finite field is absolutely irreducible
##
## this function does the work for an absolute irreducibility test but does
## not actually set the flags.
## The function calculates the centralizer of the module.
## The centralizer should be isomorphic to the multiplicative
## group of the field GF(q^e) for some e, or rather to the group of
## dim/e x dim/e scalar matrices over GF(q^e), or equivalently,
## dim x dim matrices composed of identical e x e blocks along the diagonal.
##  e = 1 <=> the module is absolutely irreducible.
## The .fieldExtDeg component is set to e during the function call.
## The function shouldn't be called if the module has not already been
## shown to be irreducible, using IsIrreducible.
##
SMTX.AbsoluteIrreducibilityTest:=function( module )
local dim, ndim, gcd, div, e, ct, F, q, ok,
      M, v, M0, v0, C, C0, centmat, one, zero,
      pow, matrices, newmatrices, looking,
      basisN, basisB, basisBN, P, Pinv, i, offset, nblocks;

   if not SMTX.IsMTXModule(module) then
      Error("Argument of IsAbsoluteIrreducible is not a module");
   fi;

   if not SMTX.IsIrreducible(module) then
      Error("Argument of iIsAbsoluteIrreducible s not an irreducible module");
   fi;

   if not module.IsOverFiniteField then
      return Error("Argument of IsAbsoluteIrreducible is not over a finite field.");
   fi;
   dim:=SMTX.Dimension(module);
   F:=SMTX.Field(module);
   q:=Size(F);
   matrices:=module.generators;

   # M acts irreducibly on N, which is canonically defined with respect to M
   # as the nullspace of fac(M), where fac is a factor of the char poly of M.
   # ndim is the dimension of N, and v is a vector of N. All these come from
   # the irreducibility test for the module.
   # An element of the centralizer must centralize every element, and
   # therefore M, and so must preserve N, since N is canonically defined
   # wrt M. Our plan is therefore first to find an element which centralizes
   # the restriction of M to N, and then extend it to the whole space.

   M:=SMTX.AlgElMat(module);
   ndim:=SMTX.AlgElNullspaceDimension(module);
   v:=SMTX.AlgElNullspaceVec(module);

   # e will have to divide both dim and ndim, and hence their gcd.
   gcd:=GcdInt(dim, ndim);
   Info(InfoMeatAxe,4,"GCD of module and nullspace dimensions = ", gcd, ".");
   if gcd = 1 then
      SMTX.SetDegreeFieldExt(module,1);
      return true;
   fi;
   div:=DivisorsInt(gcd);

   # It's easy to find elements  in the centralizer of an element in Frobenius
   # (=rational canonical) form (centralizing elements are defined by their
   # action on the first basis element).
   # M0  is the Frobenius form for the action of M on N.
   # basisN is set by the function SMTX.FrobeniusAction to be the
   # basis v, vM, vM^2, .. for N

   basisN:=[];
   Info(InfoMeatAxe,4,
     "Calc. Frobenius action of element from group algebra on nullspace.");
   M0:=SMTX.FrobeniusAction(F,M,v,basisN);

   zero:=Zero(F);
   one:= One(F);
   v0:=ListWithIdenticalEntries(Length(M0[1]),zero);
   v0[1]:=one;
   ConvertToVectorRep(v0, F);

   # v0 is just the vector (1, 0, 0....0) of length ndim. It has nothing
   # in particular to do with M0[1], but multiplying a vector that happens to be
   # around by 0 is a good way to get a zero vector of the right length.

   # we try all possible divisors of gcd (biggest first) as possibilities for e
   # We're looking for a centralizing element with order dividing q^e - 1, and
   # blocks size e on N.
   for ct in Reversed([2..Length(div)]) do
      e:=div[ct];
      Info(InfoMeatAxe,4,"Trying dimension ",e," for centralising field.");
      # if ndim = e, M0 will do.
      if ndim > e then
         C:=M0;
         # Take the smallest power of C guaranteed to have order dividing
         # q^e - 1, and try that.
         pow:=(q^ndim - 1) / (q^e - 1);
         Info(InfoMeatAxe,4,"Looking for a suitable centralising element.");
         repeat
            # The first time through the loop C is M0, otherwise we choose C
            # at random from the centralizer of M0. Since M0 is in Frobenius
            # form any centralising element is determined by its top row
            # (which may be anything but the zero vector).

            if Length(C)=0 then
               C[1]:=[];
               repeat
                  ok:=0;
                  for i in [1..ndim] do
                     C[1][i]:=Random(F);
                     if C[1][i] <> zero then  ok:=1; fi;
                  od;
               until ok=1;
               ConvertToVectorRep(C[1], F);
               for i in [2..ndim] do C[i]:=C[i - 1] * M0; od;
               C:=ImmutableMatrix(F,C);
            fi;
            # C0 is the Frobenius form for the action of this power on one
            # of its blocks, B (all blocks have the same size). basisBN will
            # be set to be a basis for B, in terms of the elements of basisN.
            # A matrix product gives us the basis for B in terms of the
            # original basis for the module.
            basisBN:=[];
            C0:=SMTX.FrobeniusAction(F,C^pow,v0,basisBN);
            C:=[];
         until Length(C0) = e;
         Info(InfoMeatAxe,5,"Found one.");
         basisB:=List(
           ImmutableMatrix(F,basisBN) *
          ImmutableMatrix(F,basisN));
      else
         C0:=M0;
         basisB:=ShallowCopy(basisN);
      fi;
      C0:=ImmutableMatrix(F,C0);
      # Now try to extend basisB to a basis for the whole module, by
      # translating it by the generating matrices.
      Info(InfoMeatAxe,4,"Trying to extend basis to whole module.");
      if SMTX.CompleteBasis(matrices,basisB) then
         # We succeeded in extending the basis (might not have done).
         # So now we have a full basis, which we think of now as a base
         # change matrix.
         Info(InfoMeatAxe,4,"Succeeded. Calculating centralising matrix.");
         newmatrices:=[];
         P:=ImmutableMatrix(F,basisB);
         Pinv:=P^-1;
         for i in [1..Length(matrices)] do
            newmatrices[i]:=P * matrices[i] * Pinv;
         od;
         # Make the sum of copies of C0 as centmat
         centmat:=NullMat(dim, dim, F);
         ConvertToMatrixRep(centmat, F);
         nblocks:=dim/e;
         Assert(0, NrRows(C0) = e); # ???
         for i in [1..nblocks] do
            offset := (i - 1) * e;
            CopySubMatrix(C0, centmat,
                          [1..e], [offset+1..offset+e],
                          [1..e], [offset+1..offset+e]);
         od;
         centmat := ImmutableMatrix(F, centmat);
         Info(InfoMeatAxe,2,"Checking that it centralises the generators.");
         # Check centralizing.
         looking:=true;
         i:=1;
         while looking and i <= Length(newmatrices) do
            if newmatrices[i] * centmat <> centmat * newmatrices[i] then
               looking:=false;
            fi;
            i:=i + 1;
         od;
         if looking then
            Info(InfoMeatAxe,2,"It did!");
            SMTX.SetDegreeFieldExt(module, e);
            SMTX.SetCentMat(module, Pinv * centmat * P); # get the base right
            # We will also record the minimal polynomial of C0 (and hence of
            # centmat) in case we need it at some future date.
            SMTX.SetCentMatMinPoly(module, MinimalPolynomial(F,C0,1));
            return false;
         fi;
         Info(InfoMeatAxe,2,"But it didn't.");
      else
         Info(InfoMeatAxe,2,"Failed!");
      fi;
   od;

   Info(InfoMeatAxe,2,
     "Tried all divisors. Must be absolutely irreducible.");
   SMTX.SetDegreeFieldExt(module, 1);
   return true;
end;

SMTX.DegreeFieldExt:=function(module)
  if not IsBound(module.smashMeataxe.degreeFieldExt) then
    SMTX.AbsoluteIrreducibilityTest( module );
  fi;
  return module.smashMeataxe.degreeFieldExt;
end;

SMTX.DegreeSplittingField:=function(module)
  return DegreeOverPrimeField(SMTX.Field(module))
         *SMTX.DegreeFieldExt(module);
end;

#############################################################################
##
#F  FieldGenCentMat( module ) . . find a centralizing matrix that generates
##                                the centralizing field of an irred. module
##
## FieldGenCentMat( ) should only be applied to modules that have already
## been proved irreducible using IsIrreducible. It then tests for absolute
## irreducibility (if not already known) and does nothing if module is
## absolutely irreducible. Otherwise, it returns a matrix that generates
## (multiplicatively) the centralizing field (i.e. its multiplicative order
## is q^e - 1, where e is the degree of the centralizing field. This is not
## yet used, but maybe in future, if we wish to reduce the group to matrices
## over the larger field.
SMTX.FieldGenCentMat:=function( module )
   local e, F, R, q, qe, minpol, pp,
         centmat, newcentmat, genpol, looking,
         okd;

  if SMTX.FGCentMat(module)=fail then
    if SMTX.IsMTXModule(module) = false then
      Error("Argument of IsIrreducible is not a module.");
    fi;

    if not SMTX.IsIrreducible(module) then
      Error("GModule is not irreducible.");
    fi;

    # enforce absirred knowledge as well.
    #if not SMTX.IsAbsolutelyIrreducible(module) then
    #  Error("GModule is not absolutely irreducible.");
    #fi;

    if SMTX.CentMat(module)=fail then
      Error("No CentMat component!");
    fi;

    F:=SMTX.Field(module);
    R:=PolynomialRing(F);
    q:=Size(F);
    e :=SMTX.DegreeFieldExt(module);
    qe:=q^e - 1;
    minpol:=SMTX.CentMatMinPoly(module);
    # Factorise q^e - 1
    pp:=PrimePowersInt(qe);
    # We seek a generator of the field of order q^e - 1. In other words, a
    # polynomial genpol of degree e, which has multiplicative order q^e - 1
    # modulo minpol. We first try the polynomial x, which is the element we
    # have already. If this does not work, then we try random nonconstant
    # polynomials until we find one with the right order.

    genpol:=Indeterminate(F);

    looking:=true;
    while looking do
      if genpol <> minpol then
      okd:=FFPOrderKnownDividend(R, genpol, minpol, pp);
      if okd[1] * Order(One(F)*okd[2]) = qe then
          looking:=false;
      fi;
      fi;
      if looking then
          repeat
            genpol:=RandomPol(F, e,1);
          until DegreeOfUnivariateLaurentPolynomial(genpol) > 0;
          genpol:=StandardAssociate(R, genpol);
      fi;
    od;
    # Finally recalculate centmat and its minimal polynomial.
    centmat:=SMTX.CentMat(module);
    newcentmat:=SMTX_Value(genpol, centmat,centmat^0);
    SMTX.SetFGCentMat(module, newcentmat);
    SMTX.SetFGCentMatMinPoly(module,MinimalPolynomialMatrixNC(F,newcentmat,1));
    # Ugh! That was very inefficient - should work out the min poly using
    # polynomials, but will sort that out if its ever needed.
  fi;
  return SMTX.FGCentMat(module);
end;

###############################################################################
##
#F  SMTX.CollectedFactors( module ) . . find composition factors of a module
##
## 01/01/01 Try to deal more efficiently with large numbers of repeated
## small factors by using SMTX.Homomorphisms
##
## SMTX.CollectedFactors calls IsIrreducible repeatedly to find the
## composition factors of the GModule `module'. It also calls
## IsomorphismGModule to determine which are isomorphic.
## It returns a list [f1, f2, ..fr], where each fi is a list [m, n],
## where m is an irreducible composition factor of module, and n is the
## number of times it occurs in module.
##
SMTX.CollectedFactors:= function( module )
  local field,dim, factors, factorsout, queue, cmod, new,
      d, i, j, l, lf, q, smod, ds, homs, mat;
   if SMTX.IsMTXModule(module) = false then
      return Error("Argument is not a module.");
   fi;

   dim:=SMTX.Dimension(module);
   field:= SMTX.Field(module);
   factors:=[];
   for i in [1..dim] do
      factors[i]:=[];
   od;
   # factors[i] will contain a list [f1, f2, ..., fr] of the composition factors
   # of module of dimension i. Each fi will have the form [m, n], where m is
   # the module, and n its multiplicity.

   queue:=[module];
   # queue is the list of modules awaiting processing.

   while Length(queue) > 0 do
      cmod:=Remove(queue);
      d:=SMTX.Dimension(cmod);
      Info(InfoMeatAxe,3,"Length of queue = ", Length(queue)+1, ", dim = ", d, ".");

      if SMTX.IsIrreducible(cmod) then
         Info(InfoMeatAxe,2,"Irreducible: ");
         # module is irreducible. See if it is already on the list.
         new:=true;
         lf:=Length(factors[d]);
         i:=1;
--> --------------------

--> maximum size reached

--> --------------------

[ zur Elbe Produktseite wechseln0.57Quellennavigators  Analyse erneut starten  ]