# GAP script to implement UOP algorithm
# Revision 3: using Seress' routine for finding minimal blocks.
# 21/4/1998
# main4.g: major surgery to improve the code and presentation for publication

UOPcore := function(G, n)

# Global variables
local m,      # degree of solution, n = {m\choose 2}
      omega,  # base point (in fact omega = 1 throughout)
      ss,     # G_omega-orbits on [1..n]
      ll,     # lengths of the elements of ss
      trans,  # factored inverse transversal for G_omega in G

# Miscellaneous variables with only localised usage

     ind, inds, i, L, sqrt, ppresA, ppresB, ppres, res, q, x, d 
     , i, j, k, is, js, ks, LL, LLL, img, gg, gen, img1, orbit, pnt

     ,Q, nxt, nxtfree, g


# Function names
 
     ,TestAffine, MaximalBlocksAff, Index2Op, SLOverMoreThanF2, TestIndex2Aff
     , TestIndex2AS
     ,TestIsomorphism, SLOverF2Only, AffFindIndex2,TestSet, TestSolnSubset

     ,TestSympA, TestSympB, TestPSL, TestThreeTrans, TestLR1, TestFnLR1 
     , PrimePower;

# Structure of the code for UOPcore:
# Functions:
#   Functions for affine groups (classes A_o, A_e)
#   Functions for groups of Lie rank 1 (class L and PSL(2,q))
#   Functions for class Z groups
#   Multipurpose utility functions
# Main program

# NB read annotations in conjunction with the LMS JCM paper 
# 'Algorithmic Recognition of Actions of $2$-Homogeneous Groups on Pairs',
# Graham R. Sharp, 1998 

####### Functions ###########################################################
#############################################################################

####### Affine Groups #######################################################
# Structure of the affine groups section
#
# MaximalBlocksAff : adaptation of GAP implementation of Sch\"onert-Seress
#   to find a maximal block whose size is not divisible by p.
# TestAffine : main function called to check an action to see whether it is 
#   the action on pairs of an affine $2$-hom group of degree m=p^d
#   distinguishes the odd case from the even, class A_e from the Z and
#   finds estar if applicable (A_e); handles exceptional cases where
#   estar is wrong or non-existant
# SLOverMoreThanF2 : find a largest suborbit X in K^{2} (ie size divisible by
#   2*estar), pass it to AffFindIndex2 to split X into a block system and
#   thence to TestIndex2Aff (NB in title SL=SemiLinear not Special Linear) 
# TestIndex2Aff : given a list of subgroups of G_omega of index 2 (given by 
#   G_omega blocks), form the action on cosets, call MaximalBlocksAff, 
#   check for an isomorphism
# TestIsomorphism : only called for odd characteristic case where projecting
#   to get a solution subset doesn't work, this tests the action on pairs of
#   a putative solution for isomorphism with the input action
# SLOverF2Only : handles those affine groups in the class Z
# AffFindIndex2 : takes a suborbit on which G_omega should act in a certain,
#   uncomplicated way, and calculates up to three block systems whose 
#   kernels have index 2 in G_omega and might be the stabilizer of an ordered
#   pair
# (TestFnNotDivisibleBy_p) : not used

#############################################################################
##
#F  PermGroupOps.Blocks( <G>, <D> [, <seed>] [, <operation>] )
##
#PermGroupOps.BlocksNoSeed := function ( G, D )
MaximalBlocksAff := function(G, D, d, p, stabgen)
# This is the code for BlocksNoSeed (which is based on the Sch\"onert-Seress
# primitivity-testing algorithm) adapted to 
# (i) look at each minimal block in turn until it finds a suitable one
#      (suitable == size not divisible by p)
# (ii) having found a suitable block, factor out by the corresponding 
#      equivalence relation and continue until a maximal suitable block
#      is found.
# We assume the generators of G_1 are known (and therefore assume omega=1
#   elsewhere in the program) and given in stabgen.
# It should be easy to change the definition of suitable (to 'having p-power
#   size') to adapt this for the Lie rank one case (latest work shows that 
#   such an algorithm would work).
    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, tentative result
            g,          # group element
            i,          # loop variable
            GG,         # action of G on tmpblocks
#            stabgen,    # generators for GG_1
            tmpblocks,  # minimal block system from last iteration
            poss,       # flag to indicate whether we might have a block
            iter,       # which points we've checked when
            stind,      # index in orbit of <start>
            start;      # first cur for this iteration

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

    
    blocks := [];
    for pnt in D do
        blocks[pnt] := [pnt];
    od;

    # repeat to ascend a structure tree. Start with known generators of G_1 
#    stabgen := Stabilizer(G, D[1]).generators;
    repeat
