#############################################################################
##
##  This file is part of GAP, a system for computational discrete algebra.
##  This file's authors include Ákos Seress.
##
##  Copyright of GAP belongs to its developers, whose names are too numerous
##  to list here. Please refer to the COPYRIGHT file for details.
##
##  SPDX-License-Identifier: GPL-2.0-or-later
##
##  This  file contains  the  functions for  a random Schreier-Sims algorithm
##  with verification.
##

#############################################################################
##
#F  StabChainRandomPermGroup( <gens>, <id>, <opts> )  .  random Schreier-Sims
##
##  The  method  consists of  2  phases:  a  heuristic construction  and
##  either a deterministic or two random checking phases.
##  In the random checking  phases,  we  take random  elements of
##  created set, multiply by  random subproduct of  generators, and sift down
##  by dividing with the  appropriate coset representatives. The stabchain is
##  correct when  all  siftees are  ().    In the first  checking  phase, all
##  computations are carried out with words  in strong generators and siftees
##  are checked by plugging in random points in words; in the second checking
##  phase, permutations are multiplied.
##
##  During the construction, we create new records for stabilizers:
##  S.aux = set of strong generators for S, temporarily created during the
##          construction. On the other hand, S.generators contains a strong
##          generating set for S which would have been created by the
##          deterministic method
##  S.treegen = elements of S.aux used in S.transversal
##  S.treegeninv = inverses of S.treegen; also used in S.transversal,
##                 and they are also elements of S.aux
##  S.treedepth = depth of Schreier tree of S
##  S.diam = sum of treedepths in stabilizer chain of S
##
##  After the  stabilizer chain  is  ready,  the  extra records are  deleted.
##  Transversals are rebuilt using the generators in the .generators records.
##
InstallGlobalFunction( StabChainRandomPermGroup, function( gens, id, options)
    local S,        # stabilizer chain
          degree,   # degree of S
          givenbase,# list of points from which first base points should come
          correct,  # boolean; true if a correct base is given
          size,     # size of <G> as constructed
          order,    # size of G if given in input
          limit,    # upper bound on Size(G) given in input
          orbits,   # list of orbits of G
          orbits2,  # list of orbits of G
          k,        # number of pairs of generators checked
          param,    # list of parameters guiding number of repetitions
                    # in random constructions
          where,    # list indicating which orbit contains points in domain
          basesize, # list; i^th entry = number of base points in orbits[i]
          i,j,
          ready,    # boolean; true if stabilizer chain ready
          warning,  # used at warning if given and computed size differ
          new,      # list of permutations to be added to stab. chain
          result,   # output of checking phase; nontrivial if stabilizer
                    # chain is incorrect
          base,     # ordering of domain from which base points are taken
          missing,  # if a correct base was provided by input, missing
                    # contains those points of it which are not in
                    # constructed base
          T;        # a stabilizer chain containing the usual records

    S:= rec( generators := ShallowCopy( gens ), identity := id );
    if options.random = 1000 then
       #case of deterministic computation with known size
       k := 1;
    else
       k:=First([1..14],x->(3/5)^x<1-options.random/1000);
    fi;

    degree := LargestMovedPoint( S.generators );

    if IsBound( options.knownBase) and
      Length(options.knownBase)<4+LogInt(degree,10)  then
        param:=[k,4,0,0,0,0];
    else
        param:=[QuoInt(k,2),4,QuoInt(k+1,2),4,50,5];
        options:=ShallowCopy(options);
        Unbind(options.knownBase);
    fi;
    if options.random <= 200 then
       param[2] := 2;
       param[4] := 2;
    fi;

    #param[1] = number of pairs of random subproducts from generators in
    #           first checking phase
    #param[2] = (number of random elements from created set)/S.diam
    #           in first checking phase
    #param[3] = number of pairs of random subproducts from generators in
    #           second checking phase
    #param[4] = (number of random elements from created set)/S.diam
    #           in second checking phase
    #param[5] = maximum size of orbits in  which we evaluate words on all
    #           points of orbit
    #param[6] = minimum number of random points from orbit to plug in to check
    #           whether given word is identity on orbit


    # prepare input of construction
    if IsBound(options.base)  then
        givenbase := options.base;
    else
        givenbase := [];
    fi;

    if IsBound(options.size) then
        order := options.size;
        warning := 0;
        limit := 0;
    else
        order := 0;
        if IsBound(options.limit) then
            limit := options.limit;
        else
            limit := 0;
        fi;
    fi;

    if IsBound( options.knownBase )  then
        correct := true;
    else
        correct := false;
    fi;

    if correct then
        # if correct  base was given as input, no need for orbit information
        base:=Concatenation(givenbase,Difference(options.knownBase,givenbase));
        missing := Set(options.knownBase);
        basesize := [];
        where := [];
        orbits := [];
    else
        # create ordering of domain used in choosing base points and
        # compute orbit information
        base:=Concatenation(givenbase,Difference([1..degree],givenbase));
        missing:=[];
        orbits2:=OrbitsPerms(S.generators,[1..degree]);
        #throw away one-element orbits
        orbits:=[];
        j:=0;
        for i in [1..Length(orbits2)] do
            if Length(orbits2[i]) >1 then
               j:=j+1; orbits[j]:= orbits2[i];
            fi;
        od;
        basesize:=[];
        where:=[];
        for i in [1..Length(orbits)] do
            basesize[i]:=0;
            for j in [1..Length(orbits[i])] do
                where[orbits[i][j]]:=i;
            od;
        od;
        # temporary solution to speed up of handling of lots of small orbits
        # until compiler
        if Length(orbits) > degree/40 then
           param[1] := 0;
           param[3] := k;
        fi;
    fi;

    ready:=false;
    new:=S.generators;

    while not ready do
        SCRMakeStabStrong
           (S,new,param,orbits,where,basesize,base,correct,missing,true);
        # last parameter of input is true if function called for original G
        # in recursive calls on stabilizers, it is false
        # reason: on top level,
        #         we do not want to add anything to generating set

        # start checking
        size := 1;  T := S;
        while Length( T.generators ) <> 0  do
            size := size * Length( T.orbit );
            T := T.stabilizer;
        od;
        if size = order  or  size = limit  then
            ready := true;
        elif size < order then
            # we have an incorrect stabilizer chain
            # repeat checking until a new element is discovered
            result := id;
            if options.random = 1000 then
                if correct then
                   result := SCRStrongGenTest(S,[1,10/S.diam,0,0,0,0],orbits,
                                      basesize,base,correct,missing);
                else
                   result := SCRStrongGenTest2(S,[0,0,1,10/S.diam,0,0]);
                fi;
                if result = id then
                   T := SCRRestoredRecord(S);
                   result := VerifySGS(T,missing,correct);
                fi;
                if result = id then
                   Print("Warning, computed and given size differ","\n");
                   ready := true;
                   S := T;
                else
                   if not IsPerm(result) then
                      repeat
                         result := SCRStrongGenTest2(S,[0,0,1,10/S.diam,0,0]);
                      until result <> id;
                   fi;
                   new := [result];
                fi;
            else
               warning := 0;
               if correct then
                 # if correct base was provided, it is enough to check
                 # images of base points to check whether a word is trivial
                 # no need for second checking phase
                 while result = id do
                    warning := warning + 1;
                    if warning > 5 then
                       Print("Warning, computed and given size differ","\n");
                    fi;
                    result := SCRStrongGenTest(S,param,orbits,
                                      basesize,base,correct,missing);
                 od;
               else
                 while result = id do
                    warning := warning + 1;
                    if warning > 5 then
                       Print("Warning, computed and given size differ","\n");
                    fi;
                    result:=SCRStrongGenTest(S,param,orbits,
                                    basesize,base,correct,missing);
                    if result = id then
                        # Print("entering SGT2","\n");
                        result:=SCRStrongGenTest2(S,param);
                    fi;
                 od;
               fi;  # correct or not
               new:=[result];
            fi;   # end of random checking or not, when size is known
        else
            #no information or only upper bound about Size(S)
            if options.random = 1000 then
              if correct then
                   result := SCRStrongGenTest(S,[1,10/S.diam,0,0,0,0],orbits,
                                      basesize,base,correct,missing);
               else
                   result := SCRStrongGenTest2(S,[0,0,1,10/S.diam,0,0]);
               fi;
               if result = id then
                  T := SCRRestoredRecord(S);
                  result := VerifySGS(T,missing,correct);
               fi;
               if result = id then
                  S := T;
                  ready := true;
               else
                  if not IsPerm(result) then
                     repeat
                        result := SCRStrongGenTest2(S,[0,0,1,10/S.diam,0,0]);
                     until result <> id;
                  fi;
                  new := [result];
               fi;
            else
               # Print("entering SGT", "\n");
               result:=SCRStrongGenTest(S,param,orbits,basesize,
                                  base,correct,missing);
               if result <> id then
                  new:=[result];
               elif correct then
                  # no need for second checking phase
                  ready:=true;
               else
                   # Print("entering SGT2","\n");
                   result:=SCRStrongGenTest2(S,param);
                   if result = id then
                      ready:=true;
                   else
                      new:=[result];
                   fi;
               fi;
            fi;   # random checking or not, when size is not known
        fi;   # size known
    od;

    #restore usual record elements
    if not IsBound(S.labels) then
       S := SCRRestoredRecord(S);
    fi;
    return S;
end );


