// Functions used in paper 11:
// _Computing the primitive permutation groups of degree less than 1000_
// by Colva M. Roney-Dougal and William R. Unger

// Section 3: Determining conjugacy

// Subsection 3.1: Proving that groups are not conjugate

NotConjTest := function(H, K)
  if #H ne #K or Transitivity(H) ne Transitivity(K) then
    return true, false;
  end if;
  if H eq K then
    return true, true;
  end if;
  k := Transitivity(K);
  Hs := {* t[1] : t in OrbitRepresentatives(BasicStabilizer(H, k + 1)) *};
  Ks := {* t[1] : t in OrbitRepresentatives(BasicStabilizer(K, k + 1)) *};
  if Hs ne Ks then
    return true, false;
  end if;
  h_solv := IsSolvable(H);
  k_solv := IsSolvable(K);
  if not h_solv eq k_solv then
    return true, false;
  end if;
  if h_solv then
    h_classes := Classes(H); k_classes := Classes(K);
    if not #h_classes eq #k_classes then
      return true, false;
    end if;
    classlength := #h_classes;
    for i in [1..classlength] do
      if not exists(t){ cl : cl in k_classes |
                    h_classes[i][2] eq cl[2] and
                    CycleStructure(h_classes[i][3]) eq CycleStructure(cl[3]) }
      then
        return true, false;
      end if;
    end for;
    if #H le 1000 and not IsIsomorphic(PCGroup(H), PCGroup(K)) then
      return true, false;
    end if;
  end if;
  return false, _;
end function;

// Subsection 3.2: Proving that groups are conjugate

ConjSemilin := function(H, K, phi)
  n := Degree(H); assert n eq Degree(K);
  p := #BaseRing(H); assert p eq #BaseRing(K);
  glp := GL(n, p) @ phi;
  is_abs_irred_h, C1, ext_h := IsAbsolutelyIrreducible(H);
  is_abs_irred_k, C2, ext_k := IsAbsolutelyIrreducible(K);
  if not (is_abs_irred_h cmpeq is_abs_irred_k) then
    return false, _;
  end if;
  if not is_abs_irred_h then
    if not ext_h eq ext_k then return false, _; end if;
    C1 := GL(n, p) ! C1; C2 := GL(n, p) ! C2;
  end if;
  if is_abs_irred_h then
    C1 := GL(n, p) ! CentralisingMatrix(H);
    C2 := GL(n, p) ! CentralisingMatrix(K);
    if not DegreeOfFieldExtension(H) eq DegreeOfFieldExtension(K) then
      return false, _;
    end if;
  end if;
  o1 := Order(C1); o2 := Order(C2);
  gcd := Gcd(o1, o2);
  C1 := C1^(o1 div gcd); C2 := C2^(o2 div gcd);
  c1grp := sub< GL(n, p) | C1 >;
  c2grp := sub< GL(n, p) | C2 >;
  centralisers_conj := false;
  for elt in c2grp do
    if Order(elt) eq Order(C1) then
      centralisers_conj, x := IsSimilar(elt, C1);
      if centralisers_conj then
        x := GL(n, p) ! x;
        break;
      end if;
    end if;
  end for;
  if not centralisers_conj then
    return false, _;
  end if;
  n := Normaliser(glp, c1grp @ phi);
  conj_in_gaml, y := IsConjugate(n, H @ phi, ((K^x) @ phi));
  if not conj_in_gaml then
    return false, _;
  else
    x2 := (y @@ phi)*x^-1;
    return true, x2;
  end if;
end function;

// The next three functions (MakeImprim, StandardiseBasis, ConjImprim)
// were not given in the text, although the algorithm used by ConjImprim
// is sketched there.  The other two are used by it.

MakeImprim := function(group, a, process)
  i := 0;
  for i in [1..100] do
    if IsPrimitive(group^a) cmpeq false then
      return true, a;
    end if;
    a := Random(process);
  end for;
  /* If we reach here, it has failed */
  return false, _;
end function;

StandardiseBasis := function(listofblocks, d, p)
  nmr_of_blocks := #listofblocks;

  /* The following line produces a single list of vectors */
  list := &cat[ Basis(listofblocks[i]) : i in [1..nmr_of_blocks] ];

  /* The following line turns the list of vectors into list of rows
     in a matrix */
  mat := GL(d, p) ! &cat[ [ list[i][j] : j in [1..d] ] : i in [1..d] ];

  /* We will need  to conjugate by the inverse of that matrix */
  return mat^-1;
end function;

