#############################################################################
##
#F  PermGroupOps.Blocks( <G>, <D> [, <seed>] [, <operation>] )
##
#PermGroupOps.BlocksNoSeed := function ( G, D )

# Adaptation of the code for BlocksNoSeed to return _all_ minimal blocks 
# containing D[1].
# Graham Sharp (Oxford)
# August 1997
# Revised: Sept 23 1997 for GAP 4.
# (NB will be included in GAP 4b2 library, expected release date Oct 1997)

MinimalBlocks := function(G, D)
# <D> should be a list of positive integers, <G> should act by OnPoints
    local   blocks,     # block system of <G>, result
            orbit,      # orbit of 1 under <G>
            trans,      # factored inverse transversal for <orbit>
            eql,        # '<i> = <eql>[<k>]' means $\beta(i)  = \beta(k)$,
            next,       # the points that are equivalent are linked
            last,       # last point on the list linked through 'next'
            leq,        # '<i> = <leq>[<k>]' means $\beta(i) <= \beta(k)$
            gen,        # one generator of <G> or 'Stab(<G>,1)'
            rnd,        # random element of <G>
            pnt,        # one point in an orbit
            img,        # the image of <pnt> under <gen>
            cur,        # the current representative of an orbit
            rep,        # the representative of a block in the block system
            block,      # the block, result
            changed,    # number of random Schreier generators
            nrorbs,     # number of orbits of subgroup $H$ of $G_1$
            i,          # loop variable
            minblocks,  # set of minimal blocks, result
            poss,       # flag to indicate whether we might have a block
            iter,       # which points we've checked when
            start;      # index of first cur for this iteration (non-dec)

    # handle trivial domain
    if Length( D ) = 1  or IsPrime( Length( D ) )  then
        return [ D ];
    fi;
 
    # handle trivial group
    if IsTrivial( G )  then
        Error("<G> must operate transitively on <D>");
    fi;

    # compute the orbit of $G$ and a factored transversal
    orbit := [ D[1] ];
    trans := [];
    trans[ D[1] ] := ();
    for pnt  in orbit  do
        for gen  in GeneratorsOfGroup(G)  do
            if not IsBound( trans[ pnt / gen ] )  then
                Add( orbit, pnt / gen );
                trans[ pnt / gen ] := gen;
            fi;
        od;
    od;

    # check that the group is transitive
    if Length( orbit ) <> Length( D )  then
        Error("<G> must operate transitively on <D>");
    fi;
    #InfoOperation1("#I  BlocksNoSeed transversal computed.\n");
    nrorbs := Length( orbit );

    # since $i \in k^{G_1}$ implies $\beta(i)=\beta(k)$,  we initialize <eql>
    # so that the connected components are orbits of some subgroup  $H < G_1$
    eql := [];
    leq := [];
    next := [];
    last := [];
    iter := [];
    for pnt  in orbit  do
        eql[pnt]  := pnt;
        leq[pnt]  := pnt;
        next[pnt] := 0;
        last[pnt] := pnt;
        iter[pnt] := 0;
    od;

    # repeat until we run out of points
    minblocks := [];
    changed := 0;
    rnd := ();

    for start in orbit{[2..Length(D)]} do
  
      # repeat until we have a block system
      cur := start;
      # unless this is a new point, ignore and go on to the next
      # -we could do this by a linked list to avoid these checks but the 
      #  O(n) overheads involved in setting it up are greater than those saved
	if iter[cur] = 0 then
    
	repeat
    
	    # compute such an $H$ by taking random  Schreier generators  of $G_1$
	    # and stop if 2 successive generators dont change the orbits any more
	    while changed < 2  do

		# compute a random Schreier generator of $G_1$
		i := Length( orbit );
		while 1 <= i  do
		    rnd := rnd * Random( GeneratorsOfGroup(G));
		    i   := QuoInt( i, 2 );
		od;
		gen := rnd;
		while D[1] ^ gen <> D[1]  do
		    gen := gen * trans[ D[1] ^ gen ];
		od;
		changed := changed + 1;

		# compute the image of every point under <gen>
		for pnt  in orbit  do
		    img := pnt ^ gen;

		    # find the representative of the orbit of <pnt>
		    while eql[pnt] <> pnt  do
			pnt := eql[pnt];
		    od;

		    # find the representative of the orbit of <img>
		    while eql[img] <> img  do
			img := eql[img];
		    od;

		    # if the don't agree merge their orbits
		    if   pnt < img  then
			eql[img] := pnt;
			next[ last[pnt] ] := img;
			last[pnt] := last[img];
			nrorbs := nrorbs - 1;
			changed := 0;
		    elif img < pnt  then
			eql[pnt] := img;
			next[ last[img] ] := pnt;
			last[img] := last[pnt];
			nrorbs := nrorbs - 1;
			changed := 0;
		    fi;

		od;

	    od;
	    #InfoOperation1("#I  BlocksNoSeed ",
		#	   "number of orbits of <H> < <G>_1 is ",nrorbs,"\n");

	    # take arbitrary point <cur>,  and an element <gen> taking 1 to <cur>
	    while eql[cur] <> cur  do
		cur := eql[cur];
	    od;
	    # Mark the points in this new H-orbit as visited
            if iter[cur] <> start then
  	      img := cur;
	      while img <> 0 do
		iter[img] := start;
		img := next[img];
	      od;
            fi;
	    gen := [];
	    img := cur;
	    while img <> D[1]  do
		Add( gen, trans[img] );
		img := img ^ trans[img];
	    od;
	    gen := Reversed( gen );

	    # compute an alleged block as orbit of 1 under $< H, gen >$
	    pnt := cur;
	    poss := true;
	    while pnt <> 0  do

		# compute the representative of the block containing the image
		img := pnt;
		for i  in gen  do
		    img := img / i;
		od;
		while eql[img] <> img  do
		    img := eql[img];
		od;

		# if its not our current block but a new block
		if   img <> D[1]  and img <> cur and leq[img] = img
			and (iter[img] = 0 or iter[img] = start) then

		    # then try <img> as a new start
		    leq[cur] := img;
		    cur := img;
                    if iter[cur] <> start then
  		      img := cur;
		      while img <> 0 do
			iter[img] := start;
			img := next[img];
		      od;
                    fi;
		    gen := [];
		    img := cur;
		    while img <> D[1]  do
			Add( gen, trans[img] );
			img := img ^ trans[img];
		    od;
		    gen := Reversed( gen );
		    pnt := cur;

		# otherwise if its not our current block but contains it
		# by construction a nonminimal block contains the current block
		# - not any more it doesn't! Now we also have to check whether 
		# the block appeared this time or earlier. 
		elif img <> D[1]  and img <> cur  
			   and leq[img] <> img  and iter[img] = start then

		    # then merge all blocks it contains with <cur>
		    while img <> cur  do
			eql[img] := cur;
			next[ last[cur] ] := img;
			last[ cur ] := last[ img ];
			img := leq[img];
			while img <> eql[img]  do
			    img := eql[img];
			od;
		    od;
		    pnt := next[pnt];

		# else if the block appeared in a previous iteration
		elif iter[img] <> start and iter[img] <> 0 then

		   # then end this iteration as this is not a minimal block
		    pnt := 0;
		    poss := false;

		# otherwise go on to the next point in the orbit
		else

		    pnt := next[pnt];

		fi;

	    od;

	    # Skip this bit if we know we haven't got a block
	    if poss = true then
	    # make the alleged block
	    block := [ D[1] ];
	    pnt := cur;
	    while pnt <> 0  do
		Add( block, pnt );
		pnt := next[pnt];
	    od;
	    block := Set( block );
	    blocks := [ block ];
	    #InfoOperation1("#I  BlocksNoSeed ",
			#   "length of alleged block is ",Length(block),"\n");

	    # quick test to see if the group is primitive
	    if Length( block ) = Length( orbit )  then
		#InfoOperation1("#I  BlocksNoSeed <G> is primitive\n");
		return [ D ];
	    fi;

	    # quick test to see if the orbit can be a block
	    if Length( orbit ) mod Length( block ) <> 0  then
		#InfoOperation1("#I  BlocksNoSeed ",
			#       "alleged block is clearly not a block\n");
		changed := -1000;
	    fi;

	    # '<rep>[<i>]' is the representative of the block containing <i>
	    rep := [];
	    for pnt  in orbit  do
		rep[pnt] := 0;
	    od;
	    for pnt  in block  do
		rep[pnt] := 1;
	    od;

	    # compute the block system with an orbit algorithm
	    i := 1;
	    while 0 <= changed  and i <= Length( blocks )  do

		# loop over the generators
		for gen  in GeneratorsOfGroup(G) do

		    # compute the image of the block under the generator
		    img := OnSets( blocks[i], gen );

		    # if this block is new
		    if rep[ img[1] ] = 0  then

			# add the new block to the list of blocks
			Add( blocks, img );

			# check that all points in the image are new
			for pnt  in img  do
			    if rep[pnt] <> 0  then
				#InfoOperation1("#I  BlocksNoSeed ",
					#       "alleged block is not a block\n");
				changed := -1000;
			    fi;
			    rep[pnt] := img[1];
			od;

		    # if this block is old
		    else

			# check that all points in the image lie in the block
			for pnt  in img  do
			    if rep[pnt] <> rep[img[1]]  then
			       #InfoOperation1("#I  BlocksNoSeed ",
					#       "alleged block is not a block\n");
				changed := -1000;
			    fi;
			od;

		    fi;

		od;

		# on to the next block in the orbit
		i := i + 1;
	    od;
	    fi;

	until 0 <= changed;
	if poss = true then AddSet(minblocks, block); fi;
  
      # loop back to get another minimal block
      fi;
    od;

    # return the block system
    return minblocks;
end;