Impressum genss.gi   Interaktion und
Portierbarkeitunbekannt

#############################################################################
##
##  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

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