|
|
|
|
Quelle genss.gi
Sprache: unbekannt
|
|
#############################################################################
##
## genss.gi genss package
## Max Neunhoeffer
## Felix Noeske
##
## Copyright 2006 Lehrstuhl D für Mathematik, RWTH Aachen
##
## Implementation stuff for generic Schreier-Sims
##
#############################################################################
#############################################################################
# The following global record contains default values for options for the
# main function "StabilizerChain":
#############################################################################
# Initial hash size for orbits:
GENSS.InitialHashSize := NextPrimeInt(1000);
# Number of points to process before reporting:
GENSS.Report := 30000;
# Number of random elements to consider for the determination of short orbits:
GENSS.ShortOrbitsNrRandoms := 10;
# Absolute limit for orbit length in search for short orbits:
GENSS.ShortOrbitsOrbLimit := 80000;
# First limit for the parallel enumeration of orbits looking for short orbits:
GENSS.ShortOrbitsInitialLimit := 400;
# Absolute limit for single orbit length:
GENSS.OrbitLengthLimit := 10000000;
# Number of points in the previous orbit to consider for the next base point:
GENSS.NumberPrevOrbitPoints := 10;
# Number of (evenly distributed) random generators for the stabilizer:
GENSS.RandomStabGens := 3;
# Product replacement parameters for the stabilizer element generation:
# Now actually used for the generation of all random elements on top
# level. Random elements further down are created on top and then
# sifted.
GENSS.StabGenScramble := 30;
GENSS.StabGenScrambleFactor := 6;
GENSS.StabGenAddSlots := 3;
GENSS.StabGenMaxDepth := 400;
# Number of random elements used for verification,
# note that this is changed by the "random" and "ErrorBound" options!
GENSS.VerifyElements := 10; # this amounts to an error bound of 1/1024
GENSS.DeterministicVerification := false;
# Number of random elements used to do immediate verification:
GENSS.ImmediateVerificationElements := 3;
# Are we working projectively?
GENSS.Projective := false;
# To find a very short orbit we try two basis vectors and a random vector:
GENSS.VeryShortOrbLimit := 500;
# Never consider more than this number of candidates for short orbits:
GENSS.LimitShortOrbCandidates := 50;
# Do not throw Errors but return fail:
GENSS.FailInsteadOfError := false;
# Number of Schreier generators to create in TC verification:
GENSS.NumberSchreierGens := 20;
# Maximal number of Schreier generators to create in TC verification:
GENSS.MaxNumberSchreierGens := 100;
# By default do 3 rounds of birthday paradox method:
GENSS.TryBirthdayParadox := 3;
# By default do 1 short orbit tries:
GENSS.TryShortOrbit := 1;
# Limit for orbit length during orbit estimation using birthday paradox:
GENSS.OrbitLimitBirthdayParadox := 1000000;
# We immediately take an orbit if its estimate is lower than this limit:
GENSS.OrbitLimitImmediatelyTake := 10000;
# Limit for number of random elements during orbit estimation:
GENSS.NrRandElsBirthdayParadox := 6000;
# Set this to false to allow orbits of length 1:
GENSS.Reduced := true;
# Product replacement parameters for Stab:
GENSS.StabScramble := 10;
GENSS.StabScrambleFactor := 1;
GENSS.StabAddSlots := 1;
GENSS.StabMaxDepth := 400;
# Initial limit for orbit length for Stab computation:
GENSS.StabInitialLimit := 1000;
# Patience to find random elements mapping the orbit into itself:
GENSS.StabInitialPatience := 10;
# Maximal amount of memory used for the orbit:
GENSS.StabOrbitLimit := 1000000;
# If the probability of a wrong stabiliser is smaller than this, do no
# longer try to create stabiliser elements:
GENSS.StabAssumeCompleteLimit := 1/(10^7);
# The following switches off the log used in all orbits for stabilizer
# chains if set to false. Do not do this unless you know what you are
# doing since since could prevent Schreier trees from being shallow.
GENSS.OrbitsWithLog := true;
#############################################################################
# A few helper functions needed elsewhere:
#############################################################################
InstallGlobalFunction( GENSS_CopyDefaultOptions,
function( defopt, opt )
local n;
for n in RecNames(defopt) do
if not(IsBound(opt.(n))) then
opt.(n) := defopt.(n);
fi;
od;
end );
InstallGlobalFunction( GENSS_MapBaseImage,
function( bi, el, S )
# bi must be a base image belonging to S, el a group element
local i,l,res;
l := Length(bi);
res := 0*[1..l];
i := 1;
while true do
res[i] := S!.orb!.op(bi[i],el);
i := i + 1;
if i = l+1 then break; fi;
S := S!.stab;
od;
return res;
end );
InstallGlobalFunction( GENSS_FindVectorsWithShortOrbit,
# This implements Murray/O'Brien-like heuristics.
# It produces new random elements using GENSS_RandomElementFromAbove
# and stores them in opt.FindBasePointCandidatesData.randpool for
# later usage.
function(g,opt,parentS)
# Needs opt components "ShortOrbitsNrRandoms"
local l, f, data, x, c, onlydegs, v, vv, w, i, nw, inters, sortfun,
wb, j, ww;
Info(InfoGenSS,4,"Trying Murray/O'Brien heuristics...");
l := ShallowCopy(GeneratorsOfGroup(g));
f := DefaultFieldOfMatrixGroup(g);
data := opt.FindBasePointCandidatesData;
for i in [1..opt.ShortOrbitsNrRandoms] do
if parentS = false then
x := GENSS_RandomElementFromAbove(opt,0);
else
x := GENSS_RandomElementFromAbove(parentS,parentS!.layer);
fi;
Add(l,x);
Add(data.randpool,x);
od;
ForgetMemory(l);
c := List(l,x->Set(Factors(CharacteristicPolynomial(x,1):
onlydegs := [1..3])));
v := [];
for i in [1..Length(l)] do
for j in [1..Length(c[i])] do
vv := EmptyPlist(Length(c[i]));
Add(vv,[NullspaceMat(Value(c[i][j],l[i])),
Degree(c[i][j]),
WeightVecFFE(CoefficientsOfLaurentPolynomial(c[i][j])[1]),
1]);
od;
Add(v,vv);
od;
Info(InfoGenSS,4,"Have eigenspaces.");
# Now collect a list of all those spaces together with all
# possible intersects:
w := [];
i := 1;
while i <= Length(l) and Length(w) < GENSS.LimitShortOrbCandidates do
nw := [];
for j in [1..Length(v[i])] do
for ww in w do
inters := SumIntersectionMat(ww[1],v[i][j][1])[2];
if Length(inters) > 0 then
Add(nw,[inters,Minimum(ww[2],v[i][j][2]),
Minimum(ww[3],v[i][j][3]),ww[4]+v[i][j][4]]);
fi;
od;
Add(nw,v[i][j]);
od;
Append(w,nw);
i := i + 1;
od;
sortfun := function(a,b)
if a[2] < b[2] then return true;
elif a[2] > b[2] then return false;
elif a[3] < b[3] then return true;
elif a[3] > b[3] then return false;
elif a[4] < b[4] then return true;
elif a[4] > b[4] then return false;
elif Length(a[1]) < Length(b[1]) then return true;
else return false;
fi;
end;
Sort(w,sortfun);
wb := List(w,ww->ww[1][1]);
Info(InfoGenSS,3,"Have ",Length(wb)," vectors for possibly short orbits.");
for ww in [1..Length(wb)] do
if not(IsMutable(wb[ww])) then
wb[ww] := ShallowCopy(wb[ww]);
fi;
ORB_NormalizeVector(wb[ww]);
od;
return wb;
end );
InstallGlobalFunction( GENSS_FindShortOrbit,
function( g, opt, parentS )
# Needs opt components:
# "ShortOrbitsNrRandoms" (because it uses GENSS_FindVectorsWithShortOrbit)
# "ShortOrbitsOrbLimit"
# "ShortOrbitsInitialartLimit"
local ThrowAwayOrbit,found,gens,hashlen,i,j,limit,newnrorbs,nrorbs,o,wb;
wb := GENSS_FindVectorsWithShortOrbit(g,opt,parentS);
if Length(wb) = 0 then return fail; fi;
# Now we have a list of vectors with (hopefully) short orbits.
# We start enumerating all those orbits, but first only 50 elements:
nrorbs := Minimum(Length(wb),64); # take only the 64 first
gens := GeneratorsOfGroup(g);
o := [];
hashlen := NextPrimeInt(QuoInt(opt.ShortOrbitsOrbLimit,2));
for i in [1..nrorbs] do
Add(o,Orb(gens,ShallowCopy(wb[i]),OnLines,
rec( treehashsize := hashlen )));
od;
limit := opt.ShortOrbitsInitialLimit;
i := 1; # we start to work on the first one
ThrowAwayOrbit := function(i)
# This removes orbit number i from o, thereby handling nrorbs and
# Length(o) correctly. If you want to use o[i] further, please
# make a copy (of the reference) before calling this function.
if Length(o) > nrorbs then
o[i] := o[nrorbs+1];
o{[nrorbs+1..Length(o)-1]} := o{[nrorbs+2..Length(o)]};
Unbind(o[Length(o)]);
else
o{[i..nrorbs-1]} := o{[i+1..nrorbs]};
Unbind(o[nrorbs]);
nrorbs := nrorbs-1;
fi;
end;
repeat
Enumerate(o[i],limit);
found := IsClosed(o[i]);
if Length(o[i]) = 1 then
Info(InfoGenSS,3,"Orbit Number ",i," has length 1.");
found := false;
# Now throw away this orbit:
ThrowAwayOrbit(i);
# we intentionally do not increase i here!
elif not(found) then
i := i + 1;
fi;
if i > nrorbs then
Info(InfoGenSS,3,"Done ",nrorbs,
" orbit(s) to limit ",limit,".");
limit := limit * 2;
if limit > opt.ShortOrbitsOrbLimit then
Info(InfoGenSS,3,"Limit reached, giving up.");
return fail;
fi;
i := 1;
if nrorbs < i then
Info(InfoGenSS,3,"No orbits left, giving up.");
return fail;
fi;
if nrorbs > 1 then
newnrorbs := QuoInt((nrorbs+1),2);
for j in [newnrorbs+1..nrorbs] do
Unbind(o[j]);
od;
nrorbs := newnrorbs;
fi;
fi;
until found;
Info(InfoGenSS,2,"Found short orbit of length ",Length(o[i])," (#",i,").");
return o[i];
end );
InstallGlobalFunction( GENSS_IsOneProjective,
function(el)
local s;
s := el[1][1];
if IsZero(s) then return false; fi;
if not(IsOne(s)) then
s := s^-1;
el := s * el;
fi;
return IsOne( el );
end );
#############################################################################
# Now to the heart of the method, the Schreier-Sims:
#############################################################################
InstallMethod( FindBasePointCandidates,
"for a scalar matrix group over a FF",
[ IsGroup and IsMatrixGroup and IsFinite, IsRecord, IsInt, IsObject ],
50, # highest weight,
function( grp, opt, i, parentS )
local F, q, d, gens, v;
Info( InfoGenSS, 3, "Finding nice base points (scalar)..." );
F := DefaultFieldOfMatrixGroup(grp);
q := Size(F);
d := DimensionOfMatrixGroup(grp);
gens := GeneratorsOfGroup(grp);
if IsObjWithMemory(gens[1]) then
gens := StripMemory(gens);
fi;
if ForAny(gens,x->not(GENSS_IsOneProjective(x))) then
opt.FindBasePointCandidatesData := rec( randpool := [], vecs := [] );
# This is needed for communication between different methods
TryNextMethod();
fi;
v := ZeroMutable(gens[1][1]);
v[1] := PrimitiveRoot(F); # to indicate the field!
return rec( points := [v], ops := [OnRight], used := 0 );
end );
InstallMethod( FindBasePointCandidates,
"for a matrix group over a FF, very short orbit",
[ IsGroup and IsMatrixGroup and IsFinite, IsRecord, IsInt, IsObject ],
40, # highest weight,
function( grp, opt, i, parentS )
local F, q, d, gens, op, v, vv, k, kk, o, cand, j;
Info( InfoGenSS, 3, "Finding nice base points (very short)..." );
F := DefaultFieldOfMatrixGroup(grp);
q := Size(F);
d := DimensionOfMatrixGroup(grp);
gens := GeneratorsOfGroup(grp);
if IsObjWithMemory(gens[1]) then
gens := StripMemory(gens);
fi;
# Try two standard basis vectors and a random vector to find a very
# short orbit:
if q = 2 then
op := OnPoints;
else
op := OnLines;
fi;
v := [];
# Find first standard basis vector that is moved:
vv := ZeroMutable( gens[1][1] );
k := 1;
while k <= d do
vv[k] := One(F);
if ForAny(gens,x->op(vv,x) <> vv) then
Add(v,vv);
break;
fi;
vv[k] := Zero(F);
k := k + 1;
od;
# Find last standard basis vector that is moved:
vv := ZeroMutable( gens[1][1] );
kk := d;
while kk >= k do
vv[kk] := One(F);
if ForAny(gens,x->op(vv,x) <> vv) then
Add(v,vv);
break;
fi;
vv[kk] := Zero(F);
kk := kk - 1;
od;
# Pick a random vector:
vv := ZeroMutable( gens[1][1] );
if IsPlistRep(vv) then
for j in [1..Length(vv)] do
vv[j] := Random(F);
od;
else
Randomize(vv);
fi;
ORB_NormalizeVector(vv);
Add(v,vv);
# Now investigate these up to a certain limit:
for j in [1..Length(v)] do
o := Orb(gens,v[j],op,
rec(treehashsize := QuoInt(opt.VeryShortOrbLimit,2)+1));
Enumerate(o,opt.VeryShortOrbLimit);
if Length(o) > 1 and Length(o) < opt.VeryShortOrbLimit then
Info( InfoGenSS, 3, "Found orbit of length ",Length(o) );
cand := rec( points := [v[j]], ops := [op], used := 0 );
# Note that if we work non-projectively, then the same
# point will be taken next with OnRight automatically!
return cand;
fi;
od;
Info( InfoGenSS, 3, "Found no very short orbit up to limit ",
opt.VeryShortOrbLimit );
Append(opt.FindBasePointCandidatesData.vecs,v); # hand on vectors
TryNextMethod();
end );
# Method with rank 30 taking some commutators and invariant spaces
InstallMethod( FindBasePointCandidates,
"for a matrix group over a FF, using birthday paradox method",
[ IsGroup and IsMatrixGroup and IsFinite, IsRecord, IsInt, IsObject ], 20,
function( grp, opt, mode, parentS )
local F, q, d, randels, immorblimit, orblimit, data, op, v, l, c, e, ht,
val, x, w, cand, minest, minpos, round, i, j, gens;
F := DefaultFieldOfMatrixGroup(grp);
q := Size(F);
d := DimensionOfMatrixGroup(grp);
randels := opt.NrRandElsBirthdayParadox;
immorblimit := opt.OrbitLimitImmediatelyTake;
orblimit := opt.OrbitLimitBirthdayParadox;
Info( InfoGenSS, 3, "Finding base points (birthday paradox, limit=",
orblimit,", randels=",randels,")..." );
data := opt.FindBasePointCandidatesData; # this we get from earlier methods
if q = 2 then
op := OnPoints;
else
op := OnLines;
fi;
gens := GeneratorsOfGroup(grp);
if IsObjWithMemory(gens[1]) then
gens := StripMemory(gens);
fi;
for round in [1..opt.TryBirthdayParadox] do
v := Set(GENSS_FindVectorsWithShortOrbit(grp,opt,parentS));
if round = 1 then
Append(v,data.vecs); # take previously tried ones as well
fi;
v := Filtered(v,vv->ForAny(gens,x-> vv <> op(vv,x)));
l := Length(v);
if l = 0 then
# Find a vector on which grp acts non-trivially.
# At least one basis vector is guaranteed to be moved,
# unless grp is diagonal and op is OnLines, in which case
# they might all be fixed. However, in that case, the
# vector [1,1,...,1] is moved.
v := OneMutable(gens[1]); # List of basis vectors
Add(v, Sum(v)); # add vector [1,1,...,1]
v := Filtered(v,vv->ForAny(gens,x-> vv <> op(vv,x)));
l := Length(v);
fi;
c := 0*[1..l]; # the number of coincidences
e := ListWithIdenticalEntries(l,infinity); # the current estimates
ht := HTCreate(v[1]*PrimitiveRoot(F),
rec(hashlen := NextPrimeInt(l * randels * 4)));
for i in [1..l] do
val := HTValue(ht,v[i]);
if val = fail then
HTAdd(ht,v[i],[i]);
else
AddSet(val,i);
fi;
od;
for i in [1..randels] do
if parentS = false then
x := GENSS_RandomElementFromAbove(opt,0);
else
x := GENSS_RandomElementFromAbove(parentS,parentS!.layer);
fi;
Add(data.randpool,x);
for j in [1..l] do
if IsObjWithMemory(x) then
w := op(v[j],x!.el);
else
w := op(v[j],x);
fi;
val := HTValue(ht,w);
if val <> fail then # we know this point!
if j in val then # a coincidence!
c[j] := c[j] + 1;
e[j] := QuoInt(i^2,2*c[j]);
if (c[j] >= 3 and e[j] <= immorblimit) or
(c[j] >= 15 and e[j] <= orblimit) then
Info( InfoGenSS, 2, "Found orbit with estimated ",
"length ",e[j]," (coinc=",c[j],")" );
cand := rec(points := [v[j]], ops := [op],
used := 0);
for i in [1..l] do
if i <> j and c[i] >= 10 and
e[i] <= orblimit then
Add(cand.points,v[i]);
Add(cand.ops,op);
fi;
od;
if Length(cand.points) > 1 then
Info( InfoGenSS, 2, "Adding ",
Length(cand.points)-1, " more vectors",
" to candidates.");
fi;
return cand;
fi;
else
AddSet(val,j);
fi;
else
HTAdd(ht,w,[j]);
fi;
od;
od;
minest := Minimum(e);
minpos := Position(e,minest);
Info( InfoGenSS,2,"Birthday #", round, ": no small enough estimate. ",
"MinEst=",minest," Coinc=",c[j] );
randels := randels * 2;
orblimit := orblimit * 4;
od;
TryNextMethod();
end );
InstallMethod( FindBasePointCandidates,
"for a matrix group over a FF, original try short orbit method",
[ IsGroup and IsMatrixGroup and IsFinite, IsRecord, IsInt, IsObject ], 10,
function( grp, opt, i, parentS )
local F, d, data, cand, res;
F := DefaultFieldOfMatrixGroup(grp);
d := DimensionOfMatrixGroup(grp);
data := opt.FindBasePointCandidatesData;
Info( InfoGenSS, 3, "Finding nice base points (TryShortOrbit)..." );
# Next possibly "TryShortOrbit":
cand := rec( points := [], used := 0 );
if IsBound(opt.TryShortOrbit) and opt.TryShortOrbit > 0 then
repeat
opt.TryShortOrbit := opt.TryShortOrbit - 1;
Info(InfoGenSS,1,"Looking for short orbit (",opt.TryShortOrbit,
")...");
res := GENSS_FindShortOrbit(grp,opt,parentS);
until res <> fail or opt.TryShortOrbit = 0;
if res <> fail then
if Size(F) > 2 then
cand.points := [res[1]];
cand.ops := [OnLines];
# Note that if we work non-projectively, then the same
# point will be taken next with OnRight automatically!
else
cand.points := [res[1]];
cand.ops := [OnRight];
fi;
return cand;
fi;
fi;
TryNextMethod();
end );
InstallMethod( FindBasePointCandidates,
"for a matrix group over a FF, traditional Murray/O'Brien",
[ IsGroup and IsMatrixGroup and IsFinite, IsRecord, IsInt, IsObject ],
function( grp, opt, i, parentS )
local F, d, bv, cand, w, v;
F := DefaultFieldOfMatrixGroup(grp);
d := DimensionOfMatrixGroup(grp);
Info( InfoGenSS, 3, "Finding nice base points (Murray/O'Brien)..." );
# Standard Murray/O'Brien heuristics:
if i = 0 and
((opt!.Projective = false and Size(F)^d > 300000) or
(opt!.Projective = true and Size(F)^(d-1) > 300000)) then
bv := GENSS_FindVectorsWithShortOrbit(grp,opt,parentS);
if Length(bv) < 3 then
bv := MutableCopyMat(One(grp));
fi;
bv := bv{[1..3]}; # just take 3 of them
else
bv := One(grp);
fi;
cand := rec( points := [], ops := [], used := 0 );
for v in bv do
w := ORB_NormalizeVector(ShallowCopy(v));
if Size(F) = 2 then
Add(cand.points,w);
Add(cand.ops,OnRight);
else
Add(cand.points,w);
Add(cand.ops,OnLines);
# Note that if we work non-projectively, then the same
# point will be taken next with OnRight automatically!
fi;
od;
return cand;
end );
InstallMethod( FindBasePointCandidates, "for a permutation group",
[ IsGroup and IsPermGroup, IsRecord, IsInt, IsObject ],
function( grp, opt, i, parentS )
local ops,points;
if i = 0 then
points := [1..Minimum(20,LargestMovedPoint(grp))];
else
points := [1..LargestMovedPoint(grp)];
fi;
ops := List([1..Length(points)],x->OnPoints);
return rec( points := points, ops := ops, used := 0 );
end );
GENSS_HACK := OnPoints;
InstallGlobalFunction( GENSS_OpFunctionMaker, function(op,index)
local name,s,f;
name := NAME_FUNC(op);
if name = "unknown" or name = "GENSS_HACK" then return fail; fi;
s := Concatenation( "GENSS_HACK := function(x,el) return ",
name, "(x,el[",String(index),"]); end;" );
f := InputTextString(s);
Read(f);
return GENSS_HACK;
end);
InstallMethod( FindBasePointCandidates, "for a direct product",
[ IsGroup, IsRecord, IsInt, IsObject ],
function( grp, opt, i, parentS )
local gens,l,cand,j,factgens,fac,cand2,k,op,S2,op2;
gens := GeneratorsOfGroup(grp);
if not(ForAll(gens,IsDirectProductElement)) or Length(gens) = 0 then
TryNextMethod();
fi;
l := Length(gens[1]);
cand := rec( points := [], ops := [] );
for j in [1..l] do
factgens := List(gens,x->x[j]);
fac := Group(factgens);
S2 := StabilizerChain(fac);
cand2 := BaseStabilizerChain(S2);
for k in [1..Length(cand2.ops)] do
op := cand2.ops[k];
op2 := GENSS_OpFunctionMaker(op,j);
if op2 <> fail then
Add(cand.points,cand2.points[k]);
Add(cand.ops,op2);
fi;
od;
od;
cand.used := 0;
return cand;
end );
InstallGlobalFunction( GENSS_NextBasePoint,
function( gens, cand, opt, S )
local NotFixedUnderAllGens,i,notfixed;
if IsBound(opt.FindBasePointCandidatesData) then
Unbind(opt.FindBasePointCandidatesData);
fi;
NotFixedUnderAllGens := function( gens, x, op )
if IsObjWithMemory(gens[1]) then
return ForAny( gens, g->op(x,g!.el) <> x );
else
return ForAny( gens, g->op(x,g) <> x );
fi;
end;
# S can be false or a stabilizer chain record
if S <> false and not(opt.StrictlyUseCandidates) then
# try points in previous orbit
for i in [2..Minimum(Length(S!.orb),opt.NumberPrevOrbitPoints)] do
if NotFixedUnderAllGens(gens,S!.orb[i],S!.orb!.op) then
Info(InfoGenSS,3,"Taking another point in same orbit.");
return rec( point := S!.orb[i], op := S!.orb!.op,
cand := cand );
fi;
od;
# Maybe we can take the last base point now acting non-projectively?
if opt.Projective = false and IsIdenticalObj(S!.orb!.op,OnLines) and
NotFixedUnderAllGens(gens,S!.orb[1],OnRight) then
Info(InfoGenSS,3,"Taking same point in non-projective action.");
return rec( point := S!.orb[1], op := OnRight, cand := cand );
fi;
fi;
# Now use the candidates we already have:
repeat
if cand.used >= Length(cand.points) then
if IsBound(cand.points2) and Length(cand.points2) > 0 then
cand := rec( points := cand.points2, ops := cand.ops2,
used := 0 );
else
cand := FindBasePointCandidates(Group(gens),opt,1,S);
opt.StrictlyUseCandidates := false;
opt.Reduced := true;
fi;
fi;
cand.used := cand.used + 1;
notfixed := NotFixedUnderAllGens(gens,
cand.points[cand.used],cand.ops[cand.used]);
if notfixed = false and opt.Reduced = true then # everything is fixed!
if IsBound(cand.points2) then
# Maybe our current stabilizer is too small, we still
# might want to use these points again:
Add(cand.points2,cand.points[cand.used]);
Add(cand.ops2,cand.ops[cand.used]);
fi;
fi;
until opt.Reduced = false or notfixed;
return rec( point := cand.points[cand.used], op := cand.ops[cand.used],
cand := cand );
end );
InstallGlobalFunction( GENSS_CreateStabChainRecord,
function( parentS, gens, size, nextpoint, nextop, cand, opt )
# parentS can be false or the parent in the stabiliser chain
local base, layer, stronggens, layergens, nr, orb, S, hashsize;
Info( InfoGenSS, 4, "Creating new stab chain record..." );
if parentS = false then
base := [];
layer := 1;
stronggens := ShallowCopy(gens);
layergens := [1..Length(gens)]; # indices in stronggens
else
base := parentS!.base;
layer := parentS!.layer + 1;
stronggens := parentS!.stronggens;
nr := Length(stronggens);
Append(stronggens,gens);
layergens := [nr+1..Length(stronggens)];
fi;
gens := ShallowCopy(gens);
# Note that we do ShallowCopy such that the original list and the
# one in the orbit record are different from each other.
if IsInt(size) then
hashsize := NextPrimeInt(Minimum(size,opt.InitialHashSize));
else
hashsize := opt.InitialHashSize;
fi;
orb := Orb( gens, nextpoint, nextop,
rec( treehashsize := hashsize, schreier := true,
log := opt.OrbitsWithLog,
report := opt.Report ) );
S := rec( stab := false, orb := orb, cand := cand, base := base,
opt := opt, layer := layer, parentS := parentS,
stronggens := stronggens, layergens := layergens,
size := size, randpool := [], IsOne := opt.IsOne );
if parentS = false then
S!.topS := S;
else
S!.topS := parentS!.topS;
fi;
if IsBound(opt.FindBasePointCandidatesData) then
S.randpool := opt.FindBasePointCandidatesData.randpool;
fi;
Add(base,nextpoint);
S.orb!.stabilizerchain := S;
Objectify( StabChainByOrbType, S );
return S;
end );
InstallGlobalFunction( GENSS_RandomElementFromAbove,
function( S, i )
# This function provides a random element for the i-th stabiliser
# in the chain "from above". These elements eventually come from
# the one product replacer object in the opt record for the whole
# group. They are sifted through i layers of the stabiliser chain
# infrastructure until they stabilise the first i points (i=0 simply
# means a random element in the whole group). If we find out underways
# that an orbit is too small (that is, a stabiliser was not complete),
# we fix the stabiliser chain as we go. Random elements that have
# been generated in a layer j < i previously and are not yet used
# as strong generators can be taken from a pool that was kept in
# the stabiliser chain. S must either be layer i of the stabiliser
# or (for i=0) it can be equal to the options record opt chain.
local SS, x, topS, j, o, p, po;
if i = 0 then
return Next(S.pr); # in this case we got the options record
fi;
if S!.layer <> i then
Error("i must be equal to the layer");
fi;
SS := S;
while true do # will be left by break
if Length(SS!.randpool) > 0 then
x := Remove(SS!.randpool,Length(SS!.randpool));
topS := SS!.topS;
break;
fi;
if SS!.parentS = false then # we have reached the top
if IsBound(SS!.opt.RandomElmFunc) then
x := SS!.opt.RandomElmFunc();
else
x := Next(SS!.opt.pr);
fi;
topS := SS; # remember top
break;
fi;
SS := SS!.parentS;
od;
# We now have a random element x in layer SS!.layer, that is,
# it is contained in the SS!.layer-1-th stabiliser, sift it down:
j := SS!.layer-1;
while j < i do
o := SS!.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
Info(InfoGenSS,3,"Random element from top found error in layer ",
SS!.layer);
AddGeneratorToStabilizerChain(topS,x);
# We add it at the top to add it in every orbit above except
# the first one!
return GENSS_RandomElementFromAbove(S,i);
fi;
# Now sift through Schreier tree:
while po > 1 do
x := x * SS!.orb!.gensi[o!.schreiergen[po]];
po := o!.schreierpos[po];
od;
SS := SS!.stab;
j := j + 1;
od;
# After this, we have successfully reached i-th stabilizer
return x;
end );
InstallGlobalFunction( GENSS_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 GENSS_StabilizerChainInner layer=",layer);
next := GENSS_NextBasePoint(gens,cand,opt,parentS);
cand := next.cand; # This could have changed
S := GENSS_CreateStabChainRecord(parentS,gens,size,
next.point,next.op,next.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
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);
# Are we done?
if size <> false and Length(S!.orb) = size then
S!.proof := true;
Info(InfoGenSS,4,"Leaving GENSS_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 := GENSS_StabilizerChainInner(stabgens,size/Length(S!.orb),
layer+1,cand,opt,S);
else
S!.stab := GENSS_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;
#Error("foobar");
while i < opt.ImmediateVerificationElements do
i := i + 1;
x := GENSS_RandomElementFromAbove(S,layer);
if 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 GENSS_StabilizerChainInner layer=",layer);
return S;
end );
InstallGlobalFunction( GENSS_DeriveCandidatesFromStabChain,
function( S )
local cand;
cand := rec( points := [], ops := [], used := 0 );
while S <> false do
Add( cand.points, S!.orb[1] );
Add( cand.ops, S!.orb!.op );
S := S!.stab;
od;
return cand;
end );
InstallGlobalFunction( GENSS_TrivialOp,
function( p, x )
return p;
end );
InstallMethod( StabilizerChain, "for a group object", [ IsGroup ],
function( grp )
return StabilizerChain( grp, rec() );
end );
InstallMethod(VerifyStabilizerChainMC,
"for a stabilizer chain and a positive integer",
[ IsStabilizerChainByOrb, IsInt ],
function( S, nrels )
local i,x;
Info(InfoGenSS,2,"Doing randomized verification...");
i := 0;
while i < nrels do
i := i + 1;
if IsBound(S!.opt.RandomElmFunc) then
x := S!.opt.RandomElmFunc();
else
x := Next(S!.opt.pr);
fi;
if AddGeneratorToStabilizerChain(S,x) then
Info( InfoGenSS, 2, "Verification found error ... ",
"new size ", Size(S) );
i := 0;
fi;
od;
end );
InstallMethod( StabilizerChain, "for a group object and a record",
[ IsGroup, IsRecord ],
function( grp, opt )
# Computes a stabilizer chain for the group grp
# Possible options:
# random: compatibility to StabChain, works as there
# ErrorBound: rational giving the error bound, probability of an
# error will be proven to be lower than this
# if given, this sets VerifyElements and
# DeterministicVerification
# VerifyElements: number of random elements for verification
# DeterministicVerification: flag, whether to do or not
# (not yet implemented)
# Base: specify a known base, this can either be the
# StabilizerChain of a known supergroup or a list of
# points, if "BaseOps" is not given, it is either
# set to OnPoints if Projective is not set or false,
# or is set to OnLines if Projective is set to true,
# if Base is set
# BaseOps: a list of operation functions for the base points
# Projective: if set to true indicates that the matrices given
# as generators are to be thought as projective
# elements, that is, the whole StabilizerChain will
# be one for a projective group!
# StrictlyUseCandidates: if set to true exactly the points in
# opt.cand will be used as base points and no other
# points from previous orbits are used first
# this is set when Base was specified
# Reduced: if set to true no orbits of length 1 will occur
# (except for the trivial group), is true by default
# Cand: initializer for cand, which are base points
# candidates
# must be a record with components "points", "ops",
# "used", the latter indicates the largest index
# that has already been used. "used" is automatically
# set to 0 if not set.
# TryShortOrbit: Number of tries for the short orbit finding alg.
# StabGenScramble,
# StabGenScrambleFactor,
# StabGenAddSlots,
# StabGenMaxDepth: parameters for product replacer for generating
# random elements
# VeryShortOrbLimit: when looking for short orbits try a random
# vector and enumerate its orbit until this limit
#
# ... to be continued
#
local S,cand,i,pr,prob,x,gens,SS;
# First a few preparations, then we delegate to GENSS_StabilizerChainInner:
# Add some default options:
if (HasSize(grp) and not(IsBound(opt.Projective) and opt.Projective))
or IsBound(opt.Size) then
if not(IsBound(opt.ImmediateVerificationElements)) then
opt.ImmediateVerificationElements := 0;
fi;
fi;
if IsBound(opt.Projective) and opt.Projective then
opt.IsOne := GENSS_IsOneProjective;
elif not(IsBound(opt.IsOne)) then
opt.IsOne := IsOne;
fi;
# Now opt.IsOne is set to a function to check whether or not a group
# element is equal to the identity.
GENSS_CopyDefaultOptions(GENSS,opt);
# Check for the identity group:
gens := GeneratorsOfGroup(grp);
if Length(gens) = 0 or
ForAll(gens,opt.IsOne) then
# Set up a trivial stabilizer chain record:
S := GENSS_CreateStabChainRecord(false,gens,1,1,GENSS_TrivialOp,
rec( points := [], ops := [],
used := 0),opt);
Enumerate(S!.orb);
S!.orb!.gensi := List(S!.orb!.gens,x->x^-1);
S!.proof := true;
S!.trivialgroup := true;
if (not(IsBound(opt.Projective)) or
opt.Projective = false) and
IsIdenticalObj(opt.IsOne,IsOne) then
SetStoredStabilizerChain(grp,S);
fi;
return S;
fi;
# Setup a random element generator for the whole group and store
# it in the opt record:
opt.pr := ProductReplacer( gens,
rec( scramble := opt.StabGenScramble,
scramblefactor := opt.StabGenScrambleFactor,
addslots := opt.StabGenAddSlots,
maxdepth := opt.StabGenMaxDepth ));
# Old style error probability for compatibility:
if IsBound(opt.random) then
if opt.random = 0 then
opt.VerifyElements := 0;
elif opt.random = 1000 then
opt.DeterministicVerification := true;
Info(InfoGenSS,1,"Warning: Deterministic verification not yet ",
"implemented!");
else
prob := 1/2;
opt.VerifyElements := 1;
while prob * (1000-opt.random) >= 1 do
prob := prob / 2;
opt.VerifyElements := opt.VerifyElements + 1;
od;
fi;
fi;
# The new style error bounds:
if IsBound(opt.ErrorBound) then
prob := 1/2;
opt.VerifyElements := 1;
while prob >= opt.ErrorBound do
prob := prob / 2;
opt.VerifyElements := opt.VerifyElements + 1;
od;
fi;
# Find base point candidates:
if IsBound(opt.Base) then
if IsStabilizerChain(opt.Base) then
cand := GENSS_DeriveCandidatesFromStabChain(opt.Base);
else # directly take the base points:
cand := rec( points := opt.Base, used := 0, points2 := [],
ops2 := [] );
if not(IsBound(opt.BaseOps)) then
# Let's guess the ops:
if IsBound(opt.Projective) and opt.Projective = true then
cand.ops := ListWithIdenticalEntries(Length(opt.Base),
OnLines);
else
cand.ops := ListWithIdenticalEntries(Length(opt.Base),
OnPoints);
fi;
else
cand.ops := opt.BaseOps;
fi;
fi;
opt.StrictlyUseCandidates := true;
elif IsBound(opt.Cand) then
cand := opt.Cand;
if not(IsBound(cand.used)) then
cand.used := 0;
fi;
cand.points2 := [];
cand.ops2 := [];
else
# Otherwise try different things later using generic methods:
cand := rec( points := [], ops := [], used := 0 );
fi;
if not(IsBound(opt.StrictlyUseCandidates)) then
opt.StrictlyUseCandidates := false;
fi;
if HasSize(grp) and not(opt.Projective) then
S := GENSS_StabilizerChainInner(GeneratorsOfGroup(grp), Size(grp), 1,
cand, opt, false);
elif IsBound(opt.Size) then
S := GENSS_StabilizerChainInner(GeneratorsOfGroup(grp), opt.Size, 1,
cand, opt, false);
else
S := GENSS_StabilizerChainInner(GeneratorsOfGroup(grp), false, 1,
cand, opt, false);
fi;
if IsString(S) then return S; fi;
# Do we already have a proof?
if S!.proof then
if not(IsBound(opt.Projective)) or opt.Projective = false then
SetStoredStabilizerChain(grp,S);
fi;
return S;
fi;
if opt.VerifyElements = 0 then return S; fi;
Info(InfoGenSS,2,"Current size found: ",Size(S));
# Now a possible verification phase:
if S!.size <> false then # we knew the size in advance
Info(InfoGenSS,2,"Doing verification via known size...");
while Size(S) < S!.size do
Info(InfoGenSS,2,"Known size not reached, throwing in a random ",
"element...");
if IsBound(opt.RandomElmFunc) then
x := opt.RandomElmFunc();
else
x := Next(opt.pr);
fi;
if AddGeneratorToStabilizerChain(S,x) then
Info( InfoGenSS, 2, "Increased size to ",Size(S) );
fi;
od;
S!.proof := true;
if (not(IsBound(opt.Projective)) or opt.Projective = false) and
IsIdenticalObj(opt.IsOne,IsOne) then
SetStoredStabilizerChain(grp,S);
fi;
else
# Do some verification here:
VerifyStabilizerChainMC(S,opt.VerifyElements);
fi;
# Now clean up the random element pools to save memory:
SS := S;
while SS <> false do
if IsBound(SS!.randpool) and Length(SS!.randpool) > 0 then
SS!.randpool := [];
fi;
SS := SS!.stab;
od;
return S;
end );
InstallMethod( AddGeneratorToStabilizerChain,
"for a stabilizer chain and a new generator",
[ IsStabilizerChain and IsStabilizerChainByOrb, IsObject ],
function( S, el )
# 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);
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 stabilizer chain ",
"in layer ", SS!.layer, "..." );
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 := GENSS_StabilizerChainInner([r.rem],false,
SS!.layer+1,SS!.cand, 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 );
InstallMethod( SiftGroupElement, "for a stabilizer chain and a group element",
[ IsStabilizerChain and IsStabilizerChainByOrb, IsObject ],
function( S, x )
local o,p,po,preS,r,isone;
isone := S!.IsOne;
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
return rec( isone := false, rem := x, S := S, preS := preS );
fi;
# Now sift through Schreier tree:
while po > 1 do
x := x * S!.orb!.gensi[o!.schreiergen[po]];
po := o!.schreierpos[po];
od;
preS := S;
S := S!.stab;
od;
r := rec( rem := x, S := false, preS := preS, isone := isone(x) );
return r;
end );
InstallMethod( SiftGroupElementSLP,
"for a stabilizer chain and a group element",
[ IsStabilizerChain and IsStabilizerChainByOrb, IsObject ],
function( S, x )
local preS, nrstrong, slp, o, p, po, r, isone;
preS := false;
isone := S!.IsOne;
nrstrong := Length(S!.stronggens);
slp := []; # will be reversed in the end
while S <> false do
o := S!.orb;
p := o!.op(o[1],x);
po := Position(o,p);
if po = fail then # not in current stabilizer
return rec( isone := false, rem := x, S := S, preS := preS,
slp := fail );
fi;
# Now sift through Schreier tree:
while po > 1 do
x := x * S!.orb!.gensi[o!.schreiergen[po]];
Add(slp,1);
Add(slp,S!.layergens[o!.schreiergen[po]]);
po := o!.schreierpos[po];
od;
preS := S;
S := S!.stab;
od;
r := rec( rem := x, S := false, preS := preS, isone := isone(x) );
if r.isone then
if Length(slp) = 0 then # element was the identity!
r.slp := StraightLineProgramNC([[1,0]],nrstrong);
else
r.slp := StraightLineProgramNC(
[slp{[Length(slp),Length(slp)-1..1]}],nrstrong);
fi;
else
r.slp := fail;
fi;
return r;
end );
InstallMethod( StrongGenerators, "for a stabilizer chain",
[ IsStabilizerChain and IsStabilizerChainByOrb ],
function( S )
return S!.stronggens;
end );
InstallMethod( NrStrongGenerators, "for a stabilizer chain",
[ IsStabilizerChain and IsStabilizerChainByOrb ],
function( S )
return Length(S!.stronggens);
end );
InstallMethod( ForgetMemory, "for a stabilizer chain",
[ IsStabilizerChain and IsStabilizerChainByOrb ],
function( S )
ForgetMemory(S!.stronggens);
while S <> false do
ForgetMemory(S!.orb);
ForgetMemory(S!.orb!.gensi);
S := S!.stab;
od;
end );
InstallMethod( AddNormalizingGenToLayer,
"for a stabilizer chain, a group element, and a prime number",
[ IsStabilizerChain and IsStabilizerChainByOrb, IsObject, IsPosInt ],
function( S, x, p )
# This is a very specialised function which is not officially documented.
# It is used for the matrix PCGS code where a stabiliser chain is
# built up using a known base in a generator by generator fashion
# where the generators are the generators of the PCGS. Thus a lot
# is known about them.
# This function assumes that x is a new generator normalising the
# group generated by the previous generators in the layer described
# by S. So in particular, S is not necessarily the top of the
# stabiliser chain but the layer where actually something happens.
# Namely, since it is known that the new generator generates
# a group which is p times as big as the previous one (p a prime)
# it is known that the orbit in the given layer will become
# p times as big and we can enumerate it relatively cheaply.
local lmp,o,oldnrgens;
o := S!.orb;
oldnrgens := Length(o!.gens);
Add(o!.gens,x);
Add(o!.gensi,x^-1);
if IsPermOnIntOrbitRep(o) then
lmp := LargestMovedPoint(o!.gens);
if lmp > Length(o!.tab) then
Append(o!.tab,ListWithIdenticalEntries(lmp-Length(o!.tab),0));
fi;
fi;
ResetFilterObj(o,IsClosed);
o!.pos := 1;
o!.genstoapply := [oldnrgens+1..Length(o!.gens)];
repeat
Enumerate(o,S!.opt.OrbitLengthLimit);
if not(IsClosed(o)) 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(S!.orb);
o!.genstoapply := [1..Length(o!.gens)];
# Now fix up the stabilizer chain record:
if not x in S!.stronggens then
Add(S!.stronggens,x);
Add(S!.layergens,Length(S!.stronggens));
fi;
Unbind(S!.size); # whatever we knew before, we know no longer
end );
InstallMethod( GENSS_CreateSchreierGenerator,
"for a stabilizer chain, a position in the first orbit and a gen number",
[ IsStabilizerChain and IsStabilizerChainByOrb, IsPosInt, IsPosInt ],
function( S, i, j )
local o,p,w1,w2;
o := S!.orb;
p := Position(o,o!.op(o[i],o!.gens[j]));
if o!.schreierpos[p] = i and o!.schreiergen[p] = j then
return fail;
fi;
w1 := TraceSchreierTreeForward(o,i);
w2 := TraceSchreierTreeBack(o,p);
return [w1,j,w2];
end );
InstallGlobalFunction( GENSS_Prod,
function( gens, word )
if Length(word) = 0 then
return gens[1]^0;
else
return Product(gens{word});
fi;
end );
InstallGlobalFunction( VerifyStabilizerChainTC,
function( S ) Error("Currently not functional!"); return fail; end );
## function( S )
## local Grels,Hrels,Prels,MakeSchreierGens,ct,f,gens,gensi,i,j,k,l,li,max,
## newpres,nrgens,nrschr,o,ords,pres,sb,sgs,slp,subgens,v,w,x,
## cosetnrlimitfactor;
##
## allstrong := function(S)
## st := Set(S!.layergens);
## while S!.stab <> false do
## S := S!.stab;
## UniteSet(st,S!.layergens);
## od;
## return st;
## end;
## if S!.stab <> false then
## pres := VerifyStabilizerChainTC(S!.stab);
## if IsList(pres) then return pres; fi;
## strongbelow := allstrong(S!.stab);
## else
## pres := StraightLineProgram([[]],0);
## strongbelow := [];
## fi;
## strong := allstrong(S);
## Info(InfoGenSS,1,"Verifying stabilizer chain in layer ",S!.layer);
##
## # First create a few Schreier generators:
## sgs := [];
## i := 1;
## j := 1;
## o := S!.orb;
## nrgens := Length(o!.gens);
## MakeSchreierGens := function(n)
## local sg;
## Info(InfoGenSS,3,"Creating ",n," Schreier generators...");
## while Length(sgs) < n and
## i <= Length(o) do
## sg := GENSS_CreateSchreierGenerator(S,i,j);
## j := j + 1;
## if j > nrgens then
## j := 1;
## i := i + 1;
## fi;
## if sg <> fail then
## Add(sgs,sg);
## fi;
## od;
## end;
##
## nrschr := S!.opt.NumberSchreierGens;
## MakeSchreierGens(nrschr);
## f := FreeGroup(Length(strong));
## gens := GeneratorsOfGroup(f);
## gensi := List(gens,x->x^-1);
## subgens := gens{List(strongbelow,i->Position(strong,i))};
## Hrels := ResultOfStraightLineProgram(pres,subgens);
## if S!.opt.Projective then
## ords := List([1..nrgens],i->ProjectiveOrder(o!.gens[i]));
## else
## ords := List([1..nrgens],i->Order(o!.gens[i]));
## fi;
## Prels := List([1..nrgens],i->gens[i+sb]^ords[i]);
## Grels := [];
## cosetnrlimitfactor := 4;
## while true do # will be left by return eventually
## for k in [Length(Grels)+1..Length(sgs)] do
## Grels[k] := GENSS_Prod(gens,sgs[k][1]+sb) * gens[sgs[k][2]+sb] *
## GENSS_Prod(gensi,sgs[k][3]+sb);
## x := GENSS_Prod(o!.gens,sgs[k][1]) * o!.gens[sgs[k][2]] *
## GENSS_Prod(o!.gensi,sgs[k][3]);
## if S!.stab <> false then
## slp := SiftGroupElementSLP(S!.stab,x);
## if not(slp.isone) then
## return [fail,S!.layer];
## fi;
## Grels[k] := Grels[k] / ResultOfStraightLineProgram(slp.slp,
## subgens);
## sgs[k][4] := slp.slp;
## else
## if not(S!.IsOne(x)) then
## return [fail,S!.layer];
## fi;
## sgs[k][4] := false;
## fi;
## od;
## Info(InfoGenSS,2,"Doing coset enumeration with limit ",
## cosetnrlimitfactor*Length(o));
## ct := CosetTableFromGensAndRels(gens,Concatenation(Prels,Hrels,Grels),
## subgens:max := cosetnrlimitfactor*Length(o),silent);
## if ct = fail then # did not close!
## cosetnrlimitfactor := QuoInt(cosetnrlimitfactor*3,2);
## Info(InfoGenSS,2,"Coset enumeration did not finish!");
## #Error(1);
## if nrschr > Length(sgs) # or
## # nrschr > S!.opt.MaxNumberSchreierGens
## then # we are done!
## # Something is wrong!
## return [fail, S!.layer];
## fi;
## else
## Info(InfoGenSS,2,"Coset enumeration found ",Length(ct[1]),
## " cosets.");
## #Error(2);
## if Length(ct[1]) = Length(o) then
## # Verification is OK, now build a presentation:
## l := GeneratorsWithMemory(
## ListWithIdenticalEntries(S!.nrstrong,()));
## li := List(l,x->x^-1);
## newpres := ResultOfStraightLineProgram(pres,
## l{[1..S!.strongbelow]});
## for k in [1..nrgens] do
## Add(newpres,l[k+sb]^ords[k]);
## od;
## for k in [1..Length(sgs)] do
## if sgs[k][4] <> false then
## Add(newpres,
## GENSS_Prod(l,sgs[k][1]+sb)*l[sgs[k][2]+sb]*
## GENSS_Prod(li,sgs[k][3]+sb)*
## ResultOfStraightLineProgram(sgs[k][4],l)^-1);
## else
## Add(newpres,GENSS_Prod(l,sgs[k][1]+sb)*l[sgs[k][2]+sb]*
## GENSS_Prod(li,sgs[k][3]+sb));
## fi;
## od;
## Info(InfoGenSS,2,"Found presentation for layer ",S!.layer,
## " using ",Length(newpres)," relators.");
## return SLPOfElms(newpres);
## fi;
## fi;
## nrschr := QuoInt(nrschr*4,2);
## MakeSchreierGens(nrschr);
## od;
## end);
GENSS_CosetTableFromGensAndRelsInit :=
function(fgens,grels,fsgens,limit)
local next, prev, # next and previous coset on lists
firstFree, lastFree, # first and last free coset
firstDef, lastDef, # first and last defined coset
table, # columns in the table for gens
rels, # representatives of the relators
relsGen, # relators sorted by start generator
subgroup, # rows for the subgroup gens
i, gen, inv, # loop variables for generator
g, # loop variable for generator col
rel, # loop variables for relation
p, p1, p2, # generator position numbers
app, # arguments list for 'MakeConsequences'
j, # integer variable
length, length2, # length of relator (times 2)
cols,
nums,
l,
nrdef, # number of defined cosets
nrmax, # maximal value of the above
nrdel, # number of deleted cosets
nrinf; # number for next information message
# give some information
Info( InfoFpGroup, 3, " defined deleted alive maximal");
nrdef := 1;
nrmax := 1;
nrdel := 0;
nrinf := 1000;
# define one coset (1)
firstDef := 1; lastDef := 1;
firstFree := 2; lastFree := limit;
# make the lists that link together all the cosets
next := [ 2 .. limit + 1 ]; next[1] := 0; next[limit] := 0;
prev := [ 0 .. limit - 1 ]; prev[2] := 0;
# compute the representatives for the relators
rels := RelatorRepresentatives( grels );
# make the columns for the generators
table := [];
for gen in fgens do
g := ListWithIdenticalEntries( limit, 0 );
Add( table, g );
if not ( gen^2 in rels or gen^-2 in rels ) then
g := ListWithIdenticalEntries( limit, 0 );
fi;
Add( table, g );
od;
# make the rows for the relators and distribute over relsGen
relsGen := RelsSortedByStartGen( fgens, rels, table, true );
# make the rows for the subgroup generators
subgroup := [];
for rel in fsgens do
#T this code should use ExtRepOfObj -- its faster
# cope with SLP elms
if IsStraightLineProgElm(rel) then
rel:=EvalStraightLineProgElm(rel);
fi;
length := Length( rel );
length2 := 2 * length;
nums := [ ]; nums[length2] := 0;
cols := [ ]; cols[length2] := 0;
# compute the lists.
i := 0; j := 0;
while i < length do
i := i + 1; j := j + 2;
gen := Subword( rel, i, i );
p := Position( fgens, gen );
if p = fail then
p := Position( fgens, gen^-1 );
p1 := 2 * p;
p2 := 2 * p - 1;
else
p1 := 2 * p - 1;
p2 := 2 * p;
fi;
nums[j] := p1; cols[j] := table[p1];
nums[j-1] := p2; cols[j-1] := table[p2];
od;
Add( subgroup, [ nums, cols ] );
od;
# make the structure that is passed to 'MakeConsequences'
app := [ table, next, prev, relsGen, subgroup,
firstFree, lastFree, firstDef, lastDef, 0, 0 ];
# we do not want minimal gaps to be marked in the coset table
app[12] := 0;
return rec( nrdef := nrdef, nrmax := nrmax, nrdel := nrdel, nrinf := nrinf,
app := app, closed := false, table := table, limit := limit,
fgens := fgens );
end;
GENSS_DoToddCoxeter := function( r )
local app,fgens,firstDef,firstFree,gen,i,inv,lastDef,lastFree,next,nrdef,
nrdel,nrinf,nrmax,prev,relsGen,subgroup,table;
fgens := r.fgens;
nrdef := r.nrdef;
nrmax := r.nrmax;
nrdel := r.nrdel;
nrinf := r.nrinf;
app := r.app;
table := app[1];
next := app[2];
prev := app[3];
relsGen := app[4];
subgroup := app[5];
firstFree := app[6];
lastFree := app[7];
firstDef := app[8];
lastDef := app[9];
# run over all the cosets:
while firstDef <> 0 do
# run through all the rows and look for undefined entries
for i in [ 1 .. Length( table ) ] do
gen := table[i];
if gen[firstDef] <= 0 then
inv := table[i + 2*(i mod 2) - 1];
# if necessary expand the table
if firstFree = 0 then
app[6] := firstFree;
app[7] := lastFree;
app[8] := firstDef;
app[9] := lastDef;
r.nrdef := nrdef;
r.nrmax := nrmax;
r.nrdel := nrdel;
r.nrinf := nrinf;
return fail;
fi;
# update the debugging information
nrdef := nrdef + 1;
if nrmax <= firstFree then
nrmax := firstFree;
fi;
# define a new coset
gen[firstDef] := firstFree;
inv[firstFree] := firstDef;
next[lastDef] := firstFree;
prev[firstFree] := lastDef;
lastDef := firstFree;
firstFree := next[firstFree];
next[lastDef] := 0;
# set up the deduction queue and run over it until it's empty
app[6] := firstFree;
app[7] := lastFree;
app[8] := firstDef;
app[9] := lastDef;
app[10] := i;
app[11] := firstDef;
--> --------------------
--> maximum size reached
--> --------------------
[ Dauer der Verarbeitung: 0.46 Sekunden
(vorverarbeitet)
]
|
2026-04-02
|
|
|
|
|