#############################################################################
##
#F  SCRMakeStabStrong( ... )  . . . . . . . . . . . . . . . . . . . . . local
##
##  heuristic stabilizer  chain  construction, with one random  subproduct on
##  each  level,  and one or two  (defined  by mlimit) random cosets  to make
##  Schreier generators
##
InstallGlobalFunction( SCRMakeStabStrong,
    function ( S, new, param, orbits, where, basesize, base,
                                correct, missing, top )
    local   x,m,j,l,      # loop variables
            ran1,         # random permutation
            string,       # random 0-1 string
            w,            # random subproduct of generators
            len,          # number of generators of S
            mlimit,       # number of random elements to be tested
            coset,        # word representing coset of S
            residue,      # first component: remainder of Schreier generator
                          # after factorization; second component > 0
                          # if factorization unsuccessful
            jlimit,       # number of random points to plug into residue[1]
            ran,          # index of random point in an orbit of S
            g,            # permutation to be added to S.stabilizer
            gen,          # permutations used in S.transversal
            inv,          # their inverses
            firstmove;    # first point of base moved by an element of new

    if new <> [] then
        firstmove := First( base, x->ForAny( new, gen->x^gen<>x ) );
        # if necessary add a new stabilizer to the stabchain
        if not IsBound(S.stabilizer) then
            S.orbit                    := [firstmove];
            S.transversal              := [];
            S.transversal[S.orbit[1]]  := S.identity;
            S.generators               := [];
            S.treegen                  := [];
            S.treegeninv               := [];
            S.stabilizer               := rec();
            S.stabilizer.identity      := S.identity;
            S.stabilizer.aux           := [];
            S.stabilizer.generators    := [];
            S.stabilizer.diam          := 0;
            if not correct then
                basesize[where[S.orbit[1]]]
                    := basesize[where[S.orbit[1]]] + 1;
            fi;
            missing := Difference( missing, [ firstmove ] );
        else
            if Position(base,firstmove) < Position(base,S.orbit[1]) then
                S.stabilizer               := ShallowCopy(S);
                S.orbit                    := [firstmove];
                S.transversal              := [];
                S.transversal[S.orbit[1]]  := S.identity;
                S.generators := ShallowCopy(S.stabilizer.generators);
                S.treegen                  := [];
                S.treegeninv               := [];
                if not correct then
                    basesize[where[S.orbit[1]]]
                      := basesize[where[S.orbit[1]]] + 1;
                fi;
                missing := Difference( missing, [ firstmove ] );
            fi;
        fi;

        # on top level, we want to keep the original generators
        if not top or Length(S.generators) = 0  then
          for j in new do
            StretchImportantSLPElement(j);
          od;
          Append(S.generators,new);
        fi;
        #construct orbit of basepoint
        SCRSchTree(S,new);

        #check whether new elements are really new in the system
        while new <> [] do
            g := SCRSift( S, Remove(new) );
            if g <> S.identity then
                SCRMakeStabStrong(S.stabilizer,[g],param,orbits,
                                  where,basesize,base,correct,missing,false);
                S.diam:=S.treedepth+S.stabilizer.diam;
                S.aux:=Concatenation(S.treegen,
                                     S.treegeninv,S.stabilizer.aux);
            fi;
        od;
    fi;

    #check random Schreier generators
    gen := Concatenation(S.treegen,S.treegeninv,[S.identity]);
    inv := Concatenation(S.treegeninv,S.treegen,[S.identity]);
    len := Length(S.aux);
    #in case of more than one generator for S, form a random subproduct
    #otherwise, use the generator
    if len > 1 then
        ran1 := SCRRandomPerm(len);
        string := SCRRandomString(len);
        w := S.identity;
        for x in [1..len] do
            w := w*(S.aux[x^ran1]^string[x]);
        od;
    else
        w:=S.aux[1];
    fi;

    # take random coset(s)
    mlimit:=1;
    m:=0;
    while m < mlimit do
        m := m+1;
        ran := Random(1, Length(S.orbit));
        coset := CosetRepAsWord(S.orbit[1],S.orbit[ran],S.transversal);
        coset := InverseAsWord(coset,gen,inv);
        if w <> S.identity then
            # form Schreier generator and factorize
            Add(coset,w);
            residue := SiftAsWord(S,coset);
            # check whether factorization is successful
            if residue[2] > 0  then
                # factorization is unsuccessful; use remainder for
                # construction in stabilizer
                g := Product(residue[1]);
                SCRMakeStabStrong(S.stabilizer,[g],param,orbits,where,
                                  basesize,base,correct,missing,false);
                S.diam := S.treedepth+S.stabilizer.diam;
                S.aux := Concatenation(S.treegen,S.treegeninv,
                                       S.stabilizer.aux);
                # get out of current loop
                m := 0;
            elif correct then
                # enough to check images of points in given base
                l := 0;
                while l < Length(missing) do
                    l := l+1;
                    if ImageInWord(missing[l],residue[1]) <> missing[l] then
                        # factorization is unsuccessful;
                        # use remainder for construction in stabilizer
                        g := Product(residue[1]);
                        SCRMakeStabStrong(S.stabilizer,[g],param,
                                          orbits,where,basesize,
                                          base,correct,missing,false);
                        S.diam := S.treedepth+S.stabilizer.diam;
                        S.aux := Concatenation(S.treegen,
                                               S.treegeninv,S.stabilizer.aux);
                        # get out of current loop
                        m := 0;
                        l := Length(missing);
                    fi;
                od;
            else
                l:=0;
                while l < Length(orbits) do
                    l:=l+1;
                    if Length(orbits[l]) > param[5] then
                        # in large orbits, plug in random points
                        j:=0;
                        jlimit:=Maximum(param[6],basesize[l]);
                        while j < jlimit do
                            j:=j+1;
                            ran:=Random(1, Length(orbits[l]));
                            if ImageInWord(orbits[l][ran],residue[1])
                              <> orbits[l][ran]
                            then
                                # factorization is unsuccessful;
                                # use remainder for construction in stabilizer
                                g := Product(residue[1]);
                                SCRMakeStabStrong(S.stabilizer,[g],param,
                                                  orbits,where,basesize,
                                                  base,correct,missing,false);
                                S.diam := S.treedepth+S.stabilizer.diam;
                                S.aux := Concatenation(S.treegen,S.treegeninv,
                                                       S.stabilizer.aux);
                                # get out of current loop
                                m := 0;
                                j := jlimit;
                                l := Length(orbits);
                            fi;
                        od; #j loop
                    else
                        # in small orbits, check images of all points
                        j := 0;
                        while j < Length(orbits[l]) do
                            j := j+1;
                            if ImageInWord(orbits[l][j],residue[1])
                              <> orbits[l][j]
                            then
                                # factorization is unsuccessful;
                                # use remainder for construction in stabilizer
                                g := Product(residue[1]);
                                SCRMakeStabStrong(S.stabilizer,[g],param,
                                                  orbits,where,basesize,
                                                  base,correct,missing,false);
                                S.diam := S.treedepth+S.stabilizer.diam;
                                S.aux := Concatenation(S.treegen,S.treegeninv,
                                                       S.stabilizer.aux);
                                # get out of current loop
                                m := 0;
                                j := Length(orbits[l]);
                                l := Length(orbits);
                            fi;
                        od; #j loop
                    fi;
                od; #l loop
            fi;
        fi;
    od; #m loop

end );