#Print(blocks[1],"\n");

    # operate on [1..Length(blocks)] in the same way as on blocks
    if Size(blocks[1]) > 1 then
        GG := Operation(GG, tmpblocks, OnSets);

        # Now recalculate stabgen
        stabgen := List( stabgen, g -> Permutation(g, tmpblocks, OnSets));

    else
        GG := Operation(G, D);
    fi;


    # compute the orbit of $G$ and a factored transversal
    orbit := [ 1 ];
    trans := [];
    trans[ 1 ] := ();
    for pnt  in orbit  do
        for gen  in GG.generators  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");

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

    # Deterministically calculate suborbits
    for gen in stabgen do

	# 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];
	    elif img < pnt  then
		eql[pnt] := img;
		next[ last[img] ] := pnt;
		last[img] := last[pnt];
	    fi;

	od;

    od;

    # Repeat until we run out of points
    stind := 2;
    block := [1];
    repeat
        start := orbit[stind];
        cur := start;
        # unless this is not 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
    
	    # take arbitrary point <cur>, and 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 <> 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 ..
		if   img <> 1  and img <> cur then
                  # .. but a new block
                  if 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 <> 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 current block
		  # Not any more it doesn't! Now we also have to check whether 
		  # the block appeared this time or earlier. 
		  elif 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;
                  # That should be all the possibilities if it's not the 
                  # current block, so this is just to avoid accidents:
                  else
                    pnt := next[pnt];
                  fi;

		# else it is the current block,
                # 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 := [1];
		 pnt := cur;
		 while pnt <> 0  do
		     Add( block, pnt );
		     pnt := next[pnt];
		 od;
		 block := Set( block );

		 # Test to see if this is the right size
		 if  RemInt(Size(block),p) <> 0 then
		   # If so then add gen to generators of stabilizer
                   g := ();
                   for i in gen do
                       g := g / i;
                   od;
                   Add(stabgen, g);
                   # and get out of this loop
                   stind := Length(orbit)+1;

		 # quick test to see if the action is primitive
		 elif Length( block ) = Length( orbit )  then
                   # if so, end iteration
                   stind := Length(orbit)+1; 
     
                 else
                   stind := stind + 1;
		 fi;
            else
              stind := stind + 1;
	    fi;

        else
          stind := stind + 1;
        fi;
        # loop back to get another minimal block
    until stind > Length(orbit);

    # If a block of suitable size has been found then blocks needs updating
    if Size(block) > 1 and RemInt(Size(block), p) <> 0 then

	# compute the block system in the original domain 
        # with an orbit algorithm
	tmpblocks := [block];
	# save space:
	rep := iter;
	for pnt  in orbit  do
		rep[pnt] := 0;
	od;
	for pnt  in block  do
		rep[pnt] := 1;
	od;
	    
	i := 1;
	while i <= Length( tmpblocks )  do
    
	    # loop over the generators
	    for gen  in GG.generators  do
    
		# compute the image of the block under the generator
		img := OnSets( tmpblocks[i], gen );
    
		# if this block is new
		if rep[ img[1] ] = 0  then
    
		    # add the new block to the list of blocks
		    Add( tmpblocks, img );
		    for pnt  in img  do
			rep[pnt] := img[1];
		    od;
    
		fi;
    
	    od;
    
	    # on to the next block in the orbit
	    i := i + 1;
	od;
        # Adjust blocks in line with the new system
        blocks := List(tmpblocks, blk -> Union( blocks{blk} ) );

    fi;

#Print(block,"\n");

    until Size(block) = Length(orbit) or Length(blocks) = p^d or
            RemInt(Size(block), p) = 0;

    # return the block system
    return blocks;
end;

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

TestAffine := function(p,d)

local divs, a, l, estar, estars, blocks, i, j, ind, X, Zs;

# Suborbits(1); # Sets up global variables ss, ls, omega

if p <> 2 then 
  blocks := MaximalBlocksAff(G, [1..n], d, p, Stabilizer(G, 1).generators);
  if Length(blocks) = m then # need to test for isomorphism.
    return TestIsomorphism( blocks );
  else
    return [false];
  fi;
else

  # set estars to the set of integers st Q_estar has size 2^estar choose 2.
  divs := DivisorsInt(d);
  estars := Filtered(divs, function(a)
      local ds, targ;
      ds := DivisorsInt(2*a);
      targ :=  2^(a-1)*((2^a)-1);
      return ( Sum( List( Filtered( Collected( ll ), l -> l[1] in ds ), 
                       l -> l[1]*l[2] ) ) = targ );
      end);

  if d <> 4 then
    if estars=[] then 
      if d = 8 then
        # ASL(2,16) \leq G \leq ASigmaL(2,16)
        if Collected( ll )  =  [[1,8],[2,56],[8,225],[32,960]] then
           # ASL(2,16)
           Zs := AffFindIndex2( ss[ PositionProperty( [1..1249],
                                      i -> ll[i] = 2 ) ] );
           return TestIndex2Aff( Zs, d, p );
        elif Collected( ll ) = [[1,4],[2,10],[4,24],[8,9],[16,108],
                    [32,48],[64,456]] then
           # ASL(2,16).2
           Zs := AffFindIndex2( ss[ PositionProperty( [1..659],
                                      i -> ll[i] = 4 ) ] );
           return TestIndex2Aff( Zs, d, p  );
        elif Collected( ll ) = [[1,2],[2,1],[4,5],[8,13],[16,4],[32,58],
                    [64,22],[128,228]] then
           # AZL(2,16)
           i := PositionProperty( [1..333], j -> ll[j] = 8 );
           Zs := Concatenation([ AffFindIndex2( ss[i] ), 
                   AffFindIndex2( ss[ PositionProperty( [(i+1)..333],
                                         j -> ll[ j ] = 8) ] ) ] );
           return TestIndex2Aff( Zs, d, p  );
        else
          return [false];
        fi;
      else 
        return [false];
      fi;
    else 
      estar := Maximum(estars);
      if estar = 1 then 
        # Type 4,4a, 5, 5a, 6, 6a
        return SLOverF2Only(d);
      else 
        # G semilinear over F_q, q= 2^{estar}>2
        return SLOverMoreThanF2(estar, p, d);
      fi;
    fi;
  else
    if estars = [] then
      # G = ASigmaL(2,4)
      if Collected( ll )  =  [[1,2],[2,1],[4,5],[8,4],[16,4]] then
        Zs := Concatenation( List( Filtered( [1..16], i -> ll[i] = 4 ),
                            i -> AffFindIndex2( ss[i] ) ) );
        return TestIndex2Aff( Zs, d, p );
      else
        return [false];
      fi;
    else
      estar := Maximum(estars);
      if estar = 1 then
        # Type 4,4a,5,5a
        return SLOverF2Only(d);
      elif estar = 2 then
        #  G semilinear over F_q, q= 2^{estar}=4, in fact AGL(2,4)\leq G
        return SLOverMoreThanF2(estar, p, d);
      else # estar = 4; all suborbits have lengths dividing 8 
        if (4 in ll) or not (4 in ll or 8 in ll) then 
          # G semilinear over F_16 (ie G \leq AGL(1,16) )
          return SLOverMoreThanF2(estar, p, d);
        else
          # G = ASL(2,4). Apply Prop 5.4 to each suborbit of length 2,
          # checking first there are the right number of each length
          if Collected( ll ) = [[1,2], [2,11], [8,12]] then 
            Zs := Concatenation( List( Filtered( [1..25], i -> ll[i] = 2 ),
                            i -> AffFindIndex2( ss[i] ) ) );
            return TestIndex2Aff( Zs, d, p  );
          else
            return [false];
          fi;
        fi;
      fi;
    fi;
  fi;
fi;
end;

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

SLOverMoreThanF2 := function(estar, p, d)

