|
#############################################################################
##
#W recograt.gi matgrp package Alexander Hulpke
##
##
#Y Copyright (C) 2014-18, Alexander Hulpke
##
## basic setup for matrix fitting free.
##
SetInfoLevel(InfoRecog,0); # recog will print status messages otherwise
# mod to recog -- rings
# newer recog has some changes to these internal APIs
if not IsBound(RIFac) and IsBound(ImageRecogNode) then
RIFac := ImageRecogNode;
RIKer := KernelRecogNode;
fi;
if IsBound(BindRecogMethod) then
BindRecogMethod(FindHomMethodsMatrix, "Nonfield",
"catch matrix groups defined over nonfield rings",
function(ri, G)
if IsBound(ri!.ring) and not IsBound(ri!.field) then
Error("hereIAm");
fi;
return false;
end);
AddMethod( FindHomDbMatrix, FindHomMethodsMatrix.Nonfield, 5100);
else
FindHomMethodsMatrix.Nonfield := function(ri, G)
if IsBound(ri!.ring) and not IsBound(ri!.field) then
Error("hereIAm");
fi;
return false;
end;
AddMethod( FindHomDbMatrix, FindHomMethodsMatrix.Nonfield,
5100, "Nonfield",
"catch matrix groups defined over nonfield rings" );
fi;
OnSubmoduleCosets:=function(cset,g)
local v;
return [SiftedVector(cset[2],cset[1]*g),cset[2]];
end;
MakeSubmoduleCosetAction:=function(basis)
return function(vec,g)
return SiftedVector(basis,vec*g);
end;
end;
MakeSubmoduleColineAction:=function(basis)
return function(vec,g)
local c;
vec:=SiftedVector(basis,vec*g);
c:=PositionNonZero(vec);
if c<=Length(vec) then
vec := Inverse( vec[c] ) * vec;
fi;
return vec;
end;
end;
FUNCSPACEHASH:=[];
# mod to genss -- use submodules and quotients for base points
MSSFBPC:=function( grp, opt, ii, parentS ) # owf fct to call easily
local F,dim,orb,orbs,i,fct,
mo,cs,j,k,dims,bas,basc,basinv,nb,lastdim,cand,fcand,sel,limit,trysel,submodule;
trysel:=function(recsub,recfac)
local lgens;
# nor the trivial action
if ForAny(mo,x->Order(x{sel}{sel})>1) then
Info(InfoFFMat,2,"range ",sel," have ",Length(cand.points));
lgens:=List(mo,x->x{sel}{sel});
fcand:=FindBasePointCandidates(Group(lgens),opt,ii,
false:Subrecurse:=recsub,Facrecurse:=recfac);
Info( InfoGenSS, 3, "Subfactor module of range ",sel,", ",Length(fcand.ops),
" candidates");
for k in [1..Length(fcand.ops)] do
if ForAll(lgens,x->fcand.ops[k](fcand.points[k],x)=fcand.points[k]) then
Info(InfoFFMat,2,"Ignoring fixed element for base");
elif fcand.ops[k]=OnRight or fcand.ops[k]=OnPoints or fcand.ops[k]=OnLines then
nb:=fcand.points[k]*bas{sel};
if lastdim=0 then
# proper subspace -- just vectors
Add(cand.points,nb);
Add(cand.ops,fcand.ops[k]);
elif true then
# # action on cosets
submodule:=SemiEchelonBasis(VectorSpace(F,bas{[1..lastdim]},Zero(bas[1])));
nb:=SiftedVector(submodule,nb);
if fcand.ops[k]=OnLines then
fct:=MakeSubmoduleColineAction(submodule);
Add(FUNCSPACEHASH,[fct,submodule]);
else
fct:=MakeSubmoduleCosetAction(submodule);
Add(FUNCSPACEHASH,[fct,submodule]);
fi;
Add(cand.points,nb);
Add(cand.ops,fct);
elif Length(sel)=1 then
# TODO: 1-dim factor -- need to do cosets
Info(InfoWarning,1,"Case not yet implemented");
else
# subfactor -- take subspace preimage
nb:=OnSubspacesByCanonicalBasis(Concatenation(bas{[1..lastdim]},[nb]),
One(grp));
Add(cand.points,nb);
Add(cand.ops,OnSubspacesByCanonicalBasis);
fi;
elif ForAny(FUNCSPACEHASH,x->x[1]=fcand.ops[k]) then
fct:=First(FUNCSPACEHASH,x->x[1]=fcand.ops[k]);
submodule:=SemiEchelonBasis(VectorSpace(F,
Concatenation(bas{[1..lastdim]},
BasisVectors(fct[2])*bas{sel})));
Add(cand.points,fcand.points[k]*bas{sel});
fct:=MakeSubmoduleCosetAction(submodule);
Add(FUNCSPACEHASH,[fct,submodule]);
Add(cand.ops,fct);
Info(InfoFFMat,2,"ACTPOP");
elif fcand.ops[k]=OnSubspacesByCanonicalBasis then
nb:=fcand.points[k]*bas{sel};
nb:=OnSubspacesByCanonicalBasis(Concatenation(bas{[1..lastdim]},nb),
One(grp));
Add(cand.points,nb);
Add(cand.ops,OnSubspacesByCanonicalBasis);
else
Info(InfoWarning,1,"Action not recognized");
fi;
od;
return true;
else
return false;
fi;
end;
if IsBound(opt.VeryShortOrbLimit) then
limit:=2*opt.VeryShortOrbLimit;
else
limit:=10^4;
fi;
F := DefaultFieldOfMatrixGroup(grp);
# don't bother in small cases
if Size(F)^Length(One(grp))<=opt.ShortOrbitsOrbLimit then
TryNextMethod();
fi;
cs:=GeneratorsOfGroup(grp);
if ForAny(cs,IsObjWithMemory) then
cs:=Concatenation(Filtered(cs,x->not IsObjWithMemory(x)),
List(Filtered(cs,x->IsObjWithMemory(x)),
x->x!.el));
fi;
cand:=rec(ops:=[],points:=[],used:=0);
mo:=GModuleByMats(cs,F);
dim:=mo.dimension;
if MTX.IsIrreducible(mo) or Size(F)^dim<=limit then
TryNextMethod();
fi;
# build new basis corresponding to comp ser.
cs:=MTX.BasesCompositionSeries(mo);
Info(InfoFFMat,2,"dims=",List(cs,Length));
dims:=List(cs,Length);
bas:=[];
basc:=[];
for j in [2..Length(cs)] do
nb:=BaseSteinitzVectors(cs[j],basc);
Append(bas,nb.factorspace);
basc:=Concatenation(nb.subspace,nb.factorspace);
Sort(basc,function(a,b) return PositionNonZero(a)<PositionNonZero(b);end);
od;
basinv:=bas^-1;
mo:=List(mo.generators,x->bas*x*basinv);
# now step up in sizes indicated by short orbit lengths
lastdim:=0;
j:=2;
if ValueOption("Subrecurse")<>false then
while j<=Length(dims) and Length(cand.points)<=5 do
# don't bother is the space is too small
if (j=Length(dims) and lastdim>0) or
(j<Length(dims) and Size(F)^(dims[j+1]-lastdim)>limit) then
sel:=[lastdim+1..dims[j]];
if trysel(false,true) then
# we tried a space
lastdim:=dims[j];
fi;
fi;
j:=j+1;
od;
fi;
if lastdim=0 and ValueOption("Facrecurse")<>false then
# all but last submodule was trivial -- step down on factors
j:=Length(dims)-1;
while j>0 and Length(cand.points)<=5 do
lastdim:=dims[j]-1;
if (j>1 and Size(F)^(dim-dims[j-1])>limit) then
sel:=[lastdim+1..dim];
if trysel(false,false) then
j:=0;
fi;
fi;
j:=j-1;
od;
fi;
if Length(cand.points)>0 then return cand; fi;
# refer to parent
if IsBound(opt.PCand) then
# try which ones are short
sel:=[];
for i in [1..Length(opt.PCand.points)] do
orb:=[opt.PCand.points[i]];
mo:=opt.PCand.ops[i];
orbs:=Set(orb); # as short, set is fine.
j:=1;
while j<=Length(orb) and Length(orb)<=limit do
for k in GeneratorsOfGroup(grp) do
cs:=mo(orb[j],k);
if not cs in orbs then
Add(orb,cs);
AddSet(orbs,cs);
fi;
od;
j:=j+1;
od;
if Length(orb)<=limit then
Add(sel,i);
Add(cand.points,orb[1]);
Add(cand.ops,mo);
fi;
od;
Info(InfoFFMat,2,"Selected ",Length(sel)," of ",Length(opt.PCand.points)," group basis vectors");
if Length(cand.points)>0 then return cand; fi;
fi;
TryNextMethod();
end;
InstallMethod(FindBasePointCandidates,
"for reducible matrix group over a FF, use submodules and quotients",
[ IsGroup and IsMatrixGroup and IsFinite, IsRecord, IsInt, IsObject ], 21,
MSSFBPC);
GetInformationFromRecog:=function(recog)
local treerecurse,n,factors,homs,leafgens,niceranges,genum,sz,leafs,g,
minranges,mingens,permap,furtherhom,fn,extragens,extranum,extra,allgens,
nicegp,map,stc,orbit,gd,cnt;
# the main worker function -- traverse the tree.
treerecurse:=function(r,h)
local
hom,f,k,sel,u,i,j,x,cs,nsol,xf,bigger,nicehom,ngi,stbc,opt,act,ng,dom,v;
if IsLeaf(r) then
f:=Length(NiceGens(r));
Info(InfoFFMat,2,"Leaf",Size(r)," ",genum," ",f);
k:=[genum+1..genum+f];
if Size(r)=1 then
Info(InfoFFMat,2,"ignoring trivial factor");
elif Length(Set(Factors(Size(r))))<3 then
Info(InfoFFMat,2,"ignoring solvable factor of order ",Size(r));
# store info
n:=n+Length(Factors(Size(r)));
SetIsSolvableGroup(Grp(r),true);
permap[n]:=IdentityMapping(Grp(r));
leafs[n]:=r;
homs[n]:=h;
factors[n]:=false;
sz[n]:=Size(r);
leafgens[n]:=f;
sel:=[1..Length(NiceGens(r))];
mingens[n]:=sel;
niceranges[n]:=k;
minranges[n]:=k{sel};
else
nicehom:=false;
if IsPermGroup(Grp(r)) then
nicehom:=IdentityMapping(Grp(r));
elif IsBound(r!.stabilizerchain) then
#use ActionOnOrbit?
stbc:=BaseStabilizerChain(r!.stabilizerchain);
act:=stbc.ops[1];
if ForAll(stbc.ops,x->x=act) then
# all the same action -- union
orbit:=[];
for j in [1..Length(stbc.points)] do
if not stbc.points[j] in orbit then
orbit:=Union(orbit,Orbit(Grp(r),stbc.points[j],stbc.ops[j]));
fi;
od;
if Length(orbit)>1 and Length(orbit)<50000 then
nicehom:=ActionHomomorphism(Grp(r),orbit,stbc.ops[1],"surjective");
Info(InfoFFMat,2,"Got degree ",Length(orbit));
fi;
fi;
fi;
if nicehom=false then
# Hasfhmethsel, success method: StabilizerChain
if IsBound(r!.projective) and r!.projective then
g:=Grp(r);
v:=OnLines(g.1[1],One(g));
dom:=ShallowCopy(Orbit(g,v,OnLines));
nicehom:=ActionHomomorphism(g,dom,OnLines,"surjective");
while Size(Image(nicehom))<Size(r) do
cnt:=0;
repeat
v:=OnLines(Random(DefaultFieldOfMatrixGroup(g)^Length(One(g))),
One(g));
cnt:=cnt+1;
if cnt>1000 then Error("no vector found");fi;
until not v in dom;
Append(dom,Orbit(g,v,OnLines));
nicehom:=ActionHomomorphism(g,dom,OnLines,"surjective");
od;
else
nicehom:=IsomorphismPermGroup(Grp(r));
fi;
fi;
g:=Image(nicehom);
if Size(g)<>Size(r) then
Error("some discrepancy happened");
fi;
#test!!
# .isknownsimple, .isknownalmostsimple
if IsBound(r!.comment) and r!.comment[1]<>'_' then
SetIsSimpleGroup(g,true);
gd:=g;
elif not IsSolvableGroup(g) then
gd:=PerfectResiduum(g);
else
gd:=fail;
fi;
if gd<>fail then
if NrMovedPoints(g)> SufficientlySmallDegreeSimpleGroupOrder(Size(gd))
then
nicehom:=nicehom*SmallerDegreePermutationRepresentation(g);
gd:=Image(nicehom);
Info(InfoFFMat,2,"Improved degree ",NrMovedPoints(g),"->",
NrMovedPoints(gd));
if HasIsSimpleGroup(g) and IsSimpleGroup(g) then
SetIsSimpleGroup(gd,true);
SetIsSolvableGroup(gd,false);
fi;
g:=gd;
fi;
fi;
## indicate unsolvable
#nsol:=IsSimpleGroup(g) or
# (Length(Set(Factors(Size(r))))>2 and not IsSolvableGroup(g));
if IsSolvableGroup(g) then
Info(InfoFFMat,2,"Ignoring solvable factor of order ",Size(g));
n:=n+Length(Factors(Size(g)));
elif not IsSimpleGroup(g) then
Info(InfoFFMat,2,"doing size",Size(g)," ",IsSimpleGroup(g));
# not simple -- split
cs:=CompositionSeries(g);
ng:=List(NiceGens(r),x->ImagesRepresentative(nicehom,x));
for i in [2..Length(cs)] do
extra:=[];
n:=n+1;
# test generators
u:=cs[i];
sel:=[];
for j in [1..f] do
x:=ng[j];
if x in cs[i-1] and not x in u then
Add(sel,j);
u:=ClosureSubgroup(u,x);
fi;
od;
nicegp:=Group(ng);
map:=EpimorphismFromFreeGroup(nicegp);
if Size(u)<Size(cs[i-1]) then
Info(InfoFFMat,2,"cannot compatibilize generators -- add extras");
while Size(u)<Size(cs[i-1]) do
x:=First(GeneratorsOfGroup(cs[i-1]),x->not x in u);
u:=ClosureSubgroup(u,x);
# decompose into word
xf:=Factorization(nicegp,x);
if xf=fail then Error("factorization error");fi;
x:=MappedWord(xf,MappingGeneratorsImages(map)[1],allgens{k});
extragens[extranum]:=x;
Add(extra,extranum);
extranum:=extranum+1;
od;
fi;
hom:=NaturalHomomorphismByNormalSubgroup(cs[i-1],cs[i]);
if not IsAbelian(Image(hom)) then
permap[n]:=IsomorphismPermGroup(Image(hom));
fi;
fn:=fn+1;
if Size(cs[i])=1 then
furtherhom[fn]:=nicehom;
else
furtherhom[fn]:=nicehom*hom;
fi;
leafs[n]:=false;
homs[n]:=Concatenation(h,[fn]);
factors[n]:=Image(hom);
sz[n]:=Size(Image(hom));
leafgens[n]:=false;
niceranges[n]:=false;
minranges[n]:=Concatenation(k{sel},extra);
od;
else
# found a simple case
n:=n+1;
permap[n]:=nicehom;
# get smaller generating set
sel:=[];
u:=TrivialSubgroup(Image(nicehom));
for i in [1..Length(NiceGens(r))] do
ngi:=ImagesRepresentative(nicehom,NiceGens(r)[i]);
if not ngi in u then
u:=ClosureGroup(u,ngi);
Add(sel,i);
fi;
od;
leafs[n]:=r;
homs[n]:=h;
factors[n]:=g;
sz[n]:=Size(r);
leafgens[n]:=f;
mingens[n]:=sel;
niceranges[n]:=k;
minranges[n]:=k{sel};
fi;
fi;
genum:=genum+f;
else
#hom:=Homom(r);
f:=RIFac(r);
treerecurse(f,Concatenation(h,[1]));
k:=RIKer(r);
if k<>fail then
treerecurse(k,Concatenation(h,[0]));
fi;
fi;
end;
extragens:=[];
furtherhom:=[];
fn:=2;
n:=0;
genum:=0;
factors:=[];
homs:=[];
leafgens:=[];
leafs:=[];
niceranges:=[];
minranges:=[];
mingens:=[];
sz:=[];
permap:=[];
allgens:=NiceGens(recog);
extranum:=Length(allgens)+1;
treerecurse(recog,[]);
return rec(
group:=Grp(recog),
n:=n,
genum:=Length(allgens),
recog:=recog,
leafs:=leafs,
factors:=factors,
furtherhom:=furtherhom,
sz:=sz,
homs:=homs,
leafgens:=leafgens,
minranges:=minranges,
mingens:=mingens,
extragens:=extragens,
permap:=permap,
niceranges:=niceranges);
end;
# assume that hom i can be applied
CSIImageHomNr:=function(csi,n,x)
local r,h,i;
r:=csi.recog;
h:=csi.homs[n];
for i in h do
if i=0 then
r:=RIKer(r);
elif i=1 then
x:=ImageElm(Homom(r),x);
r:=RIFac(r);
else
x:=ImageElm(csi.furtherhom[i],x);
fi;
od;
return x;
end;
# since nice ranges can contain negatives, we cannot simply sublist
CSINiceGens:=function(csi,a)
if IsList(a) then
return List(a,x->CSINiceGens(csi,x));
elif a<=csi.genum then
return NiceGens(csi.recog)[a];
else
return csi.extragens[a];
fi;
end;
CSIDepthElm:=function(csi,x)
local i;
i:=1;
while i<=csi.n and
(not IsBound(csi.leafs[i]) or # early skipped multiple solvable factor
(
(IsBool(csi.leafs[i]) or not (IsBound(csi.leafs[i]!.projective) and csi.leafs[i]!.projective)
or csi.leafs[i]!.projective=false)
and Order(CSIImageHomNr(csi,i,x))=1)
or
(not IsBool(csi.leafs[i]) and (IsBound(csi.leafs[i]!.projective) and csi.leafs[i]!.projective) and
Order(ImagesRepresentative(csi.permap[i],CSIImageHomNr(csi,i,x)))=1)) do
i:=i+1;
od;
return i;
end;
# to test whether an element act trivially projectively, it is not enough to
# test on base, but also need to have sum of base
# elm must be an element that is already image
CSIProjectiveBases:=function(csi,i,elm)
local m;
if not IsBound(csi.projbases) then
csi.projbases:=[];
fi;
if not IsBound(csi.projbases[i]) then
m:=ShallowCopy(One(elm));
Add(m,Sum(m));
csi.projbases[i]:=m;
fi;
return csi.projbases[i];
end;
CSIDepthElm:=function(csi,x)
local i,ximg;
i:=1;
while i<=csi.n do
if not IsBound(csi.leafs[i]) or
((IsBool(csi.leafs[i]) or not (IsBound(csi.leafs[i]!.projective) and csi.leafs[i]!.projective)
or csi.leafs[i]!.projective=false)
and Order(CSIImageHomNr(csi,i,x))=1) then
i:=i+1;
elif (not IsBool(csi.leafs[i]) and IsBound(csi.leafs[i]!.projective) and
csi.leafs[i]!.projective) then
ximg:=CSIImageHomNr(csi,i,x);
if ForAll(CSIProjectiveBases(csi,i,ximg),z->OnLines(z,ximg)=z) then
i:=i+1;
else
return i;
fi;
else
return i;
fi;
od;
return i;
end;
CSIAelement:=function(a,localgens,l)
local limg,rep;
limg:=List(l,x->Image(a[2],x));
rep:= RepresentativeAction(a[1],localgens,limg,OnTuples);
if rep=fail then
Error("no rep");
fi;
return rep;
end;
FindAct:=function(csi)
local csinice,dci,pool,
act,n,i,j,k,l,gens,lgens,c,d,genims,gp,hom,auts,isoms,pools,poolnum,a,x,
isom,diso,conj,dgp,imgdepth,perms,kn,process,ii,goodgens,conjgens,genimgs,
doesaut,biggens,wrimages,m,w,e,poolimggens,img,localgens,dfgens,wrs,dfimgs,b,perm,lims,map,aels,wrsoc,dfs;
# get a list of generator images and construct the appropriate element of
# a[1] inducing these generator images by conjugation
aels:=[];
# dci is decomposition information needed for restriction to proper
# subgroups
dci:=rec(csi:=csi);
# isoms[i] An isomorphism from the first group in its pool to group i
goodgens:=[];
poolimggens:=[];
genimgs:=List([1..Maximum(Union(csi.minranges))],x->[]);
n:=csi.n;
csinice:=NiceGens(csi.recog);
act:=[];
auts:=[];
pools:=[];
perms:=[];
isoms:=[];
n:=csi.n;
doesaut:=[];
process:=[n];
ii:=1;
while ii<=n do
# do we need to feed in another layer?
if Length(process)<ii then
Add(process,First([n,n-1..1],x->not x in process));
fi;
i:=process[ii];
if IsBound(csi.permap[i]) and not IsSolvableGroup(Image(csi.permap[i])) then
if not IsBound(goodgens[i]) then
# no generators set yet -- just take new ones
#lgens:=csinice{csi.minranges[i]};
lgens:=CSINiceGens(csi,csi.minranges[i]);
goodgens[i]:=lgens;
else
lgens:=goodgens[i];
fi;
gp:=Image(csi.permap[i]);
act[i]:=[];
gens:=List(lgens,x->ImagesRepresentative(csi.permap[i],CSIImageHomNr(csi,i,x)));
# for x in [1..Length(gens)] do
# j:=CSIImageHomNr(csi,i,csinice[csi.minranges[i][x]]);
# k:=Image(csi.permap[i],j);
# if k<>gens[x] then
# Error("GENS");
# fi;
# od;
# nonsolvable factor -- Is it in pools?
poolnum:=First([1..Length(pools)],x->i in pools[x]);
# pools will always be joined TO the current one, so isom will not
# have to change on the way.
if poolnum=fail then
Add(pools,[i]);
poolnum:=Length(pools);
auts[poolnum]:=[];
isoms[i]:=false; # representative
isom:=IdentityMapping(gp);
Add(poolimggens,gens);
elif i<>pools[poolnum][1] then
isom:=isoms[i]; # iso
else
isom:=IdentityMapping(gp); # first one -- iso is trivial
fi;
#find out what the previous factors do on it
for j in Filtered([1..i-1],x->IsBound(csi.minranges[x])) do
Info(InfoFFMat,2,"Try ",i," mapped by ",j);
for kn in csi.minranges[j] do
#k:=csinice[kn];
k:=CSINiceGens(csi,kn);
if not IsBound(perms[kn]) then
perms[kn]:=[];
genimgs[kn]:=[];
doesaut[kn]:=[];
fi;
conjgens:=[];
genims:=[];
imgdepth:=fail;
# calculate images
for l in lgens do
c:=l^k; # conjugate
Add(conjgens,c);
d:=CSIDepthElm(csi,c);
if imgdepth=fail then
# new image -- isomorphism
imgdepth:=d;
perms[kn][i]:=d; # record permutations
elif imgdepth<>d then
# this cannot happen
Error("incompatible depths");
fi;
Add(genims,ImagesRepresentative(csi.permap[d],CSIImageHomNr(csi,d,c)));
od;
dgp:=Image(csi.permap[d]);
if AssertionLevel()>2 then
hom:=GroupHomomorphismByImages(gp,dgp,gens,genims);
else
hom:=GroupHomomorphismByImagesNC(gp,dgp,gens,genims);
fi;
if hom=fail then Error("should not happen");fi;
if d=i then
# pull generator images back in original group
genimgs[kn][i]:=List(genims,x->PreImagesRepresentativeNC(isom,x));
#if ForAny(Flat(genimgs),x->not x in Image(csi.permap[pools[poolnum][1]])) then
# Error("imerrD");
#fi;
# component is not moved we just need to conjugate with isom
hom:=isom*hom*InverseGeneralMapping(isom);
SetIsBijective(hom,true);
if not IsInnerAutomorphism(hom) then
Add(auts[poolnum],hom);
AddSet(doesaut[kn],i);
elif not IsOne(hom) then
# not automorphism, but still acting
AddSet(doesaut[kn],i);
fi;
elif d in pools[poolnum] then
# we know already that the groups are isomorphic -- just store
# the new automorphism
# isomorphism is always to the first element in the pool
if d=pools[poolnum][1] then
# the group is the normal form -- we do not need to translate
diso:=IdentityMapping(dgp);
else
# isomorphism to canonical form
diso:=InverseGeneralMapping(isoms[d]);
fi;
genimgs[kn][i]:=List(genims,x->ImagesRepresentative(diso,x));
# if ForAny(Flat(genimgs),x->not x in Image(csi.permap[pools[poolnum][1]])) then
# Error("imerrC");
# fi;
hom:=isom*hom*diso;
SetIsBijective(hom,true);
if not IsInnerAutomorphism(hom) then
Add(auts[poolnum],hom);
AddSet(doesaut[kn],i);
fi;
else
# isomorphism wasn't known yet -- we need to join pools
# hom is the joiner, isom*hom*isoms[d]^-1 the conjugator of the
# canonical map
# the pool we add. We always join the image pool to the current
# one.
a:=First([1..Length(pools)],x->d in pools[x]);
if a=fail then
# we did not yet know layer d -- add it
Info(InfoFFMat,2,"Found Layer ",d);
Add(pools[poolnum],d);
isoms[d]:=isom*hom;
Assert(3,IsBijective(isoms[d]));
# newly included component -- the generators *are* simply the
# images
genimgs[kn][i]:=List(genims,
x->PreImagesRepresentativeNC(isoms[d],x));
#if ForAny(Flat(genimgs),x->not x in Image(csi.permap[pools[poolnum][1]])) then
#Error("imerrB");
#fi;
goodgens[d]:=conjgens;
Add(process,d);
else
Error("This cannot happen??");
# we knew the layer and join pools
Info(InfoFFMat,2,"Join ",pools[a]," to ",pools[poolnum]);
conj:=isom*hom;
if not d=pools[a][1] then
# translate to normal form
conj:=conj*InverseGeneralMapping(isoms[d]);
fi;
# join pools
for x in pools[a] do
Add(pools[poolnum],x);
isoms[x]:=conj*isoms[x]; # reroot isomorphism
od;
# translate automorphisms
for x in auts[a] do
Add(auts[poolnum],conj*x*InverseGeneralMapping(conj));
od;
# delete old
pools[a]:=[];
auts[a]:=[];
fi;
fi;
od;
od;
fi;
ii:=ii+1;
od;
wrs:=[];
wrsoc:=[];
dfgens:=[];
dfimgs:=[];
#if ForAny(Flat(genimgs),x->not x in Image(csi.permap[pools[poolnum][1]])) then
# Error("imerrA");
#fi;
if Length(pools)=0 then
# solvable
d:=Group(());
a:=GeneratorsOfGroup(csi.group);
b:=List(a,x->One(d));
RUN_IN_GGMBI:=true; # hack to skip Nice treatment
hom:=GroupHomomorphismByImagesNC(csi.group,d,a,b);
RUN_IN_GGMBI:=false;
dci:=rec(isTrivial:=true);
SetRecogDecompinfoHomomorphism(hom,dci);
SetImagesSource(hom,Group(()));
return hom;
fi;
dci.pools:=pools;
dci.wreathemb:=[];
dci.goodgens:=goodgens;
dci.isoms:=isoms;
dci.poollocalgens:=[];
# now build each group
for i in [1..Length(pools)] do
pool:=pools[i];
m:=Length(pools[i]);
a:=Image(csi.permap[pools[i][1]]);
a:=[a,Concatenation(auts[i],
List(GeneratorsOfGroup(a),x->InnerAutomorphism(a,x)))];
a:=AutomorphismRepresentingGroup(a[1],a[2]);
aels[i]:=a;
localgens:=List(poolimggens[i],x->Image(a[2],x));
dci.poollocalgens[i]:=localgens;
b:=SymmetricGroup(m);
w:=WreathProduct(a[1],b);
e:=List([1..m+1],x->Embedding(w,x));
dci.wreathemb[i]:=e;
biggens:=[];
wrimages:=[];
for j in Filtered([1..n],x->IsBound(csi.minranges[x])) do
c:=csi.minranges[j];
Info(InfoFFMat,2,"Images for ",j," ",c);
# does any of the elements actually do something -- if so, what?
# if true or ForAny(c,x->IsBound(doesaut[x])
# and ForAny(pools[i],y->y in doesaut[x])) then
Info(InfoFFMat,2,c," acts ",pools[i],j);
for kn in c do
#Add(biggens,csinice[kn]);
Add(biggens,CSINiceGens(csi,kn));
img:=One(w);
for l in [1..Length(pools[i])] do
if IsBound(genimgs[kn][pools[i][l]]) then
d:=Image(e[l],CSIAelement(a,localgens,genimgs[kn][pools[i][l]]));
else
if pools[i][l]=j then
d:=List(goodgens[pools[i][l]],
x->ImagesRepresentative(csi.permap[pools[i][l]],
CSIImageHomNr(csi,pools[i][l],x^CSINiceGens(csi,kn))));
if isoms[j]<>false then
d:=List(d,x->PreImagesRepresentativeNC(isoms[j],x));
fi;
d:=Image(e[l],CSIAelement(a,localgens,d));
else
# component is not acted on as it lies higher
d:=One(a[1]);
fi;
fi;
img:=img*d;
od;
# is it a nontrivial permutation
if IsBound(perms[kn]) and ForAny([1..Length(perms[kn])],
x->IsBound(perms[kn][x]) and perms[kn][x]<>x) then
# fill up fixed points
for d in pool do
if not IsBound(perms[kn][d]) then perms[kn][d]:=d;fi;
od;
# permuting part
d:=PermList(List(perms[kn]{pool},x->Position(pool,x)));
img:=img*Image(e[Length(pool)+1],d);
fi;
Add(wrimages,img);
od;
od;
Add(wrs,w);
a:=List([1..m],x->PerfectResiduum(Image(Embedding(w,x))));
b:=TrivialSubgroup(w);
for l in a do
b:=ClosureGroup(b,l);
od;
SetDirectFactorsFittingFreeSocle(w,a);
Add(wrsoc,b); # socle
Add(dfgens,biggens);
Add(dfimgs,wrimages);
od;
d:=DirectProduct(wrs);
dci.dirprod:=d;
dci.embeddings:=List([1..Length(pools)],x->Embedding(d,x));
dci.aels:=aels;
a:=List([1..Length(pools)],x->List(DirectFactorsFittingFreeSocle(wrs[x]),
y->Image(dci.embeddings[x],y)));
dfs:=Concatenation(a);
SetDirectFactorsFittingFreeSocle(d,dfs);
a:=dfgens[1];
b:=List(a,x->One(d));
w:=TrivialSubgroup(d);
for i in [1..Length(pools)] do
e:=dci.embeddings[i];
w:=ClosureGroup(w,Image(e,GeneratorsOfGroup(wrsoc[i])));
for j in [1..Length(a)] do
Assert(1,dfgens[i][j]=a[j]);
b[j]:=b[j]*Image(e,dfimgs[i][j]);
od;
od;
if IsPermGroup(csi.group) and AssertionLevel()>2 then
hom:=GroupHomomorphismByImages(csi.group,d,a,b);
else
RUN_IN_GGMBI:=true; # hack to skip Nice treatment
hom:=GroupHomomorphismByImagesNC(csi.group,d,a,b);
RUN_IN_GGMBI:=false;
fi;
if hom=fail then
Error("fail!");
fi;
b:=Subgroup(d,b);
#SetSocle(b,w);
dci.socle:=w;
dfs:=Filtered(dfs,x->IsSubset(b,x));
SetDirectFactorsFittingFreeSocle(w,dfs);
SetDirectFactorsFittingFreeSocle(b,dfs);
SetImagesSource(hom,b);
SetRecogDecompinfoHomomorphism(hom,dci);
return hom;
end;
CSIElementAct:=function(dci,elm)
local csi,pools,result,i,e,a,img,b,dp,d,kn,perm,l;
csi:=dci.csi;
pools:=dci.pools;
result:=One(dci.dirprod);
for i in [1..Length(pools)] do
e:=dci.wreathemb[i];
a:=dci.aels[i];
img:=One(Image(e[1]));
perm:=[];
for l in [1..Length(pools[i])] do
b:=List(dci.goodgens[pools[i][l]],x->x^elm);
dp:=CSIDepthElm(csi,b[1]);
Assert(2,ForAll(b,x->CSIDepthElm(csi,x)=dp));
perm[l]:=Position(pools[i],dp);
d:=List(dci.goodgens[pools[i][l]],
x->ImagesRepresentative(csi.permap[dp],
CSIImageHomNr(csi,dp,x^elm)));
if dci.isoms[dp]<>false then
d:=List(d,x->PreImagesRepresentativeNC(dci.isoms[dp],x));
fi;
d:=Image(e[l],CSIAelement(a,dci.poollocalgens[i],d));
img:=img*d;
od;
# permuting part
d:=PermList(perm);
img:=img*Image(e[Length(pools[i])+1],d);
result:=result*Image(dci.embeddings[i],img);
od;
return result;
end;
#############################################################################
##
#M ImagesRepresentative(<hom>,<x>)
##
InstallMethod(ImagesRepresentative,"for recognition mappings",
FamSourceEqFamElm,
[ IsGroupGeneralMapping and
HasRecogDecompinfoHomomorphism,IsMultiplicativeElementWithInverse], 0,
function(hom, x)
local d;
d:=RecogDecompinfoHomomorphism(hom);
if IsBound(d.isTrivial) then return ();fi;
if IsBound(d.fct) then
return d.fct(x);
else
return CSIElementAct(RecogDecompinfoHomomorphism(hom),x);
fi;
end);
#############################################################################
##
#M PreImagesSetNC(<hom>,<x>)
##
InstallMethod(PreImagesSetNC,"for recognition mappings", CollFamRangeEqFamElms,
[ IsGroupGeneralMapping and
HasRecogDecompinfoHomomorphism,IsGroup], 0,
function(hom, U)
local gens,pre;
gens:=SmallGeneratingSet(U);
pre:=List(gens,x->PreImagesRepresentativeNC(hom,x));
U:=RecogDecompinfoHomomorphism(hom).LiftSetup;
U:=SubgroupByFittingFreeData(Source(hom),pre,gens,U.pcgs);
return U;
end);
#TODO: Detect Nonsolvable permuters (via perms) and leave out of pool
BasePointsActionsOrbitLengthsStabilizerChain:=function(c)
local l,o;
l:=[];
while c<>false do
o:=c!.orb;
Add(l,[o!.orbit[1],o!.op,Length(o!.orbit)]);
c:=c!.stab;
od;
return l;
end;
FactorspaceActfun:=function(field,bas)
local heads;
bas:=TriangulizedMat(bas);
heads:=HeadsInfoOfSemiEchelonizedMat(bas,Length(bas[1]));
return function(v,g)
local c;
v:=v*g;
v:=SiftedVectorForGaussianRowSpace(field,bas, heads, v );
c:=PositionNonZero(v);
if c<=Length(v) then
v:=Inverse(v[c])*v;
fi;
MakeImmutable(v);
return v;
end;
end;
INVTRANSPCACHE:=[];
AsInverseTranspose:=function(x,g)
local a,p;
a:=fail;
p:=0;
while a=fail and p<Length(INVTRANSPCACHE) do
p:=p+1;
if IsIdenticalObj(INVTRANSPCACHE[p][1],g) then
a:=INVTRANSPCACHE[p][2];
if p>50 then
p:=Concatenation([p],[1..p-1],[p+1..Length(INVTRANSPCACHE)]);
INVTRANSPCACHE:=INVTRANSPCACHE{p};
fi;
fi;
od;
if a=fail then
a:=TransposedMat(Inverse(g));
if Length(INVTRANSPCACHE)>100 then
# clean out old
INVTRANSPCACHE:=Concatenation([[g,a]],INVTRANSPCACHE{[1..90]});
else
INVTRANSPCACHE:=Concatenation([[g,a]],INVTRANSPCACHE);
fi;
fi;
return OnRight(x,a);
end;
ModuleStructureBase:=function(mats)
local orbtranslimit,f,total,gens,mo,bas,basn,dims,a,p,vec,orb,t,dict,use,fct,
new,alt,g,pnt,limit,whole,siftchain,sub,b,trymultipleorbits;
orbtranslimit:=function(vec,fct,limit)
local dict,orb,t,pnt,g,img,p,coinc;
dict:=NewDictionary(vec,true,f^Length(vec));
orb:=[vec];
AddDictionary(dict,vec,1);
t:=[One(gens[1])];
pnt:=1;
coinc:=true;
while pnt<=Length(orb) do
for g in gens do
img:=fct(orb[pnt],g);
p:=LookupDictionary(dict,img);
if p=fail then
Add(orb,img);
if Length(orb)>limit then
return fail;
fi;
AddDictionary(dict,img,Length(orb));
Add(t,t[pnt]*g);
elif coinc then
coinc:=false;
fi;
od;
pnt:=pnt+1;
od;
return [dict,orb,t];
end;
trymultipleorbits:=function(seeds,fct,limit)
local best,bestl,bestc,i,orb;
best:=fail;
bestl:=0;
bestc:=0;
for i in seeds do
# form orbit with possible initial canonization
orb:=orbtranslimit(fct(i,One(gens[1])),fct,limit);
if orb<>fail and Length(orb[2])>1 then
if Length(orb[2])>bestl then
best:=orb;
bestl:=Length(orb[2]);
bestc:=1;
if bestl>100 then
return best;
fi;
else
bestc:=bestc+1;
# if we got 30 times the same length then use it.
if bestc>30 then
return best;
fi;
fi;
elif bestl>0 then
# we found at least one and others failed
return best;
fi;
od;
return best;
end;
siftchain:=function(use,x)
local i,img,u;
i:=1;
for u in use do
img:=u.fct(u.vec,x);
img:=LookupDictionary(u.dict,img);
if img=fail then
return [fail,i,x];
fi;
x:=x/u.t[img];
i:=i+1;
od;
return x;
end;
limit:=1000;
whole:=Group(mats);
f:=FieldOfMatrixList(mats);
total:=1;
gens:=mats;
mo:=GModuleByMats(gens,f);
bas:=MTX.BasesCompositionSeries(mo);
dims:=List(bas,Length);
if ForAll([2..Length(bas)],x->IsSubset(bas[x],bas[x-1])) then
basn:=List([2..Length(bas)],x->Filtered(bas[x],y->not y in bas[x-1]));
bas:=Concatenation(basn);
else
a:=ShallowCopy(bas[2]); # 1 is empty
p:=3;
while p<=Length(bas) do
t:=Difference(bas[p],a);
t:=Difference(t,bas[p-1]);
repeat
vec:=List([1..Length(bas[p])-Length(bas[p-1])],x->Random(t));
until RankMat(Concatenation(a,vec))=Length(bas[p]);
a:=Concatenation(a,vec);
p:=p+1;
od;
bas:=a;
fi;
use:=[];
while Length(gens)>0 do
basn:=List(bas,x->OnLines(x,One(gens[1])));
p:=PositionProperty(basn,x->ForAny(gens,y->OnLines(x,y)<>x));
if p=fail then
p:=PositionProperty(bas,x->ForAny(gens,y->x*y<>x));
vec:=bas{[p..Minimum(p+10,Length(bas))]};
fct:=OnRight;
elif Size(f)^p<limit then
p:=PositionProperty(dims,x->Size(f)^x>limit);
if p=fail then
p:=Length(dims);
else
p:=p-1;
fi;
orb:=OrbitsDomain(Group(gens),NormedRowVectors(VectorSpace(f,bas{[1..dims[p]]},Zero(bas[1]))),OnLines);
orb:=Filtered(orb,x->Length(x)>1);
Sort(orb,function(a,b) return Length(a)>Length(b);end);
vec:=List(orb,x->x[1]);
fct:=OnLines;
else
vec:=bas{[p..Minimum(p+10,Length(bas))]};
fct:=OnLines;
fi;
# nontrivial orbit
orb:=trymultipleorbits(vec,fct,limit);
if orb=fail then
b:=p;
a:=Reversed(Filtered(dims,x->x<b));
p:=fail;
while p=fail and a[1]<>0 do
sub:=a[1];
fct:=FactorspaceActfun(f,bas{[1..sub]});
basn:=List(bas,x->fct(x,One(gens[1])));
p:=Filtered([a[1]+1..Minimum(a[1]+10,Length(bas))],
x->ForAny(gens,y->fct(basn[x],y)<>basn[x]));
if p=[] then
p:=Filtered([a[1]+10..Length(bas)],
x->ForAny(gens,y->fct(basn[x],y)<>basn[x]));
fi;
if p=[] then
p:=fail;
else
# try to find one with orbit not just length p but also not too
# large
p:=Reversed(p);
while Length(p)>1 and Size(f)^(p[1]-a[1])/(Size(f)-1)>limit do
p:=p{[2..Length(p)]};
od;
p:=p[1];
fi;
a:=a{[2..Length(a)]};
od;
if p<>fail then
orb:=trymultipleorbits(bas{[p..Minimum(p+10,Length(bas))]},fct,limit);
if orb=fail then
# try dual action
fct:=AsInverseTranspose;
p:=First([b..Length(bas)],x->ForAny(gens,y->fct(bas[x],y)<>bas[x]));
if p<>fail then
p:=[p];
for orb in [1..5] do
Add(p,p[1]+orb);
Add(p,p[1]-orb);
od;
p:=Filtered(p,x->x>0 and x<Length(bas));
orb:=trymultipleorbits(bas{p},fct,limit);
fi;
if orb=fail then
Info(InfoFFMat,2,"even too long in dual");
return rec(points:=[],ops:=[]);
else
Info(InfoFFMat,2,"dual ");
fi;
else
Info(InfoFFMat,2,"factor space ");
fi;
else
return rec(points:=[],ops:=[]);
Error("p=fail -- should not happen");
fi;
fi;
dict:=orb[1];
t:=orb[3];
orb:=orb[2];
vec:=orb[1];
Info(InfoFFMat,2,"got Length ",Length(orb),":",
Collected(Factors(total*Length(orb))));
Add(use,rec(vec:=vec,fct:=fct,orb:=orb,dict:=dict,t:=t,gens:=gens));
#sift 20 random elements
new:=[];
p:=Maximum(Minimum(1000,Length(orb)*Length(gens)),10);
alt:=0;
a:=fail;
while a=fail and Length(new)<20 and alt<p do
a:=PseudoRandom(whole);
a:=siftchain(use,a);
alt:=alt+1;
if a[1]<>fail then
if fct(vec,a)<>vec then
Error("not stab");
fi;
if not IsOne(a) and not a in new then
Add(new,a);
alt:=0;
fi;
a:=fail; # sifted through
fi;
od;
if a<>fail then
gens:=ShallowCopy(use[a[2]].gens);
Add(gens,a[3]);
use:=use{[1..a[2]-1]};
total:=Product(List(use,x->Length(x.orb)));
Info(InfoFFMat,2,"construction fail on level ",a[2]," of ",Length(use),
" -- redo ",total);
else
total:=total*Length(orb);
gens:=new;
fi;
od;
return rec(points:=List(use,x->x.vec),ops:=List(use,x->x.fct),
order:=total);
end;
MATGRP_AddGeneratorToStabilizerChain:="2bdefined";
MATGRP_StabilizerChainInner:=
function( gens, size, layer, cand, opt, parentS )
# Computes a stabilizer chain for the group generated by gens
# with known size size (can be false if not known). This will be
# layer layer in the final stabilizer chain. cand is a (shared)
# record for base point candidates and opt the (shared) option
# record. This is called in StabilizerChain and calls itself.
# It also can be called if a new layer is needed.
local base,gen,S,i,merk,merk2,next,pr,r,stabgens,x;
Info(InfoGenSS,4,"Entering MATGRP_StabilizerChainInner layer=",layer);
next:=rec(point:=cand.points[1],op:=cand.ops[1]);
#next := GENSS_NextBasePoint(gens,cand,opt,parentS);
#cand := next.cand; # This could have changed
S := GENSS_CreateStabChainRecord(parentS,gens,size,
next.point,next.op,cand,opt);
base := S!.base;
Info( InfoGenSS, 3, "Entering orbit enumeration layer ",layer,"..." );
repeat
Enumerate(S!.orb,opt.OrbitLengthLimit);
if not(IsClosed(S!.orb)) then
if opt.FailInsteadOfError then
return "Orbit too long, increase opt.OrbitLengthLimit";
else
Error("Orbit too long, increase opt.OrbitLengthLimit");
fi;
fi;
until IsClosed(S!.orb);
Info(InfoGenSS, 2, "Layer ", layer, ": Orbit length is ", Length(S!.orb));
if layer > 1 then
parentS!.stab := S; # such that from now on random element
# generation works!
else
if (Length(S!.orb) > 50 or S!.orb!.depth > 5) and
S!.opt.OrbitsWithLog then
Error("ARGH!");
Info(InfoGenSS, 3, "Trying to make Schreier tree shallower (depth=",
S!.orb!.depth,")...");
merk := Length(S!.orb!.gens);
merk2 := Length(S!.stronggens);
MakeSchreierTreeShallow(S!.orb);
Append(S!.stronggens,S!.orb!.gens{[merk+1..Length(S!.orb!.gens)]});
Append(S!.layergens,[merk2+1..Length(S!.stronggens)]);
Info(InfoGenSS, 3, "Depth is now ",S!.orb!.depth);
fi;
fi;
S!.orb!.gensi := List(S!.orb!.gens,x->x^-1);
# as we add normalizing, we are done
return S;
# Are we done?
if size <> false and Length(S!.orb) = size then
S!.proof := true;
Info(InfoGenSS,4,"Leaving MATGRP_StabilizerChainInner layer=",layer);
return S;
fi;
# Now create a few random stabilizer elements:
stabgens := EmptyPlist(opt.RandomStabGens);
for i in [1..opt.RandomStabGens] do
x := GENSS_RandomElementFromAbove(S,layer);
if not(S!.IsOne(x)) then
Add(stabgens,x);
fi;
od;
Info(InfoGenSS,3,"Created ",opt.RandomStabGens,
" random stab els, ",
Length(stabgens)," non-trivial.");
if Length(stabgens) > 0 then # there is a non-trivial stabiliser
Info(InfoGenSS,3,"Found ",Length(stabgens)," non-trivial ones.");
if size <> false then
S!.stab := MATGRP_StabilizerChainInner(stabgens,size/Length(S!.orb),
layer+1,cand,opt,S);
else
S!.stab := MATGRP_StabilizerChainInner(stabgens,false,
layer+1,cand,opt,S);
fi;
if IsString(S!.stab) then return S!.stab; fi;
if opt.ImmediateVerificationElements > 0 then
Info(InfoGenSS,2,"Doing immediate verification in layer ",
S!.layer," (",opt.ImmediateVerificationElements,
" elements)...");
i := 0;
while i < opt.ImmediateVerificationElements do
i := i + 1;
x := GENSS_RandomElementFromAbove(S,layer);
if MATGRP_AddGeneratorToStabilizerChain(S!.topS,x) then
Info( InfoGenSS, 2, "Immediate verification found error ",
"(layer ",S!.layer,")..." );
i := 0;
fi;
od;
fi;
S!.proof := S!.stab!.proof; # hand up information
else
# We are not sure that the next stabiliser is trivial, but we believe!
Info(InfoGenSS,3,"Found no non-trivial ones.");
S!.proof := false;
fi;
Info(InfoGenSS,4,"Leaving MATGRP_StabilizerChainInner layer=",layer);
return S;
end;
MATGRP_AddGeneratorToStabilizerChain:=
function( S, el,pt,op )
# Increases the set represented by S by the generator el.
local SS, r, n, pr, i, newstrongnr;
if IsBound(S!.trivialgroup) and S!.trivialgroup then
if S!.IsOne(el) then
return false;
fi;
#SS := StabilizerChain(Group(el),S!.opt);
SS:=StabilizerChain(Group(el),rec(Cand:=rec(ops:=[op],points:=[pt]),Reduced:=false,StrictlyUseCandidates:=true));
if IsString(SS) then return SS; fi;
for n in NamesOfComponents(SS) do
S!.(n) := SS!.(n);
od;
Unbind(S!.trivialgroup);
return true;
fi;
r := SiftGroupElement( S, el );
# if this is one then el is already contained in the stabilizer chain
if r.isone then # already in the group!
return false;
fi;
# Now there remain two cases:
# (1) the sift stopped somewhere and we have to add a generator there
# (2) the sift ran all through the chain and the element still was not
# the identity, then we have to prolong the chain
if r.S <> false then # case (1)
SS := r.S;
Info( InfoGenSS, 2, "Adding new generator to stab. chain ",
"in layer ", SS!.layer, " from ",Length(SS!.stronggens) );
Add(SS!.stronggens,r.rem);
Add(SS!.layergens,Length(SS!.stronggens));
AddGeneratorsToOrbit(SS!.orb,[r.rem]);
Add(SS!.orb!.gensi,r.rem^-1);
newstrongnr := Length(SS!.stronggens);
Info( InfoGenSS, 4, "Entering orbit enumeration layer ",SS!.layer,
"..." );
repeat
Enumerate(SS!.orb,S!.opt.OrbitLengthLimit);
if not(IsClosed(SS!.orb)) then
if S!.opt.FailInsteadOfError then
return "Orbit too long, increase S!.opt.OrbitLengthLimit";
else
Error("Orbit too long, increase S!.opt.OrbitLengthLimit!");
fi;
fi;
until IsClosed(SS!.orb);
Info( InfoGenSS, 4, "Done orbit enumeration layer ",SS!.layer,"..." );
SS!.proof := false;
else # case (2)
# Note that we do not create a pr instance here for one
# generator, this will be done later on as needed...
SS := r.preS;
newstrongnr := Length(SS!.stronggens)+1; # r.rem will end up there !
SS!.stab := MATGRP_StabilizerChainInner([r.rem],false,
SS!.layer+1,rec(points:=[pt],ops:=[op]), SS!.opt, SS );
if IsString(SS!.stab) then return SS!.stab; fi;
SS := SS!.stab;
fi;
# Now we have added a new generator (or a new layer) at layer SS,
# the new gen came from layer S (we were called here, after all),
# thus we have to check, whether all the orbits between S (inclusively)
# and SS (exclusively) are also closed under the new generator r.rem,
# we add it to all these orbits, thereby also making the Schreier trees
# shallower:
while S!.layer < SS!.layer do
Info(InfoGenSS,2,"Adding new generator to orbit in layer ",S!.layer);
Add(S!.layergens,newstrongnr);
AddGeneratorsToOrbit(S!.orb,[r.rem]);
Add(S!.orb!.gensi,r.rem^-1);
S := S!.stab;
od;
# Finally, we have to add it to the product replacer!
AddGeneratorToProductReplacer(S!.opt!.pr,r.rem);
return true;
end;
SolvableBSGS:=function(arg)
local CBase, normalizingGenerator,df,ops,firstmoved,i,
solvNC,S,pcgs,x,r,c,w,a,bound,U,xp,depths,oldsz,prime,relord,gens,acter,ogens,stabs,n,strongs,stronglevs,laynums,slvec,layerzero,p,laynum,layers,sel,vals,stronglayers,layervecs,slpval,slp,baspts,levp,blocksz,lstrongs,lstrongsinv,bl,opt,check,primes,shortlim,orb,orbs,j,sz,goodbase,CHAINTEST,orblens;
goodbase:=[];
CHAINTEST:=function(X,str)
while X<>false do
#if IsBound(X!.stronggens) and
# Length(X!.layergens)>Length(Factors(Size(X))) then
# Error("UGH");
#fi;
if not X!.orb!.orbit[1] in goodbase.points then
Error("new point!");
fi;
#if Length(X!.orb!.gens)<>Length(Set(X!.orb!.gens)) then
# Error("duplicate!");
#fi;
#if IsBound(X!.opt) and IsBound(X!.opt.StrictlyUseCandidates) and
# X!.opt.StrictlyUseCandidates=false then Error("eh5!"); fi;
if IsBound(X!.orb!.gens) and Length(X!.orb!.gens)<>Length(Factors(Size(X))) then
Error("length!");
fi;
if IsBound(X!.orb!.gens) and Length(Difference(X!.orb!.gens,strongs))>0 then
Error("XTRA!");
fi;
if IsBound(X!.orb!.gensi) and List(X!.orb!.gens,Inverse)<>X!.orb!.gensi then
Error(str,"inverse!");
fi;
if IsBound(X!.orb!.gensi) and ForAny(X!.orb!.schreiergen,IsInt) and
Maximum(Filtered(X!.orb!.schreiergen,IsInt))>Length(X!.orb!.gensi)
then
Error("length!");
fi;
X:=X!.stab;
od;
end;
gens:=arg[Minimum(2,Length(arg))];
if IsGroup(gens) then
a:=gens;
gens:=GeneratorsOfGroup(gens);
else
a:=Group(gens);
fi;
if false and Length(arg)>2 and Length(gens)>5+Length(Factors(arg[3])) then
gens:=List([1..5+Length(Factors(arg[3]))],x->PseudoRandom(a));
gens:=Set(gens);
a:=Group(gens);
check:=true;
#SetSize(a,sz);
else
check:=false;
fi;
if IsBound(arg[3]) then
sz:=arg[3];
#SetSize(a,sz);
else
sz:=fail;
fi;
SetIsSolvableGroup(a,true);
if Length(arg)=1 then
acter:=gens;
else
acter:=arg[1];
if IsGroup(acter) then
acter:=GeneratorsOfGroup(acter);
fi;
fi;
shortlim:=Size(DefaultFieldOfMatrixGroup(a))*3000;
df:=DefaultFieldOfMatrixGroup(a);
if false then
# try vectors from big group
w:=BasisVectorsForMatrixAction(Group(acter));
w:=ImmutableMatrix(df,w);
if Size(df)=2 then
ops:=List(w,x->OnRight);
else
ops:=List(w,x->OnLines);
fi;
fi;
w:=ModuleStructureBase(gens);
if IsBound(arg[3]) and IsBound(w.order) and w.order<arg[3] then
for i in gens[1] do
Add(w.points,i);
Add(w.ops,OnRight);
od;
fi;
opt:=rec(RandomStabGens:=10,
Cand:=rec(points:=w.points,ops:=w.ops),
#TODO XXX Find better strategy for limit iof field is larger
VeryShortOrbLimit:=shortlim
);
FUNCSPACEHASH:=[]; # needed in base point selection code
CBase:=StabilizerChain(a,opt);
if sz<>fail and Size(CBase)<sz then
Info(InfoWarning,1,"Wrong size -- redo with `doall' option");
return fail;
fi;
primes:=Set(Factors(Size(CBase)));
FUNCSPACEHASH:=[];
w:=BasePointsActionsOrbitLengthsStabilizerChain(CBase);
Info(InfoGenSS,2,"Suggested Base Points",List(w,x->Position(opt.Cand.points,x[1])),
" Lengths ",List(w,x->x[3]));
goodbase:=BaseStabilizerChain(CBase);
w:=1;
opt:=1;
# Dixon's (1968 - the solvable length of a solvable linear group) bounds
if IsPerm(gens[1]) then
a:=Length(MovedPoints(gens));
bound := Int( LogInt( Maximum(1,a ^ 5), 3 ) / 2 );
elif IsMatrixOrMatrixObj(gens[1]) then
a:=NrRows(gens[1]); # dimension
# we would need proper log
#bound := Int( (8.55*Log(a+1,10)+0.36 );
# since it does not exist, replace 8.55 by 9 and simplify rounded up.
# this anyhow only matters for catching wrong input, so its harmless
bound := Log((a+1)^9,10)+1;
else
# no bound known
bound:=infinity;
fi;
normalizingGenerator:=function(g)
local oldstrong,oldsz,phq,pbp,pba;
oldsz:=Size(S);
oldstrong:=ShallowCopy(StrongGenerators(S));
# work around the ommission of prior redundant base points
pba:=Filtered([1..Length(goodbase.points)],x->goodbase.ops[x](goodbase.points[x],g)<>goodbase.points[x]);
pbp:=pba[1];
#Print("\n");
#View(S);
#Print("\n use ",pba,"\n");
MATGRP_AddGeneratorToStabilizerChain(S,g,goodbase.points[pbp],goodbase.ops[pbp]);
while Size(S)=oldsz do
Error("orderr BUG !!!\n");
MATGRP_AddGeneratorToStabilizerChain(S,g);
od;
return Difference(StrongGenerators(S),oldstrong);
end;
solvNC:=function()
local N,process,layergens,U,g,h,r,V,x,comm,pow,wp,sift;
N:=StabilizerChain(Group(StrongGenerators(S)),rec(Base:=CBase,Size:=Size(S),Reduced:=false,StrictlyUseCandidates:=true)); # really should be copy
# determine relative order
pow:=2;
wp:=w*w;
while not SiftGroupElement(N,wp).isone do
wp:=wp*w;
pow:=pow+1;
od;
if not IsPrimeInt(pow) then
prime:=Factors(pow)[1];
w:=w^(pow/prime);
else
prime:=pow;
fi;
Info(InfoFFMat,2,"Relative order ",pow," prime ",prime);
process:=[w];
layergens:=[];
U:=[];
for g in process do
sift:=SiftGroupElement(S,g);
if not sift.isone then
Info(InfoFFMat,2,"SS=",Size(S)," ",Length(U));
for h in layergens do
comm:=Comm(g,h);
#Print(Order(comm),"\n");
if not SiftGroupElement(N,comm).isone then
# wrong element -- get new generator for layer below
a:=false;
w:=comm;
S:=N; # don't we want to forget what we added?
return false;
fi;
od;
# the sifted element will be as good as generator, but allows for
# cleaner decomposition.
g:=sift.rem;
V:=normalizingGenerator(g);
#g:=V[1];
if g in U then Error("duplicate");fi;
U:=Concatenation([g],U);
Add(layergens,g);
for x in acter do
Add(process,g^x);
od;
fi;
od;
a:=true;
return U;
end;
S:=StabilizerChain(Group(One(gens[1])),rec(Base:=CBase,Reduced:=false,StrictlyUseCandidates:=true));
pcgs:=[];
depths:=[];
relord:=[];
xp:=1;
while xp<=Length(gens) do
Info(InfoFFMat,2,"Processing ",xp);
x:=gens[xp];
# feed through derived series if necessary
while not SiftGroupElement(S,x).isone do
c:=0;
w:=x;
a:=false;
oldsz:=Size(S);
while not a do
c:=c+1;
if c>bound then
# not solvable
Error("not solvable");
fi;
U:=solvNC(); # will change a,w
od;
Info(InfoFFMat,2,"Found Layer ",Size(S)/oldsz);
#View(S);
Info(InfoFFMat,2,"n");
if prime^Length(U)<>Size(S)/oldsz then Error("layerlen"); fi;
pcgs:=Concatenation(U,pcgs);
Add(depths,Length(pcgs));
Append(relord,ListWithIdenticalEntries(Length(U),prime));
od;
xp:=xp+1;
od;
n:=Length(pcgs);
depths:=Reversed(List(depths,x->n+1-x));
relord:=Reversed(relord);
Add(depths,n+1);
w:=BasePointsActionsOrbitLengthsStabilizerChain(S);
Info(InfoGenSS,2,"USED Base Points",List(w,x->Position(goodbase.points,x[1])),
" Lengths ",List(w,x->x[3]));
# now build the chain once more with memory so we can decompose
blocksz:=[];
Info(InfoFFMat,2,"NOW",Size(S)," ",Length(pcgs));
if Length(Factors(Size(S)))<>Length(pcgs) then Error("redundancies");fi;
# now sift the generators through so that the strong generators *ARE* a pcgs
gens:=Reversed(pcgs); # start with low level gens again
pcgs:=[];
acter:=[];
#S:=StabilizerChain(Group(One(gens[1])),rec(Base:=CBase,Reduced:=false,StrictlyUseCandidates:=true));
#S := GENSS_CreateStabChainRecord(false,
# [],1,
# goodbase.points[1],
# goodbase.ops[1],goodbase,
# rec(Base:=CBase,Reduced:=false,StrictlyUseCandidates:=true,InitialHashSize:=10));
S:=false; # work around special treatment for trivial group.
oldsz:=1;
stabs:=[];
xp:=1;
strongs:=[];
stronglevs:=[];
while xp<=Length(gens) do
if S<>false then
#CHAINTEST(S,"E");
oldsz:=Size(S);
Info(InfoFFMat,2,"ProcessiNg ",xp," ",Size(S));
#if Length(StrongGenerators(S))>Length(Factors(Size(S))) then Error("more2\n");fi;
x:=gens[xp];
x:=SiftGroupElement(S,x).rem;
#SiftGroupElement(S,x);
if not IsOne(x) then
normalizingGenerator(x);
fi;
else
x:=gens[1];
# OrbitsWithLog will add powers as strong generators to make the tree shallower.
# This messes with the correspondence of strong generators and pcgs elements
S:=StabilizerChain(Group(x),rec(Base:=CBase,Reduced:=false,StrictlyUseCandidates:=true,
OrbitsWithLog := false,RandomStabGens:=1,Size:=Order(x)));
# Remove those ~!@#$ extra generators (powers of x) that are added
# to make the tree shallower, but cause problems here.
strongs:=[x];
a:=S;
while a<>false do
a!.stronggens:=ShallowCopy(strongs);
a!.layergens:=[1];
if Length(a!.orb!.gens)>0 then
if a!.orb!.gens<>strongs then
a!.orb:=Orb(ShallowCopy(strongs),a!.orb[1],a!.orb!.op,rec(schreier:=true));
a!.orb!.gensi:=List(a!.orb!.gens,Inverse);
repeat
Enumerate(a!.orb);
until IsClosed(a!.orb);
a!.orb!.depth:=Size(a!.orb)-1;
#if Length(a!.orb)=1 then
# a!.orb!.schreiergen:=ShallowCopy(strongs);
#else
# a!.orb!.schreiergen:=[];
#fi;
fi;
fi;
a:=a!.stab;
od;
#CHAINTEST(S,"F2");
strongs:=[];
fi;
if Size(S)=oldsz then Error("no change!");fi;
Add(pcgs,x);
a:=Set(Filtered(StrongGenerators(S),x->not x in strongs));
a:=Difference(a,List([0..Order(x)],y->x^y));
if Length(a)>0 then
# redo
Info(InfoFFMat,1,"nanu-redo");
#Error("nanu");
return SolvableBSGS(arg[1],arg[2],arg[3]);
fi;
Append(strongs,[x]);
#CHAINTEST(S,"F");
Append(stronglevs,ListWithIdenticalEntries(1,n+1-xp));
#if n+1-xp in depths then
# Add(stabs,StructuralCopy(S));
#fi;
xp:=xp+1;
od;
pcgs:=Reversed(pcgs);
CBase:=BaseStabilizerChain(S);
w:=BasePointsActionsOrbitLengthsStabilizerChain(S);
Info(InfoGenSS,2,"Used Base Points",List(w,x->Position(goodbase.points,x[1])),
" Lengths ",List(w,x->x[3]));
baspts:=[];
levp:=[];
xp:=[1..Length(pcgs)];
a:=S;
for i in [1..Length(CBase.points)] do
p:=List(xp,x->Position(a!.orb!.orbit,CBase.ops[i](CBase.points[i],pcgs[x])));
sel:=Filtered([1..Length(xp)],x->p[x]<>1);
baspts{xp{sel}}:=ListWithIdenticalEntries(Length(sel),i);
levp[i]:=xp{sel};
blocksz{xp{sel}}:=p{sel}-1;
xp:=Difference(xp,xp{sel});
a:=a!.stab;
od;
slpval:=function(w)
local a,i;
a:=w[2]*vals[w[1]];
for i in [3,5..Length(w)-1] do
a:=(a+w[i+1]*vals[w[i]]) mod p;
od;
return a;
end;
layerzero:=[];
laynum:=Length(depths)-1;
layers:=List([1..n],x->First([laynum,laynum-1..1],y->x>=depths[y]));
stronglayers:=layers{stronglevs};
# get vectors for strong generators (b/c we decompose into these)
slvec:=[];
for i in [1..laynum] do
p:=relord[depths[i]];
a:=IdentityMat(depths[i+1]-depths[i]);
layerzero[i]:=Zero(a[1]);
sel:=Positions(stronglayers,i);
slvec{sel}:=a{List(strongs{sel},x->Position(pcgs,x))-depths[i]+1};
# layervecs:=ListWithIdenticalEntries(n,fail); # this will trap use of
# # higher gens
# layervecs{[depths[i]..depths[i+1]-1]}:=a;
#
# for j in [depths[i+1]..n] do
# layervecs[j]:=Zero(a[1]);
# od;
# slp:=LinesOfStraightLineProgram(SLPOfElms(strongs{sel}));
# vals:=Reversed(layervecs);
# for j in slp do
# if ForAll(j,IsList) then
# #last line
# slvec{sel}:=List(j,x->slpval(x));
# elif IsList(j[1]) then
# Error("reassign syntax?");
# else
# Add(vals,slpval(j));
# fi;
# od;
od;
ForgetMemory(strongs);
ForgetMemory(gens);
ForgetMemory(S);
S!.pcgs:=pcgs;
S!.layerzero:=layerzero;
S!.layerprimes:=relord{depths{[1..laynum]}};
S!.strongvecs:=slvec;
S!.layranges:=List([1..laynum],x->[depths[x]..depths[x+1]-1]);
S!.revranges:=List([1..laynum],x->Reversed([1..Length(S!.layranges[x])]));
# reindexing
a:=S;
while a<>false do
p:=List(a!.orb!.gensi,x->Position(strongs,x^-1));
if fail in p then
Error("not found!");
fi;
a!.gensistrongpos:=p;
a!.gensilayers:=stronglayers{p};
a:=a!.stab;
od;
#TODO: Use chain for decomposition -- at the moment it uses subgroups
a:=pcgs;
pcgs:=PcgsByPcSequenceCons(IsPcgsDefaultRep,
IsPcgsMatGroupByStabChainRep and IsPcgs and IsPrimeOrdersPcgs,
FamilyObj(gens[1]),pcgs,[]);
pcgs!.stabilizerChain:=S;
#pcgs!.stabs:=Reversed(stabs);
pcgs!.laynums:=[1..laynum];
pcgs!.layranges:=S!.layranges;
pcgs!.baspts:=baspts;
pcgs!.blocksz:=blocksz;
pcgs!.levp:=levp;
pcgs!.invpows:=
List([1..Length(pcgs)],x->List([1..relord[x]-1],y->a[x]^(-y)));
SetRelativeOrders(pcgs,relord);
SetIndicesEANormalSteps(pcgs,depths);
return rec(
pcgs:=pcgs,
depths:=depths,
relord:=relord,
baspts:=baspts,
blocksz:=blocksz);
end;
# old ersion of exponent routines, allowing for arbitrary orbit arrangemnts
MatPcgsSift:=function( S, x,l )
local o,p,po,preS,r,v,vecs,prime;
v:=S!.layerzero[l];
prime:=S!.layerprimes[l];
vecs:=S!.strongvecs;
preS := false;
while S <> false do
o := S!.orb;
if IsObjWithMemory(x) then
p := o!.op(o[1],x!.el);
else
p := o!.op(o[1],x);
fi;
po := Position(o,p);
if po = fail then # not in current stabilizer
Error("element not in group");
return rec( v:=v,rem := x, S := S );
fi;
# Now sift through Schreier tree:
while po > 1 do
r:=o!.schreiergen[po];
x := x * S!.orb!.gensi[r];
po := o!.schreierpos[po];
if S!.gensilayers[r]=l then
# generator on the desired layer -- add up vectors
v:=(v+vecs[S!.gensistrongpos[r]]) mod prime;
fi;
od;
S := S!.stab;
od;
return v;
end ;
MatPcgsExponentsOld:=function(S,laynums,x)
local a,j,i,v;
v:=[];
for i in [1..Length(laynums)] do
#S:=stabs[laynums[i]];
a:=MatPcgsSift(S,x,laynums[i]);
Append(v,a);
if i<Length(laynums) then
# need to divide off what is given by the vector
for j in [1..Length(a)] do
if a[j]<>0 then
x:=LeftQuotient(S!.pcgs[S!.layranges[laynums[i]][j]]^a[j],x);
fi;
od;
fi;
od;
return v;
end;
# new routine, using orbit positions
#decomposition as product, but not in canonical order, thus in multi stages
#TODO compare cost matrix multiplication vs. collection
MatPcgsExponents:=function(arg)
local pcgs,laynums,ox,o,p,po,preS,r,isone,ind,i,prd,S,q,rem,bs,pS,x,dep,e,layer,delta,
z,prime,curran,seli,md,pos,sel,deponly;
pcgs:=arg[1];
laynums:=arg[2];
ox:=arg[3];
deponly:=Length(arg)>3 and arg[4];
dep:=IndicesEANormalSteps(pcgs);
md:=laynums[Length(laynums)]; #the layer depth which we ignore.
z:=ListWithIdenticalEntries(dep[Length(dep)]-1,0);
e:=ShallowCopy(z);
pS:=pcgs!.stabilizerChain;
repeat
# factor the element using orbit positions -- this is correct only for
# the topmost layer used
x:=ox;
S:=pS;
prd:=[];
pos:=[];
preS := false;
ind:=1;
layer:=md; # Maximal layer we do
delta:=ShallowCopy(z); # we need to forget all in the previous layer
curran:=pcgs!.layranges[layer];
prime:=RelativeOrders(pcgs)[curran[1]];
while S <> false do
sel:=pcgs!.levp[ind];
bs:=pcgs!.blocksz{sel};
o := S!.orb;
if IsObjWithMemory(x) then
p := o!.op(o[1],x!.el);
else
p := o!.op(o[1],x);
fi;
po := Position(o,p)-1;
# decompose
for i in [1..Length(sel)] do
seli:=sel[i];
q:=QuoInt(po,bs[i]);
rem:=po-q*bs[i];
if q>0 then
if seli<dep[layer] then
# decrease layer in which we decompose
layer:=layer-1;
while seli<dep[layer] do
layer:=layer-1;
od;
delta:=ShallowCopy(z); # we need to forget all in the previous layer
prime:=RelativeOrders(pcgs)[seli];
curran:=pcgs!.layranges[layer];
fi;
#Add(prd,[sel[i],q]);;
# remember exponent entry if right layer
if seli<dep[layer+1] then
delta[seli]:=delta[seli]+q mod prime;
fi;
#x:=x/(pcgs[sel[i]]^q);
x:=x*pcgs!.invpows[seli][q]; # use stored inverse powers
fi;
po:=rem;
od;
#if o[1]<>o!.op(o[1],x) then Error("not all off!");fi;
S := S!.stab;
ind:=ind+1;
od;
if delta<>fail then
# now the affected (topmost) layer is correct
e:=e+delta;
# we don't need to correct the last layer we use
if deponly and not IsZero(delta) then
delta:=fail; # pop out.
elif layer<laynums[Length(laynums)] then
# we have d describe a the highest level. So we want x=d*newx
for i in curran do
if delta[i]>0 then
ox:=pcgs!.invpows[i][delta[i]]*ox;
fi;
od;
fi;
fi;
# stop when we've set the lowest layer or if the remainder is the identity
until layer=laynums[Length(laynums)] or delta=fail;
# do we only want some?
if laynums=pcgs!.laynums then
return e;
else
return e{Concatenation(pcgs!.layranges{laynums})};
fi;
end;
InstallMethod(ExponentsOfPcElement,"matrix pcgs",IsCollsElms,
[IsPcgsMatGroupByStabChainRep and IsPcgs and IsPrimeOrdersPcgs,
IsMatrixOrMatrixObj],
function(pcgs,x)
return MatPcgsExponents(pcgs,pcgs!.laynums,x);
end);
InstallOtherMethod(ExponentsOfPcElement,"matrix pcgs,indices",IsCollsElmsX,
[IsPcgsMatGroupByStabChainRep and IsPcgs and IsPrimeOrdersPcgs,
IsMatrixOrMatrixObj,IsList],
function(pcgs,x,inds)
local i,j,sel,lay,ip,sl,r,use;
# is it a layer?
for i in pcgs!.laynums do
if inds=pcgs!.layranges[i] then
# only this layer
return MatPcgsExponents(pcgs,[i],x);
fi;
od;
# which layers do we need?
if not IsSortedList(inds) then
# perverse case -- old
return MatPcgsExponents(pcgs,pcgs!.laynums,x){inds};
fi;
sel:=[];
lay:=[];
ip:=1;
sl:=0;
for i in pcgs!.laynums do
r:=pcgs!.layranges[i];
use:=false;
for j in [1..Length(r)] do
if ip<=Length(inds) and r[j]=inds[ip] then
Add(sel,j+sl);
use:=true;
ip:=ip+1;
fi;
od;
if use then
AddSet(lay,i);
sl:=sl+Length(r);
fi;
od;
return MatPcgsExponents(pcgs,lay,x){sel};
end);
InstallMethod(DepthOfPcElement,"matrix pcgs",IsCollsElms,
[IsPcgsMatGroupByStabChainRep and IsPcgs and IsPrimeOrdersPcgs,
IsMatrixOrMatrixObj],
function(pcgs,x)
local e;
e:=MatPcgsExponents(pcgs,pcgs!.laynums,x,true);
return PositionNonZero(e);
end);
BindGlobal("TFDepthLeadExp",function(pcgs,x)
local p;
x:=MatPcgsExponents(pcgs,pcgs!.laynums,x,true); # first nonzero
p:=PositionNonZero(x);
if p>Length(x) then
# identity
return [p,0];
else
return [p,x[p]];
fi;
end);
InstallMethod(FittingFreeLiftSetup,"fields, using recognition",true,
[IsMatrixGroup],0,
function(G)
local csi,r,factorhom,sbs,k,pc,hom,rad,it,i,sz,x,stop;
r:=DefaultFieldOfMatrixGroup(G);
if not IsField(r) then
TryNextMethod();
fi;
r:=RecognizeGroup(G);
SetSize(G,Size(r));
csi:=GetInformationFromRecog(r);
G!.storedrecog:=csi;
factorhom:=FindAct(csi);
#TODO: Better kernel gens by random selection
sz:=Size(G)/Size(Image(factorhom));
if Size(Image(factorhom))=1 then
# solvable
k:=GeneratorsOfGroup(G);
else
it:=CoKernelGensIterator(InverseGeneralMapping(factorhom));
k:=Filtered(List([1..csi.genum],x->CSINiceGens(csi,x)),x->IsOne(ImagesRepresentative(factorhom,x)));
if sz>1 then
stop:=true;
repeat
for i in [1..3*Length(Factors(sz))] do
if not IsDoneIterator(it) then
x:=NextIterator(it);
if not IsOne(x) and not x in k then
Add(k,x);
fi;
fi;
od;
if sz>1 and Length(k)<Length(Factors(sz)) then stop:=false; fi; # work around issue
if ValueOption("doall")=true then stop:=false;fi;
until stop or IsDoneIterator(it);
fi;
fi;
Info(InfoGenSS,3,"|k|=",Length(k));
# TODO: Proper test
if ForAll(k,IsOne) then
rad:=TrivialSubgroup(G);
sbs:=rec(pcgs:=[],depths:=[1],relord:=[],pcisom:=IsomorphismPcGroup(rad),radical:=rad);
else
sbs:=SolvableBSGS(G,k,sz);
while sbs=fail do
sbs:=SolvableBSGS(G,k,sz:doall);
od;
rad:=SubgroupNC(G,sbs.pcgs);
SetSize(rad,Product(RelativeOrders(sbs.pcgs)));
sbs.radical:=rad;
pc:=PcGroupWithPcgs(sbs.pcgs);
RUN_IN_GGMBI:=true; # hack to skip Nice treatment
hom:=GroupHomomorphismByImagesNC(rad,pc,sbs.pcgs,FamilyPcgs(pc));
RUN_IN_GGMBI:=false;
SetIsBijective(hom,true);
sbs.pcisom:=hom;
fi;
sbs.csi:=csi;
SetKernelOfMultiplicativeGeneralMapping(factorhom,rad);
sbs.factorhom:=factorhom;
RecogDecompinfoHomomorphism(factorhom)!.LiftSetup:=sbs;
return sbs;
end);
FFStats:=function(g)
local start,f;
start:=Runtime();
f:=FittingFreeLiftSetup(g);
Print("Time:",Runtime()-start,"\n");
Print("Factordegree ",NrMovedPoints(Range(f.factorhom)),"\n");
Print("PcgsDegrees ",Maximum(List(
BasePointsActionsOrbitLengthsStabilizerChain(f.pcgs!.stabilizerChain),
x->x[3])),"\n");
end;
# work over residue class rings
MyZmodnZObj:=function(a,b)
if IsPrimeInt(b) and b<65536 then
return Z(b)^0*a;
else
return ZmodnZObj(a,b);
fi;
end;
ReduceModM:=function(a,m)
local b,r,i,j;
if IsList(a) then
if IsList(a[1]) then
# matrix
b:=[];
for i in a do
r:=[];
for j in i do
if IsZmodnZObjNonprime(j) then
Add(r,MyZmodnZObj(j![1],m));
else
Add(r,MyZmodnZObj(j,m));
fi;
od;
Add(b,r);
od;
if IsPrimeInt(m) then
b:=ImmutableMatrix(m,b);
fi;
return b;
else
# vector
b:=[];
for j in a do
--> --------------------
--> maximum size reached
--> --------------------
[ Verzeichnis aufwärts0.73unsichere Verbindung
Übersetzung europäischer Sprachen durch Browser
]
|