#############################################################################
##
#F  SCRStrongGenTest( ... ) . . . . . . . . . . . . . . . . . . . . . . local
##
##  tests whether product of a random element of S and a random subproduct of
##  the strong generators  of S is in S.   Computations are  carried out with
##  words in generators   representing  group  elements; random   points  are
##  plugged in to test whether a word represents  the identity.  If SGS for S
##  not complete, returns a permutation not in S.
##
InstallGlobalFunction( SCRStrongGenTest,
    function ( S, param, orbits, basesize,
                               base, correct, missing)
    local   x,i,k,m,j,l,  # loop variables
            ran1,ran2,    # random permutations
            string,       # random 0-1 string
            w,            # list containing random subproducts of generators
            len,          # number of generators of S
            len2,         # length of short random subproduct
            mlimit,       # number of random elements to be tested
            ranword,      # random element of S as a word in generators
            residue,      # first component: remainder of ranword
                          # after factorization; second component > 0
                          # if factorization unsuccessful
            jlimit,       # number of random points to plug into residue[1]
            ran,          # index of random point in an orbit of S
            g;            # product of residue[1]

    k := 0;
    while k < param[1] do
        k := k+1;
        len := Length(S.aux);
        #in case of large S.aux, form random subproducts
        #otherwise, try all of them
        if len > 2*param[3] then
            ran1 := SCRRandomPerm(len);
            ran2 := SCRRandomPerm(len);
            len2 := Random(1, QuoInt(len,2));
            string := SCRRandomString(len+len2);

            # we form two random subproducts:
            # w[1] in a random ordering of all generators
            # w[2] in a random ordering of a random subset of them
            w:=[];
            w[1] := S.identity;
            for x in [1 .. len] do
                w[1] := w[1]*(S.aux[x^ran1]^string[x]);
            od;
            w[2] := S.identity;
            for x in [1 .. len2] do
                w[2] := w[2]*(S.aux[x^ran2]^string[x+len]);
            od;

        else
            # take next two generators of S (unless only one is left)
            w := [];
            w[1] := S.aux[2*k-1];
            if len > 2*k-1 then
                w[2] := S.aux[2*k];
            else
                w[2] := S.identity;
            fi;
        fi;

        # take random elements of S as words
        m := 0;
        mlimit := param[2]*S.diam;
        while m < mlimit do
            m:=m+1;
            ranword := RandomElmAsWord(S);
            i := 0;
            while i < 2 do
                i := i+1;
                if w[i] <> S.identity then
                    Append(ranword,[w[i]]);
                    residue := SiftAsWord(S,ranword);
                    if residue[2]>0 then
                        # factorization is unsuccessful;
                        # remainder is witness that SGS for S is not complete
                        g := Product(residue[1]);
                        # Print("k=",k," i=",i," m=",m," mlimit=",mlimit,"\n");
                        return g;
                    elif correct then
                        # enough to check whether base points are fixed
                        l:=0;
                        while l < Length(missing) do
                            l:=l+1;
                            if ImageInWord(missing[l],residue[1])
                              <> missing[l]
                            then
                                # remainder is not in S
                                g := Product(residue[1]);
                                return g;
                            fi;
                        od;
                    else
                        # plug in points from each orbit to check whether
                        # action on orbit is trivial
                        l:=0;
                        while l < Length(orbits) do
                            l:=l+1;
                            if Length(orbits[l]) > param[5] then
                                # on large orbits, plug in random points
                                j := 0;
                                jlimit := Maximum(param[6],basesize[l]);
                                while j < jlimit do
                                    j := j+1;
                                    ran := Random(1, Length(orbits[l]));
                                    if ImageInWord(orbits[l][ran],residue[1])
                                      <> orbits[l][ran]
                                    then
                                        #remainder is not in S
                                        g := Product(residue[1]);
                                        return g;
                                    fi;
                                od; #j loop
                            else
                                # on small orbits, plug in all points
                                j := 0;
                                while j < Length(orbits[l]) do
                                    j := j+1;
                                    if ImageInWord( orbits[l][j],residue[1] )
                                      <> orbits[l][j]
                                    then
                                        # remainder is not in S
                                        g:=Product(residue[1]);
                                        return g;
                                    fi;
                                od; #j loop
                            fi;
                        od; #l loop
                    fi;
                fi;
            od; #i loop
        od; #m loop
        if len <= 2*k then
            #finished making Schr. generators with all in S.aux
            k := param[1];
        fi;
    od; #k loop

    return S.identity;
end );


