Quelle search.gi
Sprache: unbekannt
|
|
#############################################################################
##
## orb package
## search.gi
## Juergen Mueller
## Max Neunhoeffer
## Felix Noeske
##
## Copyright 2005-2008 by the authors.
## This file is free software, see license information at the end.
##
## Implementation stuff for searching in groups.
##
#############################################################################
# Random vectors and matrices:
# Methods for Randomize for gf2 and 8bit vectors have been taken out
# they will be in GAP >= 4.5 in the library. In the meantime use
# cmats.
InstallGlobalFunction( MakeRandomVectors,
function( arg )
local i,l,number,randomsource,sample;
if Length(arg) <= 1 or Length(arg) > 3 then
ErrorNoReturn("Usage: MakeRandomVectors( sample, number [,randomsource] )");
fi;
sample := arg[1];
number := arg[2];
if Length(arg) = 3 then
randomsource := arg[3];
else
randomsource := GlobalRandomSource;
fi;
l := [];
for i in [number,number-1..1] do
sample := ShallowCopy(sample);
Randomize(sample,randomsource);
l[i] := sample;
od;
return l;
end );
InstallGlobalFunction( MakeRandomLines,
function( arg )
local i,l,number,pos,randomsource,sample;
if Length(arg) <= 1 or Length(arg) > 3 then
ErrorNoReturn("Usage: MakeRandomLines( sample, number [,randomsource] )");
fi;
sample := arg[1];
number := arg[2];
if Length(arg) = 3 then
randomsource := arg[3];
else
randomsource := GlobalRandomSource;
fi;
l := [];
for i in [number,number-1..1] do
sample := ShallowCopy(sample);
repeat
Randomize(sample,randomsource);
pos := PositionNonZero(sample);
until pos <= Length(sample);
MultVector(sample,sample[pos]^-1);
l[i] := sample;
od;
return l;
end );
# Product replacers:
BindGlobal( "ProductReplacersType",
NewType( ProductReplacersFamily, IsProductReplacer and IsMutable) );
InstallMethod( ProductReplacer, "for a list of group generators",
[IsList], function( gens )
return ProductReplacer( gens, rec( ) );
end );
InstallMethod( ProductReplacer,
"for a group object and an options record",
[IsGroup, IsRecord],
function(g,r) return ProductReplacer(GeneratorsOfGroup(g),r); end );
InstallMethod( ProductReplacer,
"for a group object",
[IsGroup],
function(g) return ProductReplacer(GeneratorsOfGroup(g),rec()); end );
InstallMethod( ProductReplacer,
"for a list of group generators and an options record",
[IsList, IsRecord],
function( gens, opt )
# First add some default options if not set:
local pr;
if Length(gens) = 0 then
Error("ProductReplacer: need at least one generator");
return fail;
fi;
pr := ShallowCopy(opt);
if CanEasilySortElementsFamily(FamilyObj(gens[1])) then
pr.gens := Set(gens);
else
pr.gens := ShallowCopy(gens);
fi;
pr.nrgens := Length(pr.gens);
if not IsBound(pr.addslots) then
pr.addslots := 5;
fi;
if not IsBound(pr.minslots) then
pr.minslots := 0;
fi;
if Length(gens) = 1 and pr.addslots = 0 then
# Make sure that the case of one generator is OK!
pr.minslots := Maximum(2,pr.minslots);
fi;
if not IsBound(pr.randomsource) then
pr.randomsource := GlobalRandomSource;
fi;
if not IsBound(pr.scramble) then
pr.scramble := 30;
fi;
if not IsBound(pr.scramblefactor) then
pr.scramblefactor := 4;
fi;
pr.scramble := Maximum(pr.scramble,pr.scramblefactor*pr.nrgens);
if not IsBound(pr.retirecaptain) then
pr.retirecaptain := 2 * pr.scramble;
fi;
if not IsBound(pr.maxdepth) then
pr.maxdepth := infinity;
fi;
if not IsBound(pr.noaccu) then
pr.noaccu := false;
fi;
if not IsBound(pr.accus) or pr.accus < 0 then
pr.accus := 5;
elif pr.accus = 0 then
pr.noaccu := true;
fi;
if not IsBound(pr.accelerator) then
pr.accelerator := true;
fi;
if not IsBound(pr.normalin) then
pr.normalin := false;
else
if pr.normalin <> false and
not(IsProductReplacer(pr.normalin)) then
Error("normalin option must be a product replacer");
fi;
fi;
pr.slots := pr.nrgens + pr.addslots;
if pr.slots < pr.minslots then
pr.slots := pr.minslots;
fi;
pr.initialized := false;
pr.steps := 0;
if pr.noaccu = false then
pr.accu := ListWithIdenticalEntries(pr.accus,pr.gens[1]^0);
pr.accupos := 0;
else
pr.accu := fail;
pr.accupos := 0;
fi;
Objectify(ProductReplacersType,pr);
pr!.team := ShallowCopy(pr!.gens);
while Length(pr!.team) < pr!.slots do
Add(pr!.team,pr!.gens[Random(pr!.randomsource,1,pr!.nrgens)]);
od;
return pr;
end );
InstallMethod( Reset, "for a product replacer", [IsProductReplacer],
function(pr)
if not(pr!.initialized) then return; fi;
pr!.team := ShallowCopy(pr!.resetteam);
pr!.accu := ShallowCopy(pr!.resetaccu);
pr!.accupos := pr!.resetaccupos;
pr!.steps := pr!.resetsteps;
if pr!.normalin <> false then
Reset(pr!.normalin);
fi;
end );
InstallMethod( Next, "for a product replacer", [IsProductReplacer],
function(pr)
local OneStep,i;
OneStep := function(pr)
local a, b, c, l, x, result, i;
pr!.steps := pr!.steps + 1;
c := Random(pr!.randomsource,1,2);
if pr!.accelerator and pr!.steps <= pr!.retirecaptain then
a := Random(pr!.randomsource,2,pr!.slots);
b := Random(pr!.randomsource,2,pr!.slots);
if c = 1 then
if pr!.normalin <> false then
pr!.team[1] := pr!.team[1]*pr!.team[b]^Next(pr!.normalin);
pr!.team[a] := pr!.team[a]*pr!.team[1];
else
pr!.team[1] := pr!.team[1]*pr!.team[b];
pr!.team[a] := pr!.team[a]*pr!.team[1];
fi;
else
if pr!.normalin <> false then
pr!.team[1] := pr!.team[b]^Next(pr!.normalin)*pr!.team[1];
pr!.team[a] := pr!.team[1] * pr!.team[a];
else
pr!.team[1] := pr!.team[b] * pr!.team[1];
pr!.team[a] := pr!.team[1] * pr!.team[a];
fi;
fi;
result := pr!.team[a];
else
a := Random(pr!.randomsource,1,pr!.slots);
b := Random(pr!.randomsource,1,pr!.slots-1);
if b >= a then b := b + 1; fi;
if c = 1 then
if pr!.normalin <> false then
pr!.team[a]:=pr!.team[a]*pr!.team[b]^Next(pr!.normalin);
else
pr!.team[a]:=pr!.team[a]*pr!.team[b];
fi;
else
if pr!.normalin <> false then
pr!.team[a]:=pr!.team[b]^Next(pr!.normalin)*pr!.team[a];
else
pr!.team[a]:=pr!.team[b]*pr!.team[a];
fi;
fi;
result := pr!.team[a];
fi;
if pr!.noaccu or
pr!.steps <= Maximum(pr!.nrgens * pr!.scramblefactor, pr!.scramble)
- pr!.accus then
return result;
else
pr!.accupos := pr!.accupos + 1;
if pr!.accupos > pr!.accus then pr!.accupos := 1; fi;
pr!.accu[pr!.accupos] := pr!.accu[pr!.accupos] * result; # Rattle
return pr!.accu[pr!.accupos];
fi;
end;
if pr!.steps > pr!.maxdepth then Reset(pr); fi;
if not(pr!.initialized) then
for i in [1..Maximum(pr!.nrgens * pr!.scramblefactor, pr!.scramble)] do
OneStep(pr);
od;
pr!.resetteam := ShallowCopy(pr!.team); # Remember for reset
pr!.resetaccu := ShallowCopy(pr!.accu);
pr!.resetaccupos := pr!.accupos;
pr!.resetsteps := pr!.steps;
pr!.initialized := true;
fi;
return OneStep(pr);
end );
InstallMethod( AddGeneratorToProductReplacer,
"for a product replacer and a new generator",
[ IsProductReplacer, IsObject ],
function(pr,el)
local i;
Add(pr!.gens,el);
pr!.nrgens := pr!.nrgens + 1;
for i in [1..Length(pr!.team)] do
pr!.team[i] := pr!.team[i] * el^Random(0,7);
od;
Add(pr!.team,el);
if pr!.initialized then
for i in [1..Length(pr!.resetteam)] do
pr!.resetteam[i] := pr!.resetteam[i] * el^Random(0,7);
od;
Add(pr!.resetteam,el);
fi;
pr!.slots := pr!.slots + 1;
end );
InstallMethod( ViewObj, "for a product replacer", [IsProductReplacer],
function(pr)
Print("<product replacer gens=",pr!.nrgens," slots=",pr!.slots,
" scramble=",Maximum(pr!.nrgens*pr!.scramblefactor,pr!.scramble),
" maxdepth=",pr!.maxdepth,"\n steps=",pr!.steps);
if pr!.noaccu then Print(" (without accu");
else Print(" (rattle accus=",pr!.accus); fi;
if pr!.accelerator then Print(", with accelerator)");
else Print(")"); fi;
if pr!.normalin <> false then
Print("\n normalin=");
ViewObj(pr!.normalin);
fi;
Print(">");
end );
# Finding things in groups, random searchers:
BindGlobal( "RandomSearchersType",
NewType( RandomSearchersFamily,
IsRandomSearcher and IsMutable ) );
InstallMethod( RandomSearcher,
"for a list of generators and a test function",
[IsList, IsFunction],
function( gens, testfunc )
return RandomSearcher( gens, testfunc, rec() );
end );
InstallMethod( RandomSearcher,
"for a list of generators, a test function, and an options record",
[IsList, IsFunction, IsRecord],
function( gens, testfunc, opt )
# Set some default options:
local r;
r := ShallowCopy(opt);
if not IsBound(r.exceptions) then
r.exceptions := [];
fi;
if not IsBound(r.maxdepth) then
r.maxdepth := 4 * Maximum(100,10*Length(gens));
fi;
if not IsBound(r.addslots) then
r.addslots := 10;
fi;
if not IsBound(r.limit) then
r.limit := infinity;
fi;
r.stops := 0;
r.tries := 0;
r.testfunc := testfunc;
if IsBound( r.scramble ) then
r.productreplacer := ProductReplacer( gens,
rec( maxdepth := r.maxdepth,
scramble := 100,
scramblefactor := 10,
addslots := r.addslots ) );
# this is the standard PseudoRandom behaviour
else # the standard case using all random elements:
r.productreplacer := ProductReplacer( gens,
rec( maxdepth := r.maxdepth,
scramble := 0,
scramblefactor := 0,
addslots := r.addslots ) );
# this uses all random elements from the first on
fi;
Objectify( RandomSearchersType,r );
return r;
end );
InstallMethod( ViewObj, "for a random searcher", [IsRandomSearcher],
function( rs )
Print("<random searcher");
if IsBound(rs!.target) then
Print(" for ",rs!.target);
fi;
Print(", tries: ",rs!.tries," stops: ",rs!.stops,", exceptions: ",
Length(rs!.exceptions),">");
end );
InstallMethod( Search, "for a random searcher", [IsRandomSearcher],
function( rs )
local x,y;
repeat
repeat
x := Next( rs!.productreplacer );
rs!.tries := rs!.tries + 1;
if rs!.tries mod 10000 = 0 then
Info(InfoOrb,2,"Tried ",rs!.tries," elements.");
fi;
if IsObjWithMemory(x) then
y := x!.el;
else
y := x;
fi;
if rs!.tries > rs!.limit then
return fail;
fi;
until not( y in rs!.exceptions );
until rs!.testfunc(y);
rs!.stops := rs!.stops + 1;
Add(rs!.exceptions,y);
return x;
end );
###################################################
# Involution centralisers and the dihedral trick: #
###################################################
InstallGlobalFunction( FindInvolution,
function(pr)
# g a matrix group
local i,o,x;
for i in [1..1000] do
x := Next(pr);
o := Order(StripMemory(x));
if o mod 2 = 0 then
return x^(o/2);
fi;
od;
return fail;
end );
InstallGlobalFunction( FindCentralisingElementOfInvolution,
function(pr,x)
# pr a product replacer for G, x an involution in G
local o,r,y,z;
r := Next(pr);
y := x^r;
# Now x and y generate a dihedral group
if x=y then return r; fi;
z := x*y;
o := Order(StripMemory(z));
if IsEvenInt(o) then
return z^(o/2);
else
return z^((o+1)/2)*r^(-1);
fi;
end );
InstallGlobalFunction( FindInvolutionCentraliser,
function(pr,x,n)
# pr a product replacer for G, x an involution in G, n an integer
local i,l,y;
l := [];
for i in [1..n] do # find 20 generators of the centraliser
y := FindCentralisingElementOfInvolution(pr,x);
AddSet(l,y);
od;
return l;
end );
InstallGlobalFunction( ReduceNumberOfGeneratorsUsingRecog,
function(l,Sizerecogniser)
local drop,i,ll,ri,ri2,si,wholegroup;
wholegroup := Group(StripMemory(l));
ri := Sizerecogniser(wholegroup);
si := Size(ri);
Print("Whole size: ",si,"\n");
i := Length(l);
l := ShallowCopy(l);
while i >= 1 do
drop := false;
ll := Concatenation(l{[1..i-1]},l{[i+1..Length(l)]});
ri := Sizerecogniser(Group(StripMemory(ll)));
if Size(ri) <> si then
ri2 := Sizerecogniser(Group(StripMemory(ll)));
if Size(ri2) <> Size(ri) then
if Size(ri2) = si then
# we do not need this generator
drop := true;
fi;
# otherwise non-conclusive, leave the generator where it is
fi;
else
# we most probably do not need this generator
drop := true;
fi;
if drop then
Remove(l,i);
i := i - 1;
Print("-\c");
else
i := i - 1;
Print("+\c");
fi;
od;
Print("\n");
return l;
end );
###############################################
# Find class representatives using powermaps: #
###############################################
InstallGlobalFunction( ClassMaker,
function(ct,cyc)
# ct an ordinary character table and cyc a list of class names we have
local cln,done,erg,i,l,m,p,po,pow,todo,start;
pow := ComputedPowerMaps(ct);
cln := ClassNames(ct);
cyc := List(cyc,LowercaseString);
for i in [1..Length(cyc)] do
p := Position(cyc[i],'-');
if p <> fail then
Info(InfoOrb,1,"Warning: Class name containing \"-\"!");
cyc[i] := cyc[i]{[1..p-1]};
fi;
od;
erg := [];
todo := [];
for i in [1..Length(cyc)] do
po := Position(cln,cyc[i]);
erg[po] := [i,1];
Add(todo,po);
od;
done := Length(cyc);
while done < Length(cln) and Length(todo) > 0 do
i := todo[1];
todo := todo{[2..Length(todo)]};
for p in [1..Length(pow)] do
if IsBound(pow[p]) then
m := pow[p];
if not(IsBound(erg[m[i]])) then
erg[m[i]] := [erg[i][1],erg[i][2]*p];
Add(todo,m[i]);
done := done + 1;
fi;
fi;
od;
od;
if done < Length(cln) then
Error("Unvorhergesehener Fall Nummer 1!");
fi;
if cln[1] = "1a" then
l := [[1,0]];
start := 2;
else
l := [];
start := 1;
fi;
for i in [start..Length(cln)] do
Add(l,erg[i]);
od;
return StraightLineProgramNC([l],Length(cyc));
end);
###########################
# Finding nice quotients: #
###########################
InstallGlobalFunction( OrbitStatisticOnVectorSpace,
function(gens,size,ti)
# gens must be a list of compressed matrices, size the order of the group
local len,nr,o,t,total,v;
v := ShallowCopy(gens[1][1]);
t := Runtime();
total := 0;
nr := 0;
while Runtime() < t + ti*1000 do
Randomize(v);
o := Orb(gens,v,OnRight,rec(orbsizebound := size,
hashlen := 3*size, report := 0));
Enumerate(o);
len := Length(o!.orbit);
total := total + len;
nr := nr + 1;
Print("Found length ",String(Length(o!.orbit),9),", have now ",
String(nr,4)," orbits, average length: ",
QuoInt(total+QuoInt(nr,2),nr)," \r");
od;
Print("\n");
end);
InstallGlobalFunction( OrbitStatisticOnVectorSpaceLines,
function(gens,size,ti)
# gens must be a list of compressed matrices, size the order of the group
local len,nr,o,t,total,v,c;
v := ShallowCopy(gens[1][1]);
t := Runtime();
total := 0;
nr := 0;
while Runtime() < t + ti*1000 do
Randomize(v);
c := PositionNonZero(v);
if c <= Length(v) then
v := v / v[c];
fi;
o := Orb(gens,v,OnLines,rec(orbsizebound := size,
hashlen := 3*size, report := 0));
Enumerate(o);
len := Length(o!.orbit);
total := total + len;
nr := nr + 1;
Print("Found length ",String(Length(o!.orbit),9),", have now ",
String(nr,4)," orbits, average length: ",
QuoInt(total+QuoInt(nr,2),nr)," \r");
od;
Print("\n");
end);
###############################################
# Finding nice generating sets for subgroups: #
###############################################
InstallGlobalFunction( ORB_PowerSet,
function(s)
local i,l,le,ll;
if Length(s) = 0 then
return [[]];
elif Length(s) = 1 then
return [[],[s[1]]];
else
l := ORB_PowerSet(s{[1..Length(s)-1]});
le := Length(l);
for i in [1..le] do
ll := ShallowCopy(l[i]);
Add(ll,s[Length(s)]);
Add(l,ll);
od;
return l;
fi;
end );
InstallGlobalFunction( ORB_SLPLineFromWord,
function(wo)
local li,i,j;
li := [];
i := 1;
while i <= Length(wo) do
j := i+1;
while j <= Length(wo) and wo[j] = wo[i] do
j := j + 1;
od;
if wo[i] > 0 then
Add(li,wo[i]);
Add(li,j-i);
else
Add(li,-wo[i]);
Add(li,-(j-i));
fi;
i := j;
od;
return li;
end );
InstallMethod( FindShortGeneratorsOfSubgroup, "without option rec or func",
[ IsGroup, IsGroup ],
function(G,U)
return FindShortGeneratorsOfSubgroup(G,U,
rec( membershiptest := \in, sizetester := Size ) );
end );
InstallMethod( FindShortGeneratorsOfSubgroup, "with option rec or func",
[ IsGroup, IsGroup, IsObject ],
function(G,U,opt)
local su,o,si,s,ps,min,minsi,subgens,subwords,i,l,membershiptest,sizetester;
if IsFunction(opt) then
membershiptest := opt;
sizetester := Size;
su := Size(U);
elif IsRecord(opt) then
if IsBound(opt.membershiptest) then
membershiptest := opt.membershiptest;
else
membershiptest := \in;
fi;
if IsBound(opt.sizetester) then
sizetester := opt.sizetester;
else
sizetester := Size;
fi;
if IsBound(opt.sizeU) then
su := opt.sizeU;
else
if HasSize(U) then
su := Size(U);
else
su := sizetester(U);
fi;
fi;
fi;
if su = 1 then # the trivial subgroup is easy to generate:
return rec( gens := [One(U)],
slp := StraightLineProgram( [[[1,0]]],
Length(GeneratorsOfGroup(G)) ) );
fi;
o := Orb(GeneratorsOfGroup(G),One(G),OnRight,
rec( lookingfor := function(o,x) return membershiptest(x,U); end,
schreier := true ) );
Enumerate(o);
subgens := [o!.orbit[o!.found]];
subwords := [TraceSchreierTreeForward(o,o!.found)];
l := 1; # will always be the length of subgens and subwords
si := sizetester(Group(ShallowCopy(subgens)));
Info(InfoOrb,2,"Found subgroup of size ",si,":",subwords);
if si = su then
# Cyclic subgroup
s := StraightLineProgram([[ORB_SLPLineFromWord(subwords[1])]],
Length(GeneratorsOfGroup(G)));
return rec( gens := subgens, slp := s );
fi;
while true do
Enumerate(o);
Add(subgens,o!.orbit[o!.found]);
Add(subwords,TraceSchreierTreeForward(o,o!.found));
l := l + 1;
if l <> Length(subgens) or l <> Length(subwords) then
Error();
fi;
s := sizetester(Group(ShallowCopy(subgens)));
if s = su then
# OK, we have got a generating set:
# Now try shortening generating set:
Info(InfoOrb,2,"Found full subgroup of size ",si,":",subwords);
Info(InfoOrb,2,"Need ",l," generators.");
if l <= 8 then
ps := ORB_PowerSet([1..l]);
min := Length(ps);
minsi := l;
for i in [2..Length(ps)-1] do
s := sizetester(Group(subgens{ps[i]}));
if s = su and Length(ps[i]) < minsi then
min := i;
minsi := Length(ps[i]);
Info(InfoOrb,2,"Need ",minsi," generators.");
fi;
od;
subgens := subgens{ps[min]};
subwords := subwords{ps[min]};
else
i := 1;
while i <= Length(subgens) do
s := sizetester(
Group(subgens{Concatenation([1..i-1],[i+1..l])}));
if s = su then # we do not need generator i
Remove(subgens,i);
Remove(subwords,i);
l := l - 1;
else
i := i + 1;
fi;
od;
fi;
l := List(subwords,ORB_SLPLineFromWord);
s := StraightLineProgram([l],Length(GeneratorsOfGroup(G)));
return rec( gens := subgens, slp := s );
fi;
if s=si then
Unbind(subgens[l]);
Unbind(subwords[l]);
l := l - 1;
else
si := s;
Info(InfoOrb,2,"Found subgroup of size ",si,":",subwords);
fi;
od;
end );
##############################################################
# Helpers for permutation characters for certain operations: #
##############################################################
InstallGlobalFunction( NumberFixedVectors,
function( mat )
local f,n,q;
n := NullspaceMat(mat - mat^0);
f := BaseField(mat);
q := Size(f);
return q^Length(n);
end );
InstallGlobalFunction( NumberFixedLines,
function( mat )
local el,f,n,nr,q;
f := BaseField(mat);
q := Size(f);
nr := 0;
for el in f do
if not(IsZero(el)) then
n := NullspaceMat(mat - el*mat^0);
nr := nr + (q^Length(n)-1)/(q-1);
fi;
od;
return nr;
end );
InstallGlobalFunction( SpacesOfFixedLines,
function( mats )
local el,f,i,n,q,spcs;
f := BaseField(mats);
q := Size(f);
spcs := List(mats,i->[]);
for i in [1..Length(mats)] do
for el in f do
if not(IsZero(el)) then
n := NullspaceMat(mats[i] - el*mats[i]^0);
Add(spcs[i],n);
fi;
od;
od;
# to be continued...
return spcs;
end );
##################################################
# Helpers for making short SLPs from word lists: #
##################################################
InstallGlobalFunction( SLPForWordList,
function( wordlist, nrgens )
local bestlen,bestn,bestword,havewords,i,j,l,line,n,p,slp,w,where,wo,
wordset,writepos,ww;
havewords := [[]];
where := [0];
for i in [1..nrgens] do
Add(havewords,[i]);
Add(where,i);
Add(havewords,[-i]);
Add(where,-i);
od;
wordset := Set(wordlist);
slp := [];
writepos := nrgens+1;
for w in wordset do
if not(w in havewords) then
wo := ShallowCopy(w);
line := [];
while Length(wo) > 0 do
# Look through all havewords and see which one helps most:
bestlen := 0;
bestword := 0;
bestn := 0;
for j in [2..Length(havewords)] do
ww := havewords[j];
l := Length(ww);
if l <= Length(wo) and wo{[1..l]} = ww then
# How often does it fit here?
n := 1;
while (n+1) * l <= Length(wo) and
ww = wo{[n*l+1..(n+1)*l]} do
n := n + 1;
od;
if n*l > bestlen then
bestlen := n*l;
bestword := j;
bestn := n;
fi;
fi;
od;
if where[bestword] > 0 then
Add(line,where[bestword]);
Add(line,bestn);
else
Add(line,-where[bestword]);
Add(line,-bestn);
fi;
l := Length(havewords[bestword]);
wo := wo{[l*bestn+1..Length(wo)]};
od;
Add(slp,line);
Add(havewords,w);
Add(where,writepos);
Add(havewords,-Reversed(w));
Add(where,-writepos);
writepos := writepos+1;
fi;
od;
line := [];
for w in wordlist do
p := Position(havewords,w);
if where[p] > 0 then
Add(line,[where[p],1]);
elif where[p] = 0 then # The empty word
Add(line,[1,0]);
else
Add(line,[-where[p],-1]);
fi;
od;
Add(slp,line);
return [StraightLineProgram(slp,nrgens),havewords,where];
end );
#############################################################################
# A generic way to find stabilizers:
#############################################################################
InstallGlobalFunction( ORB_EstimatePermGroupSize,
function( gens )
local g,s;
g := Group(gens);
s := StabChain(g,rec( random := 900 ));
return SizeStabChain(s);
end );
InstallGlobalFunction( ORB_FindStabilizerMC,
function( gens, pt, op, memory, limit, errorbound, estimatefunc, opt )
local counter,done,est,gensi,mem,o,oldest,orbpart,p,pr,prob,stab,tries,
x,y,z;
mem := Memory(pt);
if mem = fail then
mem := SIZE_OBJ(pt);
if mem = 0 then mem := GAPInfo.BytesPerVariable; fi;
fi;
orbpart := QuoInt(memory,mem);
o := Orb( gens, pt, op, rec( report := QuoInt(orbpart,10),
schreier := true,
hashlen := NextPrimeInt( orbpart*2 ) ) );
Info( InfoOrb, 2, "Enumerating up to ",orbpart," elements of orbit..." );
Enumerate(o,orbpart);
# Now find stabilizer generators:
stab := [];
done := false;
tries := 0;
pr := ProductReplacer(gens,opt);
est := 1;
prob := 1;
gensi := List(gens,x->x^-1);
while prob > errorbound do
counter := 0;
repeat
counter := counter + 1;
x := Next(pr);
y := op(pt,x);
p := Position(o,y);
Print(".\c");
until p <> fail or counter > limit;
Print("\n");
if p = fail then
Info( InfoOrb, 1, "Giving up..." );
return rec(stab := stab, prob := fail);
fi;
z := x * Product(gensi{TraceSchreierTreeBack(o,p)});
if IsOne(z) then
Info( InfoOrb, 2, "Found identity element." );
else
AddSet(stab,z);
oldest := est;
if IsFunction(estimatefunc) then
Info( InfoOrb, 2, "Starting estimation (old est=", est,")..." );
est := estimatefunc(stab);
Info( InfoOrb, 2, "Finished estimation: ", est );
fi;
if est = oldest then
prob := prob / 2;
Info( InfoOrb, 2, "New estimate is the same, prob=",prob );
else
prob := 1;
Info( InfoOrb, 2, "New estimate is ",est,", prob=",prob );
fi;
fi;
od;
return rec(stab := stab, prob := prob);
end );
InstallGlobalFunction( ORB_FindNeedleMappers,
function( gens, pt, op, needles, memory, mappers, opt )
local foundmappers,gensi,i,j,lookfunc,mem,o,p,r,size,x;
size := fail;
if IsGroup(gens) then
if HasSize(gens) then
size := Size(gens);
fi;
gens := GeneratorsOfGroup(gens);
fi;
if not IsBound(opt.orblen) then
mem := Memory(pt);
if mem = fail then
mem := SIZE_OBJ(pt);
if mem = 0 then
mem := GAPInfo.BytesPerVariable;
fi;
fi;
opt.orblen := QuoInt(memory,Length(needles)*mem);
if size <> fail and opt.orblen * opt.orblen > size then
opt.orblen := 2^(QuoInt(Log2Int(size),2)+1);
fi;
fi;
o := [];
Info( InfoOrb, 1, "Enumerating ",Length(needles)," orbits up to length ",
opt.orblen, "..." );
for i in [1..Length(needles)] do
o[i] := Orb( gens, needles[i], op, rec( report := QuoInt(opt.orblen,5),
schreier := true ) );
Enumerate(o[i],opt.orblen);
Info( InfoOrb, 2, "Finished one orbit." );
od;
j := fail;
p := fail;
lookfunc := function( x )
# Uses o from outer function!
# Writes j and p in outer function!
# This uses the feature of GAP, that j and p here always refer to
# the outer function ORB_FindNeedleMappers
# Reads pt and op from outer function.
j := 1;
while j <= Length(o) do
p := Position(o[j],op(pt,x));
if p <> fail then return true; fi;
j := j + 1;
od;
j := fail;
p := fail;
return false;
end;
gensi := List(gens,x->x^-1);
r := RandomSearcher( gens, lookfunc, opt );
# we just hand down the options record to the random searcher!
foundmappers := [];
i := 1;
while Length(foundmappers) <= mappers and i <= 15 * mappers do
x := Search(r);
if x = fail then return foundmappers; fi;
AddSet(foundmappers,x * Product(gensi{TraceSchreierTreeBack(o[j],p)}));
Info( InfoOrb, 2, "Found a needle mapper, have now ",
Length(foundmappers), ".");
i := i + 1;
od;
return foundmappers;
end );
############################################################################
# A method to find transversals in matrix groups:
############################################################################
InstallGlobalFunction( FindWordsForRightTransversal,
function( Ggens, Hgens, index )
# Ggens and Hgens must be lists of invertible matrices over a
# finite field. The group H generated by Hgens must be a subgroup
# of G, which is generated by Ggens, index must be the index.
# This function returns words in Ggens that are a complete set
# of coset representatives of the right cosets {Hg}.
local Check,i,inv,l,n,res,v;
# Find invariant vectors:
Info(InfoOrb,2,"Computing invariant spaces of generators...");
l := List(Hgens,x->NullspaceMat(x-One(x)));
Info(InfoOrb,2,"Found dimensions: ",List(l,Length));
if Length(l) = 1 then
inv := l[1];
else
inv := SumIntersectionMat(l[1],l[2]);
inv := inv[2];
Info(InfoOrb,2,"Intersection has dim ",NrRows(inv));
i := 3;
while i <= Length(l) and NrRows(inv) > 0 do
inv := SumIntersectionMat(inv,l[i]);
inv := inv[2];
Info(InfoOrb,2,"Intersection has dim ",NrRows(inv));
i := i + 1;
od;
fi;
# Now inv contains vectors that span the invariant space.
if NrRows(inv) = 0 then
Info(InfoOrb,1,"no H-invariants except 0");
return fail;
fi;
Unbind(l);
Check := function( pt, op )
local i,o,words;
o := Orb(Ggens,pt,op,rec( report := 2000, schreier := true,
hashlen := NextPrimeInt( index*3 ) ));
Enumerate(o);
if Length(o) = index then
# Hurra!
words := [[]];
for i in [2..Length(o)] do
Add(words,TraceSchreierTreeForward(o,i));
od;
return rec( words := words, orbit := o );
fi;
Info(InfoOrb,2,"Found orbit length ",Length(o));
return fail;
end;
# First try whether one basis vector provides an orbit of the
# right length:
Info(InfoOrb,2,"Trying to find a nice orbit on vectors...");
for i in [1..NrRows(inv)] do
Info(InfoOrb,2,"Using vector #",i,"(",NrRows(inv),")");
res := Check(ShallowCopy(inv[i]),OnRight);
if res <> fail then return res; fi;
od;
# Now try tuples of vectors:
n := 2; # first pairs, then maybe more
repeat
l := Matrix([],NrCols(inv),inv);
for i in [1..n] do
v := ZeroVector(NrRows(inv),inv[1]);
Randomize(v);
Add(l,v*inv);
od;
Info(InfoOrb,2,"Using ",n,"-tuple...");
res := Check(l,OnRight);
if res <> fail then return res; fi;
n := n + 1;
until n = NrRows(inv);
# Now try the last resort:
Info(InfoOrb,2,"Trying full invariant space...");
return Check(inv,OnRight);
end );
InstallGlobalFunction( FindWordsForLeftTransversal,
function( Ggens, Hgens, index )
# The same as above, but gives words for coset reps of the left
# cosets {gH}. This is done by using the transposition
# anti-automorphism. This maps a right transversal of the
# transposed groups to a left transversal. Since we take the
# transposed generators in both groups, the words have to be
# reversed. Note that the orbit is the one obtained for the
# transposed generators!
local res;
res := FindWordsForRightTransversal(List(Ggens,TransposedMat),
List(Hgens,TransposedMat),index);
if res = fail then return fail; fi;
res.words := List(res.words,Reversed);
return res;
end );
############################################################################
# Find transforming matrices:
############################################################################
InstallGlobalFunction( TransformingMatsLSE,
function( A, B )
local M,n;
# Returns the matrix of the linear system of equations to solve TA=BT
# such that the nullspace of that matrix gives a basis of the solution
# space in flat notation.
# "Flat notation" means Concatenation(T) instead of T.
# Each row of M will be an equation:
# M := [ [A[1][1],A[2][1],A[3][1],0,0,0,0,0,0],
# [A[1][2],A[2][2],A[3][2],0,0,0,0,0,0],
# [A[1][3],A[2][3],A[3][3],0,0,0,0,0,0],
# [0,0,0,A[1][1],A[2][1],A[3][1],0,0,0],
# [0,0,0,A[1][2],A[2][2],A[3][2],0,0,0],
# [0,0,0,A[1][3],A[2][3],A[3][3],0,0,0],
# [0,0,0,0,0,0,A[1][1],A[2][1],A[3][1]],
# [0,0,0,0,0,0,A[1][2],A[2][2],A[3][2]],
# [0,0,0,0,0,0,A[1][3],A[2][3],A[3][3]] ]
# -[ [B[1][1],0,0,B[1][2],0,0,B[1][3],0,0],
# [0,B[1][1],0,0,B[1][2],0,0,B[1][3],0],
# [0,0,B[1][1],0,0,B[1][2],0,0,B[1][3]],
# [B[2][1],0,0,B[2][2],0,0,B[2][3],0,0],
# [0,B[2][1],0,0,B[2][2],0,0,B[2][3],0],
# [0,0,B[2][1],0,0,B[2][2],0,0,B[2][3]],
# [B[3][1],0,0,B[3][2],0,0,B[3][3],0,0],
# [0,B[3][1],0,0,B[3][2],0,0,B[3][3],0],
# [0,0,B[3][1],0,0,B[3][2],0,0,B[3][3]] ];
n := NrRows(A);
if n <> NrCols(A) or n <> NrRows(B) or n <> NrCols(B) then
Error("need square matrices of same size");
fi;
M := KroneckerProduct(OneMutable(A),TransposedMat(A))
- KroneckerProduct(B,OneMutable(B));
return TransposedMat(M);
end);
InstallGlobalFunction( TransformingMats,
function(A, B)
return NullspaceMat(TransformingMatsLSE(A,B));
end );
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <https://www.gnu.org/licenses/>.
##
[ Dauer der Verarbeitung: 0.41 Sekunden
(vorverarbeitet)
]
|
2026-03-28
|
