## Gap 4 script based on paper GRS1997b: Recognition of actions on pairs
## II: with applications to groups acting on graphs.
##
## pairs3.g : incorporating attempts to avoid repeating conjugate solutions

## orbitals1.g : rearranged to remove recursion, so reducing number of times
## minblks is called and facilitating use of stronger notion of equivalence
## Index2 removed and Index1 renamed - this is now a program to 
## recognise orbitals, rather than to recognise actions on pairs.
## orbitals2a.g: using equivalence of solutions to reduce output and effort
## orbitals3a.g: more on equivalence: elimination of duplicate output
##   subgroups from the same conjugacy class of N_G(G_w).
## orbitals4a.g: elimination of unneeded values of x from orbitals3a.g
## orbitals5a.g: corrected version of orbitals4a.g
## orbitalss-p1.g: adapted version of orbitals5a.g to find self-paired
##   solutions

OrbitalActionsSelfPaired := function(G, n)
                       # G input group acting on [1..n]

    local Gcos,        # Another name for G (historical reasons)
          Gw,          # Stabilizer of 1 in G
          GwOrbs,      # Orbits of Gw on [1..n]
          DtoOrbs,     # map from [1..n] to linked lists of orbits of Gw
          next,        # next entry in linked list of orbit of Gw
          last,        # last entry in linked list of orbit of Gw
          trans,       # factored inverse Schreier transversal
          N,           # elements of 1^{N_G(Gw)} (ie normalizer of Gw)
          Nlr,         # points to generate N as a block containing 1
          Npts,        # in making Nlr: contains points still to be added
          GwNorbs,     # orbits of N on GwOrbs (NB action is 'conjugation')
          GwNorbsUn,   # contains first elt of each unprocessed Gw orbit 
          reptoGwOrb,  # maps Gw orbit reps (from DtoOrbs) to the orbits
          L,           # solutions waiting for processing in size order
                       # L[i].blk       block corresponding to subgroup J
                       # L[i].blkstab   subset of block which generates it
                       #                  as a block containing 1
                       # L[i].xs        list of values of x for this soln
                       # L[i].conjstabs[j] is the conjugate of blkstab
                       #                  by (a rep of) N[j] 
          x,           # the x from some soln, as a point in [1..n]
          gx,          # inverse of rep of x as a word in the generators
          gxr,         # reverse of gx
          blk,         # the block in a solution; same under constn
          blkstab,     # points generating blk as block containing 1
          imgblk,      # image of blk under 'conjugation' by x
          output,      # the output of the whole function
          xs,          # set of elements x
          blks,        # block system of translates of blk 
          GG,          # isomorphic to action of G on blks
          blocks,      # the minimal blocks of GG
          b,           # one member of blocks
          preblk,      # pre-image of b in [1..n] under the natural map
          z,           # element of preblk but not blk
          preblkstab,  # subset of preblk that generates it as a block
          ind,         # index in L of soln. to compare with new one,
                       #   place in L to add new soln.
          imgpreblk,   # image of preblk under 'conjugation' by x
          orbit,       # orbit of G; used as trans is made
          blkstabgens, # list of inverse reps of elts of blkstab as words
          newxs,       # those elements of the old xs which we are keeping
          preblkstabgens, # same as blkstabgens only for preblkstab
          acc,         # flag to indicate whether a (blk,x) is acceptable
          found,       # says whether an equiv subgroup etc has been found
          ni,          # index for counting through N
          conjstabs,   # used to calculate L[ind].conjstabs
          npt,         # point in N
          blks1,       # block system of translates of L[ind].blk
          conjstab,    # one element of L[ind].conjstabs during construction
          gnpt,        # inverse rep of npt as word in generators
          nptinv,      # 1^gnpt (NB 1^{gnpt^{-1}}=npt)
          nis,         # set of ni s.t. N[ni] conjugates L[ind].blk to preblk
          pr,          # an element of a particular list of pairs
          blkx,        # translate fo blk or preblk by (rep of) x
    # Localised variables with repeated usage
          g,         # group element, usu. one of the generators
          gpnt, pntg,
          elt, 
          gen, 
          gens, 
          img, 
          pre, 
          pt,        # a point 
          orb,       # general name for an orbit
          rep,       # representative
          pnt,       # point
          pntgx, 
          gpre, pnt2, img2, rep3, img3,
          used,      # flags showing membership of orbit etc as it's made
          i, j;      # index variables 

    Gcos := G;
    # Find inverse factored transversal
    orbit  := [1];
    trans := [];
    trans[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;

    # Find orbits of Gw = Stab(G, 1) and store them as linked lists;
    # set up conversion functions DtoOrbs, reptoGwOrb
    Gw := Stabilizer(G, 1);
    GwOrbs := Orbits(Gw, [1..n]);
    DtoOrbs := [];
    next := [];
    last := [];
    reptoGwOrb := [];
    for orb in GwOrbs do
      rep := orb[1];
      reptoGwOrb[rep] := orb;
      last[rep] := rep;
      for pnt in orb do
        DtoOrbs[pnt] := rep;
        next[last[rep]] := pnt;
        last[rep] := pnt;
      od;
      next[last[rep]] := 0;
    od;

    # Find normalizer <N> of <Gw> (as list of generating points).
    # Npts starts off containing all the points of 1^N(Gw)
    # N ends up containing what Npts started off with
    # Nlr is a smaller set which generates N as a block
    N := [1];
    used := [];
    used[1] := 0;
    Nlr := [];
    Npts := Set(Flat(Filtered(GwOrbs, x -> Length(x) = 1)));
    RemoveSet(Npts, 1);
    while Npts <> [] do
      Add(Nlr, Npts[1]);
      # each time we add a point to Nlr, do an orbit calc to close
      # under the block generated by Nlr. Then Nlr only has size \log\#N  
      for pnt in N do 
        for x in Nlr do
          img := x;
          img2 := pnt;
          while img <> 1 do
            img2 := img2 ^ trans[img];
            img := img ^ trans[img];
          od;
          if not IsBound(used[img2]) then
            Add(N, img2);
            RemoveSet(Npts, img2);
            used[img2] := 0;
          fi;
        od;
      od;
    od;
    N := Set(N);

    # N(Gw)/Gw acts on GwOrbs by conjugation of the double cosets
    # Form the orbits under this action. Each orbit is represented by 
    # one element, the head of its linked list.
    GwNorbsUn := Set(List(GwOrbs, x -> DtoOrbs[x[1]]));
    used := [];
    GwNorbs := [];
    # for each orbit of the action of N on GwOrbs...
    while GwNorbsUn <> [] do;
      orb := [GwNorbsUn[1]];
      used[orb[1]] := 0;
      RemoveSet(GwNorbsUn, orb[1]); 
      # ...perform an orbit calculation
      for pnt in orb do
        # set gpnt to the reverse of a word in the gens representing
        # an inverse representative for pnt.
        gpnt := [];
        img := pnt;
        while img <> 1 do
          Add(gpnt, trans[img]);
          img := img ^ trans[img];
        od;
        gpnt := Reversed(gpnt);
        for x in Nlr do
          # do the conjugation: we want 1^{g^-1}pg where 1^p=pnt and 1^g=x
          # note that x and gx are used differently from the description
          img2 := 1;
          gx := [];
          img := x;
          while img <> 1 do
            Add(gx, trans[img]);
            img2 := img2 ^trans[img];
            img := img ^ trans[img];
          od;
          gx := Reversed(gx);
          for g in gpnt do
            img2 := img2 / g;
          od;
          for g in gx do
            img2 := img2 / g;
          od;
          img2 := DtoOrbs[img2];
          if not IsBound(used[img2]) then 
            Add(orb, img2);
            used[img2] := 0;
            RemoveSet(GwNorbsUn, img2);
          fi;
        od;
      od;
      Add(GwNorbs, Concatenation(List(orb, x -> reptoGwOrb[x])));
    od;


    Info(InfoOperation, 1, "OrbitalActions: transversal and suborbits found");

    # L is the main data structure storing pending solutions
    L := [];

    # Take each possible value of x, form block J = <Gw,Gw^(x^-1)>
    for orb in GwNorbs do
      x := DtoOrbs[orb[1]];
      # Find representative of x
      gx := [];
      pnt := x;
      while pnt <> 1 do
        Add(gx, trans[pnt]);
        pnt := pnt ^ trans[pnt];
      od;
      gxr := Reversed(gx);

      blk := [ 1 ];
      blkstab := [];
      used := [ ];
      used[ 1 ] := 0;
      for pnt in blk do
        # Want to translate GwOrbs by gx^-1, find the translated orbit
        # containing pnt and ensure whole orbit is in blk; however each 
        # time we add a point (pre) to blk, we add it to blkstab and close
        # blk under being a blk: thus blkstab is a small set which 
        # generates blkas a block containing 1.
        # Equivalently to translating GwOrbs, we translate pnt, find 
        # orbit, translate orbit back.
        rep := pnt;
        for g in gxr do
          rep := rep / g;
        od;
        rep := DtoOrbs[ rep ];
        img := rep;
        while img <> 0 do
          pre := img;
          for g in gx do
            pre := pre ^ g;
          od;
          if not IsBound( used[ pre] ) then
            # Instead of adding the new point <pre>, we find a rep for <pre>
            # so <pre>=1^<gpre> and then close under action of <gpre> and 
            # under <Gw>-orbits

            Add(blkstab, pre);
            # NB elements of <gpre> will be the inverses of the ones we want
            gpre := [];
            pnt2 := pre;
            while pnt2 <> 1 do
              Add(gpre, trans[pnt2]);
              pnt2 := pnt2 ^ trans[pnt2];
            od;
            gpre := Reversed(gpre);
            for pnt2 in blk do
              # set <img2> = <pnt2>^<gpre>
              img2 := pnt2;
              for g in gpre do
                img2 := img2 / g;
              od;
              if not IsBound( used[ img2 ] ) then 
                # Add the whole <Gw>-orbit containing <img2>
                rep3 := DtoOrbs[ img2 ];
                img3 := rep3;
                while img3 <> 0 do
                  if not IsBound(used[img3]) then
                     Add(blk, img3);
                     used[img3] := 0;
                  fi;
                  img3 := next[img3];
                od;
              fi;
            od;
          fi;
          img := next[img];
        od;
      od;
      # blk should now be a block with the property that if J = stab(blk)
      # so blk = 1^J then J \meet J^x contains Gw.

      blk := Immutable(Set(blk));

      # check whether block is self-paired.
      # form translate of blk by (rep of) x
      blkx := blk;
      for g in gxr do
        blkx := List(blkx, x -> x / g);
      od;
      blkx := Set(blkx);
      # Only need look at points which normalize Gw
      blkx := Intersection(blkx, N);

      acc := false;
      for pnt in blkx do
        if not acc then
	  # calculate image of 1 under square of rep of pnt
	  img := pnt;
	  pntg := [];
	  while img <> 1 do
	    Add(pntg, trans[img]);
	    img := img ^ trans[img];
	  od;
	  img := pnt;
	  for g in Reversed(pntg) do
	    img := img / g;
	  od;
	  if img = 1 then acc := true; fi;
        fi;
      od;

      if acc then
	# Find representatives of the generating points in blkstab
	# (blkstab generates blk as a block containing 1)
	blkstabgens := [];
	for pnt in blkstab do
	  img := pnt;
	  gen := [];
	  while img <> 1 do
	    Add(gen, trans[img]);
	    img := img ^ trans[img];
	  od;
	  Add(blkstabgens, Reversed(gen));
	od;
  
	# Calculate image of blk under 'conjugation' by <x> (actually
	# conjugation of the stabilizer by a representative <gxr> of x)
  
	used := [];
	pnt := 1;
	for g in gx do
	  pnt := pnt ^ g;
	od;
  
	imgblk := [pnt];
	used[pnt] := 0;
      fi;    
      i := 1;
      # we want the orbit closed under <Gw>^<x> and under the conjugates
      # of the elements of blkstabgens by <x>. imgblk is actually the 
      # translate of what we want by <x>^-1, so it starts off containing
      # 1^{x^-1} and gets closed under Gw and under blkstabgens
      while acc and i <= Length(imgblk) do
        pnt := imgblk[i];
        # Close under the generators
        j := 1;
        while acc and j <= Length(blkstabgens) do
          gen := blkstabgens[j];
          # form img = pnt ^ ( gen )
          img := pnt;
          for g in gen do
            img := img / g;
          od;
          if not IsBound( used[img] ) then
            # each new element of imgblk gets translated by x and compared
            # with blk
            img2 := img;
            for g in gxr do
              img2 := img2 / g;
            od;
            if img2 <> 1 and img2 in blk then
              acc := false;
            else
              Add(imgblk, img);
              used[img] := 0;
            fi;
          fi;
          j := j + 1;
        od;
        # Close under <Gw>
        img := DtoOrbs[pnt];
        # Avoid doing this suborbit if we've done it before
        if not IsBound(used[img]) or used[img] <> -1 then 
          while acc and img <> 0 do
            if not IsBound( used[img] ) then
              # again new elements need translating and testing
              img2 := img;
              for g in gxr do
                img2 := img2 / g;
              od;
              if img2 <> 1 and img2 in blk then
                acc := false;
              else
                Add(imgblk, img);
                used[img] := 0;
              fi;
            fi;
            img := next[img];
          od;
          # mark that we've done this suborbit
          used[ DtoOrbs[ pnt ] ] := -1;
        fi;           
        i := i + 1;
      od;

      if acc then
        # J\cap J^x= Gw. We now update the data structure L
        # Check to see if it's a conjugate of anything we've seen before.
        # Find index of first entry in L of same size or larger.
        ind := PositionSorted(L, rec(blk := blk),
                function(x, y) return (Size(x.blk) < Size(y.blk))
                  ; end );
        found := false;
        # Check each entry in L of same size to see if they are N-conjugate
        while not found and IsBound(L[ind]) and 
                     Size(L[ind].blk) = Size(blk) do
	  nis := [];
	  for ni in [1..Length(N)] do
            # L[ind].conjstabs[ni] contains generators for L[ind].blk
            # conjugated by N[ni] (as N[ni] normalizes Gw, we can get away
            # with just conjugating the generating points of the block
            # and ignoring Gw, which we couldn't when conjugating by x) 
            if IsSubsetSet(blk, L[ind].conjstabs[ni]) then
	      Add(nis, ni);
              found := true;
            fi;
          od;
          if not found then ind := ind + 1; fi;
        od;
       if found then
          # Calculate the action of Gcos on the blocks system 
          # corresponding to <L[ind].blk>
	  i := 1;
	  blks1 := [L[ind].blk]; # NB L[ind].blk is immutable
	  rep := [];
	  for pnt in [1..n] do
	    rep[pnt] := 0;
	  od;
	  for pnt in L[ind].blk do
	    rep [pnt] := 1;
	  od;
	  while i <= Length( blks1) do
	    for gen in GeneratorsOfGroup(Gcos) do
	      img := Immutable(OnSets( blks1[i], gen ));
	      if rep[ img[1] ] = 0 then
		Add( blks1, img );
		for pnt in img do
		  rep[pnt] := Length( blks1 );
		od;
	      fi;
	    od;
	    i := i + 1;
	  od;

	  used := [];
	  # Now merge to get those orbits of Gw on this system which
          # contain a point in L[ind].xs
	  # Note that N[1] = 1 so the 1 in the next line does 
	  # correspond to no conjugation
	  for pr in Concatenation(List(L[ind].xs, x -> [x, [1]]), 
		       [[x, nis]]) do
	    x := pr[1];
	    found := false;
	    i := 1;
	    while not found and i <= Length(pr[2]) do
      
	    # for each ni in nis, N[ni] conjugates L[ind].blk to preblk.
	    # conjugate x by N[ni]^-1 to find the equivalent pair.
	    # Then decide whether it's in same Gw-orbit on cos(G:J) as 
	    # a member of L[ind].xs (where J <-> L[ind].blk ).
      
	      ni := pr[2][i]; 
	      if N[ni] <> 1 then
		gx := [];
		img := x;
		while img <> 1 do
		  Add(gx, trans[img]);
		  img := img ^ trans[img];
		od;
		gxr := Reversed(gx);
		img2 := N[ni];
		for g in gxr do
		  img2 := img2 / g;
		od;
		img := N[ni];
		while img <> 1 do
		  img2 := img2 ^ trans[img];
		  img := img ^ trans[img];
		od;
	      else
		img2 := x;
	      fi;
	      if IsBound(used[img2]) then 
		found := true;
	      else
		i := i + 1;
	      fi;
	    od;
	    if not found then
              # need to extend L[ind].xs; img2 contains the conjugate by
              # one of the N[ni] as ni ranges over nis
	      AddSet(L[ind].xs, img2);
	      orb := ShallowCopy(blks1[rep[img2]]);
	      # Find orbit of Gw on L[ind].blk images containing img2
	      for pnt in orb do
		img := DtoOrbs[pnt];
		if not IsBound(used[img]) or used[img] <> 1 then
		  while img <> 0 do
		    if not IsBound(used[img]) then
		      Append(orb, blks1[rep[img]]);
		      for pt in blks1[rep[img]] do
			used[pt] := 0;
		      od;
		    fi;
		    img := next[img];
		  od;
		  used[DtoOrbs[pnt]] := 1;
		fi;
	      od;
	    fi;
	  od;

        else
          # We have to make a new entry in L; with only one x to put in it
          # there is no need for extensive equivalence checking.
          # Calculate 'conjugates' of the elements of blkstab under N -
          # these are stored to speed up future conjugacy tests
          conjstabs := [];
          for npt in N do
            conjstab := [];
            img := npt;
            nptinv := 1;
            gnpt := [];
            while img <> 1 do
              nptinv := nptinv ^ trans[img];
              Add(gnpt, trans[img]);
              img := img ^ trans[img];
            od;
            gnpt := Reversed(gnpt);
            for gpnt in blkstabgens do
              img := nptinv;
              for g in gpnt do
                img := img / g;
              od;
              for g in gnpt do
                img := img / g;
              od;
              Add(conjstab, img);
            od;
            Add(conjstabs, Set(conjstab));
          od;
          L := Concatenation( L{[1..ind-1]}, [ 
		       rec(blk := blk,
		           blkstab := blkstab,
		           xs := [x], 
                           conjstabs := conjstabs )
		     ], L{[ind..Length(L)]} );
        fi;
      fi;
    od;

    # We've set L up. Now comes the main loop where we 'recurse' over 
    # previously found solutions (stored in L) to find all the rest.

    output := [];    
    while L <> [] do

      # Next solution to process is the smallest, ie the front of L.
      blk := L[1].blk;
      blkstab := L[1].blkstab;
      xs := L[1].xs;

      L := L{[2..Length(L)]};


      Add(output, [blk, xs]); 

      # First step is to find the blocks containing blk maximally.
      # To do this we form the action on translates of blk, then
      # call the Schonert-Seress routine

      i := 1;
      blks := [blk];
      elt := [];
      elt[i] := blk[1];
      rep := [];
      for pnt in [1..n] do
        rep[pnt] := 0;
      od;
      for pnt in blk do
        rep [pnt] := 1;
      od;
      while i <= Length( blks) do
        for gen in GeneratorsOfGroup(Gcos) do
          img := OnSets( blks[i], gen );
          if rep[ img[1] ] = 0 then
            Add( blks, img );
            elt[Length( blks )] := img[1];
            for pnt in img do
              rep[pnt] := Length( blks );
            od;
          fi;
        od;
        i := i + 1;
      od;

      gens := List(GeneratorsOfGroup(Gcos), gen ->
                PermList(List(elt, x -> rep[x^gen])));
      GG := Group(gens, ());

      orb := [1..Length(blks)];
      Info(InfoOperation,1,"OrbitalActions: Action on ",Length(orb),
                                 " points being tested");

      # Schonert-Seress called by MinimalBlocks
      # MinimalBlocks not in GAP4b1 (but will be in 4b2!)
      # We require that MinimalBlocks uses <1> as its basepoint.
      blocks := MinimalBlocks(GG, orb);

      # test blocks
      for b in blocks do
        # Calculate block in [1..n] corresponding to b 
        preblk := Set( Concatenation(blks{b}));
        # ignore blocks which are too large
        if Size(preblk)*(Size(preblk)-1) <= n then

          # Get block generators (using minimality of block)
          z := First(preblk, i -> not i in blk);
          preblkstab := Concatenation(blkstab, [z]);

          # Find representatives of the generators in preblkstab
          # stored as reverse of a word in the generators for the inverse
          preblkstabgens := [];
          for pnt in preblkstab do
            img := pnt;
            gen := [];
            while img <> 1 do
              Add(gen, trans[img]);
              img := img ^ trans[img];
            od;
            Add(preblkstabgens, Reversed(gen));
          od;

          # We test the original values of x stored as L[1].xs
          # newxs will contain those that pass, ie for which J\capJ^x=Gw
          # where J = Stab(preblk)
          newxs := [];
          for x in xs do

            # contents of this loop are similar to the section above for
            # testing x.

            # Find representative of x
	    gx := [];
	    pnt := x;
	    while pnt <> 1 do
	      Add(gx, trans[pnt]);
	      pnt := pnt ^ trans[pnt];
	    od;
	    gxr := Reversed(gx);

            # block is automatically self-paired as it arises from a 
            # self-paired block.
            acc := true;

	    # Calculate image of preblk under 'conjugation' by <x> (actually
	    # conjugation of the stabilizer by <gx>)

            # Actually <imgpreblk> will be the translate of the desired 
            # conjugate by <gx>^-1 so <imgpreblk> = 1^(<gx>^(-1)<J>) where
            # <J> = Stab(<preblk>)

            used := [];
            pnt := 1;
            for g in gx do
              pnt := pnt ^ g;
            od;

            imgpreblk := [pnt];
            used[pnt] := 0;

	    i := 1;
	    while acc and i <= Length(imgpreblk) do
	      pnt := imgpreblk[i];
	      # Close under the generators
	      j := 1;
	      while acc and j <= Length(preblkstabgens) do
		gen := preblkstabgens[j];
		# form img = pnt ^ ( gen )
		img := pnt;
		for g in gen do
		  img := img / g;
		od;
		if not IsBound( used[img] ) then
		  img2 := img;
		  for g in gxr do
		    img2 := img2 / g;
		  od;
		  if img2 <> 1 and img2 in preblk then
		    acc := false;
		  else
		    Add(imgpreblk, img);
		    used[img] := 0;
		  fi;
		fi;
		j := j + 1;
	      od;
	      # Close under <Gw>^<gx>
	      img := DtoOrbs[pnt];
	      # Avoid doing this suborbit if we've done it before
	      if not IsBound(used[img]) or used[img] <> -1 then 
		while acc and img <> 0 do
		  if not IsBound( used[img] ) then
		    img2 := img;
		    for g in gxr do
		      img2 := img2 / g;
		    od;
		    if img2 <> 1 and img2 in preblk then
		      acc := false;
		    else
		      Add(imgpreblk, img);
		      used[img] := 0;
		    fi;
		  fi;
		  img := next[img];
		od;
		used[ DtoOrbs[ pnt ] ] := -1;
	      fi;           
	      i := i + 1;
	    od;

            if acc then
              AddSet(newxs, x);
            fi;
          od;

          # Now we add those that pass the test to the data structure
          # much as before except now we are adding more than one x so
          # we need to do the full orbit calculation even if a new
          # entry is required as the different x may now be equivalent
          if newxs <> [] then

	    # Check to see if it's a conjugate of anything we've seen before.
	    ind := PositionSorted(L, rec(blk := preblk),
		    function(x, y) return (Size(x.blk) < Size(y.blk)); end );
	    found := false;
	    while not found and IsBound(L[ind]) and 
			 Size(L[ind].blk) = Size(preblk) do
              nis := [];
              for ni in [1..Length(N)] do
		if IsSubsetSet(preblk, L[ind].conjstabs[ni]) then
                  Add(nis, ni);
		  found := true;
		fi;
	      od;
	      if not found then ind := ind + 1; fi;
	    od;
            if not found then 
	      # calculate 'conjugates' of the elements of preblkstab under N
	      conjstabs := [];
	      for npt in N do
		conjstab := [];
		img := npt;
		nptinv := 1;
		gnpt := [];
		while img <> 1 do
		  nptinv := nptinv ^ trans[img];
		  Add(gnpt, trans[img]);
		  img := img ^ trans[img];
		od;
		gnpt := Reversed(gnpt);
		for gpnt in preblkstabgens do
		  img := nptinv;
		  for g in gpnt do
		    img := img / g;
		  od;
		  for g in gnpt do
		    img := img / g;
		  od;
		  Add(conjstab, img);
		od;
		Add(conjstabs, Set(conjstab));
	      od;
              L := Concatenation(  L{[1..ind-1]}, [ 
			   rec(blk := preblk,
			       blkstab := preblkstab,
			       xs := [], 
			       conjstabs := conjstabs )
			 ], L{[ind..Length(L)]} );
              nis := [];
              for ni in [1..Length(N)] do
	        if IsSubsetSet(preblk, L[ind].conjstabs[ni]) then
                  Add(nis, ni);
	        fi;
	      od;
            fi;

	    # Calculate the action of Gcos on the blocks system 
	    # corresponding to <L[ind].blk>
      
	    i := 1;
	    blks1 := [L[ind].blk]; # L[ind].blk is immutable
	    rep := [];
	    for pnt in [1..n] do
	      rep[pnt] := 0;
	    od;
	    for pnt in L[ind].blk do
	      rep [pnt] := 1;
	    od;
	    while i <= Length( blks1) do
	      for gen in GeneratorsOfGroup(Gcos) do
		img := Immutable(OnSets( blks1[i], gen ));
		if rep[ img[1] ] = 0 then
		  Add( blks1, img );
		  for pnt in img do
		    rep[pnt] := Length( blks1 );
		  od;
		fi;
	      od;
	      i := i + 1;
	    od;


	    used := [];
	    # Now merge to get orbits of Gw on this system.
            # Note that N[1] = 1 so the 1 in the next line does 
            # correspond to no conjugation
	    for pr in Concatenation(List(L[ind].xs, x -> [x, [1]]), 
			List(newxs, x -> [x, nis])) do
	      x := pr[1];
	      found := false;
	      i := 1;
	      while not found and i <= Length(pr[2]) do
      
	      # for each ni in nis, N[ni] conjugates L[ind].blk to preblk.
	      # conjugate x by N[ni]^-1 to find the equivalent pair.
	      # Then decide whether it's in same Gw-orbit on cos(G:J) as 
	      # a member of L[ind].xs (where J <-> L[ind].blk ).
      
		ni := pr[2][i]; 
		if N[ni] <> 1 then
		  gx := [];
		  img := x;
		  while img <> 1 do
		    Add(gx, trans[img]);
		    img := img ^ trans[img];
		  od;
		  gxr := Reversed(gx);
		  img2 := N[ni];
		  for g in gxr do
		    img2 := img2 / g;
		  od;
		  img := N[ni];
		  while img <> 1 do
		    img2 := img2 ^ trans[img];
		    img := img ^ trans[img];
		  od;
		else
		  img2 := x;
		fi;
		if IsBound(used[img2]) then 
		  found := true;
		else
		  i := i + 1;
		fi;
	      od;
	      if not found then
		AddSet(L[ind].xs, img2);
		orb := ShallowCopy(blks1[rep[img2]]);
		# close under Gw and L[ind].blk images.
		for pnt in orb do
		  img := DtoOrbs[pnt];
		  if not IsBound(used[img]) or used[img] <> 1 then
		    while img <> 0 do
		      if not IsBound(used[img]) then
			Append(orb, blks1[rep[img]]);
			for pt in blks1[rep[img]] do
			  used[pt] := 0;
			od;
		      fi;
		      img := next[img];
		    od;
		    used[DtoOrbs[pnt]] := 1;
		  fi;
		od;
	      fi;
	    od; 


          fi;
        fi;
      od;
      
    od;
    
    return output;

end;