#############################################################################
##
#F  SCRSift( <S>, <g> ) . . . . . . . . . . . . . . . . . . . . . . . . local
##
##  tries to factor g as product of cosetreps in S; returns remainder
##
InstallGlobalFunction( SCRSift, function(S,g)
     local result;

     if not IsInternalRep( g ) then
       return SiftedPermutation(S, g);
     fi;

     result :=  SCR_SIFT_HELPER(S, g, Maximum(LargestMovedPoint(g),LargestMovedPoint(S!.generators)));

#     Assert(2,result = SiftedPermutation(S, g));
     return result;
end);



#############################################################################
##
#F  SCRStrongGenTest2( <S>, <param> ) . . . . . . . . . . . . . . . . . local
##
##  tests whether product of a random element of S and a random subproduct of
##  the  strong generators of S is  in S.  Computations  are carried out with
##  complete permutations.
##
InstallGlobalFunction( SCRStrongGenTest2, function ( S, param )
    local   x,i,k,m,      # loop variables
            ran1,ran2,    # random permutations
            string,       # random 0-1 string
            w,            # list containing random subproducts of generators
            len,          # number of generators of S
            len2,         # length of short random subproduct
            mlimit,       # number of random elements to be tested
            ranelement,   # random element of S
            T,  p,
            residue;      # remainder of ranelement after factorization

    k := 0;
    while k < param[3] do
        k := k+1;
        len := Length(S.aux);
        #in case of large S.aux, form random subproducts
        #otherwise, try all of them
        if len > 2*param[3] then
            ran1 := SCRRandomPerm(len);
            ran2 := SCRRandomPerm(len);
            len2 := Random(1, QuoInt(len,2));
            string := SCRRandomString(len+len2);

            # we form two random subproducts:
            # w[1] in a random ordering of all generators
            # w[2] in a random ordering of a random subset of them
            w:=[];
            w[1] := S.identity;
            for x in [1 .. len] do
                w[1] := w[1]*(S.aux[x^ran1]^string[x]);
            od;
            w[2] := S.identity;
            for x in [1 .. len2] do
                w[2] := w[2]*(S.aux[x^ran2]^string[x+len]);
            od;

        else
            # take next two generators of S (unless only one is left)
            w := [];
            w[1] := S.aux[2*k-1];
            if len > 2*k-1 then
                w[2] := S.aux[2*k];
            else
                w[2] := S.identity;
            fi;
        fi;

        # take random elements of S
        m := 0;
        mlimit := param[4]*S.diam;
        while m < mlimit do
            m:=m+1;
            ranelement := S.identity;
            T := S;
            while Length( T.generators ) <> 0  do
                p := Random( T.orbit );
                while p <> T.orbit[ 1 ]  do
                    ranelement := LeftQuotient( T.transversal[ p ],
                                          ranelement );
                    p := p ^ T.transversal[ p ];
                od;
                T := T.stabilizer;
            od;
            i := 0;
            while i < 2 do
                i := i+1;
                if w[i] <> S.identity then
                    # test whether product of ranelement and w[i] in S
                    ranelement := ranelement*w[i];
                    residue := SCRSift(S,ranelement);
                    if residue <> S.identity then
                        return residue;
                    fi;
                fi;
            od; #i loop
        od; #m loop
        if len <= 2*k then
            #finished checking all in S.aux
            k := param[3];
        fi;
    od; #k loop

    return S.identity;
end );


#############################################################################
##
#F  SCRNotice( <orb>, <transversal>, <genlist> )  . . . . . . . . . . . local
##
##  checks whether   orbit  is closed  for  the   action of  permutations  in
##  genlist. If not, returns orbit point and generator witnessing.
##
InstallGlobalFunction( SCRNotice,
    function ( orb, transversal, genlist )
    local flag, #first component of output; true if orb is closed for
                #action of genlist
          i,    #second component of output, index of point in orb moving out
          j ;   #third component, index of permutation in genlist moving orb[i]

    i := 0;
    flag := true;
    while i < Length(orb) and flag  do
        i := i+1;
        j := 0;
        while j < Length(genlist) and flag do
            j := j+1;
            if not IsBound(transversal[orb[i]^genlist[j]])  then
                flag := false;
            fi;
        od;
    od;

    return [flag,i,j];
end );


#############################################################################
##
#F  SCRExtend( <list> ) . . . . . . . . . . . . . . . . . . . . . . . . local
##
##  given a partial Schreier tree of depth d,
##  SCRExtends the partial Schreier tree to depth d+1
##  input, output coded in list of length 5
##
InstallGlobalFunction( SCRExtend, function ( list )
    local orb,          #partial orbit
          transversal,  #partial transversal
          treegen,      #list of generators
          treegeninv,   #inverses of elements of treegen
                        #both treegen, treegeninv are used in transversal
          i, j,         #loop variables
          previous,     #index, showing end of level d-1 in orb
          len;          #length of orb at entering routine

    orb:=list[1];
    transversal:=list[2];
    treegen:=list[3];
    treegeninv:=list[4];
    previous:=list[5];
    len:=Length(orb);

    # for each point on level d, check whether one of the generators or
    # inverses moves it out of orb. If yes, add image to orb
    for i in [previous+1..len] do
        for j in [1..Length(treegen)] do
            if not IsBound(transversal[orb[i]^treegen[j]]) then
                transversal[orb[i]^treegen[j]] := treegeninv[j];
                Add(orb, orb[i]^treegen[j]);
            fi;
            if not IsBound(transversal[orb[i]^treegeninv[j]]) then
                transversal[orb[i]^treegeninv[j]] := treegen[j];
                Add(orb, orb[i]^treegeninv[j]);
            fi;
        od;
    od;

    # return Schreier tree of depth one larger
    return [orb,transversal,treegen,treegeninv,len];

end );


