## 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
## orbitalso-s-p1.g: adapted version of orbitals5a.g to find self-paired
##   solutions corresponding to one supplied x with x^2 in Gw
## orbitalso-s-p2.g: changes so all solns returned have x as given

OrbitalActionsOneSelfPaired := function(G, n, xx)
                       # G input group acting on [1..n]
                       # xx point in [1..n] such that for any representative
                       # g of xx, g^2 in Gw and Gw^g=Gw.

    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
          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
                       #                  stored as records: .x is the value,
                       #             and .c: if r is a rep of .c then
                       #           there is a soln (J, xx) st J^r=L[i].blk
                       #           and conj of xx by r is L[i].xs[j]
                       # 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 (as records - see above)
          xv,          # one record in xs
          conj,        # the .c component of xv
          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);


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

    # L is the main data structure storing pending solutions
    L := [ rec(blk := [1],
              blkstab := [],
              xs := [rec( x := xx, c := 1)],
              conjstabs := List(N, x -> []) ) ];

    # 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 these solns to output
#      Add(output, [blk, xs]); 
      for xv in xs do
        # conjugate blk and xv.x by inverse rep of xv.c
        # as xv.c normalizes Gw we don't need to worry about Gw orbits
        if xv.c = 1 then
          Add(output, [blk, [xv.x]]);
        else
          conj := [];
          img := xv.c;
          while img <> 1 do
            Add(conj, trans[img]);
            img := img ^ trans[img];
          od;
          imgblk := [];
          for pnt in blk do
            pntg := [];
            img := pnt;
            while img <> 1 do
              Add(pntg, trans[img]);
              img := img ^ trans[img];
            od;
            img := xv.c;
            for g in Reversed(pntg) do
              img := img / g;
            od;
            for g in conj do
              img := img ^ g;
            od;
            AddSet(imgblk, img);
          od;
          gx := [];
          img := xv.x;
          while img <> 1 do
            Add(gx, trans[img]);
            img := img ^ trans[img];
          od;
          img := xv.c;
          for g in Reversed(gx) do
            img := img / g;
          od;
          for g in conj do
            img := img ^ g;
          od;
          Add(output, [imgblk, [img]]);
        fi;
      od;

      # 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 xv in xs do
            x := xv.x;

            # 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
              Add(newxs, xv);
            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, xv -> [xv, nis])) do
	      x := pr[1].x;
              conj := pr[1].c;
	      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
                # Keep x fixed unless absolutely necessary to change it
                if 1 in N{nis} then
		  AddSet(L[ind].xs, pr[1]);
                else
                  # new value of conj is the old one 'multiplied' by N[ni]^-1
                  img3 := conj;
                  img := N[ni];
                  while img <> 1 do
                    img3 := img3 ^ trans[img];
                    img := img ^ trans[img];
                  od;
                  AddSet(L[ind].xs, rec(x := img2, c := img3));
                fi;
		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;