local max, i, ind, X, Zs;

  # Find a largest suborbit in K^{2}
  max := 0;
  for i in [1..Length( ss )] do
    if ll[i] > max and RemInt( 2 * estar, ll[i] ) = 0 then 
      max := ll[i];
      ind := i;
    fi;
  od;
  X := Set( ss[ind] );
  Zs := AffFindIndex2( X );
  return TestIndex2Aff( Zs, d, p  );
end;

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

TestIndex2Aff := function(Zs, d, p)
local Delta, Z, blocks, A, x, res, GonPairs, GonPairsHom, stabgen;

# Test putative actions on ordered pairs corresponding to the subsets in Zs
# which should be such that Z^{G_omega} has size 2 for each Z in Zs.
  Delta := Cartesian( [1..n], [false, true] );
  for Z in Zs do
    GonPairs := Operation(G, Delta, Index2Op(Z));
    GonPairsHom := OperationHomomorphism(G, GonPairs);
    stabgen := List( Stabilizer(Stabilizer(G, omega), false, 
                   function(x, g) 
                     if Z[1]^g in Z then return x;
                     else return not x; fi;
                   end ).generators, g -> Image(GonPairsHom, g));
    blocks := MaximalBlocksAff(GonPairs, [1..2*n], d, p, stabgen);
    if Length(blocks) = m then # Got to test this action for isomorphism
      # Note that if (G, blocks,OnSets) is the action we want then the
      # projection of blocks[1] onto the first coordinate is a 
      # solution subset. This is how we will test isomorphism

      A := Set(List( Delta{blocks[1]}, x -> x[1] ) );
      if Size(A) = m-1 then
        res := TestSolnSubset(A);
        if res[1] = true then return res; fi;
      fi;
    fi;
  od;
  return [false];
end;

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

TestIsomorphism := function(blocks)
# ONly works where the domain is the set of integers [1..n]
# blocks should have size m

local Gomega, lookupA, lookupB, op, fixed, test, b, d, i, j, g, pr, 
        img, genimgs, etainv, solnset;   

    Gomega := Stabilizer(G, omega);
 
    # G acts on [1..m] by operation op.
    # This is just a quick way of doing Operation(G, blocks, OnSets) but 
    # ending up with G acting by op, not an image of G acting by OnPoints
    lookupA := List( blocks, b -> b[1] );
    lookupB := [];
    for i in [1..m] do 
      for j in blocks[i] do
        lookupB[j] := i;
      od;
    od;

    op := function( d, g )
      return lookupB[ lookupA[ d ] ^ g ];
    end;

    test := Combinations( [1..m], 2 );
    for g in Gomega.generators do
      fixed := [];
      for pr in test do
        # ascertain whether the pair pr is fixed
        # by g under the operation op 
        # if not remove pr from fixed ('cos it isn't fixed!)
        img := Set( List( pr, d -> op(d,g) ) );
        if img = pr then
          AddSet(fixed, pr);
        fi;
      od;
      test := fixed;
    od;
    fixed := test;
    if fixed = [] then 
      return [false];
    else
      # isom maps identify omega and [blocks[pr[1]], blocks[pr[2]]]
      # for any pr in fixed
      genimgs := List( G.generators, g -> Permutation( g, [1..m], op) );
      etainv := List( [1..n], x -> Set( List( fixed[1],
                   function(d) 
                     local img, g, gen;
                     img := x;
                     g := [];
                     while img <> omega do
                       Add(g, trans[img]);
                       img := img ^ trans[img];
                     od;
                     img := lookupA[ d ];
                     for gen in Reversed(g) do
                       img := img / gen;
                     od;
                     return lookupB[img];
                   end 
#                   d -> lookupB[ TransImg( x, lookupA[ d ] ) ]
                       ) ) );
      solnset := Filtered( [1..n], i -> 1 in etainv[i] );
      return [true, solnset, genimgs, etainv];
    fi;
end;

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

SLOverF2Only := function (d)
# Types 4,4a, 5, 5a, 6, 6a
# Call TestSet one or more times

local num, L, LL, LLL, m, i, j, is, js, res;

  num := Length(ss);
  # Recall m = 2^d. But for safety we'll use a local variable of same name
  m := 2^d;
  if num = 4 then
    # Type 4,4a
    if Collected( ll ) = Collected( [1,2*m-4,m/2-1,(m-2)*(m-4)/2] ) then
      L :=  Set( ss[Position( ll, 2*m-4 ) ] );
      AddSet( L, omega );
      return TestSet( omega, L);
    else
      return [false];
    fi;
  elif num in [9,10] and RemInt( d, 2 ) = 0 then
    # Type 5,5a
    if num = 9 and d = 4 then 
      if Collected( ll ) = [[1,1],[3,1],[4,1],[12,4],[16,1],[48,1]] then
        L := Set( ss[Position( ll, 16 )] );
        AddSet( L, omega );
        is := Filtered( [1..9], i -> ll[i] = 12 );
        for i in is do
          LL := ShallowCopy(L);
          UniteSet( LL, ss[i] );
          res := TestSet( omega, LL );
          if res[1] = true then return res; fi;
        od;
      fi;
      return [false];
    elif num = 10 and d > 4 then
      if Collected( ll ) = Collected( [1,m-4,m/4-1,(m-4)*(m/16-1),
         m*(m/4-1)/4, m, m/4, m*(m/4-1)/4, m*(m/4-1)/4, m*(m/4-1) ] ) then
        L := Set( ss[Position( ll, m ) ] );
        UniteSet( L, ss[Position( ll, m-4 ) ] );
        AddSet( L, omega );
        return TestSet( omega, L );
      fi;
    else
      return [false];
    fi;
  elif num = 25 and d = 6 then
    # Type 6
    if Collected( ll ) = [[1,1],[3,1],[12,5],[16,1],[48,5],[64,1],
                             [96,7],[192,3],[384,1]] then
      L := Set( ss[Position( ll, 64 )] );
      AddSet( L, omega );
      is := Filtered( [1..25], i -> ll[i] = 48 );
      js := Filtered( [1..25], i -> ll[i] = 12 );
      for i in is do
        LL := ShallowCopy(L);
        UniteSet( LL, ss[i] );
        for j in js do 
          LLL := ShallowCopy(LL);
          UniteSet( LLL, ss[j] );
          res := TestSet( omega, LLL );
          if res[1] = true then return res; fi;
        od;
      od;
    fi;
    return [false];
  elif num = 31 and d = 6 then
    # Type 6a  
    if Collected( ll ) = [[1,1],[3,1],[12,5],[16,1],[32,2],[48,7],[96,12],
                               [192,2]] then
      is := Filtered( [1..31], i -> ll[i] = 32 );
      L := [omega];
      for i in is do UniteSet( L, ss[i] ); od;
      is := Filtered( [1..31], i -> ll[i] = 48 );
      js := Filtered( [1..31], i -> ll[i] = 12 );
      for i in is do
        LL := ShallowCopy(L);
        UniteSet( LL, ss[i] );
        for j in js do 
          LLL := ShallowCopy(LL);
          UniteSet( LLL, ss[j] );
          res := TestSet( omega, LLL );
          if res[1] = true then return res; fi;
        od;
      od;
    fi;
    return [false];
  else
    return [false];
  fi;