#############################################################################
##
#F  SCRSchTree( <S>, <new> )  . . . . . . . . . . . . . . . . . . . . . local
##
##  creates Schreier tree for the group generated by S.generators \cup new
##
InstallGlobalFunction( SCRSchTree, function ( S, new )
    local l,        #output of notice
          flag,     #first component of output
          i,        #second component of output
          j,        #third component of output
          word,     #list of permutations coding a coset representative
          g,        #the coset representative coded by word
          witness,  #a permutation moving a point out of S.orbit
          list;     #list coding input and output of 'extend'

    l := SCRNotice(S.orbit,S.transversal,new);
    flag := l[1];
    if flag then
        #do nothing; the orbit did not change
        return;
    else
        i := l[2];
        j := l[3];
        witness := new[j];
    fi;

    while not flag do
        word := CosetRepAsWord(S.orbit[1],S.orbit[i],S.transversal);
        g := Product(word);
        #add a new generator to treegen which moves S.orbit[1] out of S.orbit
        Add(S.treegen, g^(-1)*witness);
        Add(S.treegeninv, witness^(-1)*g);

        #recompute Schreier tree to new depth
        S.orbit := [S.orbit[1]];
        S.transversal := [];
        S.transversal[S.orbit[1]] := S.identity;
        S.treedepth := 0;
        list := [S.orbit,S.transversal,S.treegen,S.treegeninv,0];
        flag := false;

        #with k generators, we build only a tree of depth 2k
        while not flag and S.treedepth < 2*Length(S.treegen) do
            list := SCRExtend(list);
            S.orbit := list[1];
            S.transversal := list[2];
            if Length(S.orbit) = list[5] then
                #the tree did not extend; orbit is closed for treegen
                flag := true;
            else
                S.treedepth := S.treedepth + 1;
            fi;
        od;
        #increased S.orbit may not be closed for all generators of S
        l := SCRNotice(S.orbit,S.transversal,S.generators);
        flag := l[1];
        if not flag then
            i := l[2];
            j := l[3];
            witness := S.generators[j];
        fi;
    od;

    #update record components aux, diam
    S.aux := Concatenation(S.treegen,S.treegeninv,S.stabilizer.aux);
    S.diam := S.treedepth+S.stabilizer.diam;

end );


#############################################################################
##
#F  SCRRandomPerm( <d> )  . . . . . . . . . . . . . . . . . . . . . . . local
##
##  constructs random permutation in Sym(d)
##  without creating group record of Sym(d)
##
InstallGlobalFunction( SCRRandomPerm, function ( d )
    local   rnd,        # random permutation, result
            tmp,        # temporary variable for swapping
            i,  k;      # loop variables

    # use Floyd\'s algorithm
    rnd := [ 1 .. d ];
    for i  in [ 1 .. d-1 ]  do
        k := Random( 1, d+1-i );
        tmp := rnd[d+1-i];
        rnd[d+1-i] := rnd[k];
        rnd[k] := tmp;
    od;

    # return the permutation
    return PermList( rnd );
end );


#############################################################################
##
#F  SCRRandomString( <n> )  . . . . . . . . . . . . . . . . . . . . . . local
##
##  constructs random 0-1 string of length n
##  same steps as Random, but uses created random number for 28 bits
##
InstallGlobalFunction( SCRRandomString, function ( n )
    local i, j,     # loop variables
          k,        # number of 28 long substrings
          rnd,      # the random number which would be created by Random
          string,   # the random string constructed
          range;    # Upper value of range used for getting ints

    range := 2^28-1;

    string:=[];
    k:=QuoInt(n-1,28);
    for i in [0..k-1] do
        rnd := Random(0,range);
        # use each bit of rnd
        for j in [1 .. 28] do
            string[28*i+j] := rnd mod 2;
            rnd := QuoInt(rnd,2);
        od;
    od;

    # construct last <= 28 bits  of string
    rnd := Random(0, range);
    for j in [28*k+1 .. n] do
        string[j] := rnd mod 2;
        rnd := QuoInt(rnd,2);
    od;

    return string;
end );

#############################################################################
##
#F  SCRRandomSubproduct( <list>, <id> ) .  random subproduct of perms in list
##
InstallGlobalFunction( SCRRandomSubproduct, function( list, id )
    local string,  # 0-1 string containing the exponents of elements of list
          random,  # the random subproduct
          i;       # loop variable

    string := SCRRandomString(Length(list));
    random := id;
    for i in [1 .. Length(list)] do
        if string[i] = 1 then
            random := random*list[i];
        fi;
    od;

    return random;
end );

#############################################################################
##
#F  SCRExtendRecord( <G> )  . . . . . . . . . . . . . . . . . . . . . . local
##
##  defines record elements used at random stabilizer chain construction
##
InstallGlobalFunction( SCRExtendRecord, function(G)
    local list,       # list of stabilizer subgroups
          len,        # length of list
          i;          # loop variable

    list := ListStabChain(G);
    len := Length(list);
    list[len].diam := 0;
    list[len].aux := [];
    for i in [1 .. len - 1] do
        # list[len - i].real := list[len - i].generators;

        if Length(list[len - i].orbit) = 1 then
            # in this case, SCRSchTree will not do anything;
            # we have to define records treedepth, aux, and diameter
            list[len - i].treedepth := 0;
            list[len - i].aux := list[len - i + 1].aux;
            list[len - i].diam := list[len - i + 1].diam;
        fi;

        list[len - i].orbit := [ list[len - i].orbit[1] ];
        list[len - i].transversal := [];
        list[len - i].transversal[list[len - i].orbit[1]]
            := list[len - i].identity;
        list[len - i].treegen := [];
        list[len - i].treegeninv := [];
        SCRSchTree( list[len - i], list[len - i].generators );

    od;

end );


#############################################################################
##
#F  SCRRestoredRecord( <G> )  . . . . . . . . . . . . . . . . . . . . . local
##
##  restores usual group records after random stabilizer chain construction
##
InstallGlobalFunction( SCRRestoredRecord, function( G )
    local   sgs,  T,  S,  l,  pnt,o,ind;

    S := G;
    sgs := [ S.identity ];
    while IsBound( S.stabilizer )  do
        UniteSet( sgs, S.treegen );
        UniteSet( sgs, S.treegeninv );
        S := S.stabilizer;
    od;
    T := EmptyStabChain( sgs, G.identity );
    sgs := [ 2 .. Length( sgs ) ];
    S := T;
    while IsBound( G.stabilizer )  do
        InsertTrivialStabilizer( S, G.orbit[ 1 ] );
        S.genlabels   := sgs;
        S.generators  := G.generators;
        S.orbit       := G.orbit;
        S.transversal := G.transversal;
        o:=ShallowCopy(S.orbit);
        # check identity of transversal elements first: Most are
        # identical and thus element comparisons are relatively
        # infrequent
        while Length(o)>0 do

          ind:=S.transversal[o[1]];
          ind:=Filtered([1..Length(o)],
                  i->IsIdenticalObj(S.transversal[o[i]],ind));
          for l in sgs do
            if S.transversal[o[1]]=S.labels[l] then
              for pnt in o{ind} do # all these transv. elements are same
                S.translabels[ pnt ] := l;
              od;
            fi;
          od;
          o:=o{Difference([1..Length(o)],ind)}; # the rest
        od;

            #was:
            # (The element comparisons could be expensive)
            #for pnt  in S.orbit  do
            #    if S.transversal[ pnt ] = lab  then
            #        S.translabels[ pnt ] := l;
            #    fi;
            #od;

        sgs := Filtered( sgs, l ->
                       S.orbit[ 1 ] ^ S.labels[ l ] = S.orbit[ 1 ] );
        S := S.stabilizer;
        G := G.stabilizer;
    od;
    return T;
end );

