1#############################################################################
2##
3##  This file is part of GAP, a system for computational discrete algebra.
4##  This file's authors include Derek Holt, Sarah Rees, Alexander Hulpke.
5##
6##  Copyright of GAP belongs to its developers, whose names are too numerous
7##  to list here. Please refer to the COPYRIGHT file for details.
8##
9##  SPDX-License-Identifier: GPL-2.0-or-later
10##
11##  This file contains the 'Smash'-MeatAxe modified for GAP4 and using the
12##  standard MeatAxe interface.  It defines the MeatAxe SMTX.
13##
14
15InstallGlobalFunction(GModuleByMats,function(arg)
16local l,f,dim,m;
17  l:=arg[1];
18  if Length(arg)=1 then
19    Error("Usage: GModuleByMats(<mats>,[<id>,]<field>)");
20  fi;
21  f:=arg[Length(arg)];
22  if Length(l)>0 and Characteristic(l[1])<>Characteristic(f) then
23      Error("matrices and field do not fit together");
24  fi;
25  l:=List(l,i->ImmutableMatrix(f,i));
26
27  if ForAny(l,i->Length(i)<>Length(i[1])) or
28    Length(Set(List(l,Length)))>1 then
29    Error("<l> must be a list of square matrices of the same dimension");
30  fi;
31  m:=rec(field:=f,
32         isMTXModule:=true);
33  if Length(l)>0 then
34    dim:=Length(l[1][1]);
35  elif Length(arg)=2 then
36    Error("if no generators are given the dimension must be given explicitly");
37  else
38    dim:=arg[2];
39    l:=[ ImmutableMatrix(f, IdentityMat(dim,f) ) ];
40    m.smashMeataxe:=rec(isZeroGens:=true);
41  fi;
42  m.dimension:=dim;
43  m.generators:=MakeImmutable(l);
44  m.IsOverFiniteField:= Size(f)<>infinity and IsFFECollCollColl(l);
45  return m;
46end);
47
48#############################################################################
49##
50#F  TrivialGModule( g, F ) . . . trivial G-module
51##
52##  g is a finite group, F a field, trivial smash G-module computed.
53InstallGlobalFunction(TrivialGModule,function(g, F)
54local mats;
55  mats:=List(GeneratorsOfGroup(g),i->[[One(F)]]);
56  return GModuleByMats(mats,F);
57end);
58
59#############################################################################
60##
61#F  InducedGModule( g, h, m ) . . . calculate an induced G-module
62##
63## h should be a subgroup of a finite group g, and m a smash
64## GModule for h.
65## The induced module for g is calculated.
66InstallGlobalFunction(InducedGModule,function(g, h, m)
67   local  gensh, mats, ghom, gdim, hdim, F, index, gen, genim,
68         gensim, r, i, j, k, l, elt, im;
69
70   if IsGroup(g) = false then
71      return Error("First argument is not a group.");
72   fi;
73   if SMTX.IsMTXModule(m) = false then
74      return Error("Second argument is not a meataxe module.");
75   fi;
76
77   gensh:=GeneratorsOfGroup(h);
78   mats:=SMTX.Generators(m);
79   if Length(gensh) <> Length(mats) then
80      Error("m does not have same number of generators as h = G1");
81   fi;
82
83   hdim:=SMTX.Dimension(m);
84   F:=SMTX.Field(m);
85   if Characteristic(F)=0 then
86       ghom:=GroupHomomorphismByImagesNC(h,Group(mats),gensh,mats);
87   else
88       ghom:=GroupHomomorphismByImages(h,GL(hdim,F),gensh,mats);
89   fi;
90
91   # set up transveral
92   r:=RightTransversal(g, h);
93   index:=Length(r);
94
95   gdim:=index*hdim;
96
97   # Now calculate images of generators.
98   gensim:=[];
99   for gen in GeneratorsOfGroup(g) do
100      genim:=NullMat(gdim, gdim, F);
101      for i in [1..index] do
102         j:=PositionCanonical(r, r[i]*gen);
103         elt:=r[i]*gen/r[j];
104         im:=Image(ghom, elt);
105         # Now insert hdim x hdim matrix im in the correct place in the genim.
106         for k in [1..hdim] do
107            for l in [1..hdim] do
108               genim[ (i-1)*hdim+k][ (j-1)*hdim+l]:=im[k][l];
109            od;
110         od;
111      od;
112      Add(gensim, genim);
113   od;
114
115   return GModuleByMats(gensim, F);
116
117end);
118
119#############################################################################
120##
121#F PermutationGModule( g, F) . permutation module
122##
123## g is a permutation group, F a field.
124## The corresponding permutation module is output.
125InstallGlobalFunction(PermutationGModule,function(g, F)
126   local gens, deg;
127   gens:=GeneratorsOfGroup(g);
128   deg:=LargestMovedPoint(gens);
129   return GModuleByMats(List(gens,g->PermutationMat(g,deg,F)),F);
130end);
131
132###############################################################################
133##
134#F  TensorProductGModule( m1, m2 )  . . tensor product of two G-modules
135##
136## TensorProductGModule calculates the tensor product of smash
137## modules m1 and m2.
138## They are assumed to be modules over the same algebra so, in particular,
139## they  should have the same number of generators.
140##
141InstallGlobalFunction(TensorProductGModule,function( m1, m2)
142   local mat1, mat2, F1, F2,  gens, i, l;
143
144   mat1:=SMTX.Generators(m1); mat2:=SMTX.Generators(m2);
145   F1:=SMTX.Field(m1); F2:=SMTX.Field(m2);
146   if F1 <> F2 then
147      Error("GModules are defined over different fields.\n");
148   fi;
149   l:=Length(mat1);
150   if l <> Length(mat2) then
151      Error("GModules have different numbers of generators.");
152   fi;
153
154   gens:=[];
155   for i in [1..l] do
156      gens[i]:=KroneckerProduct(mat1[i], mat2[i]);
157   od;
158
159   return GModuleByMats(gens, F1);
160end);
161
162###############################################################################
163##
164#F  WedgeGModule( module ) . . . . . wedge product of a G-module
165##
166## WedgeGModule calculates the wedge product of a G-module.
167## That is the action on antisymmetrix tensors.
168##
169InstallGlobalFunction(WedgeGModule,function( module)
170   local mats, mat, newmat, row, F, gens, dim, nmats, i, j, k, m, n, x;
171
172   mats:=SMTX.Generators(module);
173   F:=SMTX.Field(module);
174   nmats:=Length(mats);
175   dim:=SMTX.Dimension(module);
176
177   gens:=[];
178   for i in [1..nmats] do
179      mat:=mats[i];
180      newmat:=[];
181      for j in [1..dim] do
182         for k in [1..j - 1] do
183            row:=[];
184            for m in [1..dim] do
185               for n in [1..m - 1] do
186                  x:=mat[j][m] * mat[k][n] - mat[j][n] * mat[k][m];
187                  Add(row, x);
188               od;
189            od;
190            Add(newmat, row);
191         od;
192      od;
193      Add(gens, newmat);
194   od;
195
196   return GModuleByMats(gens, F);
197end);
198
199SMTX.Setter:=function(string)
200  MakeImmutable(string);
201  return function(module,obj)
202    if not IsBound(module.smashMeataxe) then
203      module.smashMeataxe:=rec();
204    fi;
205    module.smashMeataxe.(string):=obj;
206  end;
207end;
208
209SMTX.IsMTXModule:=function(module)
210  return IsBound(module.isMTXModule) and
211         IsBound(module.field) and
212         IsBound(module.generators) and
213         IsBound(module.dimension);
214end;
215
216SMTX.IsZeroGens:=function(module)
217  return IsBound(module.smashMeataxe)
218     and IsBound(module.smashMeataxe.isZeroGens)
219     and module.smashMeataxe.isZeroGens=true;
220end;
221
222SMTX.Dimension:=function(module)
223  return module.dimension;
224end;
225
226SMTX.Field:=function(module)
227  return module.field;
228end;
229
230SMTX.Generators:=function(module)
231  if SMTX.IsZeroGens(module) then
232    return [];
233  else
234    return module.generators;
235  fi;
236end;
237
238SMTX.SetIsIrreducible:=function(module,b)
239  module.IsIrreducible:=b;
240end;
241
242SMTX.HasIsIrreducible:=function(module)
243  return IsBound(module.IsIrreducible);
244end;
245
246
247SMTX.IsAbsolutelyIrreducible:=function(module)
248  if not IsBound(module.IsAbsolutelyIrreducible) then
249    if not SMTX.IsIrreducible(module) then
250      return false;
251    fi;
252    module.IsAbsolutelyIrreducible:=SMTX.AbsoluteIrreducibilityTest(module);
253  fi;
254  return module.IsAbsolutelyIrreducible;
255end;
256
257SMTX.SetIsAbsolutelyIrreducible:=function(module,b)
258  module.IsAbsolutelyIrreducible:=b;
259end;
260
261SMTX.HasIsAbsolutelyIrreducible:=function(module)
262  return IsBound(module.IsAbsolutelyIrreducible);
263end;
264
265SMTX.SetSmashRecord:=SMTX.Setter("dummy");
266SMTX.Subbasis:=SMTX.Getter("subbasis");
267SMTX.SetSubbasis:=SMTX.Setter("subbasis");
268SMTX.AlgEl:=SMTX.Getter("algebraElement");
269SMTX.SetAlgEl:=SMTX.Setter("algebraElement");
270SMTX.AlgElMat:=SMTX.Getter("algebraElementMatrix");
271SMTX.SetAlgElMat:=SMTX.Setter("algebraElementMatrix");
272SMTX.AlgElCharPol:=SMTX.Getter("characteristicPolynomial");
273SMTX.SetAlgElCharPol:=SMTX.Setter("characteristicPolynomial");
274SMTX.AlgElCharPolFac:=SMTX.Getter("charpolFactors");
275SMTX.SetAlgElCharPolFac:=SMTX.Setter("charpolFactors");
276SMTX.AlgElNullspaceVec:=SMTX.Getter("nullspaceVector");
277SMTX.SetAlgElNullspaceVec:=SMTX.Setter("nullspaceVector");
278SMTX.AlgElNullspaceDimension:=SMTX.Getter("ndimFlag");
279SMTX.SetAlgElNullspaceDimension:=SMTX.Setter("ndimFlag");
280
281SMTX.CentMat:=SMTX.Getter("centMat");
282SMTX.SetCentMat:=SMTX.Setter("centMat");
283SMTX.CentMatMinPoly:=SMTX.Getter("centMatMinPoly");
284SMTX.SetCentMatMinPoly:=SMTX.Setter("centMatMinPoly");
285
286SMTX.FGCentMat:=SMTX.Getter("fieldGenCentMat");
287SMTX.SetFGCentMat:=SMTX.Setter("fieldGenCentMat");
288SMTX.FGCentMatMinPoly:=SMTX.Getter("fieldGenCentMatMinPoly");
289SMTX.SetFGCentMatMinPoly:=SMTX.Setter("fieldGenCentMatMinPoly");
290
291SMTX.SetDegreeFieldExt:=SMTX.Setter("degreeFieldExt");
292
293
294#############################################################################
295##
296#F  SMTX.OrthogonalVector( subbasis ) single vector othogonal to a submodule,
297##  N.B. subbasis is assumed to consist of normed vectors,
298##  submodule is assumed proper.
299##
300SMTX_OrthogonalVector:=function( subbasis )
301   local zero, one, v, i, j, k, x, dim, len;
302   subbasis:=ShallowCopy(subbasis);
303   Sort(subbasis);
304   subbasis:=Reversed(subbasis);
305   # Now subbasis is in order so that the vector whose leading coefficient
306   # comes furthest to the left comes first.
307   len:=Length(subbasis);
308   dim:=Length(subbasis[1]);
309   i:= 1;
310   v:=[];
311   one:=One(subbasis[1][1]);
312   zero:=Zero(one);
313   for i in [1..dim] do
314      v[i]:=zero;
315   od;
316   i:=1;
317   while i <= len and subbasis[i][i] = one do
318      i:= i + 1;
319   od;
320   v[i]:=one;
321   for j in Reversed([1..i-1]) do
322      x:=zero;
323      for k in [j + 1..i] do
324         x:=x + v[k] * subbasis[j][k];
325      od;
326      v[j]:=-x;
327   od;
328
329   return v;
330end;
331SMTX.OrthogonalVector:=SMTX_OrthogonalVector;
332
333SubGModLeadPos:=function(sub,dim,subdim,zero)
334local leadpos,i,j,k;
335   ## As in SpinnedBasis, leadpos[i] gives the position of the first nonzero
336   ## entry (which will always be 1) of sub[i].
337
338   leadpos:=[];
339   for i in [1..subdim] do
340      j:=1;
341      while j <= dim and sub[i][j]=zero do j:=j + 1; od;
342      leadpos[i]:=j;
343      for k in [1..i - 1] do
344         if leadpos[k] = j then
345            Error("Subbasis isn't normed.");
346         fi;
347      od;
348   od;
349  return leadpos;
350end;
351
352#############################################################################
353##
354#F  SpinnedBasis( v, matrices, F, [ngens] ) . . . .
355##
356## The first argument v  can either be a vector over the module on
357## which matrices act or a subspace.
358##
359## SpinnedBasis computes a basis for the submodule defined by the action of the
360## matrix group generated by the list matrices on v.
361## F is the field over which we act.
362## It is returned as a list of normed vectors.
363## If the optional third argument is present, then only the first ngens
364## matrices in the list are used.
365SMTX_SpinnedBasis:=function( arg  )
366   local   v, matrices, ngens, zero,
367           ans, dim, subdim, leadpos, w, i, j, k, l, m,F;
368
369   if Length(arg) < 3 or Length(arg) > 4 then
370      Error("Usage:  SpinnedBasis( v, matrices, F, [ngens] )");
371   fi;
372   v:=arg[1];
373   matrices:=arg[2];
374   F:=arg[3];
375   if Length(arg) = 4 then
376      ngens:=arg[4];
377      if ngens <= 0 or ngens > Length(matrices) then
378         ngens:=Length(matrices);
379      fi;
380   else
381      ngens:=Length(matrices);
382   fi;
383   ans:=[];
384   if Length(v)=0 then
385     return [];
386   fi;
387   if not IsList(v[1]) then
388     v:=[v];
389   fi;
390   zero:=Zero(matrices[1][1][1]);
391   ans:=ShallowCopy(Basis(VectorSpace(F,v)));
392   if Length(ans)=0 then
393     return ans;
394   fi;
395   ans:=List(ans, v -> ImmutableVector(F,v));
396   dim:=Length(ans[1]);
397   subdim:=Length(ans);
398   leadpos:=SubGModLeadPos(ans,dim,subdim,zero);
399
400   i:=1;
401   while i <= subdim do
402      for l in [1..ngens] do
403         m:=matrices[l];
404         # apply generator m to submodule generator i
405         w:=ShallowCopy(ans[i] * m);
406         # try to express w in terms of existing submodule generators
407         j:=1;
408         for  j in [1..subdim] do
409            k:=w[leadpos[j]];
410            if k <> zero then
411               #w:=w - k * ans[j];
412               AddRowVector(w,ans[j],-k);
413            fi;
414         od;
415
416         j:=1;
417         while j <= dim and w[j] = zero do j:=j + 1; od;
418         if j <= dim then
419            # we have found a new generator of the submodule
420            subdim:=subdim + 1;
421            leadpos[subdim]:=j;
422            #w:=(w[j]^-1) * w;
423            MultVector(w,w[j]^-1);
424            Add( ans, w );
425            if subdim = dim then
426               ans:=ImmutableMatrix(F,ans);
427               return ans;
428            fi;
429         fi;
430      od;
431      i:=i + 1;
432   od;
433
434   Sort(ans);
435   ans:=Reversed(ans); # To bring it into semi-echelonised form.
436   ans:=ImmutableMatrix(F,ans);
437   return ans;
438end;
439SMTX.SpinnedBasis:=SMTX_SpinnedBasis;
440
441SMTX_SubGModule:=function(module, subspace)
442## The submodule of module generated by <subspace>.
443  return SMTX.SpinnedBasis(subspace, SMTX.Generators(module),
444                                    SMTX.Field(module));
445end;
446
447SMTX.SubGModule:=SMTX_SubGModule;
448SMTX.SubmoduleGModule:=SMTX_SubGModule;
449
450#############################################################################
451##
452#F  SMTX.SubQuotActionsModule(matrices,sub,dim,subdim,field,typ) . . .
453##  generators of sub- and quotient-module and original module wrt new basis
454##
455##  IT IS ASSUMED THAT THE GENERATORS OF SUB ARE NORMED.
456##
457##  this function is used to compute all submodule/quotient stuff, as
458##  indicated by  typ: 1=Sub, 2=Quotient, 4=Common
459##  The function returns a record with components 'smatrices', 'qmatrices',
460##  'nmatrices' and 'nbasis' if applicable.
461##
462##  See the description for 'SMTX.InducedAction' for
463##  description of the matrices
464##
465SMTX_SubQuotActions:=function(matrices,sub,dim,subdim,F,typ)
466local s, c, q, leadpos, zero, zerov, smatrices, newg, im, newim, k, subi,
467      qmats, smats, nmats, sr, qr, g, h, erg, i, j;
468
469  s:=(typ mod 2)=1; # subspace indicator
470  typ:=QuoInt(typ,2);
471  q:=(typ mod 2)=1; # quotient indicator
472  c:=typ>1; # common indicator
473
474  zero:=Zero(F);
475  leadpos:=SubGModLeadPos(sub,dim,subdim,zero);
476
477  if subdim*2<dim and not (q or c) then
478    # the subspace dimension is small and we only want the subspace action:
479    # performing a base change is too expensive
480
481    zerov:=ListWithIdenticalEntries(subdim,zero);
482    ConvertToVectorRep(zerov,F);
483
484    smatrices:=[];
485    for g in matrices do
486      newg:=[];
487      for i in [1..subdim] do
488        im:=ShallowCopy(sub[i] * g);
489        newim:=ShallowCopy(zerov);
490        for j in [1..subdim] do
491          k:=im[leadpos[j]];
492          if k<> zero then
493            newim[j]:=k;
494            AddRowVector(im,sub[j],-k);
495          fi;
496        od;
497
498        # Check that the vector is now zero - if not, then sub was
499        # not the basis of a submodule
500        if im <> Zero(im) then return fail; fi;
501        Add(newg, newim);
502      od;
503      Add(smatrices,ImmutableMatrix(F,newg));
504    od;
505    return rec(smatrices:=smatrices);
506  else
507    # we want the quotient or all or the subspace dimension is big enough to
508    # merit a basechange
509
510    # first extend the basis
511    sub:=ShallowCopy(sub);
512    Append(sub,One(matrices[1]){Difference([1..dim],leadpos)});
513    sub:=ImmutableMatrix(F,sub);
514    subi:=sub^-1;
515    qmats:=[];
516    smats:=[];
517    nmats:=[];
518    sr:=[1..subdim];qr:=[subdim+1..dim];
519    for g in matrices do
520      g:=sub*g*subi;
521      if s then
522        h:=g{sr}{sr};
523        h:=ImmutableMatrix(F,h);
524        Add(smats,h);
525      fi;
526      if q then
527        h:=g{qr}{qr};
528        h:=ImmutableMatrix(F,h);
529        Add(qmats,h);
530      fi;
531      if c then Add(nmats,g);fi;
532    od;
533    erg:=rec();
534    if s then
535      erg.smatrices:=smats;
536    fi;
537    if q then
538      erg.qmatrices:=qmats;
539    fi;
540    if c then
541      erg.nmatrices:=nmats;
542    fi;
543    if q or c then
544      erg.nbasis:=sub;
545    fi;
546    return erg;
547  fi;
548end;
549
550
551SMTX.SubQuotActions:=SMTX_SubQuotActions;
552
553#############################################################################
554##
555##  SMTX.NormedBasisAndBaseChange(sub)
556##
557##  returns a list [bas,change] where bas is a normed basis for <sub> and
558##  change is the base change from bas to sub (the basis vectors of bas
559##  expressed in coefficients for sub)
560SMTX.NormedBasisAndBaseChange:=function(sub)
561local l,m,d;
562  l:=Length(sub);
563  d:=Length(sub[1]);
564  m:= IdentityMat(d,One(sub[1][1]));
565  sub:=List([1..l],i->Concatenation(ShallowCopy(sub[i]),m[i]));
566  TriangulizeMat(sub);
567  m:=d+l;
568  return [sub{[1..l]}{[1..d]},sub{[1..l]}{[d+1..m]}];
569end;
570
571#############################################################################
572##
573#F  SMTX.InducedActionSubmoduleNB( module, sub ) . . . . construct submodule
574##
575## module is a module record, and sub is a list of generators of a submodule.
576## IT IS ASSUMED THAT THE GENERATORS OF SUB ARE NORMED.
577## (i.e. each has leading coefficient 1 in a unique place).
578## SMTX.InducedActionSubmoduleNB( module, sub ) computes the submodule of
579## module for which sub is the basis.
580## If sub does not generate a submodule then fail is returned.
581SMTX.InducedActionSubmoduleNB:=function( module, sub )
582   local   ans, dim, subdim, smodule,F;
583
584   subdim:=Length(sub);
585   if subdim = 0 then
586      return List(module.generators,i->[[]]);
587   fi;
588   dim:=SMTX.Dimension(module);
589   F:=SMTX.Field(module);
590
591   ans:=SMTX.SubQuotActions(module.generators,sub,dim,subdim,F,1);
592
593   if ans=fail then
594     return fail;
595   fi;
596
597   if SMTX.IsZeroGens(module) then
598     smodule:=GModuleByMats([],Length(ans.smatrices[1]),F);
599   else
600     smodule:=GModuleByMats(ans.smatrices,F);
601   fi;
602   return smodule;
603end;
604
605# Ditto, but allowing also unnormed modules
606SMTX.InducedActionSubmodule:=function(module,sub)
607local nb,ans,dim,subdim,smodule,F;
608  nb:=SMTX.NormedBasisAndBaseChange(sub);
609  sub:=nb[1];
610  nb:=nb[2];
611
612   subdim:=Length(sub);
613   if subdim = 0 then
614      return List(module.generators,i->[[]]);
615   fi;
616   dim:=SMTX.Dimension(module);
617   F:=SMTX.Field(module);
618
619   ans:=SMTX.SubQuotActions(module.generators,
620                                sub,dim,subdim,F,1);
621
622   if ans=fail then
623     return fail;
624   fi;
625
626   # conjugate the matrices to correspond to given sub
627   if SMTX.IsZeroGens(module) then
628     smodule:=GModuleByMats([],Length(ans.smatrices[1]),F);
629   else
630    smodule:=GModuleByMats(List(ans.smatrices,i->i^nb),F);
631   fi;
632   return smodule;
633end;
634
635SMTX.ProperSubmoduleBasis:=function(module)
636  if SMTX.IsIrreducible(module) then
637    return fail;
638  fi;
639  return SMTX.Subbasis(module);
640end;
641
642
643#############################################################################
644##
645#F  SMTX.InducedActionFactorModule( module, sub [,compl] )
646##
647## module is a module record, and sub is a list of generators of a submodule.
648## (i.e. each has leading coefficient 1 in a unique place).
649## Qmodule is returned, where qmodule
650## is the quotient module.
651##
652SMTX.InducedActionFactorModule:=function(arg)
653local module,sub,  ans, dim, subdim, F,qmodule;
654
655   module:=arg[1];
656   sub:=arg[2];
657
658   sub:=TriangulizedMat(sub);
659
660   subdim:=Length(sub);
661   dim:=SMTX.Dimension(module);
662   if subdim = dim then
663      return List(module.generators,i->[[]]);
664   fi;
665
666   F:=SMTX.Field(module);
667
668   ans:=SMTX.SubQuotActions(module.generators,
669                                sub,dim,subdim,F,2);
670
671   if ans=fail then
672     return fail;
673   fi;
674
675   if Length(arg)=3 then
676     # compute basechange
677     sub:=Concatenation(sub,arg[3]);
678     sub:=sub*Inverse(ans.nbasis);
679     ans.qmatrices:=List(ans.qmatrices,i->i^sub);
680   fi;
681
682   if SMTX.IsZeroGens(module) then
683     qmodule:=GModuleByMats([],Length(ans.qmatrices[1]),F);
684   else
685    qmodule:=GModuleByMats(ans.qmatrices, F);
686   fi;
687   return qmodule;
688
689end;
690
691#############################################################################
692##
693#F  SMTX.InducedActionFactorModuleWithBasis( module, sub )
694##
695# FIXME: this function is never used and documented. Keep it or remove it?
696SMTX.InducedActionFactorModuleWithBasis:=function(module,sub)
697local ans, dim, subdim, F,qmodule;
698
699   sub:=TriangulizedMat(sub);
700
701   subdim:=Length(sub);
702   dim:=SMTX.Dimension(module);
703   if subdim = dim then
704      return List(module.generators,i->[[]]);
705   fi;
706
707   F:=SMTX.Field(module);
708
709   ans:=SMTX.SubQuotActions(module.generators,
710                                sub,dim,subdim,F,2);
711
712   if ans=fail then
713     return fail;
714   fi;
715
716   # fetch new basis
717   sub:=ans.nbasis{[Length(sub)+1..module.dimension]};
718
719   if SMTX.IsZeroGens(module) then
720     qmodule:=GModuleByMats([],Length(ans.qmatrices[1]),F);
721   else
722    qmodule:=GModuleByMats(ans.qmatrices, F);
723   fi;
724   return [qmodule,sub];
725
726end;
727
728#############################################################################
729##
730#F  SMTX.InducedAction( module, sub, typ )
731##  generators of sub- and quotient-module and original module wrt new basis
732##  and new basis
733##
734## module is a module record, and sub is a list of generators of a submodule.
735## IT IS ASSUMED THAT THE GENERATORS OF SUB ARE NORMED.
736## (i.e. each has leading coefficient 1 in a unique place).
737## SMTX.InducedAction computes the submodule and quotient
738## and the original module with its matrices written wrt to the basis used
739## to compute smodule and qmodule.
740## [smodule, qmodule, nmodule] is returned,
741## where smodule is the submodule and qmodule the quotient module.
742## The matrices of nmodule have the form  A  0  where  A  and  B  are the
743##                                        C  B
744## corresponding matrices of smodule and qmodule resepctively.
745## If sub is not the basis of a submodule then fail is returned.
746SMTX.InducedAction:=function( arg )
747local module,sub,typ,ans,dim,subdim,F,one,erg;
748
749   module:=arg[1];
750   sub:=arg[2];
751   if Length(arg)>2 then
752     typ:=arg[3];
753   else
754     typ:=7;
755   fi;
756   subdim:=Length(sub);
757   dim:=SMTX.Dimension(module);
758   F:=SMTX.Field(module); one:=One(F);
759
760   erg:=SMTX.SubQuotActions(module.generators,
761                                sub,dim,subdim,F,typ);
762
763   if erg=fail then
764     return fail;
765   fi;
766
767   ans:=[];
768
769   if IsBound(erg.smatrices) then
770     if SMTX.IsZeroGens(module) then
771       Add(ans,GModuleByMats([],Length(erg.smatrices[1]), F));
772     else
773       Add(ans,GModuleByMats(erg.smatrices, F));
774     fi;
775   fi;
776   if IsBound(erg.qmatrices) then
777     if SMTX.IsZeroGens(module) then
778       Add(ans,GModuleByMats([],Length(erg.qmatrices[1]), F));
779     else
780       Add(ans,GModuleByMats(erg.qmatrices, F));
781     fi;
782   fi;
783   if IsBound(erg.nmatrices) then
784     if SMTX.IsZeroGens(module) then
785       Add(ans,GModuleByMats([],Length(erg.nmatrices[1]), F));
786     else
787       Add(ans,GModuleByMats(erg.nmatrices, F));
788     fi;
789   fi;
790   if IsBound(erg.nbasis) then
791     Add(ans,erg.nbasis);
792   fi;
793
794   return ans;
795
796end;
797
798#############################################################################
799##
800#F  SMTX.InducedActionSubMatrixNB( mat, sub ) . . . . construct submodule
801##
802##  as InducedActionSubmoduleNB but for a matrix.
803# FIXME: this function is never used and documented. Keep it or remove it?
804SMTX.InducedActionSubMatrixNB:=function( mat, sub )
805local subdim, dim, F, ans;
806
807   subdim:=Length(sub);
808   if subdim = 0 then
809      return [];
810   fi;
811   dim:=Length(mat);
812   F:=DefaultFieldOfMatrix(mat);
813
814   ans:=SMTX.SubQuotActions([mat],sub,dim,subdim,F,1);
815
816   if ans=fail then
817     return fail;
818   else
819     return ans.smatrices[1];
820   fi;
821
822end;
823
824# Ditto, but allowing also unnormed modules
825# FIXME: this function is never used and documented. Keep it or remove it?
826SMTX.InducedActionSubMatrix:=function(mat,sub)
827local nb, subdim, dim, F, ans;
828  nb:=SMTX.NormedBasisAndBaseChange(sub);
829  sub:=nb[1];
830  nb:=nb[2];
831
832   subdim:=Length(sub);
833   if subdim = 0 then
834      return [];
835   fi;
836   dim:=Length(mat);
837   F:=DefaultFieldOfMatrix(mat);
838
839   ans:=SMTX.SubQuotActions([mat],sub,dim,subdim,F,1);
840
841   if ans=fail then
842     return fail;
843   else
844    # conjugate the matrices to correspond to given sub
845     return ans.smatrices[1]^nb;
846   fi;
847
848end;
849
850#############################################################################
851##
852#F  SMTX.InducedActionFactorMatrix( mat, sub [,compl] )
853##
854##  as InducedActionFactor, but for a matrix.
855##
856SMTX.InducedActionFactorMatrix:=function(arg)
857local mat, sub, subdim, dim, F, ans;
858
859   mat:=arg[1];
860   sub:=arg[2];
861
862   sub:=TriangulizedMat(sub);
863
864   subdim:=Length(sub);
865   dim:=Length(mat);
866   if subdim = dim then
867      return [];
868   fi;
869
870   F:=DefaultFieldOfMatrix(mat);
871
872   ans:=SMTX.SubQuotActions([mat],sub,dim,subdim,F,2);
873
874   if ans=fail then
875     return fail;
876   fi;
877
878   if Length(arg)=3 then
879     # compute basechange
880     sub:=Concatenation(sub,arg[3]);
881     sub:=sub*Inverse(ans.nbasis);
882     ans.qmatrices:=List(ans.qmatrices,i->i^sub);
883   fi;
884
885   return ans.qmatrices[1];
886
887end;
888
889SMTX_SMCoRaEl:=function(matrices,ngens,newgenlist,dim,F)
890local g1,g2,coefflist,M,pol;
891  g1:=Random(1, ngens);
892  g2:=g1;
893  while g2=g1 and ngens>1 do
894     g2:=Random(1, ngens);
895  od;
896  ngens:=ngens + 1;
897  matrices[ngens]:=matrices[g1] * matrices[g2];
898  Add(newgenlist, [g1, g2]);
899  # Take a random linear sum of the existing generators as new generator.
900  # Record the sum in coefflist
901  coefflist:=[Random(F)];
902  M:=coefflist[1]*matrices[1];
903  for g1 in [2..ngens] do
904     g2:=Random(F);
905     if IsOne(g2) then
906       M:=M + matrices[g1];
907     elif not IsZero(g2) then
908       M:=M + g2 * matrices[g1];
909     fi;
910     Add(coefflist, g2);
911  od;
912  Info(InfoMeatAxe,2,"Evaluated random element in algebra.");
913  pol:=CharacteristicPolynomialMatrixNC(F,M,1);
914  return [M,coefflist,pol];
915end;
916SMTX.SMCoRaEl:=SMTX_SMCoRaEl;
917
918# how many random elements should we try before (temporarily ) giving up?
919# This number is set relatively high to minimize the chance of an unlucky
920# random run in functions such as composition series computation.
921SMTX.RAND_ELM_LIMIT:=5000;
922
923#############################################################################
924##
925#F  SMTX.IrreduciblityTest( module ) try to reduce a module over a finite
926##                                      field
927##
928## 27/12/2000.
929## New version incorporating Ivanyos/Lux method of handling one difficult case
930## for proving reducibility.
931## (See G.Ivanyos and K. Lux, `Treating the exceptional cases of the meataxe',
932##  Experimental Mathematics 9, 2000, 373-381.
933##
934## module is a module record
935## IsIrreducible( ) attempts to decide whether module is irreducible.
936## When it succeeds it returns true or false.
937## We choose at random elements of the group algebra of the group.
938## If el is such an element, we define M, p, fac, N, e and v as follows:-
939## M is the matrix corresponding to el, p is its characteristic polynomial,
940## fac an irreducible factor of p, N the nullspace of the matrix fac(M),
941## ndim the dimension of N, and v a vector in N.
942## If we can find the above such that ndim = deg(fac) then we can test
943## conclusively for irreducibility. Then, in the case where irreducibility is
944## proved, we store the information as fields for the module, since it may be
945## useful later (e.g. to test for absolute irreducibility, equivalence with
946## another module).
947## These  fields are accessed by the functions
948## AlgEl()(el), AlgElMat(M), AlgElCharPol(p),
949## AlgElCharPolFac(fac), AlgElNullspaceDimension(ndim), and
950## AlgElNullspaceVec(v).
951##
952## If we cannot find such a set with ndim = deg(fac) we may nonetheless prove
953## reducibility  by finding a submodule. However we can never prove
954## irreducibility without such a set (and hence the algorithm could run
955## forever, but hopefully this will never happen!)
956## Where reducibility is proved, we set the field .subbasis
957## (a basis for the submodule, normed in the sense that the first non-zero
958## component of each basis vector is 1, and is in a different position from
959## the first non-zero component of every other basis vector).
960## The test for irreducibility is based on the meataxe method (but in the
961## meataxe, ndim is always very small, usually 1. The modification here is put
962## in to enable the method to work over modules with large centralizing fields).
963## We simply spin v. If we do not get  the whole space, we have a submodule,
964## on the other hand, if we do get the whole space, we calculate the
965## nullspace NT of the transpose of fac(M), spin that under the group
966## generated by the transposes of the generating matrices, and thus either
967## find the transpose of a submodule or conclusively prove irreducibility.
968##
969## This function can also be used to get a random submodule. Therefore it
970## is not an end-user function but only called internally
971SMTX_IrreducibilityTest:=function( module )
972   local matrices, tmatrices, ngens, ans,  M, mat, g1, g2, maxdeg,
973         newgenlist, coefflist, orig_ngens, zero,
974         N, NT, v, subbasis, fac, sfac, pol, orig_pol, q, dim, ndim, i,
975         l, trying, deg, facno, bestfacno, F, count, R, rt0,idmat,
976         pfac1, pfac2, quotRem, pfr, idemp, M2, mat2, mat3;
977
978   rt0:=Runtime();
979   Info(InfoMeatAxe,1,"Calling MeatAxe. All times will be in milliseconds");
980   if not SMTX.IsMTXModule(module) then
981      Error("Argument of IsIrreducible is not a module.");
982   fi;
983
984   if not module.IsOverFiniteField then
985      Error("Argument of IsIrreducible is not over a finite field.");
986   fi;
987   matrices:=ShallowCopy(module.generators);
988   dim:=SMTX.Dimension(module);
989   ngens:=Length(matrices);
990   orig_ngens:=ngens;
991   F:=SMTX.Field(module);
992   zero:=Zero(F);
993   R:=PolynomialRing(F);
994
995   # Now compute random elements M of the group algebra, calculate their
996   # characteristic polynomials, factorize, and apply the irreducible factors
997   # to M to get matrices with nontrivial nullspaces.
998   # tmatrices will be a list of the transposed generators if required.
999
1000   tmatrices:=[];
1001   trying:=true;
1002   # trying will become false when we have an answer
1003   maxdeg:=1;
1004   newgenlist:=[];
1005   # Do a small amount of preprocessing to increase the generator set.
1006   for i in [1..1] do
1007      g1:=Random(1, ngens);
1008      g2:=g1;
1009      while g2=g1 and Length(matrices) > 1 do
1010         g2:=Random(1, ngens);
1011      od;
1012      ngens:=ngens + 1;
1013      matrices[ngens]:=matrices[g1] * matrices[g2];
1014      Add(newgenlist, [g1, g2]);
1015   od;
1016   Info(InfoMeatAxe,1,"Done preprocessing. Time = ",Runtime()-rt0,".");
1017   count:=0;
1018
1019   # Main loop starts - choose a random element of group algebra on each pass
1020   while trying  do
1021      count:=count + 1;
1022      if count mod SMTX.RAND_ELM_LIMIT = 0 then
1023         Error("Have generated ",SMTX.RAND_ELM_LIMIT,
1024                "random elements and failed to prove\n",
1025                "or disprove irreducibility. Type return to keep trying.");
1026      fi;
1027      maxdeg:=Minimum(maxdeg * 2,dim);
1028      # On this pass, we only consider irreducible factors up to degree maxdeg.
1029      # Using higher degree factors is very time consuming, so we prefer to try
1030      # another element.
1031      # To choose random element, first add on a new generator as a product of
1032      # two randomly chosen unequal existing generators
1033      # Record the product in newgenlist.
1034      Info(InfoMeatAxe,1,"Choosing random element number ",count);
1035
1036      M:=SMTX.SMCoRaEl(matrices,ngens,newgenlist,dim,F);
1037
1038      ngens:=Length(matrices);
1039
1040      coefflist:=M[2];
1041      pol:=M[3];
1042      M:=ImmutableMatrix(F,M[1]);
1043      idmat:=M^0;
1044
1045      orig_pol:=pol;
1046      Info(InfoMeatAxe,2,"Evaluated characteristic polynomial. Time = ",
1047           Runtime()-rt0,".");
1048      # Now we extract the irreducible factors of pol starting with those
1049      # of low degree
1050      deg:=0;
1051      fac:=[];
1052      # The next loop is through the degrees of irreducible factors
1053      while DegreeOfLaurentPolynomial(pol) > 0 and deg < maxdeg and trying do
1054         repeat
1055            deg:=deg + 1;
1056            if deg > Int(DegreeOfLaurentPolynomial(pol) / 2) then
1057               fac:=[pol];
1058            else
1059               fac:=Factors(R, pol: factoroptions:=rec(onlydegs:=[deg]));
1060               fac:=Filtered(fac,i->DegreeOfLaurentPolynomial(i)=deg);
1061               Info(InfoMeatAxe,2,Length(fac)," factors of degree ",deg,
1062                    ", Time = ",Runtime()-rt0,".");
1063            fi;
1064         until fac <> [] or deg = maxdeg;
1065
1066         if fac <> [] then
1067            if DegreeOfLaurentPolynomial(fac[1]) = dim then
1068               # In this case the char poly is irreducible, so the
1069               # module is irreducible.
1070               ans:=true;
1071               trying:=false;
1072               bestfacno:=1;
1073               v:=ListWithIdenticalEntries(dim,zero);
1074               v[1]:=One(F);
1075               ndim:=dim;
1076            fi;
1077            # Otherwise, first see if there is a non-repeating factor.
1078            # If so it will be decisive, so delete the rest of the list
1079            l:=Length(fac);
1080            facno:=1;
1081            while facno <= l and trying do
1082               if facno = l  or  fac[facno] <> fac[facno + 1] then
1083                  fac:=[fac[facno]]; l:=1;
1084               else
1085                  while facno < l and fac[facno] = fac[facno + 1] do
1086                     facno:=facno + 1;
1087                  od;
1088               fi;
1089               facno:=facno + 1;
1090            od;
1091            # Now we can delete repetitions from the list fac
1092            sfac:=Set(fac);
1093
1094            if DegreeOfLaurentPolynomial(fac[1]) <> dim then
1095               # Now go through the factors and attempt to find a submodule
1096               facno:=1; l:=Length(sfac);
1097               while facno <= l and trying do
1098                  mat:=Value(sfac[facno], M,idmat);
1099                  Info(InfoMeatAxe,2,"Evaluated matrix on factor. Time = ",
1100                       Runtime()-rt0,".");
1101                  N:=NullspaceMat(mat);
1102                  v:=N[1];
1103                  ndim:=Length(N);
1104
1105                  Info(InfoMeatAxe,2,"Evaluated nullspace. Dimension = ",
1106                       ndim,". Time = ",Runtime()-rt0,".");
1107                  subbasis:=SMTX.SpinnedBasis(v, matrices, F,orig_ngens);
1108                  Info(InfoMeatAxe,2,"Spun up vector. Dimension = ",
1109                       Length(subbasis),". Time = ",Runtime()-rt0,".");
1110                  if Length(subbasis) < dim then
1111                     # Proper submodule found
1112                     trying:=false;
1113                     ans:=false;
1114                     SMTX.SetSubbasis(module, subbasis);
1115                  elif ndim = deg then
1116                     trying:=false;
1117                     # if we transpose and find no proper submodule, then the
1118                     # module is definitely irreducible.
1119                     mat:=TransposedMat(mat);
1120                     if Length(tmatrices)=0 then
1121                        for i in [1..orig_ngens] do
1122                           Add(tmatrices, TransposedMat(matrices[i]));
1123                        od;
1124                     fi;
1125                     Info(InfoMeatAxe,2,"Transposed matrices. Time = ",
1126                          Runtime()-rt0,".");
1127                     NT:=NullspaceMat(mat);
1128                     Info(InfoMeatAxe,2,"Evaluated nullspace. Dimension = ",
1129                          Length(NT),". Time = ",Runtime()-rt0, ".");
1130                     subbasis:=SMTX.SpinnedBasis(NT[1],tmatrices,F,orig_ngens);
1131                     Info(InfoMeatAxe,2,"Spun up vector. Dimension = ",
1132                          Length(subbasis),". Time = ",Runtime()-rt0, ".");
1133                     if Length(subbasis) < dim then
1134                        # subbasis is a basis for a submodule of the transposed
1135                        # module, and the orthogonal complement of this is a
1136                        # submodule of the original module. So we find a vector
1137                        # v in that, and then spin it. Of course we won't
1138                        # necessarily get the full orthogonal complement
1139                        # that way, but we'll certainly get a proper submodule.
1140                        v:=SMTX.OrthogonalVector(subbasis);
1141                        SMTX.SetSubbasis(module,
1142                          SMTX.SpinnedBasis(v,matrices,F,orig_ngens));
1143                        ans:=false;
1144                     else
1145                        ans:=true;
1146                        bestfacno:=facno;
1147                     fi;
1148                  fi;
1149                  if trying and deg>1 and count>2 then
1150                     Info(InfoMeatAxe,1,"Trying Ivanyos/Lux Method");
1151                     # first find the appropriate idempotent
1152                     pfac1:=sfac[facno];
1153                     pfac2:=orig_pol;
1154                     while true do
1155                       quotRem := QuotRemLaurpols(pfac2, sfac[facno], 3);
1156                       if quotRem[2] <> Zero(R) then
1157                         break;
1158                       fi;
1159                       pfac2 := quotRem[1];
1160                       pfac1:=pfac1*sfac[facno];
1161                     od;
1162                     pfr:=GcdRepresentation(pfac1, pfac2);
1163                     idemp:=QuotRemLaurpols(pfr[2]*pfac2, orig_pol, 2);
1164                     # Now another random element in the group algebra.
1165                     # and a random vector in the module
1166                     g2:=Random(F);
1167                     if IsOne(g2) then
1168                       M2:=matrices[1];
1169                     else
1170                       M2:=g2 * matrices[1];
1171                     fi;
1172                     for g1 in [2..ngens] do
1173                        g2:=Random(F);
1174                        if IsOne(g2) then
1175                          M2:=M2 + matrices[g1];
1176                        elif not IsZero(g2) then
1177                          M2:=M2 + g2 * matrices[g1];
1178                        fi;
1179                     od;
1180                     Info(InfoMeatAxe,2,
1181                         "Evaluated second random element in algebra.");
1182                     v:=Random(FullRowSpace(F,dim));
1183                     mat2:=Value(idemp, M,idmat);
1184                     mat3:=mat2*M2*mat2;
1185                     v:=v*(M*mat3 - mat3*M);
1186                     # This vector might lie in a proper subspace!
1187                     subbasis:=SMTX.SpinnedBasis(v, matrices, F,orig_ngens);
1188                     Info(InfoMeatAxe,2,"Spun up vector. Dimension = ",
1189                       Length(subbasis),". Time = ",Runtime()-rt0,".");
1190                     if Length(subbasis) < dim and Length(subbasis) <> 0  then
1191                       # Proper submodule found
1192                       trying:=false;
1193                       ans:=false;
1194                       SMTX.SetSubbasis(module, subbasis);
1195                    fi;
1196                  fi;
1197                  facno:=facno + 1;
1198               od; # going through irreducible factors of fixed degree.
1199               # If trying is false at this stage, then we don't have
1200               # an answer yet, so we have to go onto factors of the next degree.
1201               # Now divide p by the factors used if necessary
1202               if trying and deg < maxdeg then
1203                  for q in fac do
1204                     pol:=Quotient(R, pol, q);
1205                  od;
1206               fi;
1207            fi;           #DegreeOfLaurentPolynomial(fac[1]) <> dim
1208         fi;             #fac <> []
1209      od; # loop through degrees of irreducible factors
1210
1211      # if we have not found a submodule and trying is false, then the module
1212      # must be irreducible.
1213      if trying = false and ans = true then
1214         SMTX.SetAlgEl(module, [newgenlist, coefflist]);
1215         SMTX.SetAlgElMat(module, M);
1216         SMTX.SetAlgElCharPol(module, orig_pol);
1217         SMTX.SetAlgElCharPolFac(module, sfac[bestfacno]);
1218         SMTX.SetAlgElNullspaceVec(module, v);
1219         SMTX.SetAlgElNullspaceDimension(module, ndim);
1220      fi;
1221
1222   od;  # main loop
1223
1224   Info(InfoMeatAxe,1,"Total time = ",Runtime()-rt0," milliseconds.");
1225   return ans;
1226
1227end;
1228
1229SMTX.IrreducibilityTest:=SMTX_IrreducibilityTest;
1230
1231SMTX.IsIrreducible:=function(module)
1232  if not IsBound(module.IsIrreducible) then
1233    module.IsIrreducible:=SMTX.IrreducibilityTest(module);
1234  fi;
1235  return module.IsIrreducible;
1236end;
1237
1238#############################################################################
1239##
1240#F SMTX.RandomIrreducibleSubGModule( module ) . . .
1241## find a basis for a random irreducible
1242## submodule of module, and return that basis and the submodule, with all
1243## the irreducibility flags set.
1244## Returns false if module is irreducible.
1245SMTX_RandomIrreducibleSubGModule:=function( module )
1246   local  ranSub, subbasis, submodule, subbasis2, submodule2,
1247   F, dim, el, M, fac, N, i, matrices, ngens, genpair;
1248
1249   if not SMTX.IsMTXModule(module) then
1250      return Error("Argument of RandomIrreducibleSubGModule is not a module.");
1251   elif SMTX.HasIsIrreducible(module) and SMTX.IsIrreducible(module) then
1252      return false;
1253   fi;
1254
1255   # now call an irreducibility test that will compute a new subbasis
1256
1257   i:=SMTX.IrreducibilityTest(module);
1258
1259   if i then
1260     # we just found out it is irreducible
1261     SMTX.SetIsIrreducible(module,true);
1262     return false;
1263   elif not SMTX.HasIsIrreducible(module) then
1264     # or store reducibility
1265     SMTX.SetIsIrreducible(module,false);
1266   fi;
1267
1268   subbasis:=SMTX.Subbasis(module);
1269   submodule:=SMTX.InducedActionSubmoduleNB(module, subbasis);
1270   ranSub:=SMTX.RandomIrreducibleSubGModule(submodule);
1271   if ranSub = false then
1272      # submodule has been proved irreducible in a call to this function,
1273      # so the flags have been set.
1274      return [ subbasis, submodule];
1275   else
1276      # ranSub[1] is given in terms of the basis for the submodule,
1277      # but we want it in terms of the basis of the original module.
1278      # So we multiply it by subbasis.
1279      # Then we need our basis to be normed.
1280      # this is done by triangulization
1281      F:=SMTX.Field(module);
1282      subbasis2:=ranSub[1] * subbasis;
1283      subbasis2:=TriangulizedMat(subbasis2);
1284
1285      # But now since we've normed the basis subbasis2,
1286      # the matrices of the submodule ranSub[2] are given with respect to
1287      # the wrong basis.  So we have to recompute the submodule.
1288      submodule2:=SMTX.InducedActionSubmoduleNB(module, subbasis2);
1289      # Unfortunately, although it's clear that this submodule is
1290      # irreducible, we'll have to reset the flags that IsIrreducible sets.
1291      # AH Why can't we keep irreducibility?
1292
1293      # Some will be the same # as in ranSub[2], but some are affected by
1294      # the base change, or at least part of it, since the flags gets
1295      # screwed up by the base change.
1296      # We need to set the following flags:
1297      el:=SMTX.AlgEl(ranSub[2]);
1298      SMTX.SetAlgEl(submodule2,el);
1299      SMTX.SetAlgElCharPol(submodule2,SMTX.AlgElCharPol(ranSub[2]));
1300      fac:=SMTX.AlgElCharPolFac(ranSub[2]);
1301      SMTX.SetAlgElCharPolFac(submodule2,fac);
1302      SMTX.SetAlgElNullspaceDimension(submodule2,
1303             SMTX.AlgElNullspaceDimension(ranSub[2]));
1304
1305      # Only the actual algebra element and its nullspace have to be recomputed
1306      # This code is essentially from IsomorphismGModule
1307      dim:=SMTX.Dimension(submodule2);
1308      matrices:=ShallowCopy(submodule2.generators);
1309      ngens:=Length(matrices);
1310      for genpair in el[1] do
1311         ngens:=ngens + 1;
1312         matrices[ngens]:=matrices[genpair[1]] * matrices[genpair[2]];
1313      od;
1314      M:=ImmutableMatrix(F,Sum([1..ngens], i-> el[2][i] * matrices[i]));
1315      SMTX.SetAlgElMat(submodule2,M);
1316      N:=NullspaceMat(Value(fac,M,M^0));
1317      SMTX.SetAlgElNullspaceVec(submodule2,N[1]);
1318      return [subbasis2, submodule2];
1319   fi;
1320
1321end;
1322SMTX.RandomIrreducibleSubGModule:=SMTX_RandomIrreducibleSubGModule;
1323
1324#############################################################################
1325##
1326#F  SMTX.GoodElementGModule( module ) . .  find good group algebra element
1327##                                       in an irreducible module
1328##
1329## module is a module that is already known to be irreducible.
1330## GoodElementGModule finds a group algebra element with nullspace of
1331## minimal possible dimension. This dimension is 1 if the module is absolutely
1332## irreducible, and the degree of the relevant field extension otherwise.
1333## This is needed for testing for equivalence of modules.
1334SMTX_GoodElementGModule:=function( module )
1335local matrices, ngens, M, mat,  N, newgenlist, coefflist, orig_ngens,
1336      fac, sfac, pol, oldpol,  q, deg, i, l,
1337      trying, dim, mindim, F, R, count, rt0, idmat;
1338
1339   rt0:=Runtime();
1340   if not SMTX.IsMTXModule(module) or not SMTX.IsIrreducible(module) then
1341     ErrorNoReturn("Argument is not an irreducible module.");
1342   fi;
1343   if not SMTX.HasIsAbsolutelyIrreducible(module) then
1344      SMTX.IsAbsolutelyIrreducible(module);
1345   fi;
1346   if  SMTX.IsAbsolutelyIrreducible(module) then
1347     mindim:=1;
1348   else
1349     mindim:=SMTX.DegreeFieldExt(module);
1350   fi;
1351
1352   if SMTX.AlgElNullspaceDimension(module) = mindim then return; fi;
1353   # This is the condition that we want. If it holds already, then there is
1354   # nothing else to do.
1355
1356   dim:=SMTX.Dimension(module);
1357   matrices:=ShallowCopy(module.generators);
1358   ngens:=Length(matrices);
1359   orig_ngens:=ngens;
1360   F:=SMTX.Field(module);
1361   R:=PolynomialRing(F);
1362
1363   # Now compute random elements el of the group algebra, calculate their
1364   # characteristic polynomials, factorize, and apply the irreducible factors
1365   # to el to get matrices with nontrivial nullspaces.
1366
1367   trying:=true;
1368   count:=0;
1369   newgenlist:=[];
1370   while trying do
1371      count:=count + 1;
1372      if count mod SMTX.RAND_ELM_LIMIT = 0 then
1373         Error("Have generated ",SMTX.RAND_ELM_LIMIT,
1374                " random elements and failed ",
1375                "to find a good one. Type return to keep trying.");
1376      fi;
1377      Info(InfoMeatAxe,2,"Choosing random element number ",count,".");
1378
1379      M:=SMTX.SMCoRaEl(matrices,ngens,newgenlist,dim,F);
1380      ngens:=Length(matrices);
1381
1382      coefflist:=M[2];
1383      pol:=M[3];
1384      M:=M[1];
1385      idmat:=M^0;
1386
1387      Info(InfoMeatAxe,2,"Evaluated characteristic polynomial. Time = ",
1388           Runtime()-rt0,".");
1389      # That is necessary in case p is defined over a smaller field that F.
1390      oldpol:=pol;
1391      # Now we extract the irreducible factors of pol starting with those
1392      # of low degree
1393      deg:=0;
1394      fac:=[];
1395      while deg  <= mindim and trying do
1396         repeat
1397            deg:=deg + 1;
1398            if deg > mindim then
1399               fac:=[pol];
1400            else
1401               fac:=Factors(R, pol: factoroptions:=rec(onlydegs:=[deg]));
1402               fac:=Filtered(fac,i->DegreeOfLaurentPolynomial(i)<=deg);
1403               Info(InfoMeatAxe,2,Length(fac)," factors of degree ",deg,
1404                    ", Time = ",Runtime()-rt0,".");
1405               sfac:=Set(fac);
1406            fi;
1407         until fac <> [];
1408         l:=Length(fac);
1409         if trying and deg <= mindim then
1410            i:=1;
1411            while i <= l and trying do
1412               mat:=Value(fac[i], M,idmat);
1413               Info(InfoMeatAxe,2,"Evaluated matrix on factor. Time = ",
1414                    Runtime()-rt0,".");
1415               N:=NullspaceMat(mat);
1416               Info(InfoMeatAxe,2,"Evaluated nullspace. Dimension = ",
1417                    Length(N),". Time = ",Runtime()-rt0,".");
1418               if Length(N) = mindim then
1419                  trying:=false;
1420                  SMTX.SetAlgEl(module, [newgenlist, coefflist]);
1421                  SMTX.SetAlgElMat(module, M);
1422                  SMTX.SetAlgElCharPol(module, oldpol);
1423                  SMTX.SetAlgElCharPolFac(module, fac[i]);
1424                  SMTX.SetAlgElNullspaceVec(module, N[1]);
1425                  SMTX.SetAlgElNullspaceDimension(module, Length(N));
1426               fi;
1427               i:=i + 1;
1428            od;
1429         fi;
1430
1431         if trying then
1432            for q in fac do
1433               pol:=Quotient(R, pol, q);
1434            od;
1435         fi;
1436      od;
1437   od;
1438   Info(InfoMeatAxe,1,"Total time = ",Runtime()-rt0," milliseconds.");
1439
1440end;
1441SMTX.GoodElementGModule:=SMTX_GoodElementGModule;
1442
1443#############################################################################
1444##
1445#F  EnlargedIrreducibleGModule(module, mat) . .add a generator to a module that
1446#
1447# 2bdef!
1448
1449
1450#############################################################################
1451##
1452#F  SMTX.FrobeniusAction(A, v [, basis]) . . action of matrix A on
1453##                                  . . Frobenius block of vector v
1454##
1455## FrobeniusAction(A, v) computes the Frobenius block of the dxd matrix A
1456## generated by the length - d vector v, and returns it.
1457## It is based on code of MinPolCoeffsMat.
1458## The optional third argument is for returning the basis for this block.
1459##
1460SMTX_FrobeniusAction:=function( arg )
1461local   L, d, p, M, one, zero, R, h, v, w, i, j, nd, ans,
1462        A, basis;
1463
1464   if Length(arg) = 2  then
1465      A:=arg[1];
1466      v:=arg[2];
1467      basis:=0;
1468   elif Length(arg) = 3  then
1469      A:=arg[1];
1470      v:=arg[2];
1471      basis:=arg[3];
1472   else
1473      return Error("usage: SMTX.FrobeniusAction( <A>, <v>, [, <basis>] )");
1474   fi;
1475   one :=One(A[1][1]);
1476   zero:=Zero(one);
1477   d:=Length( A );
1478   M:=ListWithIdenticalEntries(Length(A[1]) + 1,zero);
1479   ConvertToVectorRep(M,DefaultField(v));
1480
1481   # L[i] (length d) will contain a vector with head entry 1 at position i,
1482   # which is in the current block.
1483   # R[i] (length d + 1 but (d + 1) - entry always 0) is vector expressing
1484   # L[i] in terms of the basis of the block.
1485   L:=[];
1486   R:=[];
1487
1488   # <j> - 1 gives the power of <A> we are looking at
1489   j:=1;
1490
1491   # spin vector around and construct polynomial
1492   repeat
1493
1494      # compute the head of <v>
1495      h:=1;
1496      while v[h] = zero  do
1497         h:=h + 1;
1498      od;
1499
1500      # start with appropriate polynomial x^(<j> - 1)
1501      p:=ShallowCopy( M );
1502      p[j]:=one;
1503
1504      # divide by known left sides
1505      w:=v;
1506      while h <= d and IsBound( L[h] ) do
1507         p:=p - w[h] * R[h];
1508         w:=w - w[h] * L[h];
1509         while h <= d and w[h] = zero do
1510            h:=h + 1;
1511         od;
1512      od;
1513
1514      # if <v> is not the zero vector try next power
1515      if h <= d  then
1516         #AH replaced Copy by ShallowCopy as only vector is used
1517         if basis <> 0 then basis[j]:=ShallowCopy(v); fi;
1518         R[h]:=p * w[h]^-1;
1519         L[h]:=w * w[h]^-1;
1520         j:=j + 1;
1521         v:=v * A;
1522      fi;
1523   until h > d;
1524
1525   nd:=Length(p);
1526   while 0 < nd  and p[nd] = zero  do
1527      nd:=nd - 1;
1528   od;
1529   nd:=nd - 1;
1530   ans:=[];
1531   for i in [1..nd - 1] do
1532      ans[i]:=[];
1533      for j in [1..nd] do ans[i][j]:=zero; od;
1534      ans[i][i + 1]:=one;
1535   od;
1536   ans[nd]:=[];
1537   for j in [1..nd] do
1538      ans[nd][j]:= - p[j];
1539   od;
1540
1541   return ans;
1542end;
1543SMTX.FrobeniusAction:=SMTX_FrobeniusAction;
1544
1545#############################################################################
1546##
1547#F SMTX.CompleteBasis(matrices,basis) . complete a basis under a group action
1548##
1549##  CompleteBasis( matrices, basis ) takes the partial basis 'basis' of the
1550##  underlying space of the (irreducible) module defined by matrices, and
1551##  attempts to extend it to a complete basis which is a direct sum of
1552##  translates of the original subspace under group elements. It returns
1553##  true or false according to whether it succeeds.
1554##  It is called by IsAbsolutelyIrreducible()
1555##
1556SMTX_CompleteBasis:=function( matrices, basis )
1557local  L, d, subd, subd0, zero, h, v, w, i, bno, gno, vno, newb, ngens;
1558
1559   subd:=Length(basis);
1560   subd0:=subd;
1561   d:=Length( basis[1] );
1562   if d = subd then
1563      return true;
1564   fi;
1565   # L is list of normalized generators of the subspace spanned by basis.
1566   L:=[];
1567   zero:=Zero(basis[1][1]);
1568   ngens:=Length(matrices);
1569
1570   # First find normalized generators for subspace itself.
1571   for i in [1..subd] do
1572      v:=basis[i];
1573      h:=1;
1574      while v[h] = zero  do
1575         h:=h + 1;
1576      od;
1577      w:=v;
1578      while h <= d and IsBound( L[h] )  do
1579         w:=w - w[h] * L[h];
1580         while h <= d and w[h] = zero  do
1581            h:=h + 1;
1582         od;
1583      od;
1584      if h <= d then
1585         L[h]:=w * w[h]^-1;
1586      else
1587         return Error("Initial vectors are not linearly independent.");
1588      fi;
1589   od;
1590
1591   # Now start translating
1592   bno:=1; gno:=1; vno:=1;
1593   while subd < d do
1594      # translate vector vno of block bno by generator gno
1595      v:= basis[ (bno - 1) * subd0 + vno] * matrices[gno];
1596      h:=1;
1597      while h<=d and v[h] = zero  do
1598         h:=h + 1;
1599      od;
1600      w:=v;
1601      while h <= d and IsBound( L[h] )  do
1602         w:=w - w[h] * L[h];
1603         while h <= d and w[h] = zero  do
1604            h:=h + 1;
1605         od;
1606      od;
1607      if h <= d then
1608         # new generator (and block)
1609         if vno = 1 then
1610            newb:=true;
1611         elif newb = false then
1612            return false;
1613         fi;
1614         L[h]:=w * w[h]^-1;
1615         subd:=subd + 1;
1616         basis[subd]:=v;
1617      else
1618         # in existing subspace
1619         if vno = 1 then
1620            newb:=false;
1621         elif newb = true then
1622            return false;
1623         fi;
1624      fi;
1625      vno:=vno + 1;
1626      if vno > subd0 then
1627         vno:=1;
1628         gno:=gno + 1;
1629         if gno > ngens then
1630            gno:=1;
1631            bno:=bno + 1;
1632         fi;
1633      fi;
1634   od;
1635
1636   return true;
1637end;
1638SMTX.CompleteBasis:=SMTX_CompleteBasis;
1639
1640#############################################################################
1641##
1642#F SMTX.AbsoluteIrreducibilityTest( module ) . . decide if an irreducible
1643##                    module over a  finite field is absolutely irreducible
1644##
1645## this function does the work for an absolute irreducibility test but does
1646## not actually set the flags.
1647## The function calculates the centralizer of the module.
1648## The centralizer should be isomorphic to the multiplicative
1649## group of the field GF(q^e) for some e, or rather to the group of
1650## dim/e x dim/e scalar matrices over GF(q^e), or equivalently,
1651## dim x dim matrices composed of identical e x e blocks along the diagonal.
1652##  e = 1 <=> the module is absolutely irreducible.
1653## The .fieldExtDeg component is set to e during the function call.
1654## The function shouldn't be called if the module has not already been
1655## shown to be irreducible, using IsIrreducible.
1656##
1657SMTX_AbsoluteIrreducibilityTest:=function( module )
1658local dim, ndim, gcd, div, e, ct, F, q, ok,
1659      M, v, M0, v0, C, C0, centmat, one, zero,
1660      pow, matrices, newmatrices, looking,
1661      basisN, basisB, basisBN, P, Pinv, i, j, k, nblocks;
1662
1663   if not SMTX.IsMTXModule(module) then
1664      Error("Argument of IsAbsoluteIrreducible is not a module");
1665   fi;
1666
1667   if not SMTX.IsIrreducible(module) then
1668      Error("Argument of iIsAbsoluteIrreducible s not an irreducible module");
1669   fi;
1670
1671   if not module.IsOverFiniteField then
1672      return Error("Argument of IsAbsoluteIrreducible is not over a finite field.");
1673   fi;
1674   dim:=SMTX.Dimension(module);
1675   F:=SMTX.Field(module);
1676   q:=Size(F);
1677   matrices:=module.generators;
1678
1679   # M acts irreducibly on N, which is canonically defined with respect to M
1680   # as the nullspace of fac(M), where fac is a factor of the char poly of M.
1681   # ndim is the dimension of N, and v is a vector of N. All these come from
1682   # the irreducibility test for the module.
1683   # An element of the centralizer must centralize every element, and
1684   # therefore M, and so must preserve N, since N is canonically defined
1685   # wrt M. Our plan is therefore first to find an element which centralizes
1686   # the restriction of M to N, and then extend it to the whole space.
1687
1688   M:=SMTX.AlgElMat(module);
1689   ndim:=SMTX.AlgElNullspaceDimension(module);
1690   v:=SMTX.AlgElNullspaceVec(module);
1691
1692   # e will have to divide both dim and ndim, and hence their gcd.
1693   gcd:=GcdInt(dim, ndim);
1694   Info(InfoMeatAxe,2,"GCD of module and nullspace dimensions = ", gcd, ".");
1695   if gcd = 1 then
1696      SMTX.SetDegreeFieldExt(module,1);
1697      return true;
1698   fi;
1699   div:=DivisorsInt(gcd);
1700
1701   # It's easy to find elements  in the centralizer of an element in Frobenius
1702   # (=rational canonical) form (centralizing elements are defined by their
1703   # action on the first basis element).
1704   # M0  is the Frobenius form for the action of M on N.
1705   # basisN is set by the function SMTX.FrobeniusAction to be the
1706   # basis v, vM, vM^2, .. for N
1707
1708   basisN:=[];
1709   Info(InfoMeatAxe,2,
1710     "Calc. Frobenius action of element from group algebra on nullspace.");
1711   M0:=SMTX.FrobeniusAction(M,v,basisN);
1712
1713   zero:=Zero(F);
1714   one:= One(F);
1715   v0:=ListWithIdenticalEntries(Length(M0[1]),zero);
1716   v0[1]:=one;
1717   ConvertToVectorRep(v0, F);
1718
1719   # v0 is just the vector (1, 0, 0....0) of length ndim. It has nothing
1720   # in particular to do with M0[1], but multiplying a vector that happens to be
1721   # around by 0 is a good way to get a zero vector of the right length.
1722
1723   # we try all possible divisors of gcd (biggest first) as possibilities for e
1724   # We're looking for a centralizing element with order dividing q^e - 1, and
1725   # blocks size e on N.
1726   for ct in Reversed([2..Length(div)]) do
1727      e:=div[ct];
1728      Info(InfoMeatAxe,2,"Trying dimension ",e," for centralising field.");
1729      # if ndim = e, M0 will do.
1730      if ndim > e then
1731         C:=M0;
1732         # Take the smallest power of C guaranteed to have order dividing
1733         # q^e - 1, and try that.
1734         pow:=(q^ndim - 1) / (q^e - 1);
1735         Info(InfoMeatAxe,2,"Looking for a suitable centralising element.");
1736         repeat
1737            # The first time through the loop C is M0, otherwise we choose C
1738            # at random from the centralizer of M0. Since M0 is in Frobenius
1739            # form any centralising element is determined by its top row
1740            # (which may be anything but the zero vector).
1741
1742            if Length(C)=0 then
1743               C[1]:=[];
1744               repeat
1745                  ok:=0;
1746                  for i in [1..ndim] do
1747                     C[1][i]:=Random(F);
1748                     if C[1][i] <> zero then  ok:=1; fi;
1749                  od;
1750               until ok=1;
1751               for i in [2..ndim] do C[i]:=C[i - 1] * M0; od;
1752               C:=ImmutableMatrix(F,C);
1753            fi;
1754            # C0 is the Frobenius form for the action of this power on one
1755            # of its blocks, B (all blocks have the same size). basisBN will
1756            # be set to be a basis for B, in terms of the elements of basisN.
1757            # A matrix product gives us the basis for B in terms of the
1758            # original basis for the module.
1759            basisBN:=[];
1760            C0:=SMTX.FrobeniusAction(C^pow,v0,basisBN);
1761            C:=[];
1762         until Length(C0) = e;
1763         Info(InfoMeatAxe,2,"Found one.");
1764         basisB:=ShallowCopy(basisBN * basisN);
1765      else
1766         C0:=M0;
1767         basisB:=ShallowCopy(basisN);
1768      fi;
1769      # Now try to extend basisB to a basis for the whole module, by
1770      # translating it by the generating matrices.
1771      Info(InfoMeatAxe,2,"Trying to extend basis to whole module.");
1772      if SMTX.CompleteBasis(matrices,basisB) then
1773         # We succeeded in extending the basis (might not have done).
1774         # So now we have a full basis, which we think of now as a base
1775         # change matrix.
1776         Info(InfoMeatAxe,2,"Succeeded. Calculating centralising matrix.");
1777         newmatrices:=[];
1778         P:=ImmutableMatrix(F,basisB);
1779         Pinv:=P^-1;
1780         for i in [1..Length(matrices)] do
1781            newmatrices[i]:=P * matrices[i] * Pinv;
1782         od;
1783         # Make the sum of copies of C0 as centmat
1784         centmat:=NullMat(dim, dim, F);
1785         nblocks:=dim/e;
1786         for i in [1..nblocks] do
1787            for j in [1..e] do
1788               for k in [1..e] do
1789                  centmat[ (i - 1) * e + j][ (i - 1) * e + k]:=C0[j][k];
1790               od;
1791            od;
1792         od;
1793         centmat := ImmutableMatrix(F, centmat);
1794         Info(InfoMeatAxe,2,"Checking that it centralises the generators.");
1795         # Check centralizing.
1796         looking:=true;
1797         i:=1;
1798         while looking and i <= Length(newmatrices) do
1799            if newmatrices[i] * centmat <> centmat * newmatrices[i] then
1800               looking:=false;
1801            fi;
1802            i:=i + 1;
1803         od;
1804         if looking then
1805            Info(InfoMeatAxe,2,"It did!");
1806            SMTX.SetDegreeFieldExt(module, e);
1807            SMTX.SetCentMat(module, Pinv * centmat * P); # get the base right
1808            # We will also record the minimal polynomial of C0 (and hence of
1809            # centmat) in case we need it at some future date.
1810            SMTX.SetCentMatMinPoly(module, MinimalPolynomialMatrixNC(F,C0,1));
1811            return false;
1812         fi;
1813         Info(InfoMeatAxe,2,"But it didn't.");
1814      else
1815         Info(InfoMeatAxe,2,"Failed!");
1816      fi;
1817   od;
1818
1819   Info(InfoMeatAxe,2,
1820     "Tried all divisors. Must be absolutely irreducible.");
1821   SMTX.SetDegreeFieldExt(module, 1);
1822   return true;
1823end;
1824SMTX.AbsoluteIrreducibilityTest:=SMTX_AbsoluteIrreducibilityTest;
1825
1826SMTX.DegreeFieldExt:=function(module)
1827  if not IsBound(module.smashMeataxe.degreeFieldExt) then
1828    SMTX.AbsoluteIrreducibilityTest( module );
1829  fi;
1830  return module.smashMeataxe.degreeFieldExt;
1831end;
1832
1833SMTX.DegreeSplittingField:=function(module)
1834  return DegreeOverPrimeField(SMTX.Field(module))
1835         *SMTX.DegreeFieldExt(module);
1836end;
1837
1838#############################################################################
1839##
1840#F  FieldGenCentMat( module ) . . find a centralizing matrix that generates
1841##                                the centralizing field of an irred. module
1842##
1843## FieldGenCentMat( ) should only be applied to modules that have already
1844## been proved irreducible using IsIrreducible. It then tests for absolute
1845## irreducibility (if not already known) and does nothing if module is
1846## absolutely irreducible. Otherwise, it returns a a matrix that generates
1847## (multiplicatively) the centralizing field (i.e. its multiplicative order
1848## is q^e - 1, where e is the degree of the centralizing field. This is not
1849## yet used, but maybe in future, if we wish to reduce the group to matrices
1850## over the larger field.
1851SMTX.FieldGenCentMat:=function( module )
1852   local e, F, R, q, qe, minpol, pp,
1853         M, v, M0, v0, C, C0, centmat, newcentmat, genpol, looking,
1854         i, l, okd;
1855
1856  if SMTX.FGCentMat(module)=fail then
1857    if SMTX.IsMTXModule(module) = false then
1858      Error("Argument of IsIrreducible is not a module.");
1859    fi;
1860
1861    if not SMTX.IsIrreducible(module) then
1862      Error("GModule is not irreducible.");
1863    fi;
1864
1865    # enforce absirred knowledge as well.
1866    #if not SMTX.IsAbsolutelyIrreducible(module) then
1867    #  Error("GModule is not absolutely irreducible.");
1868    #fi;
1869
1870    if SMTX.CentMat(module)=fail then
1871      Error("No CentMat component!");
1872    fi;
1873
1874    F:=SMTX.Field(module);
1875    R:=PolynomialRing(F);
1876    q:=Size(F);
1877    e :=SMTX.DegreeFieldExt(module);
1878    qe:=q^e - 1;
1879    minpol:=SMTX.CentMatMinPoly(module);
1880    # Factorise q^e - 1
1881    pp:=PrimePowersInt(qe);
1882    # We seek a generator of the field of order q^e - 1. In other words, a
1883    # polynomial genpol of degree e, which has multiplicative order q^e - 1
1884    # modulo minpol. We first try the polynomial x, which is the element we
1885    # have already. If this does not work, then we try random nonconstant
1886    # polynomials until we find one with the right order.
1887
1888    genpol:=Indeterminate(F);
1889
1890    looking:=true;
1891    while looking do
1892      if genpol <> minpol then
1893      okd:=FFPOrderKnownDividend(R, genpol, minpol, pp);
1894      if okd[1] * Order(One(F)*okd[2]) = qe then
1895          looking:=false;
1896      fi;
1897      fi;
1898      if looking then
1899          repeat
1900            genpol:=RandomPol(F, e,1);
1901          until DegreeOfUnivariateLaurentPolynomial(genpol) > 0;
1902          genpol:=StandardAssociate(R, genpol);
1903      fi;
1904    od;
1905    # Finally recalculate centmat and its minimal polynomial.
1906    centmat:=SMTX.CentMat(module);
1907    newcentmat:=Value(genpol, centmat,centmat^0);
1908    SMTX.SetFGCentMat(module, newcentmat);
1909    SMTX.SetFGCentMatMinPoly(module,MinimalPolynomialMatrixNC(F,newcentmat,1));
1910    # Ugh! That was very inefficient - should work out the min poly using
1911    # polynomials, but will sort that out if its ever needed.
1912  fi;
1913  return SMTX.FGCentMat(module);
1914end;
1915
1916###############################################################################
1917##
1918#F  SMTX.CollectedFactors( module ) . . find composition factors of a module
1919##
1920## 01/01/01 Try to deal more efficiently with large numbers of repeated
1921## small factors by using SMTX.Homomorphisms
1922##
1923## SMTX.CollectedFactors calls IsIrreducible repeatedly to find the
1924## composition factors of the GModule `module'. It also calls
1925## IsomorphismGModule to determine which are isomorphic.
1926## It returns a list [f1, f2, ..fr], where each fi is a list [m, n],
1927## where m is an irreducible composition factor of module, and n is the
1928## number of times it occurs in module.
1929##
1930SMTX_CollectedFactors:= function( module )
1931  local field,dim, factors, factorsout, queue, cmod, new,
1932      d, i, j, l, lf, q, smod, ds, homs, mat;
1933   if SMTX.IsMTXModule(module) = false then
1934      return Error("Argument is not a module.");
1935   fi;
1936
1937   dim:=SMTX.Dimension(module);
1938   field:= SMTX.Field(module);
1939   factors:=[];
1940   for i in [1..dim] do
1941      factors[i]:=[];
1942   od;
1943   # factors[i] will contain a list [f1, f2, ..., fr] of the composition factors
1944   # of module of dimension i. Each fi will have the form [m, n], where m is
1945   # the module, and n its multiplicity.
1946
1947   queue:=[module];
1948   # queue is the list of modules awaiting processing.
1949
1950   while Length(queue) > 0 do
1951      cmod:=Remove(queue);
1952      d:=SMTX.Dimension(cmod);
1953      Info(InfoMeatAxe,3,"Length of queue = ", Length(queue)+1, ", dim = ", d, ".");
1954
1955      if SMTX.IsIrreducible(cmod) then
1956         Info(InfoMeatAxe,2,"Irreducible: ");
1957         # module is irreducible. See if it is already on the list.
1958         new:=true;
1959         lf:=Length(factors[d]);
1960         i:=1;
1961         while new and i <= lf do
1962            if SMTX.IsEquivalent(factors[d][i][1], cmod) then
1963               new:=false;
1964               factors[d][i][2]:=factors[d][i][2] + 1;
1965            fi;
1966            i:=i + 1;
1967         od;
1968         if new then
1969            Info(InfoMeatAxe,2," new.");
1970            factors[d][lf + 1]:=[cmod, 1];
1971         else
1972            Info(InfoMeatAxe,2," old.");
1973         fi;
1974      else
1975         Info(InfoMeatAxe,2,"Reducible.");
1976         # module is reducible. Add sub- and quotient-modules to queue.
1977         q:=SMTX.InducedAction(cmod,
1978                  SMTX.Subbasis(cmod),3);
1979         smod:=q[1];
1980         ds:=SMTX.Dimension(smod);
1981         if ds < d/10 and SMTX.IsIrreducible(smod) then
1982           # Small dimensional submodule
1983           # test for repeated occurrences.
1984           homs:=SMTX.Homomorphisms( smod, cmod); # must have length >0
1985
1986           # build the submodule formed by their images
1987           mat:=Concatenation(homs);
1988           TriangulizeMat(mat);
1989           mat:=Filtered(mat,i->not IsZero(i));
1990           mat:=ImmutableMatrix(field,mat);
1991           if Length(mat)<cmod.dimension then
1992             # there is still some factor left
1993             Add(queue, SMTX.InducedActionFactorModule(cmod, mat));
1994           fi;
1995
1996           Info(InfoMeatAxe,2,
1997              "Small irreducible submodule X ",Length(homs),
1998              " subdim :",Length(mat)/smod.dimension,":");
1999
2000           # module is irreducible. See if it is already on the list.
2001           new:=true;
2002           lf:=Length(factors[ds]);
2003           i:=1;
2004           while new and i <= lf do
2005              if SMTX.IsEquivalent(factors[ds][i][1], smod) then
2006                 Info(InfoMeatAxe,2," old.");
2007                 new:=false;
2008                 factors[ds][i][2]:=factors[ds][i][2] +
2009                  Length(mat)/smod.dimension;
2010              fi;
2011              i:=i + 1;
2012           od;
2013           if new then
2014              Info(InfoMeatAxe,2," new.");
2015              factors[ds][lf + 1]:=[smod, Length(mat)/smod.dimension];
2016           fi;
2017
2018         else
2019           Add(queue, smod);
2020           Add(queue, q[2]);
2021         fi;
2022      fi;
2023   od;
2024
2025   # Now repack the sequence for output.
2026   l:=0;
2027   factorsout:=[];
2028   for i in [1..dim] do
2029      for j in [1..Length(factors[i])] do
2030         l:=l + 1;
2031         factorsout[l]:=factors[i][j];
2032      od;
2033   od;
2034
2035   return factorsout;
2036
2037end;
2038SMTX.CollectedFactors:=SMTX_CollectedFactors;
2039
2040SMTX.CompositionFactors:=function( module )
2041  if SMTX.IsIrreducible(module) then
2042    return [module];
2043  else
2044    module:=SMTX.InducedAction(module,
2045                  SMTX.Subbasis(module),3);
2046    return Concatenation(SMTX.CompositionFactors(module[1]),
2047                         SMTX.CompositionFactors(module[2]));
2048  fi;
2049end;
2050
2051###############################################################################
2052##
2053#F  SMTX.Distinguish( cf, i )  distinguish a composition factor of a module
2054##
2055## cf is assumed to be the output of a call to SMTX.CollectedFactors,
2056## and i is the number of one of the cf.
2057## Distinguish tries to find a group-algebra element for factor[i]
2058## which gives nullity zero when applied to all other cf.
2059## Once this is done, it is easy to find submodules containing this
2060## composition factor.
2061##
2062SMTX_Distinguish:=function( cf, i )
2063   local el, genpair, ngens, orig_ngens, mat, matsi, mats, M, idmat,
2064         dimi, dim, F, fac, sfac, p, q, oldp, found, extdeg, j, k,
2065         lcf, lf, x, y, wno, deg, trying, N, fact, R;
2066
2067   lcf:=Length(cf);
2068   ngens:=Length(cf[1][1].generators);
2069   orig_ngens:=ngens;
2070   F:=SMTX.Field(cf[1][1]);
2071   R:=PolynomialRing(F);
2072   matsi:=ShallowCopy(cf[i][1].generators);
2073   dimi:=SMTX.Dimension(cf[i][1]);
2074   idmat:=matsi[1]^0;
2075
2076   # First check that the existing nullspace has dim. 1 over centralising field.
2077   SMTX.GoodElementGModule(cf[i][1]);
2078
2079   # First see if the existing element is OK
2080   # Apply the alg. el. of factor i to every other factor and see if the
2081   # matrix is nonsingular.
2082   found:=true;
2083   el:=SMTX.AlgEl(cf[i][1]);
2084   fact:=SMTX.AlgElCharPolFac(cf[i][1]);
2085   for j in [1..lcf] do
2086      if j <> i and found then
2087         mats:=ShallowCopy(cf[j][1].generators);
2088         dim:=SMTX.Dimension(cf[j][1]);
2089         for genpair in el[1] do
2090            ngens:=ngens + 1;
2091            mats[ngens]:=mats[genpair[1]] * mats[genpair[2]];
2092         od;
2093         M:=ImmutableMatrix(F,Sum([1..ngens], k -> el[2][k] * mats[k]));
2094         ngens:=orig_ngens;
2095         mat:=Value(fact, M, M^0);
2096         if RankMat(mat) < dim then
2097            found:=false;
2098            Info(InfoMeatAxe,2,"Current element failed on factor ", j);
2099         fi;
2100      fi;
2101   od;
2102
2103   if found then
2104      Info(InfoMeatAxe,2,"Current element worked.");
2105      return;
2106   fi;
2107
2108   # That didn't work, so we have to try new random elements.
2109   wno:=0;
2110   el:=[]; el[1]:=[];
2111   extdeg:=SMTX.DegreeFieldExt(cf[i][1]);
2112
2113   while found = false do
2114      Info(InfoMeatAxe,2,"Trying new one.");
2115      wno:=wno + 1;
2116      # Add a new generator if there are less than 8 or if wno mod 10=0.
2117      if  ngens<8 or wno mod 10 = 0 then
2118         x:=Random(1, ngens);
2119         y:=x;
2120         while y = x and ngens > 1 do y:=Random(1, ngens); od;
2121         Add(el[1], [x, y]);
2122         ngens:=ngens + 1;
2123         matsi[ngens]:=matsi[x] * matsi[y];
2124      fi;
2125      # Now take the new random element
2126      el[2]:=[];
2127      for j in [1..ngens] do el[2][j]:=Random(F); od;
2128      # First evaluate on cf[i][1].
2129      M:=ImmutableMatrix(F,Sum([1..ngens], k ->  el[2][k] * matsi[k]));
2130      p:=CharacteristicPolynomialMatrixNC(F,M,1);
2131      # That is necessary in case p is defined over a smaller field that F.
2132      oldp:=p;
2133      # extract irreducible factors
2134      deg:=0;
2135      fac:=[];
2136      trying:=true;
2137      while deg <= extdeg and trying do
2138         repeat
2139            deg:=deg + 1;
2140            if deg > extdeg then
2141               fac:=[p];
2142            else
2143               fac:=Factors(R, p: factoroptions:=rec(onlydegs:=[deg]));
2144               fac:=Filtered(fac,i->DegreeOfLaurentPolynomial(i)<=deg);
2145               sfac:=Set(fac);
2146            fi;
2147         until fac <> [];
2148         lf:=Length(fac);
2149         if trying and deg <= extdeg then
2150            j:=1;
2151            while j <= lf and trying do
2152               mat:=Value(fac[j], M,idmat);
2153               N:=NullspaceMat(mat);
2154               if Length(N) = extdeg then
2155                  trying:=false;
2156                  SMTX.SetAlgEl(cf[i][1], el);
2157                  SMTX.SetAlgElMat(cf[i][1], M);
2158                  SMTX.SetAlgElCharPol(cf[i][1], oldp);
2159                  SMTX.SetAlgElCharPolFac(cf[i][1], fac[j]);
2160                  SMTX.SetAlgElNullspaceVec(cf[i][1], N[1]);
2161               fi;
2162               j:=j + 1;
2163            od;
2164         fi;
2165
2166         if trying then
2167            for q in fac do
2168               p:=Quotient(R, p, q);
2169            od;
2170         fi;
2171      od;
2172
2173      # Now see if it works against the other factors of cf
2174      if trying = false then
2175         Info(InfoMeatAxe,2,"Found one.");
2176         found:=true;
2177         fact:=SMTX.AlgElCharPolFac(cf[i][1]);
2178         # Apply the alg. el. of factor i to every other factor and
2179         # see if the matrix is nonsingular.
2180         for j in [1..lcf] do
2181            if j <> i and found then
2182               mats:=ShallowCopy(cf[j][1].generators);
2183               dim:=SMTX.Dimension(cf[j][1]);
2184               ngens:=orig_ngens;
2185               for genpair in el[1] do
2186                  ngens:=ngens + 1;
2187                  mats[ngens]:=mats[genpair[1]] * mats[genpair[2]];
2188               od;
2189               M:=ImmutableMatrix(F,Sum([1..ngens], k -> el[2][k] * mats[k]));
2190               mat:=Value(fact, M, M^0);
2191               if RankMat(mat) < dim then
2192                  found:=false;
2193                  Info(InfoMeatAxe,2,"Failed on factor ", j);
2194               fi;
2195            fi;
2196         od;
2197      fi;
2198      if found then
2199         Info(InfoMeatAxe,2,"It worked!");
2200      fi;
2201   od;
2202
2203end;
2204SMTX.Distinguish:=SMTX_Distinguish;
2205
2206###############################################################################
2207##
2208#F  SMTX.MinimalSubGModule( module, cf, i ) . .  find minimal submodule
2209##                                     containing a given composition factor.
2210##
2211## cf is assumed to be the output of a call to SMTX.CollectedFactors,
2212## and i is the number of one of the cf.
2213## It is assumed that SMTX.Distinguish(cf, i) has already been called.
2214## A basis of a minimal submodule of module containing the composition factor
2215## cf[i][1] is calculated and returned - i.e. if cf[i][2] = 1.
2216##
2217SMTX_MinimalSubGModule:=function( module, cf, i )
2218   local el, genpair, ngens, orig_ngens, mat, mats, M, dim, F,
2219         k, N, fact;
2220
2221   if SMTX.IsMTXModule(module) = false then
2222      return Error("First argument is not a module.");
2223   fi;
2224
2225   ngens:=Length(module.generators);
2226   orig_ngens:=ngens;
2227   F:=SMTX.Field(module);
2228
2229   # Apply the alg. el. of factor i to module
2230   el:=SMTX.AlgEl(cf[i][1]);
2231   mats:=ShallowCopy(module.generators);
2232   dim:=SMTX.Dimension(module);
2233   for genpair in el[1] do
2234      ngens:=ngens + 1;
2235      mats[ngens]:=mats[genpair[1]] * mats[genpair[2]];
2236   od;
2237
2238   M:=ImmutableMatrix(F,Sum([1..ngens], k -> el[2][k] * mats[k]));
2239   # Now throw away extra generators of module
2240   for k in [orig_ngens + 1..ngens] do
2241      Unbind(mats[k]);
2242   od;
2243   ngens:=orig_ngens;
2244   fact:=SMTX.AlgElCharPolFac(cf[i][1]);
2245   mat:=Value(fact, M,M^0);
2246   N:=NullspaceMat(mat);
2247   return SMTX.SpinnedBasis(N[1], mats,F, ngens);
2248
2249end;
2250SMTX.MinimalSubGModule:=SMTX_MinimalSubGModule;
2251
2252
2253#############################################################################
2254##
2255#F  SMTX.Isomomorphism(module1, module2) . . . .
2256##  decide whether two irreducible modules are isomorphic.
2257##
2258## If the 2 modules are not isomorphic, this function returns false;
2259## if they are isomorphic it returns the matrix B, whose rows form the
2260## basis of module2  which is the image of the standard basis for module1.
2261## Thus if X and Y are corresponding matrices in the generating sets
2262## for module1 and module2 respectively, Y = B^-1XB
2263## It is assumed that the same group acts on both modules.
2264## Otherwise who knows what will happen?
2265##
2266SMTX_IsomorphismComp:=function(module1, module2, action)
2267   local matrices, matrices1, matrices2, F, R, dim, swapmodule, genpair,
2268         swapped, orig_ngens, i, j, el, p, fac, ngens, M, mat, v1, v2, v,
2269         N, basis, basis1, basis2;
2270
2271   if SMTX.IsMTXModule(module1) = false then
2272      Error("Argument is not a module.");
2273   elif SMTX.IsMTXModule(module2) = false then
2274      Error("Argument is not a module.");
2275   elif SMTX.Field(module1) <> SMTX.Field(module2) then
2276      Error("GModules are defined over different fields.");
2277   fi;
2278
2279   swapped:=false;
2280   if not SMTX.HasIsIrreducible(module1) then
2281      if not SMTX.HasIsIrreducible(module2) then
2282         Error("Neither module is known to be irreducible.");
2283      else
2284         # The second module is known to be irreducible, so swap arguments.
2285         swapmodule:=module2; module2:=module1; module1:=swapmodule;
2286         swapped:=true;
2287         Info(InfoMeatAxe,2,"Second module is irreducible. Swap them round.");
2288      fi;
2289   fi;
2290
2291   # At this stage, module1 is known to be irreducible
2292   dim:=SMTX.Dimension(module1);
2293   if dim <> SMTX.Dimension(module2) then
2294      Info(InfoMeatAxe,2,"GModules have different dimensions.");
2295      return fail;
2296   fi;
2297   F:=SMTX.Field(module1);
2298   R:=PolynomialRing(F);
2299
2300   # First we must check that our nullspace is 1-dimensional over the
2301   # centralizing field.
2302
2303   Info(InfoMeatAxe,2,
2304        "Checking nullspace 1-dimensional over centralising field.");
2305   SMTX.GoodElementGModule(module1);
2306   matrices1:=module1.generators;
2307   matrices2:=ShallowCopy(module2.generators);
2308   ngens:=Length(matrices1);
2309   orig_ngens:=ngens;
2310   if ngens <> Length(matrices2) then
2311      Error("GModules have different numbers of defining matrices.");
2312   fi;
2313
2314   # Now we calculate the element in the group algebra of module2 that
2315   # corresponds to that in module1. This is done using the AlgEl flag
2316   # for module1. We first extend the generating set in the same way as
2317   # we did for module1, and then calculate the group alg. element as
2318   # a linear sum in the generators.
2319
2320   Info(InfoMeatAxe,2,"Extending generating set for second module.");
2321   el:=SMTX.AlgEl(module1);
2322   for genpair in el[1] do
2323      ngens:=ngens + 1;
2324      matrices2[ngens]:=matrices2[genpair[1]] * matrices2[genpair[2]];
2325   od;
2326   M:=ImmutableMatrix(F,Sum([1..ngens], i -> el[2][i] * matrices2[i]));
2327   # Having done that, we no longer want the extra generators of module2,
2328   # so we throw them away again.
2329   for i in [orig_ngens + 1..ngens] do
2330      Unbind(matrices2[i]);
2331   od;
2332
2333   Info(InfoMeatAxe,2,
2334        "Calculating characteristic polynomial for second module.");
2335   p:=CharacteristicPolynomialMatrixNC(F,M,1);
2336   if p <> SMTX.AlgElCharPol(module1) then
2337      Info(InfoMeatAxe,2,"Characteristic polynomial different.");
2338      return fail;
2339   fi;
2340   fac:=SMTX.AlgElCharPolFac(module1);
2341   mat:=Value(fac, M,M^0);
2342   Info(InfoMeatAxe,2,"Calculating nullspace for second module.");
2343   N:=NullspaceMat(mat);
2344   if Length(N) <> SMTX.AlgElNullspaceDimension(module1) then
2345      Info(InfoMeatAxe,2,"Null space dimensions different.");
2346      return fail;
2347   fi;
2348
2349   # That concludes the easy tests for nonisomorphism. Now we must proceed
2350   # to spin up. We first form the direct sum of the generating matrices.
2351   Info(InfoMeatAxe,2,"Spinning up in direct sum.");
2352   matrices:=SMTX.MatrixSum(matrices1, matrices2);
2353   v1:=SMTX.AlgElNullspaceVec(module1);
2354   v2:=N[1];
2355   v:=Concatenation(v1, v2);
2356   basis:=SMTX.SpinnedBasis(v, matrices, F);
2357   if Length(basis) = dim then
2358      if action<>true then
2359        return true;
2360      fi;
2361      basis1:=[]; basis2:=[];
2362      for i in [1..dim] do
2363         basis1[i]:=[]; basis2[i]:=[];
2364         for j in [1..dim] do
2365            basis1[i][j]:=basis[i][j];
2366            basis2[i][j]:=basis[i][j + dim];
2367         od;
2368      od;
2369      if swapped then
2370         return basis2^-1 * basis1;
2371      else
2372         return basis1^-1 * basis2;
2373      fi;
2374   else
2375      return fail;
2376   fi;
2377end;
2378SMTX.IsomorphismComp:=SMTX_IsomorphismComp;
2379
2380SMTX.IsomorphismIrred:=function(module1,module2)
2381  return SMTX.IsomorphismComp(module1,module2,true);
2382end;
2383
2384SMTX.Isomorphism:=SMTX.IsomorphismIrred;
2385
2386SMTX.IsEquivalent:=function(module1,module2)
2387  return SMTX.IsomorphismComp(module1,module2,false)<>fail;
2388end;
2389
2390#############################################################################
2391##
2392#F  SMTX.MatrixSum(matrices1, matrices2) direct sum of two lists of matrices
2393##
2394SMTX_MatrixSum:=function(matrices1, matrices2)
2395   local matrices, nmats, i;
2396   matrices:=[];
2397   nmats:=Length(matrices1);
2398   for i in [1..nmats] do
2399      matrices[i]:=DirectSumMat(matrices1[i],matrices2[i]);
2400   od;
2401
2402   return  matrices;
2403end;
2404SMTX.MatrixSum:=SMTX_MatrixSum;
2405
2406
2407#############################################################################
2408##
2409#F  SMTX.Homomorphisms( m1, m2) . . . . homomorphisms from an irreducible
2410##                         . . . GModule to an arbitrary GModule
2411##
2412## It is assumed that m1 is a module that has been proved irreducible
2413##  (using IsIrreducible), and m2 is an arbitrary module for the same group.
2414## A basis of the space of G-homomorphisms from m1 to m2 is returned.
2415## Each homomorphism is given as a list of base images.
2416##
2417SMTX_Homomorphisms:= function(m1, m2)
2418
2419   local F, ngens, orig_ngens, mats1, mats2, dim1, dim2, m1bas, imbases,
2420         el, genpair, fac, mat, N, imlen, subdim, leadpos, vec, imvecs,
2421         numrels, rels, leadposrels, newrels, bno, genno, colno, rowno,
2422         zero, looking, ans, i, j, k;
2423
2424   if not SMTX.IsMTXModule(m1) then
2425      return Error("First argument is not a module.");
2426   elif not SMTX.IsIrreducible(m1) then
2427      return Error("First module is not known to be irreducible.");
2428   fi;
2429
2430   if not SMTX.IsMTXModule(m2) then
2431      return Error("Second argument is not a module.");
2432   fi;
2433   mats1:=m1.generators;
2434   mats2:=ShallowCopy(m2.generators);
2435   ngens:=Length(mats1);
2436   if ngens <> Length(mats2) then
2437      return Error("GModules have different numbers of generators.");
2438   fi;
2439
2440   F:=SMTX.Field(m1);
2441   if F <> SMTX.Field(m2) then
2442      return Error("GModules are defined over different fields.");
2443   fi;
2444   zero:=Zero(F);
2445
2446   dim1:=SMTX.Dimension(m1); dim2:=SMTX.Dimension(m2);
2447
2448   if dim1=1 then
2449     # m1 is 1-dimensional -- eigenspace intersection
2450     el:=List([1..Length(m1.generators)],x->NullspaceMat(m2.generators[x]-m1.generators[x][1][1]*m2.generators[x]^0));
2451
2452     imvecs:=el[1];
2453     for j in [2..Length(el)] do
2454       imvecs:=SumIntersectionMat(imvecs,el[j])[2];
2455     od;
2456     return List(imvecs,x->ImmutableMatrix(m1.field,[x]));
2457   fi;
2458
2459   m1bas:=[];
2460   m1bas[1]:= SMTX.AlgElNullspaceVec(m1);
2461
2462   # In any homomorphism from m1 to m2, the vector in the nullspace of the
2463   # algebraic element that was used to prove irreducibility (which is now
2464   # m1bas[1]) must map onto a vector in the nullspace of the same algebraic
2465   # element evaluated in m2. We therefore calculate this nullspaces, and
2466   # store a basis in imbases.
2467
2468   Info(InfoMeatAxe,2,"Extending generating set for second module.");
2469   orig_ngens:=ngens;
2470   el:=SMTX.AlgEl(m1);
2471   for genpair in el[1] do
2472      ngens:=ngens + 1;
2473      mats2[ngens]:=mats2[genpair[1]] * mats2[genpair[2]];
2474   od;
2475   mat:=ImmutableMatrix(F,Sum([1..ngens], i -> el[2][i] * mats2[i]));
2476   # Having done that, we no longer want the extra generators of m2,
2477   # so we throw them away again.
2478   for i in [orig_ngens + 1..ngens] do
2479      Unbind(mats2[i]);
2480   od;
2481   ngens:=orig_ngens;
2482
2483   fac:=SMTX.AlgElCharPolFac(m1);
2484   mat:=Value(fac, mat,mat^0);
2485   Info(InfoMeatAxe,2,"Calculating nullspace for second module.");
2486   N:=NullspaceMat(mat);
2487   imlen:=Length(N);
2488   Info(InfoMeatAxe,2,"Dimension = ", imlen, ".");
2489   if imlen = 0 then
2490      return [];
2491   fi;
2492
2493   imbases:=List(N, vec -> [vec]);
2494
2495   # Now the main algorithm starts. We are going to spin the vectors in m1bas
2496   # under the action of the module generators, norming as we go. Every
2497   # operation that we perform on m1bas will also be performed on each of the
2498   # vectors in  imbas[1], ..., imbas[imlen].
2499   # When we find a vector that norms to zero in m1bas, then the image of this
2500   # under a homomorphism must be zero. This leads to a linear relation
2501   # amongst some vectors in imbas. We store up such relations, echelonizing as
2502   # we go. At the end, if we have numrels subch independent relations, then
2503   # there will be imlen - numrels independent homomorphisms from m1 to m2,
2504   # which we can then calculate.
2505
2506   subdim:=1; # the dimension of module spanned by m1bas
2507   numrels:=0;
2508   rels:=[];
2509
2510   # leadpos[j] will be the position of the first nonzero entry in m1bas[j]
2511   leadpos:=[];
2512   vec:=m1bas[1];
2513   j:=1;
2514   while j <= dim1 and vec[j] = zero do j:=j + 1; od;
2515   leadpos[1]:=j;
2516   k:=vec[j]^-1;
2517   m1bas[1]:=k * vec;
2518   for i in [1..imlen] do
2519      imbases[i][1]:=k * imbases[i][1];
2520   od;
2521
2522   leadposrels:=[];
2523   # This will play the same role as leadpos but for the relation matrix.
2524   Info(InfoMeatAxe,2,"Starting spinning.");
2525   bno:=1;
2526   while bno <= subdim do
2527      for genno in [1..ngens] do
2528         # apply generator no. genno to submodule generator bno
2529         vec:=m1bas[bno] * mats1[genno];
2530         # and do the same to the images
2531         imvecs:=[];
2532         for i in [1..imlen] do
2533            imvecs[i]:=imbases[i][bno] * mats2[genno];
2534         od;
2535         # try to express w in terms of existing submodule generators
2536         # make same changes to images
2537         j:=1;
2538         for  j in [1..subdim] do
2539            k:=vec[leadpos[j]];
2540            if k <> zero then
2541               vec:=vec - k * m1bas[j];
2542               for i in [1..imlen] do
2543                  imvecs[i]:=imvecs[i] - k * imbases[i][j];
2544               od;
2545            fi;
2546         od;
2547
2548         j:=1;
2549         while j <= dim1 and vec[j] = zero do j:=j + 1; od;
2550         if j <= dim1 then
2551            # we have found a new generator of the submodule
2552            subdim:=subdim + 1;
2553            leadpos[subdim]:=j;
2554            k:=vec[j]^-1;
2555            m1bas[subdim]:=k * vec;
2556            for i in [1..imlen] do
2557               imbases[i][subdim]:=k * imvecs[i];
2558            od;
2559         else
2560            # vec has reduced to zero. We get relations among the imvecs.
2561            # (these are given by the transpose of imvec)
2562            # reduce these against any existing relations.
2563            newrels:=TransposedMat(imvecs);
2564            for i in [1..Length(newrels)] do
2565               vec:=newrels[i];
2566               for j in [1..numrels] do
2567                  k:=vec[leadposrels[j]];
2568                  if k <> zero then
2569                     vec:=vec - k * rels[j];
2570                  fi;
2571               od;
2572               j:=1;
2573               while j <= imlen and vec[j] = zero do j:=j + 1; od;
2574               if j <= imlen then
2575                  # we have a new relation
2576                  numrels:=numrels + 1;
2577                  # if we have imlen relations, there can be no homomorphisms
2578                  # so we might as well give up immediately
2579                  if numrels = imlen then
2580                     return [];
2581                  fi;
2582                  k:=vec[j]^-1;
2583                  rels[numrels]:=k * vec;
2584                  leadposrels[numrels]:=j;
2585               fi;
2586            od;
2587         fi;
2588      od;
2589      bno:=bno + 1;
2590   od;
2591
2592   # That concludes the spinning. Now we do row operations on the im1bas to
2593   # make it the identity, and do the same operations to the imvecs.
2594   # Then the homomorphisms we output will be the basis images.
2595   Info(InfoMeatAxe,2,"Done. Reducing spun up basis.");
2596
2597   for colno in [1..dim1] do
2598      rowno:=colno;
2599      looking:=true;
2600      while rowno <= dim1 and looking do
2601         if m1bas[rowno][colno] <> zero then
2602            looking:=false;
2603            if rowno <> colno then
2604               # swap rows rowno and colno
2605               vec:=m1bas[rowno]; m1bas[rowno]:=m1bas[colno];
2606               m1bas[colno]:=vec;
2607               # and of course the same in the images
2608               for i in [1..imlen] do
2609                  vec:=imbases[i][rowno];
2610                  imbases[i][rowno]:=imbases[i][colno];
2611                  imbases[i][colno]:=vec;
2612               od;
2613            fi;
2614            # and then clear remainder of column
2615            for j in [1..dim1] do
2616               if j <> colno and m1bas[j][colno] <> zero then
2617                  k:=m1bas[j][colno];
2618                  m1bas[j]:=m1bas[j] - k * m1bas[colno];
2619                  for i in [1..imlen] do
2620                     imbases[i][j]:=imbases[i][j] - k * imbases[i][colno];
2621                  od;
2622               fi;
2623            od;
2624         fi;
2625         rowno:=rowno + 1;
2626      od;
2627   od;
2628
2629   # Now we are ready to compute and output the linearly independent
2630   # homomorphisms.  The coefficients for the solution are given by
2631   # the basis elements of the nullspace of the transpose of rels.
2632
2633   Info(InfoMeatAxe,2,"Done. Calculating homomorphisms.");
2634   if rels = [] then
2635      rels:=NullMat(imlen, 1, F);
2636   else
2637      rels:=TransposedMat(rels);
2638   fi;
2639   N:=NullspaceMat(rels);
2640   ans:=[];
2641   for vec in N do
2642      mat:=ImmutableMatrix(F, Sum([1..imlen], i -> vec[i] * imbases[i]));
2643      Add(ans, mat);
2644   od;
2645
2646   return ans;
2647end;
2648SMTX.Homomorphisms:=SMTX_Homomorphisms;
2649
2650#############################################################################
2651##
2652#F  SMTX.SortHomGModule( m1, m2, homs)  . . sort output of HomGModule
2653##                                           according to their images
2654##
2655## It is assumed that m1 is a module that has been proved irreducible
2656## (using IsIrreducible), and m2 is an arbitrary module for the same group,
2657## and that homs is the output of a call HomGModule(m1, m2).
2658## Let e be the degree of the centralising field of m1.
2659## If e = 1 then SMTX.SortHomGModule does nothing. If e > 1, then it replaces
2660## the basis contained in homs by a new basis arranged in the form
2661## b11, b12, ..., b1e, b21, b22, ...b2e, ..., br1, br2, ...bre,  where each
2662## block of  e  adjacent basis vectors are all equivalent under the
2663## centralising field of m1, and so they all have the same image in  m2.
2664## A complete list of the distinct images can then be obtained with a call
2665## to DistinctIms(m1, m2, homs).
2666##
2667SMTX_SortHomGModule:=function(m1, m2, homs)
2668local e, F, ngens, mats1, mats2, dim1, dim2, centmat, fullimbas, oldhoms,
2669      homno, dimhoms, newdim, subdim, leadpos, vec, nexthom,
2670      i, j, k, zero;
2671
2672   if SMTX.IsAbsolutelyIrreducible(m1) then return; fi;
2673
2674   e:=SMTX.DegreeFieldExt(m1);
2675   F:=SMTX.Field(m1);
2676   zero:=Zero(F);
2677
2678   mats1:=m1.generators;  mats2:=m2.generators;
2679   dim1:=SMTX.Dimension(m1);  dim2:=SMTX.Dimension(m2);
2680   ngens:=Length(mats1);
2681   centmat:=SMTX.CentMat(m1);
2682
2683   fullimbas:=[];
2684   subdim:=0;
2685   leadpos:=[];
2686
2687   # fullimbas will contain an echelonised basis for the submodule of m2
2688   # generated by all images of the basis vectors of hom that we have found
2689   # so far; subdim is its length.
2690
2691   # We go through the existing basis of homs.
2692   # For each hom in the basis, we first check whether the first vector in
2693   # the image  of hom is in the space spanned by fullimbas.
2694   # If so, we reject hom. If not, then hom is adjoined to the new
2695   # basis of homs, as are the other e-1 linearly independent homomorphisms
2696   # that are equivalent to hom by a multiplication by centmat. The
2697   # resulting block of e homomorphisms all have the same image in m2.
2698
2699   # first make a copy of homs.
2700
2701   oldhoms:=ShallowCopy(homs);
2702   dimhoms:=Length(homs);
2703
2704   homno:=0; newdim:=0;
2705
2706   while homno < dimhoms and newdim < dimhoms do
2707      homno:=homno + 1;
2708      nexthom:=oldhoms[homno];
2709      vec:=nexthom[1];
2710
2711      # Now check whether vec is in existing submodule spanned by fullimbas
2712      j:=1;
2713      for j in [1..subdim] do
2714         k:=vec[leadpos[j]];
2715         if k <> zero then
2716            vec:=vec - k * fullimbas[j];
2717         fi;
2718      od;
2719
2720      j:=1;
2721      while j <= dim2 and vec[j] = zero do j:=j + 1; od;
2722
2723      if j <= dim2 then
2724         # vec is not in the image, so we adjoin this homomorphism to the list;
2725         # first adjoin vec and all other basis vectors in the image to fullimbas
2726         subdim:=subdim + 1;
2727         leadpos[subdim]:=j;
2728         k:=vec[j]^-1;
2729         fullimbas[subdim]:=k * vec;
2730         for i in [2..dim1] do
2731            vec:=nexthom[i];
2732            j:=1;
2733            for  j in [1..subdim] do
2734               k:=vec[leadpos[j]];
2735               if k <> zero then
2736                  vec:=vec - k * fullimbas[j];
2737               fi;
2738            od;
2739
2740            j:=1;
2741            while j <= dim2 and vec[j] = zero do j:=j + 1; od;
2742            subdim:=subdim + 1;
2743            leadpos[subdim]:=j;
2744            k:=vec[j]^-1;
2745            fullimbas[subdim]:=k * vec;
2746         od;
2747
2748         newdim:=newdim + 1;
2749         homs[newdim]:=nexthom;
2750
2751         # Now add on the other e - 1 homomorphisms equivalent to
2752         # newhom by centmat.
2753         for k in [1..e - 1] do
2754            nexthom:=centmat * nexthom;
2755            newdim:=newdim + 1;
2756            homs[newdim]:=nexthom;
2757         od;
2758      fi;
2759   od;
2760
2761end;
2762SMTX.SortHomGModule:=SMTX_SortHomGModule;
2763
2764#############################################################################
2765##
2766#F  SMTX.Homomorphism(module1,module2,mat) . . . define a module homorphism
2767##
2768##  module1 and module2 should be meataxe modules of dimensions m and n
2769##  over the same algebra, and mat an mXn matrix over the field of
2770##  the modules that defines a homomorphism module1 -> module2, where
2771##  the i-th row of mat gives the image in module2 of the i-th basis
2772##  vector of module1.
2773##  It is checked whether mat really does define a homomorphism.
2774##  If, so then the corresponding vector space homomorphism from the underlying
2775##  row space of module1 to that of module2 is returned. This can be used
2776##  for computing images, kernel, preimages, etc.
2777
2778SMTX_Homomorphism:=function(module1, module2, mat)
2779  local F, gens1, gens2, ng, dim1, dim2, i, j;
2780  F:=SMTX.Field(module1);
2781  if F <> SMTX.Field(module2) then
2782    Error("Modules are over different fields");
2783  fi;
2784  gens1:=SMTX.Generators(module1); gens2:=SMTX.Generators(module2);
2785  dim1:=SMTX.Dimension(module1); dim2:=SMTX.Dimension(module2);
2786  ng:=Length(gens1);
2787  if ng <> Length(gens2) then
2788    Error("Modules are not over the same algebra");
2789  fi;
2790  if Length(mat) <> dim1 or Length(mat[1]) <> dim2 then
2791    Error("matrix has wrong size for a homomorphism");
2792  fi;
2793  # Check if it is a homorphism
2794  mat:=ImmutableMatrix(F,mat);
2795  for i in [1..ng] do
2796    for j in [1..dim1] do
2797      if gens1[i][j] * mat <> mat[j] * gens2[i] then
2798        Print(i,j,"\n");
2799        Error("matrix does not define a homomorphism");
2800      fi;
2801    od;
2802  od;
2803  return LeftModuleHomomorphismByImages(FullRowSpace(F,dim1),
2804                          FullRowSpace(F,dim2),IdentityMat(dim1,F),mat);
2805end;
2806SMTX.Homomorphism:=SMTX_Homomorphism;
2807
2808#############################################################################
2809##
2810#F SMTX.MinimalSubGModules(m1, m2, [max]) . .
2811## minimal submodules of m2 isomorphic to m1
2812##
2813## It is assumed that m1 is a module that has been proved irreducible
2814##  (using IsIrreducible), and m2 is an arbitrary module for the same group.
2815## MinimalSubGModules computes and outputs a list of normed bases for all of the
2816## distinct minimal submodules of m2 that are isomorphic to m1.
2817## max is an optional maximal number - if the total number of submodules
2818## exceeds max, then the procedure aborts.
2819## First HomGModule is called and then SMTX.SortHomGModule to get a basis for
2820## the homomorphisms from m1 to m2 in the correct order.
2821## It is then easy to write down the list of distinct images.
2822##
2823SMTX_MinimalSubGModules:=function(arg)
2824
2825   local m1, m2, max, e, homs, coeff,  dimhom, edimhom, F, elF, q,
2826         submodules, sub, adno, more, count, sr, er, i, j, k;
2827
2828   if Length(arg) < 2 or Length(arg) > 3 then
2829      Error("Number of arguments to MinimalSubGModules must be 2 or 3.");
2830   fi;
2831
2832   m1:=arg[1]; m2:=arg[2];
2833   if Length(arg) = 2 then max:=0; else max:=arg[3]; fi;
2834
2835   Info(InfoMeatAxe,2,"Calculating homomorphisms from m1 to m2.");
2836   homs:=SMTX.Homomorphisms(m1, m2);
2837   Info(InfoMeatAxe,2,"Sorting them.");
2838   SMTX.SortHomGModule(m1, m2, homs);
2839
2840   F:=SMTX.Field(m1);
2841   e:=SMTX.DegreeFieldExt(m1);
2842   dimhom:=Length(homs);
2843   edimhom:=dimhom / e;
2844   submodules:=[];
2845   count:=0;
2846   elF:=AsList(F);
2847   q:=Length(elF);
2848   coeff:=ListWithIdenticalEntries(dimhom,1);
2849
2850   # coeff[i] will be an integer in the range [1..q] corresponding to the
2851   # field element elF[coeff[i]].
2852   # Each submodule will be calculated as the image of the homomorphism
2853   # elF[coeff[1]] * homs[1] +...+  elF[coeff[dimhom]] * homs[dimhom]
2854   # for appropriate field elements elF[coeff[i]].
2855   # We get each distinct submodule
2856   # exactly once by making the first nonzero elF[coeff[i]] to be 1,
2857   # and all other elF[coeff[i]]'s in that block equal to zero.
2858
2859   Info(InfoMeatAxe,2,"Done. Calculating submodules.");
2860
2861   for i in Reversed([1..edimhom]) do
2862      j:=e * (i - 1) + 1;
2863      coeff[j]:=2;  # giving field element 1.
2864      for k in [j + 1..dimhom] do coeff[k]:=1; od; # field element 0.
2865      sr:=j + e; er:=dimhom;
2866      # coeff[i] for i in [sr..er] ranges over all field elements.
2867
2868      more:=true;
2869      adno:=er;
2870      while more do
2871         count:=count + 1;
2872         if max > 0 and count > max then
2873            Info(InfoMeatAxe,2,"Number of submodules exceeds ", max,
2874                 ". Aborting.");
2875            return submodules;
2876         fi;
2877
2878         # Calculate the next submodule
2879         sub:=homs[j];
2880         for k in [sr..er] do
2881            sub:=sub + elF[coeff[k]] * homs[k];
2882         od;
2883         sub:=TriangulizedMat(sub);
2884         Add(submodules, ImmutableMatrix(F,sub));
2885
2886         # Move on to next set of coefficients if any
2887         while adno >= sr and coeff[adno]=q do
2888            coeff[adno]:=1;
2889            adno:=adno - 1;
2890         od;
2891         if adno < sr then
2892            more:=false;
2893         else
2894            coeff[adno]:=coeff[adno] + 1;
2895            adno:=er;
2896         fi;
2897      od;
2898
2899   od;
2900
2901   return submodules;
2902
2903end;
2904SMTX.MinimalSubGModules:=SMTX_MinimalSubGModules;
2905
2906SMTX_BasesCompositionSeries:=function(m)
2907local q,b,s,ser,queue,F,one,mats,mo;
2908mats:=m.generators;
2909  SMTX.SetSmashRecord(m,0);
2910  F:=SMTX.Field(m);
2911  one:=One(F);
2912  b:=IdentityMat(SMTX.Dimension(m),F);
2913  b:=ImmutableMatrix(F,b);
2914  # denombasis: basis of the kernel
2915  m.smashMeataxe.denombasis:=[];
2916  # fakbasis: Urbilder der Basis, bzgl. derer csbasis angegeben wird
2917  # the first <dimension> vectors of <fakbasis> are the right ones.
2918  m.smashMeataxe.fakbasis:=b;
2919
2920  ser:=[[]];
2921  queue:=[m];
2922  while Length(queue)>0 do
2923    m:=Remove(queue);
2924    if SMTX.IsIrreducible(m) then
2925      mo:=m;
2926      Info(InfoMeatAxe,3,SMTX.Dimension(m)," ",
2927                         Length(m.smashMeataxe.denombasis));
2928      m:=Concatenation(
2929        m.smashMeataxe.denombasis,
2930        m.smashMeataxe.fakbasis{[1..SMTX.Dimension(m)]});
2931      m:=List(m,ShallowCopy);
2932      TriangulizeMat(m);
2933      m:=ImmutableMatrix(F,m);
2934      Add(ser,m);
2935    else
2936      b:=SMTX.Subbasis(m);
2937      s:=SMTX.InducedAction(m,b,3);
2938      q:=s[2];
2939      b:=s[3];
2940      s:=s[1];
2941      SMTX.SetSmashRecord(s,0);
2942      SMTX.SetSmashRecord(q,0);
2943      Info(InfoMeatAxe,1,"chopped ",SMTX.Dimension(s),"\\", SMTX.Dimension(q));
2944      s.smashMeataxe.denombasis:=m.smashMeataxe.denombasis;
2945      s.smashMeataxe.fakbasis:=b*m.smashMeataxe.fakbasis;
2946
2947      q.smashMeataxe.denombasis:=Concatenation(
2948        m.smashMeataxe.denombasis,
2949        s.smashMeataxe.fakbasis{[1..s.dimension]});
2950      q.smashMeataxe.fakbasis:=
2951        s.smashMeataxe.fakbasis{[s.dimension+1..Length(b)]};
2952
2953      Add(queue,s);
2954      Add(queue,q);
2955    fi;
2956  od;
2957  SortBy(ser,Length);
2958  return ser;
2959end;
2960SMTX.BasesCompositionSeries:=SMTX_BasesCompositionSeries;
2961
2962# composition series with small steps
2963SMTX_BasesCSSmallDimUp:=function(m)
2964local cf,F,dim,b,den,sub,i,s,q,found,qb;
2965  cf:=List(SMTX.CollectedFactors(m),x->[x[1],x[2]]); # so we can overwrite
2966  SortBy(cf,x->x[1].dimension);
2967  dim:=SMTX.Dimension(m);
2968  F:=SMTX.Field(m);
2969  b:=IdentityMat(dim,F);
2970  b:=ImmutableMatrix(F,b);
2971  den:=0;
2972  sub:=[[]];
2973  q:=m;
2974  while den<dim do
2975    i:=1;
2976    found:=false;
2977    while i<=Length(cf) and found=false do
2978      if cf[i][2]>0 then # can we still get this module?
2979        s:=MTX.Homomorphisms(cf[i][1],q);
2980        if Length(s)>0 then
2981          # found one
2982          cf[i][2]:=cf[i][2]-1;
2983          found:=true;
2984          s:=s[1];
2985          if Length(s)=q.dimension then
2986            # on top
2987            Add(sub,TriangulizedMat(b));
2988            den:=dim;
2989          else
2990            s:=TriangulizedMat(s);
2991            qb:=b{[den+1..Length(b)]}; # basis
2992            qb:=List(s,x->x*qb);
2993            qb:=Concatenation(b{[1..den]},qb);
2994            qb:=TriangulizedMat(qb);
2995            Add(sub,qb);
2996            s:=SMTX.InducedAction(q,s,3);
2997            b:=Concatenation(b{[1..den]},s[3]*b{[den+1..Length(b)]});
2998            den:=den+Length(qb);
2999            q:=s[2];
3000            den:=Length(qb);
3001          fi;
3002        fi;
3003      fi;
3004      i:=i+1;
3005    od;
3006  od;
3007  return sub;
3008end;
3009SMTX.BasesCSSmallDimUp:=SMTX_BasesCSSmallDimUp;
3010
3011SMTX_BasesCSSmallDimDown:=function(m)
3012local d,sub;
3013  d:=MTX.DualModule(m);
3014  sub:=SMTX_BasesCSSmallDimUp(d);
3015  return Concatenation([[]],
3016    Reversed(List(sub{[2..Length(sub)-1]},x->SMTX.DualizedBasis(m,x))),
3017    [IdentityMat(m.dimension,m.field)]);
3018
3019end;
3020SMTX.BasesCSSmallDimDown:=SMTX_BasesCSSmallDimDown;
3021
3022SMTX_BasesSubmodules:=function(m)
3023local cf,u,i,j,f,cl,min,neu,sq,sb,fb,k,nmin,F;
3024  F:=SMTX.Field(m);
3025  cf:=SMTX.CollectedFactors(m);
3026  cl:=Sum(cf,i->i[2]); # composition length
3027  cf:=List(cf,i->i[1]);
3028  u:=[[]];
3029  if cl>1 then
3030    min:=Concatenation(List(cf,i->SMTX.MinimalSubGModules(i,m)));
3031    u:=Concatenation(u,min);
3032  fi;
3033  for i in [2..cl-1] do
3034    neu:=[];
3035    for j in min do
3036      f:=List(j,i->List(i,i->i));
3037      sq:=SMTX.InducedAction(m,j,2);
3038      Assert(2,j=f);
3039      f:=sq[1];
3040      sb:=j;
3041      fb:=sq[2]{[Length(j)+1..Length(sq[2])]};
3042      # actually we might want to count frequencies to speed up the process,
3043      # so far I'm lazy
3044      nmin:=Concatenation(List(cf,i->SMTX.MinimalSubGModules(i,f)));
3045      Info(InfoMeatAxe,3,Length(nmin),"minimal submodules");
3046      for k in nmin do
3047        sq:=Concatenation(List(sb,ShallowCopy), # don't destroy old basis
3048                          k*fb);
3049        TriangulizeMat(sq);
3050        sq:=ImmutableMatrix(F,sq);
3051        Assert(2,SMTX.InducedAction(m,sq)<>fail);
3052        if not sq in neu then
3053          Info(InfoMeatAxe,2,"submodule dimension ",Length(sq));
3054          Add(neu,sq);
3055        fi;
3056      od;
3057    od;
3058    u:=Concatenation(u,neu);
3059    min:=neu;
3060  od;
3061  Add(u,ImmutableMatrix(SMTX.Field(m),
3062                        IdentityMat(SMTX.Dimension(m),SMTX.Field(m))));
3063  return u;
3064end;
3065SMTX.BasesSubmodules:=SMTX_BasesSubmodules;
3066
3067
3068SMTX_BasesMinimalSubmodules:=function(m)
3069local cf;
3070  cf:=SMTX.CollectedFactors(m);
3071  cf:=List(cf,i->i[1]);
3072  return Concatenation(List(cf,i->SMTX.MinimalSubGModules(i,m)));
3073end;
3074SMTX.BasesMinimalSubmodules:=SMTX_BasesMinimalSubmodules;
3075
3076SMTX.DualModule:=function(module)
3077  if SMTX.IsZeroGens(module) then
3078    return GModuleByMats([],module.dimension,SMTX.Field(module));
3079  else
3080    return GModuleByMats(List(SMTX.Generators(module),i->TransposedMat(i)^-1),
3081                        module.dimension,
3082                        SMTX.Field(module));
3083  fi;
3084end;
3085
3086###############################################################################
3087##
3088#F  DualGModule( module ) . . . . . dual of a G-module
3089##
3090## DualGModule calculates the dual of a G-module.
3091## The matrices of the module are inverted and transposed.
3092##
3093InstallGlobalFunction(DualGModule,function( module)
3094   return SMTX.DualModule(module);
3095end);
3096
3097SMTX.DualizedBasis:=function(module,sub)
3098local F,M;
3099  F:=DefaultFieldOfMatrix(sub);
3100  M:=TransposedMatMutable(sub);
3101  M:=TriangulizedNullspaceMatDestructive(M);
3102  M:=ImmutableMatrix(F,M);
3103  return M;
3104end;
3105
3106SMTX_BasesMaximalSubmodules:=function(m)
3107local d,u;
3108  d:=SMTX.DualModule(m);
3109  u:=SMTX.BasesMinimalSubmodules(d);
3110  return List(u,i->SMTX.DualizedBasis(d,i));
3111end;
3112SMTX.BasesMaximalSubmodules:=SMTX_BasesMaximalSubmodules;
3113
3114SMTX_BasesMinimalSupermodules:=function(m,sub)
3115local a,u,i,nb;
3116  a:=SMTX.InducedAction(m,sub,2);
3117  u:=SMTX.BasesMinimalSubmodules(a[1]);
3118  nb:=a[2];
3119  nb:=nb{[Length(sub)+1..Length(nb)]}; # the new basis part
3120  nb:=List(u,i->Concatenation( List( sub, ShallowCopy ),
3121                               i*nb));
3122  u:=[];
3123  for i in nb do
3124    TriangulizeMat(i);
3125    Add(u,Filtered(i,j->j<>Zero(j)));
3126  od;
3127  return u;
3128end;
3129SMTX.BasesMinimalSupermodules:=SMTX_BasesMinimalSupermodules;
3130
3131#############################################################################
3132##
3133#F SMTX.SpanOfMinimalSubGModules(m1, m2) . .
3134## span of the minimal submodules of m2 isomorphic to m1
3135##
3136## It is assumed that m1 is a module that has been proved irreducible
3137##  (using IsIrreducible), and m2 is an arbitrary module for the same group.
3138## SpanOfMinimalSubGModules computes a normed bases for the span of
3139## the minimal submodules of m2 that are isomorphic to m1,
3140## First HomGModule is called.
3141##
3142SMTX_SpanOfMinimalSubGModules:=function(m1, m2)
3143   local  homs, e, mat, i;
3144   Info(InfoMeatAxe,2,"Calculating homomorphisms from m1 to m2.");
3145   homs:=SMTX.Homomorphisms(m1, m2);
3146   if homs=[] then
3147     return [];
3148   fi;
3149   Info(InfoMeatAxe,2,"Sorting them.");
3150   SMTX.SortHomGModule(m1, m2, homs);
3151
3152   e:=SMTX.DegreeFieldExt(m1);
3153   # homs are now grouped so that each block of e have the same image.
3154   # We only want one from each block.
3155   if e > 1 then
3156     homs:=homs{Filtered([1..Length(homs)],i->(i mod e) = 1)};
3157   fi;
3158   if Length(homs) = 1 then
3159     return homs[1];
3160   fi;
3161   # The span of the images of homs is what we want!
3162   mat:=Concatenation(homs);
3163   TriangulizeMat(mat);
3164   mat:=ImmutableMatrix(m1.field,mat);
3165   return mat;
3166end;
3167SMTX.SpanOfMinimalSubGModules:=SMTX_SpanOfMinimalSubGModules;
3168
3169SMTX_BasisSocle:=function(module)
3170local cf, mat, i;
3171   cf:=SMTX.CollectedFactors(module);
3172   mat:=Concatenation(List(cf,i->SMTX_SpanOfMinimalSubGModules(i[1],module)));
3173   if Length(cf) = 1 then
3174     return ImmutableMatrix(module.field,mat);
3175   fi;
3176   TriangulizeMat(mat);
3177   mat:=ImmutableMatrix(module.field,mat);
3178   return mat;
3179end;
3180SMTX.BasisSocle:=SMTX_BasisSocle;
3181
3182SMTX_BasisRadical:=function(module)
3183local d, bs;
3184   d:=SMTX.DualModule(module);
3185   bs:=SMTX.BasisSocle(d);
3186   return SMTX.DualizedBasis(d,bs);
3187end;
3188SMTX.BasisRadical:=SMTX_BasisRadical;
3189
3190# the following assignement is for profiling
3191SMTX.funcs:=[SMTX_OrthogonalVector,SMTX_SpinnedBasis,SMTX_SubQuotActions,
3192  SMTX_SMCoRaEl,SMTX_IrreducibilityTest,SMTX_RandomIrreducibleSubGModule,
3193  SMTX_GoodElementGModule,SMTX_FrobeniusAction,SMTX_CompleteBasis,
3194  SMTX_AbsoluteIrreducibilityTest,SMTX_CollectedFactors,SMTX_Distinguish,
3195  SMTX_MinimalSubGModule,SMTX_IsomorphismComp,SMTX_MatrixSum,
3196  SMTX_Homomorphisms,SMTX_SortHomGModule,SMTX_MinimalSubGModules,
3197  SMTX_BasesCompositionSeries,SMTX_BasesSubmodules,SMTX_BasesMinimalSubmodules,
3198  SMTX_BasesMaximalSubmodules,SMTX_BasesMinimalSupermodules,SMTX_BasisSocle,
3199  SMTX_BasisRadical];
3200
3201
3202# The following functions are for finding a basis of an irreducible module
3203# that is contained in an orbit of the G-action on vectors, and for
3204# looking for G-invariant bilinear and quadratic forms of the module.
3205# The special basis is used for finding invariant quadratic forms when
3206# the characteristic of the field is 2.
3207
3208SMTX.SetBasisInOrbit:=function(module,b)
3209  module.BasisInOrbit:=b;
3210end;
3211
3212#############################################################################
3213##
3214#F  BasisInOrbit( module ) . . . .
3215##
3216## Find a basis of the irrecucible GModule module that is contained in
3217## an orbit of the action of G.
3218## The code is similar to that of SpinnedBasis.
3219SMTX_BasisInOrbit:=function( module  )
3220   local   v, matrices, ngens, zero,  ans, normedans,
3221           dim, subdim, leadpos, w, normedw, i, j, k, l, m, F;
3222
3223   if not SMTX.IsMTXModule(module) or not SMTX.IsIrreducible(module) then
3224      Error("Argument of BasisInOrbit is not an irreducible module");
3225   fi;
3226   if IsBound(module.BasisInOrbit) then return module.BasisInOrbit; fi;
3227
3228   dim:=SMTX.Dimension(module);
3229   F:=SMTX.Field(module);
3230   matrices:=module.generators;
3231   ngens:=Length(matrices);
3232
3233   zero:=Zero(F);
3234   v:=IdentityMat(dim,F)[1];
3235   ans:=[v];
3236   normedans:=[v];
3237   subdim:=1;
3238   leadpos:=SubGModLeadPos(ans,dim,subdim,zero);
3239
3240   i:=1;
3241   while i <= subdim and subdim < dim do
3242      for l in [1..ngens] do
3243         m:=matrices[l];
3244         # apply generator m to submodule generator i
3245         w:=ans[i] * m;
3246         normedw:=w;
3247         # try to express w in terms of existing submodule generators
3248         j:=1;
3249         for  j in [1..subdim] do
3250            k:=normedw[leadpos[j]];
3251            if k <> zero then
3252               normedw:=normedw - k * normedans[j];
3253            fi;
3254         od;
3255
3256         j:=1;
3257         while j <= dim and normedw[j] = zero do
3258            j:=j + 1;
3259         od;
3260         if j <= dim then
3261            # we have found a new generator of the submodule
3262            subdim:=subdim + 1;
3263            leadpos[subdim]:=j;
3264            normedw:=(normedw[j]^-1) * normedw;
3265            Add( ans, w );
3266            Add( normedans, normedw );
3267            if subdim = dim then
3268               break;
3269            fi;
3270         fi;
3271      od;
3272      i:=i + 1;
3273   od;
3274
3275   Assert(0, subdim = dim);
3276   ans:=ImmutableMatrix(F,ans);
3277   SMTX.SetBasisInOrbit(module,ans);
3278   return ans;
3279end;
3280SMTX.BasisInOrbit:=SMTX_BasisInOrbit;
3281
3282SMTX.SetInvariantBilinearForm:=function(module,b)
3283  module.InvariantBilinearForm:=b;
3284end;
3285
3286#############################################################################
3287##
3288#F  InvariantBilinearForm( module ) . . . .
3289##
3290## Look for an invariant bilinear form of the absolutely irreducible
3291## GModule module. Return fail, or the matrix of the form.
3292SMTX_InvariantBilinearForm:=function( module  )
3293   local DM, iso;
3294
3295   if not SMTX.IsMTXModule(module) or
3296                            not SMTX.IsAbsolutelyIrreducible(module) then
3297      Error(
3298 "Argument of InvariantBilinearForm is not an absolutely irreducible module");
3299   fi;
3300   if IsBound(module.InvariantBilinearForm) then
3301     return module.InvariantBilinearForm;
3302   fi;
3303   DM:=SMTX.DualModule(module);
3304   iso:=MTX.IsomorphismIrred(module,DM);
3305   if iso = fail then
3306       SMTX.SetInvariantBilinearForm(module, fail);
3307       return fail;
3308   fi;
3309   iso:=ImmutableMatrix(module.field, iso);
3310   SMTX.SetInvariantBilinearForm(module, iso);
3311   return iso;
3312end;
3313
3314SMTX.InvariantBilinearForm:=SMTX_InvariantBilinearForm;
3315
3316SMTX.MatrixUnderFieldAuto:=function(matrix, r)
3317# raise every component of matrix to r-th power
3318  local mat;
3319  mat:=List( matrix, x -> List(x, y->y^r) );
3320  mat:=ImmutableMatrix(GF(r^2), mat);
3321  return mat;
3322end;
3323
3324SMTX.TwistedDualModule:=function(module)
3325  local q, r, mats;
3326  q:=Size(module.field);
3327  r:=RootInt(q,2);
3328  if r^2 <> q then
3329    Error("Size of field of module is not a square");
3330  fi;
3331  if SMTX.IsZeroGens(module) then
3332    return GModuleByMats([],module.dimension,SMTX.Field(module));
3333  else
3334    mats:=List( SMTX.Generators(module),
3335          i->SMTX.MatrixUnderFieldAuto(TransposedMat(i)^-1,r) );
3336    return GModuleByMats( mats, module.dimension, SMTX.Field(module) );
3337  fi;
3338end;
3339
3340SMTX.SetInvariantSesquilinearForm:=function(module,b)
3341  module.InvariantSesquilinearForm:=b;
3342end;
3343
3344#############################################################################
3345##
3346#F  InvariantSesquilinearForm( module ) . . . .
3347##
3348## Look for an invariant sesquililinear form of the absolutely irreducible
3349## GModule module. Return fail, or the matrix of the form.
3350SMTX_InvariantSesquilinearForm:=function( module  )
3351   local DM, q, r, iso, isot, l;
3352
3353   if not SMTX.IsMTXModule(module) or
3354                            not SMTX.IsAbsolutelyIrreducible(module) then
3355      Error(
3356 "Argument of InvariantSesquilinearForm is not an absolutely irreducible module"
3357   );
3358   fi;
3359
3360   if IsBound(module.InvariantSesquilinearForm) then
3361     return module.InvariantSesquilinearForm;
3362   fi;
3363   DM:=SMTX.TwistedDualModule(module);
3364   iso:=MTX.IsomorphismIrred(module,DM);
3365   if iso = fail then
3366       SMTX.SetInvariantSesquilinearForm(module, fail);
3367       return fail;
3368   fi;
3369   # Replace iso by a scalar multiple to get iso twisted symmetric
3370   q:=Size(module.field);
3371   r:=RootInt(q,2);
3372   isot:=List( TransposedMat(iso), x -> List(x, y->y^r) );
3373   isot:=iso * isot^-1;
3374   if not IsDiagonalMat(isot) then
3375     Error("Form does not seem to be of the right kind (non-diagonal)!");
3376   fi;
3377   l:=LogFFE(isot[1][1],Z(q));
3378   if l mod (r-1) <> 0 then
3379     Error("Form does not seem to be of the right kind (not (q-1)st root)!");
3380   fi;
3381   iso:=Z(q)^(l/(r-1)) * iso;
3382   iso:=ImmutableMatrix(GF(q), iso);
3383   SMTX.SetInvariantSesquilinearForm(module, iso);
3384   return iso;
3385end;
3386
3387SMTX.InvariantSesquilinearForm:=SMTX_InvariantSesquilinearForm;
3388
3389SMTX.SetInvariantQuadraticForm:=function(module,b)
3390  module.InvariantQuadraticForm:=b;
3391end;
3392
3393#############################################################################
3394##
3395#F  InvariantQuadraticForm( module ) . . . .
3396##
3397## Look for an invariant quadratic form of the absolutely irreducible
3398## GModule module. Return fail, or the matrix of the form.
3399SMTX_InvariantQuadraticForm:=function( module  )
3400   local iso, bas, cgens, ciso, dim, f, z, x, i, j, qf, g, id, cqf, fix;
3401
3402   if not SMTX.IsMTXModule(module) or
3403                            not SMTX.IsAbsolutelyIrreducible(module) then
3404      Error(
3405 "Argument of InvariantQuadraticForm is not an absolutely irreducible module");
3406   fi;
3407   if IsBound(module.InvariantQuadraticForm) then
3408     return module.InvariantQuadraticForm;
3409   fi;
3410   iso:=SMTX.InvariantBilinearForm(module);
3411   if iso = fail then return fail; fi;
3412   if Characteristic(module.field) <> 2 then return iso/2; fi;
3413
3414   # In characteristic two, we change to a basis in orbit.
3415   # This makes the search for an invariant quadratic form quicker.
3416   bas:=SMTX.BasisInOrbit(module);
3417   cgens:=List(module.generators, x->bas*x*bas^-1 );
3418   ciso:=List(bas * iso * TransposedMat(bas),ShallowCopy);
3419   dim:=module.dimension;
3420   f:=module.field;
3421   z:=Zero(f);
3422
3423   # Matrix must be symplectic - perhaps it must be?
3424   for i in [1..dim] do if ciso[i][i] <> z then
3425     Print("Non-symplectic failure!\n");
3426     return fail;
3427   fi; od;
3428
3429   # If there is an invariant quadratic form, then it will be the lower
3430   # left hand part of ciso plus a scalar.
3431   for i in [1..dim-1] do for j in [i+1..dim] do ciso[i][j]:=z; od; od;
3432   id:=IdentityMat(dim, f);
3433   for x in f do
3434     qf:=ciso + x*id;
3435     fix:=true;
3436     # Form is preserved if and only if diagonal is.
3437     for g in cgens do
3438       cqf:=g * qf * TransposedMat(g);
3439       for j in [1..dim] do if cqf[j][j] <> x then
3440         fix:=false;
3441         break;
3442       fi; od;
3443       if not fix then break; fi;
3444     od;
3445     if fix then
3446       qf:=bas^-1 * qf * TransposedMat(bas^-1);
3447       # switch to lower triangular equivalent
3448       for i in [1..dim-1] do for j in [i+1..dim] do
3449         qf[j][i]:=qf[i][j] + qf[j][i];
3450         qf[i][j]:=z;
3451       od; od;
3452       qf:=ImmutableMatrix(f,qf);
3453       SMTX.SetInvariantQuadraticForm(module, qf);
3454       return qf;
3455     fi;
3456   od;
3457   SMTX.SetInvariantQuadraticForm(module, fail);
3458   return fail;
3459end;
3460
3461SMTX.InvariantQuadraticForm:=SMTX_InvariantQuadraticForm;
3462
3463#############################################################################
3464##
3465#F  OrthogonalSign( module ) . . . .
3466##
3467## When an absolutely irreducible G-module has an invariant quadratic
3468## form, this implies that it embeds in a General Orthogonal group. In
3469## even dimension there are two non-isomorphic General Orthogonal groups
3470## "plus" and "minus" type `GeneralOrthogonalGroup(+1,<n>,<q>)' and
3471## `GeneralOrthogobalGroup(-1,<n>,<q>)' in GAP terms. This function
3472## decides which one the module embeds into.
3473##
3474## It returns:
3475##  fail if the module is not absolutely irreducible, or
3476##       does not stabilize a quadratic form.
3477##  0    otherwise, if the dimension of the module is odd
3478##  +1 or -1 otherwise, according to which GO the module embeds in
3479##
3480## This is an implementation of an algorithm by Jon Thackray
3481
3482SMTX.SetOrthogonalSign:=function(module,s)
3483  module.OrthogonalSign:=s;
3484end;
3485
3486SMTX_OrthogonalSign:=function(gm)
3487    local   b,  q,  k,  n,  W,  o,  z,  lo,  lzo,  lines,  l,  w,  p,
3488            x,  y,  r,  i;
3489    if IsBound(gm.OrthogonalSign) then
3490        return gm.OrthogonalSign;
3491    fi;
3492    b:=MTX.InvariantBilinearForm(gm);
3493    q:=MTX.InvariantQuadraticForm(gm);
3494    if q = fail then
3495        return fail;
3496    fi;
3497    n:=Length(b);
3498    if n mod 2 = 1 then
3499        return 0;
3500    fi;
3501    k:=MTX.Field(gm);
3502    W:=IdentityMat(n,k);
3503
3504    #
3505    # Assemble the points of projective 3-space
3506    #
3507    o:=One(k);
3508    z:=Zero(k);
3509    lo:=[o];
3510    lzo:=[z,o];
3511    lines:=List(AsSSortedList(FullRowSpace(k,2)),x -> Concatenation(lo,x));
3512    Append(lines,List(AsSSortedList(k), x-> Concatenation(lzo,[x])));
3513    Add(lines,[z,z,o]);
3514
3515    #
3516    # Main loop of Thackray's algorithm, build up a totally isotropic
3517    # subspace and restrict it's perp until the gap between is just 2 dimensional
3518    #
3519
3520    while n > 2 do
3521
3522        #
3523        # Find an isotropic vector
3524        #
3525        for l in lines do
3526            w:=l*W;
3527            if w*q*w = z then
3528                break;
3529            fi;
3530        od;
3531        Assert(1,w*b*w = z);
3532        p:=PositionNonZero(l);
3533        #
3534        # delete it from W (add it to the subspace)
3535        #
3536        W{[p..n-1]}:=W{[p+1..n]};
3537        Unbind(W[n]);
3538        n:=n-1;
3539        #
3540        # find a vector with which it has non-zero inner product
3541        #
3542        x:=w*b;
3543        p:=PositionProperty(W, row -> x*row <> z);
3544        Assert(1, p <> fail);
3545        #
3546        # use it to find the perp of the enlarged subspace
3547        #
3548        y:=W[p];
3549        r:=x*y;
3550        for i in [p+1..n] do
3551            AddRowVector(W[i], y, - x*W[i]/r);
3552            W[i-1]:=W[i];
3553        od;
3554        Unbind(W[n]);
3555        n:=n-1;
3556        #
3557        # Now n has gone down by 2 and W is still the "gap" between the
3558        # subspace and its perp
3559        #
3560    od;
3561
3562    #
3563    # Now we need to see if the span of W contains an isotropic vector
3564    #
3565    if W[2]*q*W[2] = z then
3566        SMTX.SetOrthogonalSign(gm,1);
3567        return 1;
3568    else
3569        for x in k do
3570            w:=W[1]+x*W[2];
3571            if w*q*w = z then
3572                SMTX.SetOrthogonalSign(gm,1);
3573                return 1;
3574            fi;
3575        od;
3576        SMTX.SetOrthogonalSign(gm,-1);
3577        return -1;
3578    fi;
3579end;
3580
3581SMTX.OrthogonalSign:=SMTX_OrthogonalSign;
3582