end;

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

AffFindIndex2 := function(X)
# Run proposition 5.4 on suborbit X (bet it's got a different number now!)
# Returns a list of (1,2,or3)  subsets Z of X of size #X/2 st G_omega acts 
# non-trivially on {Z,X-Z}, under the assumptions listed in Propn 5.4
# We do extra testing: we ensure that the isomorphism type of G_omega/N is
# is C_f \times C_2

local x0, GomegaOnX, tr, lenX, cgenelts, f, g, x, sqblkseed, Ms, M, 
      Mm, MM, Zs;

  lenX := Length(X);
  GomegaOnX := Operation( Stabilizer( G, omega), X );  
  x0 := 1;
  tr := RepresentativesOperation( GomegaOnX, x0 );
  # Now tr \subset GomegaOnX = G/N where N is kernel of action of G on X
  # The hypotheses are not fulfilled unless GomegaOnX acts regularly on X 
  # and is isomorphic to C_f \times C_2. Note that since it acts faithfully
  # and transitively, testing whether it is regular is equivalent to testing
  # whether it is abelian.
  if not IsSemiRegular(GomegaOnX, [1..lenX]) then return []; fi;
 
  f := lenX/2;
  if not IsInt(f) then return []; fi;

  # Find elements of order f
  cgenelts := Set( List( Filtered(tr, g -> Order(GomegaOnX, g) = f),
                     g -> x0^g ) );

  # Form block M corresponding to subgroup generated by squares of all 
  # elements of GomegaOnX
  sqblkseed := Set(List( tr, g -> (x0^g)^g));
  if Size(sqblkseed) = 1 then 
    Ms := List( [1..lenX], x -> [x] );
  else
    Ms := Blocks( GomegaOnX, [1..lenX], sqblkseed );
  fi;
  M := First( Ms, x -> x0 in x );

  if Length(Ms) = 2 then
    if PositionProperty( M, m -> m in cgenelts ) = false then return [];
    elif RemInt(f,2) = 0 then return [];
    else
      return [X{M}];
    fi;
  elif Length(Ms) = 4 then 
    Zs := [];
    for Mm in Ms do
      if Mm <> M then
         MM := Union( M, Mm );
         if PositionProperty(MM, m -> m in cgenelts ) <> false then
            Add( Zs, X{MM} );
         fi;
      fi;
    od;
    if RemInt(f,2) = 0 then
      if Length(Zs) = 2 and RemInt(f, 4) = 0 then return Zs;
      elif Length(Zs) = 3 and RemInt(f, 4) <> 0 then return Zs;
      else
        return [];
      fi;
    fi;
  else
    return [];
  fi;

end;

####### End of Affine section ###############################################
#############################################################################
####### Functions for Lie rank 1 ############################################
#
# Outline of Lie rank 1 section
# TestLR1 : main function for recognition of PSL_2, 
#    PSU_3, Sz, R acting on pairs
# TestIndex2AS : takes subgroups of index 2 in G_omega (as G_omega blocks),
#    forms action on cosets, tests by looking for blocks of size m-1
# TestFnLR1 : function passed to MaximalBlocks (use of a different function
#    here would enable MaximalBlocks to be used in the affine case)
# MaximalBlocks : actually just used here to find blocks of size m-1



TestLR1 := function( q )
# NB m = q+1 always.

local i, L, res, lr, r, x, sr, a, Gomega, gs, M, Ms, Zs, cgens, g,
     ordH, cblks, cblk, c, sqblkseed, indM, q2, q3;

  res := TestThreeTrans();
  if res[1] = true then return res; fi;

  lr := Maximum( ll );
  r := PositionProperty( ll, x -> x = lr );
  # Can test here to see if lr is a suitable size but the list of 
  # possibilities is quite long.
  sr := ss[ r ];
  a := sr[1];
  Gomega := Stabilizer( G, omega );
  if not IsSemiRegular( Gomega, sr ) then return [false]; fi; 
  gs := RepresentativesOperation( Gomega, a );

  # Gab contains a cyclic subgroup of maximum order of Gomega (but not _all_
  # cyclic subgroups of this order, eg in PSigmaU(3,3) )
  # Gab also contains (Gw)^2 as Gab has index 2 in Gw.
  # Find elements of maximal order

  ordH := Maximum( List( gs, g -> Order( Gomega, g ) ) );
  cgens := Filtered( gs, g -> Order( Gomega, g ) = ordH );

  cblks := [];
  for c in cgens do
    if not a^c in Concatenation( cblks ) then
      cblk := [];
      x := a;
      for i in [0..(ordH-1)] do
        AddSet(cblk, x);
        x := x^c;
      od;
      Add(cblks, cblk);
    fi;
  od;

  sqblkseed := Set(List( gs, g -> (a^g)^g ) );
  
  Zs := [];
  for cblk in cblks do
    Ms := Blocks( Gomega, sr, Union( sqblkseed, cblk ) );
    indM := PositionProperty( Ms, x -> a in x );
    M := Ms[ indM ];

    # Ms should in all cases have size 2 or 4.
    if Length( Ms ) = 2 then
      AddSet( Zs, M );
    elif Length( Ms ) = 4 then
      for i in [1..4] do
        if i <> indM then
          AddSet(Zs, Union( M, Ms[i] ) );
        fi;
      od;
    else 
      return [false];
    fi;
  od;

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

TestIndex2AS := function(Zs, test)
local Delta, Z, blocks, A, x, res;