#############################################################################
##
#F  VerifyStabilizer( <S>, <z>, <missing>, <correct> )  . . . .  verification
##
InstallGlobalFunction( VerifyStabilizer, function(S,z,missing,correct)
 #z is an involution, moving first base point
 # correct is boolean telling whether a base is known
 # if yes, missing contains the base points which do not occur in the base of S
 local pt1,            # the first base point
       zinv,           # inverse of z
       pt2,            # pt1^zinv
       flag,           # Boolean; true if z exchanges pt1 and pt2
       stab,           # the stabilizer of pt1 in S
       chain,          # stab after base change, to bring pt2 first
       stabpt2,        # stabilizer of pt2 in chain
       result,         # output, witness perm if something is wrong
       g,              # generator of stabpt2
       where1,         # stores which orbit of stab pts of S.orbit belong to
       orbit1count,    # index running through orbits of stab
       leaders,        # list of orbit representatives of stab
       i, j, l,        # loop variables
       residue,        # result of sifting as word
       where2,         # boolean list to mark orbits of stabpt2 as processed
       k,              # point of an orbit of stab
       gen,            # a generator of stab
       transversal,    # transversal of stab, on all of its orbits
       orb,            # an orbit of stabpt2
       img,            # image of a point at computation of orb
       pnt,            # a point from orb
       w1, w2, w3, w4, # words/permutations coding coset representatives
       w1inv,          # inverse of w1
       schgen;         # Schreier generator

    pt1 := S.orbit[1];
    zinv := z^(-1);
    pt2 := pt1^zinv;
    result := S.identity;
    stab := S.stabilizer;
    if pt2^zinv = pt1 then
       flag := true;
       result := SCRSift(stab, z^2);
    else
       flag := false;
    fi;

    # store which orbit of stab the pts of S.orbit belong to
    # in each orbit, compute transversals from a representative
    where1 := []; # orbits of stab
    leaders := [pt2];
    orbit1count := 1;
    transversal := [];
    transversal[pt2] := S.identity;
    where1[pt2] := 1;
    orb := [pt2];
    j := 1;
    while j <= Length( orb )  do
    for gen  in stab.generators  do
       k := orb[j] / gen;
       if not IsBound( transversal[k] )  then
             transversal[k] := gen;
             Add( orb, k );
             where1[k] := orbit1count;
          fi;
       od;
       j := j + 1;
    od;
    for i in S.orbit do
        if not IsBound(where1[i]) then
           orbit1count := orbit1count + 1;
           Add(leaders, i);
           orb := [i];
           where1[i] := orbit1count;
           transversal[i] := S.identity;
           j := 1;
           while j <= Length( orb )  do
             for gen  in stab.generators  do
                 k := orb[j] / gen;
                 if not IsBound( transversal[k] )  then
                    transversal[k] := gen;
                    Add( orb, k );
                    where1[k] := orbit1count;
                 fi;
             od;
             j := j + 1;
           od;
        fi;
    od;

   # check that conjugates of point stabilizers of stab are subgroups of stab
   chain:=StructuralCopy(stab);
   for j in [1..Length(leaders)] do
       if result = S.identity then
           i := leaders[Length(leaders)+1-j];
           ChangeStabChain( chain, [i], false );
           if i = pt2 then
              w1 := z;
              w1inv := zinv;
           else
              w1 := CosetRepAsWord(pt1, i, S.transversal);
              w1 := Product(w1);
              w1inv := w1^(-1);
           fi;
           for g in chain.stabilizer.generators do
               if result = S.identity then
                  if correct then
                     residue := SiftAsWord(stab, [w1inv,g,w1]);
                     if residue[2] <> 0 then
                        result := Product(residue[1]);
                     else
                        l := 0;
                        while ( l < Length(missing) ) and ( result = S.identity ) do
                           l:=l+1;
                              if ImageInWord(missing[l],residue[1])
                                  <> missing[l]  then
                                  # remainder is not in S
                                  result := Product(residue[1]);
                              fi;
                        od;
                     fi;
                  else
                     result := SCRSift(stab, w1inv*g*w1);
                  fi;
               fi;
           od;
       fi;
   od;

 if result = S.identity then
    stabpt2 := chain.stabilizer;
    # process orbits of stabpt2
    where2:= BlistList( [1..Length(where1)],[] ) ;
    for i in S.orbit do
        if result = S.identity and (not where2[i]) then
           orb := [i];
           where2[i] := true;
           for pnt  in orb  do
               for gen  in stabpt2.generators  do
                   img := pnt^gen;
                   if not where2[img]  then
                      Add( orb, img );
                      where2[img] := true;
                   fi;
               od;
           od;

           # if we hit a new orbit of stabpt2
           # and if z exchanges pt1 and pt2 then
           # mark the orbit
           if flag and not where2[i^z] then
              orb := [i^z];
              where2[i^z] := true;
              for pnt  in orb  do
                 for gen  in stabpt2.generators  do
                    img := pnt^gen;
                    if not where2[img]  then
                       Add( orb, img );
                       where2[img] := true;
                    fi;
                 od;
              od;
           fi;

           # compute Schreier generator either as a word or perm
           # if i is not a fixed point of z,
           # we have to compute a Schreier generator
           w1 := CosetRepAsWord(pt1, leaders[where1[i]], S.transversal);
           w2 := CosetRepAsWord(leaders[where1[i]], i, transversal);
           w3 := CosetRepAsWord(leaders[where1[i^z]], i^z, transversal);
           if where1[i] <> where1[i^z] then
               w4 := CosetRepAsWord(pt1,leaders[where1[i^z]], S.transversal);
           else
               w4 := w1;
           fi;
           schgen := (Product(w1))^(-1)*(Product(w2))^(-1);
           if correct then
               schgen := Concatenation([schgen,z],w3,w4);
           else
               schgen := schgen*z*Product(w3)*Product(w4);
           fi;

           # sift Schreier generator either as a word or perm
           if correct then
              residue := SiftAsWord(stab, schgen);
              if residue[2] <> 0 then
                 result := Product(residue[1]);
              else
                  l := 0;
                  while ( l < Length(missing) ) and ( result = S.identity ) do
                     l:=l+1;
                     if ImageInWord(missing[l],residue[1])
                               <> missing[l]  then
                          # remainder is not in S
                          result := Product(residue[1]);
                     fi;
                  od;
               fi;
           else
              result := SCRSift(stab, schgen);
           fi;
        fi;  # result = S.identity and not where2[i]
     od;
 fi;

 return result;

end );