ConjImprim := function(group1, group2, hom)

  d := Degree(group1); assert d eq Degree(group2);
  p := #BaseRing(group1); assert p eq #BaseRing(group2);
  process := RandomProcess(GL(d, p));
  phi, glp:= OrbitAction(GL(d, p), RSpace(GL(d, p)).1);

  /*
   * We start by organising group1
   */

  /* First we ensure that we know that group1 is imprimitive -
   * a is the element that we will have conjugated by to
   * make IsPrimitive return false
   */

  x, a := MakeImprim(group1, Identity(GL(d, p)), process);
  if not x then          /* Unable to get IsPrimitive to return false */
    "can't make group1 imprimitive";
    return false, _;
  end if;

  /* now we find a system of impritivity and the induced action
   * on the set of blocks, and conjugate group1^a so that it
   * lies in a standard copy of the full wreath product preserving
   * this block system
   */

  new_group1 := group1^a;
  if not (IsPrimitive(new_group1) cmpeq false) then
    if not (IsPrimitive(new_group1) cmpeq false) then
      "odd thing number 1 happened";
      return false, _;
    end if;
  end if;

  assert not IsPrimitive(new_group1);
  blocks1 := Blocks(new_group1);
  number_blocks := #blocks1;
  permgroup1 := BlocksImage(new_group1);
  aprime := StandardiseBasis(blocks1, d, p);
  standardised_group1 := group1^aprime;
  overgroup := WreathProduct(GL(d div number_blocks, p), Sym(number_blocks));

  /* standardise_group1 is a
   * subgroup of GL(m, p) wr Sym(d/p) where
   * number_blocks = d/p
   */

  gens := Generators(standardised_group1);
  for g in gens do
    if not g in overgroup then
      "failed to standardise group1 correctly";
      return false, _;
    end if;
  end for;

  /********************************************************************
    We now start working on group2. I have no idea how to prove that 100
    times is likely to be enough, but small amounts of experimentation
    suggested that this really does work most of the time.
   ********************************************************************/

  b := Identity(GL(d, p));
  for i in [1..20] do
    /* This is to ensure that we
     * get a different copy of group2 each
     * time
     */
    conjelt := Random(process);
    group2_temp := group2^conjelt;

    /* This is to ensure that we have a system of imprimitivity
     */
    x, b := MakeImprim(group2_temp, b, process);

    if not x then
      /* Unable to make IsPrimitive return false on group2_temp */
      "failed to make group2 false";
      return false, _;
    end if;

    new_group2 := group2_temp^b;
    if not (IsPrimitive(new_group2) cmpeq false) then
      "IsPrimitive behaving unpredictably, or block sizes wrong";
      continue;
    end if;

    blocks2 := Blocks(new_group2);
    if not (#blocks2 eq number_blocks) then
      continue;
    end if;

    permgroup2 := BlocksImage(new_group2);
    if not IsConjugate(Sym(number_blocks), permgroup1, permgroup2) then
      "wrong action on blocks";
      continue;    /* Have found wrong block system */
    end if;

    /*******************************************************************
       if we reach this far then we hopefully have the right system of
       imprimitivity for group2_thisrun
    ********************************************************************/

    bprime := StandardiseBasis(blocks2, d, p);
    standardised_group2 := new_group2^bprime;
    gens := Generators(standardised_group2);
    for g in gens do
      if not g in overgroup then
        "failed to standardise group2 correctly";
      end if;
    end for;
    x, y := IsConjugate(overgroup @ hom, standardised_group1 @ hom,
                                          standardised_group2 @ hom);

    if x then                  /* FOUND IT!! */
      return true, (a*aprime*(y @@ hom)*bprime^-1*b^-1*conjelt^-1);
    end if;

  end for;
  // "ran 20 times - probably not conjugate";
  return false, _;
end function;

// Subsection 3.3: The main conjugacy algorithm

IsGLConjugate := function(H, K, phi)
  a,b := NotConjTest(H, K);
  if a then
    if b then
      return true, H!1;
    else
      return false, _;
    end if;
  end if;
  h_mat := H @@ phi;
  k_mat := K @@ phi;
  n := Degree(h_mat); assert n eq Degree(k_mat);
  p := #BaseRing(h_mat); assert p eq #BaseRing(k_mat);
  glp := GL(n, p) @ phi;
  hsemi_lin := IsSemiLinear(h_mat);
  ksemi_lin := IsSemiLinear(k_mat);
  if not ((hsemi_lin cmpeq ksemi_lin) or (hsemi_lin cmpeq
      "unknown") or (ksemi_lin cmpeq "unknown"))
  then
    return false, _;
  end if;
  if (hsemi_lin cmpeq true) and (ksemi_lin cmpeq true) then
    boolean, conjelt := ConjSemilin(h_mat, k_mat, phi);
    if boolean then return true, conjelt @ phi; end if;
  end if;
  h_imprim := IsPrimitive(h_mat);
  k_imprim := IsPrimitive(k_mat);
  if not ((h_imprim cmpeq k_imprim) or (h_imprim cmpeq
      "unknown") or (k_imprim cmpeq "unknown"))
  then
    return false, _;
  end if;
  if (h_imprim cmpeq false) and (k_imprim cmpeq false) then
    boolean, conjelt := ConjImprim(h_mat, k_mat, phi);
    if boolean then return true, conjelt @ phi; end if;
  end if;
  HR := Radical(H);
  KR := Radical(K);
  if #HR ne #KR then return false, _; end if;
  if #HR gt 1 and #HR lt #H then
    a, b := $$(HR, KR, phi);
    if not a then return false, _; end if;
    H := H^b;
    if H eq K then return true, b; end if;
    c, d := IsConjugate(Normaliser(glp, KR), H, K);
    if c then return true, b*d; else return false, _; end if;
  elif #HR eq 1 then
    HR := Socle(H);
    KR := Socle(K);
    if H ne HR then
      a, b := $$(HR, KR, phi);
      if not a then return false, _; end if;
      H := H^b;
      if H eq K then return true, b; end if;
      c, d := IsConjugate(Normaliser(glp, KR), H, K);
      if c then return true, b*d; else return false, _; end if;
    else
      c, d := IsConjugate(glp, H, K);
      if c then return true, d; else return false, _; end if;
    end if;
  else
    dsH := DerivedSeries(H);  dsK := DerivedSeries(K);
    len := #dsH;
    if len ne #dsK then return false, _; end if;
    if exists{ i: i in [1..len] | #dsH[i] ne #dsK[i] } then
      return false, _;
    end if;
    if len le 2 then
      c, d := IsConjugate(glp, H, K);
      if c then return c, d; else return false, _; end if;
    else
      HR := dsH[len -1];  KR := dsK[len -1];
      a, b := IsConjugate(glp, HR, KR);
      if not a then return false, _; end if;
      H := H^b;
      if H eq K then return true, b; end if;
      c, d := IsConjugate(Normaliser(glp, KR), H, K);
      if c then return true, b*d; else return false, _; end if;
    end if;
  end if;
end function;