# Test putative actions on ordered pairs corresponding to the subsets in Zs
# which should be such that Z^{G_omega} has size 2 for each Z in Zs.
  Delta := Cartesian( [1..n], [false, true] );
  for Z in Zs do
    blocks := MaximalBlocks(G, Delta, Index2Op(Z), test );
    if Length(blocks) = m then # Got to test this action for isomorphism
      # Note that if (G, blocks,OnSets) is the action we want then the
      # projection of blocks[1] onto the first coordinate is a 
      # solution subset. This is how we will test isomorphism

      A := Set(List( blocks[1], x -> x[1] ) );
      if Size(A) = m-1 then
        res := TestSolnSubset(A);
        if res[1] = true then return res; fi;
      fi;
    fi;
  od;
  return [false];
end;

  return TestIndex2AS( Zs, TestFnLR1 );
end;

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

TestFnLR1:= function( x ) 

  if x = m - 1 then 
    return 2;
  else 
    return 0;
  fi;
end;

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

MaximalBlocks := function(G, D, op, test)

# Search for a block system of G on D (assumed transitive) which is
# maximal with respect to passing condition test on block sizes. test
# takes one argument, the size of a block in a block system, and returns 0
# if that size fails the test, 1 if it passes the test, and 2 if it
# passes the test and (as a consequence of the size) must be maximal
# among blocks passing the test

  local G2, D2, S, SS, s, ss, S1, U, B, B2, X, pos, size, done, res, b, b2, i;

  if not IsPermGroup( G ) or not ForAll( D, IsInt ) or op <> OnPoints  then
    G2 := Operation( G, D, op );
    D2 := [ 1 .. Length( D ) ];
  else
    G2 := G;
    D2 := D;
  fi;

  size := 1;
  done := false;

  ## Let SS be the list of G_1 orbits
  SS := Orbits(Stabilizer(G2, D2[1]), D2, OnPoints);
  

  ## Starting positions
  S:=[ D2[1] ];
  Unbind(SS[ PositionProperty( SS, s -> D2[1] in s ) ] );
  pos := 1;

  ## For each suborbit in turn, try adding it to the seed S. If resulting 
  ## block passes test, add block to S and avoid checking it again
  while pos <= Length(SS) and not done do
    if IsBound(SS[pos]) then
      ss := SS[pos];
      S1 := Concatenation(S, ss);
      B2 := Blocks(G2, D2, S1);
      size := Length(B2[1]);
      res := test(size);
      if res > 0 then
	S := First(B2, b -> 1 in b);
	for i in [(pos+1)..Length(SS)] do
	  if IsBound(SS[i]) and SS[i][1] in S then Unbind(SS[i]);fi;
	od;
      fi;
      if res = 2 then 
        done := true;
      fi;
    fi;
    pos := pos + 1;
  od;

  if D2 <> D then 
    B := [];
    for b2 in Blocks(G2, D2, S) do
      b := [];
      for i in b2 do
        AddSet( b, D[i] );
      od;
      AddSet( B, b );
    od;
  else
    B := Blocks(G2, D2, S);
  fi;

  return B;

end;

####### End of Lie rank 1 section ###########################################
#############################################################################
####### Functions for the class Z ###########################################
#
# Outline of the class Z section
#   Much of the work is done in the main code; these are just a few 
# subroutines called to do particular tasks
# TestSympA, TestSympB : handle the symplectic groups over F_2
# TestPSL : handle PSL(d,q) for fixed d\geq 3 and q
# TestThreeTrans : handle any 3-transitive group

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

TestSympA := function( d )
# Here epsilon = 0.
# Assume d > 2 - then none of the 9 suborbits vanishes.
local z, L, x;

  z := (2^(d-2) + 1) * (2^(d-1) - 1);
  if Collected( ll ) = Collected( [ 1, 4*z, 2^(2*d-1), z, 2^(2*d-3)*z,
             4*(2^(d-2)-1)*(2^(d-3)+1)*z, 2^(2*d-3)*z, 
             2^(3*d-5)*(2^(d-1)-1), 2^(2*d-1)*z ] ) then
    # There can only be one suborbit of size 4z and one of size 2^{2d-1}
    L := Set( ss[ PositionProperty( ll, x -> x = 4*z ) ] );
    UniteSet( L, ss[ PositionProperty( ll, x -> x = 2^(2*d-1) ) ] );
    AddSet( L, omega );
    return TestSet( omega, L );
  else
    return [false];
  fi;
end;

TestSympB := function( d )
# Here epsilon = 1.
# Assume d > 2. d=3 is a special case, when there are only 8 suborbits
local z, L, x, lll;

  if d = 3 then 
    lll := Concatenation( ll, [0] );
  else
    lll := ll;
  fi;
  z := (2^(d-2) - 1) * (2^(d-1) + 1);
  if Collected( lll ) = Collected( [ 1, 4*z, 2^(2*d-1), z, 2^(2*d-3)*z,
             4*(2^(d-2)+1)*(2^(d-3)-1)*z, 2^(2*d-3)*z, 
             2^(3*d-5)*(2^(d-1)+1), 2^(2*d-1)*z ] ) then
    # There can only be one suborbit of size 4z and one of size 2^{2d-1}
    L := Set( ss[ PositionProperty( ll, x -> x = 4*z ) ] );
    UniteSet( L, ss[ PositionProperty( ll, x -> x = 2^(2*d-1) ) ] );
    AddSet( L, omega );
    return TestSet( omega, L );
  else
    return [false];
  fi;
end;

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

TestPSL := function( d, q )
# Handle PSL(d,q) (d>2) and A7 < PSL(4,2)

local is, js, i, j, L, LL, res;

  is := Filtered( [1..Length( ll )], i -> ll[ i ] = 
                      2 * q^2 * (q^(d-2) - 1) / (q - 1) );
  js := Filtered( [1..Length( ll )], j -> ll[ j ] = 
                      2 * (q - 1 )  );

  if ( ( q <> 3  and Length( is ) = 1 ) or 
       ( q = 2 and d = 4 and Length( is ) = 2 ) or
       ( q = 3 and Length( is ) = 4 )   ) and 
       1 <= Length( js ) and Length( js ) <= (q + 2)/4 then
    for i in is do
      L := Set( ss[ i ] );
      AddSet( L , omega );
      for j in js do
        LL := ShallowCopy( L );
        UniteSet( LL, ss[ j ] );
        res := TestSet( omega, LL );
        if res[1] = true then return res; fi;
      od;
    od;
    return [ false ];
  else
    return [ false ];
  fi;   