#############################################################################
##
#F  VerifySGS( <S>, <missing>, <correct> )  . . . . . . . . . .  verification
##
InstallGlobalFunction( VerifySGS, function(S,missing,correct)
 # correct is boolean telling whether a base is known
 # if yes, missing contains the base points which do not occur in the base of S
 local n,         # degree of S
       list,      # list of subgroups in stabchain
       result,    # result of the test
       len,       # length of list
       i,l,       # loop variables
       residue,   # result of sifting as word
       temp,      # subgroup of S we currently work with
       temp2,     # temp, with possible extension on blocks
       gen,       # the generator whose addition we verify
       longer,    # generator list of temp, after gen is added
       gencount,  # counts the generators to be verified
       set,       # set of orbit of temp
       orbit,     # orbit of temp after adding gen
       blks,      # images of block when gen is added
       extension, # list of length 2; first coord. is temp2, extended by
                  # the action on blks, second coord. is newgen
       newgen,    # the extension of gen
       leader,    # first point in orbit of temp2
       block,     # block containing leader
       point,     # another point from block
       pos;       # position of set in blks

 n := LargestMovedPoint(S.generators);
 list := ListStabChain(S);
 len := Length(list);
 result := S.identity;

 # verify the subgroup chain from bottom up
 i := 1;
 while i < len  and  result = S.identity  do
     temp := ShallowCopy(list[len - i + 1]);
     InsertTrivialStabilizer(temp, list[len -i].orbit[1]);
     gencount := 0;

     # add one generator a time
     while (gencount < Length(list[len - i].generators)) and (result= S.identity ) do
         gencount := gencount + 1;
         gen := list[len - i].generators[gencount];
         set := Set( temp.orbit );

         # if adding gen did not change the fundamental orbit, just sift gen
         if set = OnSets(set,gen) then
            if correct then
               residue := SiftAsWord(temp, [gen]);
               if residue[2] <> 0 then
                  result := Product(residue[1]);
               else
                  l := 0;
                  while ( l < Length(missing) ) and ( result = S.identity ) do
                     l:=l+1;
                     if ImageInWord(missing[l],residue[1])
                               <> missing[l]  then
                          # remainder is not in S
                          result := Product(residue[1]);
                     fi;
                  od;
               fi;
            else
               result := SCRSift(temp, gen);
            fi;
         # if the fundamental orbit increased, compute block system
         # obtained when adding gen
         else
            if Length(set) =1 then
                temp2 := temp;
                newgen := gen;
            # otherwise, compute the action on blocks
            else
                longer := Concatenation(temp.generators, [gen]);
                orbit := OrbitPerms( longer, temp.orbit[ 1 ] );
                blks := Blocks( GroupByGenerators( longer ), orbit, set );
                if Length(blks) * Length(set) <> Length(orbit) then
                    result := "false0";
                else
                    pos := Position( blks, set );
                    extension := ExtensionOnBlocks( temp, n, blks, [gen] );
                    temp2 := extension[1];
                    newgen := extension[2][1];
                    InsertTrivialStabilizer(temp2, n + pos);
                fi;
            fi;

            # first generator in first group can be verified easily
            if i=1 and Length(set) =1 then
                result := newgen^(CycleLengthOp(newgen,temp2.orbit[1]));
                AddGeneratorsExtendSchreierTree(temp2, [newgen]);
            elif result = S.identity then
                AddGeneratorsExtendSchreierTree( temp2, [ newgen ] );
                blks := Blocks( GroupByGenerators(temp2.generators),
                                temp2.orbit);
                if Length(blks) > 1 then
                   leader := temp2.orbit[1];
                   block := First(blks, x -> leader in x);
                   point := First(block, x -> x <> leader);
                   newgen := Product(CosetRepAsWord(leader ,point,
                                       temp2.transversal));
                   temp2 := temp2.stabilizer;
                   InsertTrivialStabilizer(temp2, leader);
                   AddGeneratorsExtendSchreierTree(temp2, [newgen]);
                   result := VerifyStabilizer(temp2,newgen,missing,correct);
                   if leader > n then
                      AddGeneratorsExtendSchreierTree(temp,[RestrictedPermNC
                          (newgen, [1..n])]);
                   else
                      temp := temp2;
                   fi;
                   gencount := gencount - 1;
                else
                   result := VerifyStabilizer(temp2,newgen,missing,correct);
                   if temp2.orbit[1] > n then
                      AddGeneratorsExtendSchreierTree(temp,[gen]);
                   fi;
                fi;
                if result <> S.identity and temp2.orbit[1] > n then
                     result := RestrictedPermNC(result, [1..n]);
                fi;

            fi;
        fi;
    od;
    i := i + 1;
 od;

 return result;

end );

#############################################################################
##
#F  ExtensionOnBlocks( <S>, <n>, <blks>, <elms> ) . . . . . . . . . extension
##
InstallGlobalFunction( ExtensionOnBlocks, function( S, n, blks, elms )
    local   where,  j,  k,  hom,  T,  newelms;

      # list which block the elements of the orbit belong to
      where := [];
      for j in [1..Length(blks)] do
          for k in blks[j] do
              where[k] := j;
          od;
      od;

      hom := function( g )
          local  perm,  j;

          perm := [1..n];
          for j in [1..Length(blks)] do
              perm[n+j] := n+where[blks[j][1]^g];
          od;
          perm := PermList(perm);
          return g * perm;
      end;

      T := EmptyStabChain( [  ], S.identity );
      ConjugateStabChain( S, T, hom, S.identity );

      # construct extensions of permutations in elms
      newelms := List( elms, hom );

      return [ T, newelms ];
end );

#############################################################################
##
#F  ClosureRandomPermGroup( <G>, <genlist>, <options> ) make closure randomly
##
InstallGlobalFunction( ClosureRandomPermGroup,
    function( G, genlist, options )
    local  k,          # number of pairs of subproducts of generators in
                       # testing result
           givenbase,  # ordering from which initial base points should
                       # be chosen
           gens,       # generators in genlist that are not in <G>
           g,          # element of gens
           degree,     # degree of closure
           orbits,     # list of orbits of closure
           orbits2,    # list of orbits of closure
           i,j,        # loop variables
           param,       # list of parameters guiding number of repetitions
                        # in random constructions
           where,       # list indicating which orbit contains points in domain
           basesize,    # list; i^th entry = number of base points in orbits[i]
           ready,       # boolean; true if stabilizer chain ready
           new,         # list of permutations to be added to stab. chain
           result,      # output of checking phase; nontrivial if stabilizer
                        # chain is incorrect
           base,        # ordering of domain from which base points are taken
           missing,     # if a correct base was provided by input, missing
                       # contains those points of it which are not in
                       # constructed base
           cnt,        # iteration counter
           correct;     # boolean; true if a correct base is given

# warning:  options.base should be compatible with BaseOfGroup(G)

    gens := Filtered( genlist, gen -> not(IsOne(SCRSift(G,gen))) );
    if Length(gens) > 0  then

        G.identity := One(gens[1]) ;
        if options.random = 1000 then
            #case of deterministic computation with known size
            k := 1;
        else
            k:=First([1..14],x->(3/5)^x<1-options.random/1000);
        fi;
        if IsBound(options.knownBase) then
            param := [k,4,0,0,0,0];
        else
            param := [QuoInt(k,2),4,QuoInt(k+1,2),4,50,5];
        fi;
        if options.random <= 200 then
            param[2] := 2;
            param[4] := 2;
        fi;

#param[1] = number of pairs of random subproducts from generators in
#           first checking phase
#param[2] = (number of random elements from created set)/S.diam
#           in first checking phase
#param[3] = number of pairs of random subproducts from generators in
#           second checking phase
#param[4] = (number of random elements from created set)/S.diam
#           in second checking phase
#param[5] = maximum size of orbits in  which we evaluate words on all
#           points of orbit
#param[6] = minimum number of random points from orbit to plug in to check
#           whether given word is identity on orbit

        degree := LargestMovedPoint( Union( G.generators, gens ) );

        # prepare input of construction
        if IsBound(options.base) then
            givenbase := options.base;
        else
            givenbase := [];
        fi;

        if IsBound(options.knownBase) then
            correct := true;
        else
            correct := false;
        fi;

        if correct then
            # if correct  base was given as input,
            # no need for orbit information
            base := Set( givenbase );
            for i in BaseStabChain(G) do
                if not i in base then
                    Add( givenbase, i );
                fi;
            od;
            base := Concatenation(givenbase,Difference(options.knownBase,
                                                       givenbase));
            missing := Difference(options.knownBase,BaseStabChain(G));
            basesize := [];
            where := [];
            orbits := [];
        else
            # create ordering of domain used in choosing base points and
            # compute orbit information
            base := Set( givenbase );
            for i in BaseStabChain(G) do
                if not i in base then
                    Add( givenbase, i );
                fi;
            od;
            base := Concatenation(givenbase,Difference([1..degree],givenbase));
            missing := [];
            orbits2 := OrbitsPerms( Union( G.generators, gens ), [1..degree] );
            #throw away one-element orbits
            orbits:=[];
            j:=0;
            for i in [1..Length(orbits2)] do
                if Length(orbits2[i]) >1 then
                    j:=j+1; orbits[j]:= orbits2[i];
                fi;
            od;
            basesize:=[];
            where:=[];
            for i in [1..Length(orbits)] do
                basesize[i]:=0;
                for j in [1..Length(orbits[i])] do
                    where[orbits[i][j]]:=i;
                od;
            od;
            # temporary solution to speed up of handling
            # of lots of small orbits until compiler
            if Length(orbits) > degree/40 then
                param[1] := 0;
                param[3] := k;
            fi;
        fi;

        if not IsBound(G.aux) then
            SCRExtendRecord(G);
        fi;
        new := gens;

        #the first call of SCRMakeStabStrong has top:=false
        #in order to add gens to the generating set of G;
        #further calls have top:=true, in order not to add
        #output of SCRStrongGenTest to generating set.
        #remark: adding gens to the generating set of G before
        #calling SCRMakeStabStrong gives a nasty error if first base
        #point changes
        for g in gens do
            if not(IsOne(SCRSift(G,g))) then
                SCRMakeStabStrong (G,[g],param,orbits,
                        where,basesize,base,correct,missing,false);
            fi;
        od;

        cnt:=0;
        ready := false;
        while not ready do
          if    IsBound(options.limit)
            and SizeStabChain(G)=options.limit
          then
              ready := true;
          else
              # we do a little random testing, to ensure heuristically a
              # correct result
              if correct then
                  result := SCRStrongGenTest(G,[1,10/G.diam,0,0,0,0],orbits,
                                      basesize,base,correct,missing);
              else
                  result := SCRStrongGenTest2(G,[0,0,1,10/G.diam,0,0]);
              fi;
              if not(IsPerm(result) and IsOne(result)) then
                  new := [result];
                  ready := false;
              elif options.random = 1000 then
                  G.restored := SCRRestoredRecord(G);
                  result := VerifySGS( G.restored, missing, correct );
                  cnt:=cnt+1;
                  if cnt>99 then
                    # in rare cases this loop iterates for a very long time.
                    # In this case, rather create a new chain, than try to
                    # fix the problematic one
                    #Error("infinite loop?");
                    return StabChainRandomPermGroup(G.generators,G.identity,
                            options);
                  fi;
              elif options.random > 0 then
                  result := SCRStrongGenTest
                          (G,param,orbits,basesize,base,correct,missing);
              fi;
              if not(IsPerm(result) and IsOne(result)) then
                  if not IsPerm(result) then
                     repeat
                         result := SCRStrongGenTest2(G,[0,0,1,10/G.diam,0,0]);
                     until not(IsPerm(result) and IsOne(result));
                  fi;
                  new := [result];
                  ready := false;
              elif correct or options.random = 0 or options.random = 1000 then
                  ready := true;
              else
                  result := SCRStrongGenTest2(G,param);
                  if IsPerm(result) and IsOne(result) then
                      ready := true;
                  else
                      new := [result];
                      ready := false;
                  fi;
              fi;
              if not ready then
                Unbind(G.restored);
                SCRMakeStabStrong (G,new,param,orbits,
                        where,basesize,base,correct,missing,true);
                #Print("D ",SizeStabChain(G),"\n");
              fi;
          fi;
        od;

        if not IsBound(options.temp) or options.temp = false then
             if IsBound( G.restored ) then
                G := G.restored;
             else
                G := SCRRestoredRecord(G);
             fi;
        else
             G.basesize := basesize;
             G.correct := correct;
             G.orbits := orbits;
             G.missing := missing;
             G.base := base;
        fi;

    fi; # if Length(gens) > 0


    # return the closure
    return G;

end );

#############################################################################