end;

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

TestThreeTrans := function()
local i, L, res, x;

  for i in Filtered( [1..Length( ll )], x -> ll[x] = 2*(m-2) ) do
    L := Set(ss[ i ]);
    AddSet( L, omega );
    res := TestSet( omega, L );
    if res[1] = true then return res; fi;
  od;
  return [false];
end;

####### End Class Z section #################################################
#############################################################################


####### Multipurpose utility functions ######################################
#
# Descriptions of the functions in this section
# Index2Op : used to construct the action on cosets of a subgroup of G_omega
#   of index 2 given a G_omega block whose stabilizer in G_omega is the
#   desired subgroup. It returns a function which can be submitted to 
#   Operation (etc) to give the action.
# PrimePower : if the argument is a prime power then returns the prime and
#   the index, otherwise returns false.
# TestSet : the TestSet of the paper, takes the point omega, and a G_omega-
#   invariant set L of size 2m-3 and decides if L is the adjacent-point set
#   of an  action on unordered pairs
# TestSolutionSubset : tests the argument to see if it is a solution subset,
#   using the algorithm in the paper, and returns a solution if it is

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

Index2Op := function(Z)
  return 
    function(delta, g)
      local res, z;
      res := [];
      res[1] := delta[1]^g;
      img := delta[1];
      gg := [];
      while img <> omega do
        Add(gg, trans[img]);
        img := img ^ trans[img];
      od;
      z := Z[1];
      for gen in Reversed(gg) do
        z := z / gen;
      od;
      z := z ^ g;
      img := res[1];
      while img <> omega do
        z := z ^ trans[img];
        img := img ^ trans[img];
      od;
#      z := TransImgQuo( res[1], ( TransImg( delta[1], Z[1] )  ) ^ g ) ;
      if z in Z then 
        res[2] := delta[2];
      else
        res[2] := not delta[2];
      fi;
      return res;
    end;
end;

PrimePower := function( x )
local p;

  p := SmallestRootInt( x );
  if not IsPrimeInt(p) then 
    return false;
  fi;
  
  return [ p, LogInt( x, p ) ];
end;

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

# GAP script for UOP problem.
# TestSet routine: transfer thesis section 3.4

# Global variables: G = input perm group = (G,[1..n], OnPoints)
# m = output degree (must be > 4)
# n = input degree = m \choose 2


TestSet := function(omega, L)

# omega in [1..n], L is a set, a subset of [1..n] of size 2m-3
# containing omega.

local omega1, omega2, A, L1, L2, I1, I2, g, i, img, x, gen;

# ASSUME that L is G_{omega}-invariant. It's easier to ensure this by
# construction than to test it now.

# CHECK that #L=2m-3 (to detect errors in programming -
# I'm not proposing that in normal operation this test will ever
# evaluate to true).

if Size(L) <> 2*m-3 then Error("TestSet: #L<>2m-3"); fi;

if not (omega in L) then Error("TestSet: omega not in L"); fi;

# Chose omega1 in L_{omega}
omega1 := L[1];
if omega1 = omega then omega1 := L[2]; fi;

# Set L1 to L_{omega1} (=L^g where omega^g = omega1)

g := [];
img := omega1;
while img <> omega do
  Add(g, trans[img]);
  img := img ^ trans[img];
od;
g := Reversed(g);

L1 := [];
for x in L do
  img := x;
  for gen in g do
    img := img / gen;
  od;
  AddSet(L1, img);
od;

# L1 := Set( List( L, x -> TransImg( omega1, x ) ) );

# Find the intersection of L1 and L
I1 := Intersection( L, L1 );
if Size(I1) <> m then return [false]; fi;

# Find omega2 in L, L1, but distinct from omega, omega1
i := 1;
repeat
  omega2 := I1[i];
  i := i+1;
until not( (omega2 = omega) or (omega2 = omega1) );

# Set L2 to L_{omega2}
g := [];
img := omega2;
while img <> omega do
  Add(g, trans[img]);
  img := img ^ trans[img];
od;
g := Reversed(g);

L2 := [];
for x in L do
  img := x;
  for gen in g do
    img := img / gen;
  od;
  AddSet(L2, img);
od;

# L2 := Set( List( L, x -> TransImg( omega2, x ) ) );

# Find the intersection I2 of L, L1, and L2
I2 := Intersection( I1, L2 );

# Find the candidate A for the solution set
# Note that if m=4 then both tests are true - but in this case if 
# I1-{omega2} is a solution set, so is I2, and vice-versa.
if Size(I2) = 3 then 
  A := ShallowCopy(I1); 
  RemoveSet( A, omega2 );
elif Size(I2) = m-1 then
  A := ShallowCopy(I2);
else
  return [false];
fi;

return TestSolnSubset( A );

end;

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

TestSolnSubset := function(A)

# Now comes the section of code that checks if A is a solution subset
# and finds a solution pair if it is.

# if TestSet returns true, then gen_images will be a list of the
# images of the corresponding elements of G.generators under an
# embedding G -> Sym(m), and eta_inv will be a G-isomorphism [1..n] ->
# [1..m]^{2} (by the isomorphism given in gen_images) (eta_inv because
# eta is meant to map the other way)

local  gen_images, eta_inv, g, i, omega1, omega2, x,
      mzero, T, W, w, found, done, q, S, U, perm;

mzero := []; mzero[m] := []; for i in [1..m] do mzero[i] := 0; od;
T := []; T[m] := [];
W := []; W[n] := [];
T[1] := A;
for w in [1..n] do W[w] := []; od;
for w in A do W[w] := [1]; od;
found := 1;
done := 0; # Set Q in paper is [(done+1) .. found], k is found

while done < found do
  done := done + 1;
  q := done;
  for g in G.generators do
    S := Set( List( T[q], x -> x^g ) );
    omega1 := S[1];
    omega2 := S[2];
    if W[omega1] = [] or W[omega2] = [] or 
        ( Length( W[omega1] ) = 1 and Length( W[omega2] ) = 1 and 
          W[omega1] <> W[omega2] ) then
      found := found + 1;
      if found > m then return [false]; fi;
      T[found] := S;
      U := ShallowCopy( mzero );
      for w in S do
        # Care to avoid accessing elements of W[w], U which aren't defined
        if Length( W[w] ) > 1 then return [false]; 
        elif Length( W[w] ) = 1 then 
          if U[ W[ w ][ 1 ] ] = 1 then return [false]; fi;
        fi;
        if Length( W[w] ) = 1 then U[ W[ w ][ 1 ] ] := 1; fi;
        Add( W[w], found );
      od;
    else 
      # By construction no two 2-element W[w]s are the same, so the
      # following is equivalent to testing whether the intersection of
      # W[w] as w ranges over S has size 1.
      # We can probably do without this test but I included it in the
      # transfer thesis and haven't proved it's not necessary.
      x := W[omega1][1];
      if not (x in W[omega2]) then 
        if Length( W[omega1] ) = 2 then x := W[omega1][2];
        else return [false];
        fi;
      fi;
      for w in S do
        if not (x in W[w]) then return [false]; fi;
      od;
      # End of test of doubtful necessity
    fi;
  od;
od;
if found < m then return [false]; fi;

# We know that we will return true; now we calculate gen_images
# and  eta_inv.

# the elements of W are sorted by construction, so W[w] is a set
eta_inv := W; 

gen_images := [];
for g in G.generators do
  perm := [];
  for i in [1..m] do
    Add( perm, Intersection( W[T[i][1]^g], W[T[i][2]^g] )[1] );
  od;
  Add( gen_images, Sortex( perm ));
od;

return [true, A, gen_images,eta_inv];
end;

####### End of utility section ##############################################

####### End of functions ####################################################
#############################################################################

####### Main program ########################################################
#
# Stucture of main program: based around working out from the degree
# what the possibilities are and calling appropriate subroutines.
#  First decide if n is m choose 2 for some m.
#  Then set up inverse factored transversal for G_omega in trans and find
#     the G_omega-suborbits
#  Check for A_m, S_m (and other 4-transitive groups);
#  Then handle small degrees and sporadic groups
#    Every value of m for which something strange happens is listed; not 
#    all require special checks so many degrees are commented out. 
#  Then the infinite families:
#    Sp(2d, 2) (two families of actions);
#    Affine groups;
#    PSL(d,q), d\geq 3;
#    Lie rank 1 groups (including PSL(2,q), as well as all of the class L).
#  If nothing recognised, give up! There are no more possibilities left.

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

# Check transitivity and that the degree is of a suitable size, find
# m such that n = m(m-1)/2

  if not IsTransitive( G, [1..n] ) then
    return [false];
  fi;

  # Find m such that n = m(m-1)/2
  sqrt := RootInt(1+8*n);
  if sqrt*sqrt <> 1+8*n then
    return [false];
  else
    m := (1 + sqrt)/2;
  fi;
# end;

#  m := InitialTests();
#  if m = false then return [false]; fi;

# Suborbits := function(omega1)

# Calculate the suborbits of (G,[1..n]) at omega= omega1.

# Set up global variables: omega, ss, ll, trans;

  omega := 1;
  # Make Schreier tree by running BFS

  orbit := [omega];
  trans := [];
  for pnt in orbit do
    for g in G.generators do
      img := pnt / g;
      if not IsBound(trans[img]) then
        Add(orbit, img);
        trans[img] := g;
      fi;
    od;
  od;
# 
#   schp := []; schg := [];
#   schp[omega] := omega;
#   schg[omega] := G.identity;
#   Q := [omega]; 
#   Q[n] := 0;
#   nxt := 1; 
#   nxtfree := 2;
#   while nxt < nxtfree do 
#     q := Q[nxt];
#     nxt := nxt + 1;
#     for g in G.generators do
#       if not IsBound( schp[q^g] ) then
#         schp[ q^g ] := q;
#         schg[ q^g ] := g;
#         Q[nxtfree] := q^g;
#         nxtfree := nxtfree + 1;
#       fi;
#     od;
#   od;
#   # Now for all w <> omega, trans[w] is given by trans[schp[w]]*schg[w]
#   # trans[omega] = schg[omega] = G.identity. Note that omega is the only 
#   # point w with schp[w] = w. Thus defined, trans has all the properties we 
#   # expect of it.

  ss := Orbits( Stabilizer( G, omega ), [1..n] );
  ll := List( ss, Length );
# end;

#  Suborbits( 1 );

  # Test for A_m \leq G. Also traps 4-trans gps: M11, M12, M23, M24
  if Length( ll ) = 3 and m > 4 then
    # Care if m = 7 since then 2*m-4 = n - 1 - 2*m-4 --- there are two 
    # suborbits of equal size.
    if m = 7 then 
      if Collected( ll ) = [[1,1], [10,2]] then
        inds := Filtered( [1..3], i -> ll[i] = 10 );
        for ind in inds do 
          L := Set( ss[ind] );
          AddSet( L, omega );
          res := TestSet( omega, L );
          if res[1] = true then return res; fi;
        od;
        return [false]; 
      else
        return [false];
      fi;
    else 
      ind := PositionProperty( [1..3], i -> ll[i] = 2*m - 4 );
      if ind <> false then 
         L := Set( ss[ind] );
         AddSet( L, omega );
         return TestSet( omega, L );
      else
         return [false];
      fi;
    fi;

  # Special cases
  elif m = 2 then return [false]; # Here n=1 and G is trivial. We don't 
                                  # count the trivial group on 2 points
                                  # as a solution. 
  elif m = 3 then # Here n = m and there is always an isomorphism. 
                  # The only transitive groups on 3 points are A(3), S(3)
    return [true, [2,3], G.generators, [[2,3],[1,3],[1,2]]];
  elif m = 4 then # Could do by affine bit but save time.
                  # By a simple argument on orders, the
                  # only possible groups are A(4) and S(4)
    if Collected( ll ) = [[1,2],[2,2]] or 
       Collected( ll ) = [[1,2],[4,1]] then
      # A(4), S(4) (resp). 
      # L is all apart from the singleton that isn't omega.
      # Note that m=4 is something of a special case for TestSet      
      L := Set( Concatenation( List( Filtered( [1..Length( ll )],
                                 i -> ll[i] <> 1 or ss[i] = [omega] ),
                                 i -> ss[i] ) ) );
      return TestSet( omega, L );
    else
      return [false];
    fi;

  # All 2-T gps of degree 6 contain PSL(2,5) or A6 as a normal subgroup
  # elif m = 6 then ;
  # All 2-T gps of deg 10 contain PSL(2,9) = A6 or A10 as a normal subgroup
  # elif m = 10 then ;

  # degree 11
  # M11 is 4-trans so caught above, so just PSL(2,11) and AGL(1,11) here
  elif m = 11 and Collected( ll ) = [[1,1],[3,2],[6,4],[12,2]] then 
    # PSL(2,11) Figures from Cameron + Cannon
    is := Filtered( [1..Length( ll )], i -> ll[ i ] = 12 );
    js := Filtered( [1..Length( ll )], j -> ll[ j ] = 6 );

    for i in is do
      L := Set( ss[ i ] );
      AddSet( L , omega );
      for j in js do
        LL := ShallowCopy( L );
        UniteSet( LL, ss[ j ] );
        res := TestSet( omega, LL );
        if res[1] = true then return res; fi;
      od;
    od;
    return [false];

  # degree 12
  # M12 is 4-trans, M11 only 3-trans on 12 pts. Also PSL(2,11)\leq G.
  elif m = 12 and Set( ll ) = [1,15,20,30] then 
    # M11
    return TestThreeTrans();

  # degree 15
  # Just PSL(4,2) and its subgroup A7 here. TestPSL adjusted to cope.
  # elif m = 15 then ;

  # Degree 22. Just M22, Aut(M22) here
  elif m = 22 then 
    if Set( ll ) = [1,30,40,160] then
      # M22, Aut(M22)
      return TestThreeTrans();
    else
      return [false];
    fi;

  # Degree 23
  # Just AGL(1,23) and the 4-trans M23 here, so no special cases
  #  elif m = 23 then ;

  # Degree 24
  # Just PSL(2,23) \leq G and the 4-trans M24 here, so no special cases 
  # elif m = 24 then ;

  # Degree 28
  # Here there are PSigmaL(2,8) and PSp(6,2) (which has one orbit fewer than 
  # standard for a PSp(2d,2) group) in addition to the 'normal' 
  # PSL(2,27)\leq G and PSU(3,3) \leq G. TestSympB adjusted accordingly.
  # Call TestSympB explicitly only when needed since otherwise TestSympB
  # will not pass
  # control on to LR1.
  elif m = 28 and Collected( ll ) = [[1,1],[5,1],[20,1],[32,1],
                     [40,2],[80,1],[160,1]] then
    return TestSympB( 3 );

  # Degree 176
  # Only A176 and HS
  elif m = 176 then 
    # HS. Figures from Cameron + Cannon
    is := Filtered( [1..Length( ll )], i -> ll[ i ] = 180);
    js := Filtered( [1..Length( ll )], j -> ll[ j ] = 144 );
    ks := Filtered( [1..Length( ll )], j -> ll[ j ] = 24 );

    for i in is do
      L := Set( ss[ i ] );
      AddSet( L , omega );
      for j in js do
        LL := ShallowCopy( L );
        UniteSet( LL, ss[ j ] );
        for k in ks do
          LLL := ShallowCopy( LL );
          UniteSet( LLL, ss[ k ] );
          res := TestSet( omega, LLL );
          if res[1] = true then return res; fi;
        od;
      od;
    od;
    return [false];

  # Degree 276
  # only A276 and Co3.
  elif m = 276 then 
     # Co3. Figures from Cameron + Cannon
    is := Filtered( [1..Length( ll )], i -> ll[ i ] = 324 );
    js := Filtered( [1..Length( ll )], j -> ll[ j ] = 224 );

    for i in is do
      L := Set( ss[ i ] );
      AddSet( L , omega );
      for j in js do
        LL := ShallowCopy( L );
        UniteSet( LL, ss[ j ] );
        res := TestSet( omega, LL );
        if res[1] = true then return res; fi;
      od;
    od;
    return [false];

  else
    # Symplectic degrees not already covered
    # Solve x^2 + x - 2m = 0 and x = 2^{d}, d>2
    # and   x^2 - x - 2m = 0 and x = 2^{d}, d>2
    sqrt := RootInt( 1 + 8 * m );
    if sqrt * sqrt = 1 + 8 * m then
      ppresA := PrimePower( (sqrt - 1) / 2 );
      ppresB := PrimePower( (sqrt + 1) / 2 );
      if ppresA <> false and ppresA[1] = 2 and ppresA[2] > 2 then
        return TestSympA( ppresA[2] );
      # If m=28 we have already called TestSympB if necessary, if execution
      # has reached here we should be calling PSL and LR1. Hence the >3
      # in the next line.
      elif ppresB <> false and ppresB[1] = 2 and ppresB[2] > 3 then
        return TestSympB( ppresB[2] );
      fi;      
    fi;

    # Affine bit
    ppres := PrimePower( m );
    if ppres <> false then
      res := TestAffine( ppres[1], ppres[2] );
      if res[1] = true then return res; fi;
    fi;

    # PSL bit
    # for all prime powers q < sqrt(m) check for d with (q^d-1)/(q-1)=m
    # Looks inefficient but remember everything is very small

    for q in Filtered( [2..( RootInt( m ) + 1 )], IsPrimePowerInt ) do
      x := q^2 + q + 1;
      d := 3;
      while x < m do
        x := q * x + 1;
        d := d + 1;
      od;
      if x = m then 
        res := TestPSL( d, q );
        if res[1] = true then return res; fi;
      fi;
    od;

    # Lie rank 1 bit

    if IsPrimePowerInt( m - 1 ) then
      return TestLR1( m - 1 );
    fi;
    
    return [false];
  fi;

end;
    


UOP := function(G, D, op)
# 'Visible' function whose purpose is to turn the action into one on [1..n]
# and call UOPcore, then translate the result back again.

local G2, n, res, w, res2;

  n := Length(D);
  if not IsPermGroup( G ) or not ForAll( D, IsInt ) or op <> OnPoints then
    G2 := Operation( G, D, op );
    res := UOPcore( G2, n );
    if res[1] then 
       res2 := [];
       res2[1] := true;
       res2[2] := [];
       for w in res[2] do 
         AddSet( res2[2], D[w] );
       od;
       return res2;
       # Rewrite to return in a different format?
    fi;
  else
    return UOPcore( G, n );
  fi;
end;