ResGroups.mag

import "Utils.mag": dedupe_conjugage_subgroups, has_root, multiplicities_to_information, dedupe_conjugage_subgroups_simple, Z, not_implemented;

declare type PGGAlg_ResGroups: PGGAlg;
declare type PGGAlgState_ResGroups: PGGAlgState;
declare attributes PGGAlgState_ResGroups: overgroup_embedding, upto;

declare type PGGAlg_ResGroups_All: PGGAlg_ResGroups;
declare attributes PGGAlg_ResGroups_All: subgroup_choice, statistic, is_subgroup_of, dedupe;

declare type PGGAlgState_ResGroups_All: PGGAlgState_ResGroups;
declare attributes PGGAlgState_ResGroups_All: galois_conjugacy, resolvent_overgroup, possible_groups, choice_state;

declare type PGGAlg_ResGroups_RootsMaximal: PGGAlg_ResGroups;
declare attributes PGGAlg_ResGroups_RootsMaximal: dedupe;

declare type PGGAlgState_ResGroups_RootsMaximal: PGGAlgState_ResGroups;
declare attributes PGGAlgState_ResGroups_RootsMaximal: galois_conjugacy, resolvent_overgroup, cur_group, possible_subgroups, conjugacy_classes;

declare type PGGAlg_ResGroups_Sequence: PGGAlg_ResGroups;
declare attributes PGGAlg_ResGroups_Sequence: seq;

declare type PGGAlgState_ResGroups_Sequence: PGGAlgState_ResGroups;
declare attributes PGGAlgState_ResGroups_Sequence: idx, cur_state, prev_state;

declare type PGGAlg_ResGroups_Null: PGGAlg_ResGroups;

declare type PGGAlgState_ResGroups_Null: PGGAlgState_ResGroups;
declare attributes PGGAlgState_ResGroups_Null: previous;

declare type PGGAlg_ResGroups_Maximal: PGGAlg_ResGroups;
declare attributes PGGAlg_ResGroups_Maximal: subgroup_choice, statistic, descend_when, descend, useful, reprocess, reset, dedupe;

declare type PGGAlgState_ResGroups_Maximal: PGGAlgState_ResGroups;
declare attributes PGGAlgState_ResGroups_Maximal: galois_conjugacy, resolvent_overgroup, choice_state, graph, pool, resolvents, useful_is_special, actually_useful_is_special;

declare type PGG_ResGroups_Graph[PGG_ResGroups_GraphNode];
declare attributes PGG_ResGroups_Graph: galois_conjugacy, dedupe, conjugacy_classes, id_to_node, class_to_id;
declare attributes PGG_ResGroups_GraphNode: parent, class, id, children, possible_children, parents, possible_parents, maybe_equal, maybe_subgroup;

declare type PGGAlgState_ResGroups_Maximal_Graph[PGGAlgState_ResGroups_Maximal_GraphNode]: PGG_ResGroups_Graph;
declare attributes PGGAlgState_ResGroups_Maximal_Graph: state;
declare type PGGAlgState_ResGroups_Maximal_GraphNode: PGG_ResGroups_GraphNode;
declare attributes PGGAlgState_ResGroups_Maximal_GraphNode: pooled_children;

declare type PGGAlg_ResGroups_Maximal2: PGGAlg_ResGroups;
declare attributes PGGAlg_ResGroups_Maximal2: subgroup_choice, statistic, reset, dedupe, descend;

declare type PGGAlgState_ResGroups_Maximal2: PGGAlgState_ResGroups;
declare attributes PGGAlgState_ResGroups_Maximal2: galois_conjugacy, resolvent_overgroup, choice_state, graph, pool;

intrinsic PGGAlg_ResGroups_All_Make(:Statistic:=false, SubgroupChoice:=false, IsSubgroupOf:=false, Dedupe:=false) -> PGGAlg_ResGroups_All
  {The "All" resolvent groups algorithm.}
  // IsSubgroupOf := OrbitImage(G, {v:v in VectorSpace(G)|v ne 0}) where G:=CSp(4,3);
  alg := New(PGGAlg_ResGroups_All);
  alg`statistic := Statistic cmpne false select Statistic else PGGStat_FactorDegrees_Make();
  alg`subgroup_choice := SubgroupChoice cmpne false select SubgroupChoice else PGGAlg_ResChoice_Tranche_Make();
  alg`is_subgroup_of := IsSubgroupOf;
  alg`dedupe := PGGAlg_SubgrpDedupe_Make(Dedupe);
  return alg;
end intrinsic;

intrinsic PGGAlg_ResGroups_RootsMaximal_Make( : Dedupe:=false) -> PGGAlg_ResGroups_RootsMaximal
  {The "RootsMaximal" resolvent groups algorithm.}
  alg := New(PGGAlg_ResGroups_RootsMaximal);
  alg`dedupe := PGGAlg_SubgrpDedupe_Make(Dedupe);
  return alg;
end intrinsic;

intrinsic PGGAlg_ResGroups_Sequence_Make(seq) -> PGGAlg_ResGroups_Sequence
  {The "Sequence" resolvent groups algorithm.}
  alg := New(PGGAlg_ResGroups_Sequence);
  alg`seq := seq;
  return alg;
end intrinsic;

intrinsic PGGAlg_ResGroups_Null_Make() -> PGGAlg_ResGroups_Null
  {The "Null" resolvent groups algorithm.}
  alg := New(PGGAlg_ResGroups_Null);
  return alg;
end intrinsic;

intrinsic PGGAlg_ResGroups_Maximal_Make(:Statistic:=false, SubgroupChoice:=false, DescendWhen:=false, Descend:=false, Useful:=false, Reprocess:=true, Reset:=true, Dedupe:=true) -> PGGAlg_ResGroups_Maximal
  {The "Maximal" resolvent groups algorithm.}
  alg := New(PGGAlg_ResGroups_Maximal);
  alg`statistic := Statistic cmpne false select Statistic else PGGStat_FactorDegrees_Make();
  alg`subgroup_choice := SubgroupChoice cmpne false select SubgroupChoice else PGGAlg_ResChoice_Tranche_Make();
  alg`descend_when := DescendWhen cmpne false select DescendWhen else "Sufficient";
  alg`descend := Descend cmpne false select Descend else "All";
  alg`useful := Useful cmpne false select Useful else "Sufficient";
  alg`reprocess := Reprocess;
  alg`reset := Reset;
  alg`dedupe := PGGAlg_SubgrpDedupe_Make(Dedupe);
  return alg;
end intrinsic;

intrinsic PGGAlg_ResGroups_Maximal2_Make(:Statistic:=false, SubgroupChoice:=false, Reset:=true, Dedupe:=true, Descend:=false) -> PGGAlg_ResGroups_Maximal2
  {The "Maximal2" resolvent groups algorithm.}
  alg := New(PGGAlg_ResGroups_Maximal2);
  alg`statistic := Statistic cmpne false select Statistic else PGGStat_FactorDegrees_Make();
  alg`subgroup_choice := SubgroupChoice cmpne false select SubgroupChoice else PGGAlg_ResChoice_Tranche_Make();
  alg`reset := Reset;
  alg`dedupe := PGGAlg_SubgrpDedupe_Make(Dedupe);
  alg`descend := Descend cmpne false select Descend else "Strict";
  return alg;
end intrinsic;

intrinsic Print(alg :: PGGAlg_ResGroups_All)
  {Print.}
  print "all";
  IndentPush();
  printf "subgroup choice = "; Print(alg`subgroup_choice); print "";
  printf "statistic = "; Print(alg`statistic); print "";
  printf "dedupe = "; Print(alg`dedupe);
  IndentPop();
end intrinsic;

intrinsic Print(alg :: PGGAlg_ResGroups_RootsMaximal)
  {"}
  print "roots maximal";
  IndentPush();
  print "dedupe = "; Print(alg`dedupe);
  IndentPop();
end intrinsic;

intrinsic Print(alg :: PGGAlg_ResGroups_Sequence)
  {"}
  if #alg`seq eq 0 then
    printf "sequence (null)";
  else
    print "sequence";
    IndentPush();
    for i in [1..#alg`seq] do
      Print(alg`seq[i]);
      if i lt #alg`seq then
        print "";
      end if;
    end for;
    IndentPop();
  end if;
end intrinsic;

intrinsic Print(alg :: PGGAlg_ResGroups_Null)
  {"}
  printf "null";
end intrinsic;

intrinsic Print(alg :: PGGAlg_ResGroups_Maximal)
  {"}
  print "maximal";
  IndentPush();
  printf "subgroup choice = "; Print(alg`subgroup_choice); print "";
  printf "statistic = "; Print(alg`statistic); print "";
  print "descend when =", alg`descend_when;
  print "descend =", alg`descend;
  print "useful =", alg`useful;
  print "reprocess =", alg`reprocess;
  print "reset =", alg`reset;
  printf "dedupe = "; Print(alg`dedupe);
  IndentPop();
end intrinsic;

intrinsic Print(alg :: PGGAlg_ResGroups_Maximal2)
  {"}
  print "maximal2";
  IndentPush();
  printf "subgroup choice = "; Print(alg`subgroup_choice); print "";
  printf "statistic = "; Print(alg`statistic); print "";
  print "descend =", alg`descend;
  print "reset =", alg`reset;
  printf "dedupe = "; Print(alg`dedupe);
  IndentPop();
end intrinsic;

intrinsic Start(alg :: PGGAlg_ResGroups_All, emb :: PGGHomGrpPerm, upto :: PGGAlg_GaloisGroup_ARM_UpTo) -> PGGAlgState_ResGroups_All
  {Starts the algorithm and returns its state.}
  s := New(PGGAlgState_ResGroups_All);
  s`algorithm := alg;
  s`overgroup_embedding := emb;
  s`upto := Start(upto, emb);
  s`galois_conjugacy := Domain(emb);
  s`resolvent_overgroup := Codomain(emb);
  // get a list of possible groups
  c := s`galois_conjugacy;
  if alg`is_subgroup_of cmpne false then
    print "subgroups ...";
    Gs := [x`subgroup : x in Subgroups(alg`is_subgroup_of : IsSolvable, IsTransitive:=GalGrpIsTransitive(c))];
    print "embeddings ...";
    Gs := [e[1] : e in Embeddings(G, c), G in Gs];
    print "consistent ...";
    Gs := [G : G in Gs | GalGrpIsConsistent(c, Gs)];
    print "dedupe ...";
    Gs := [Rep(G) : G in {Classes(alg`dedupe, c)| G : G in Gs}];
    print "done!";
    s`possible_groups := Gs;
  else
    PGG_GlobalTimer_Push("subgroups");
    Gs := [x`subgroup : x in Subgroups(GalGrpOvergroup(c) : IsSolvable, IsTransitive:=GalGrpIsTransitive(c))];
    PGG_GlobalTimer_Swap("consistent");
    Gs := [G : G in Gs | GalGrpIsConsistent(c, G)];
    if GalGrpOvergroup(c) ne Group(c) then
      PGG_GlobalTimer_Swap("dedupe");
      Gs := [Rep(G) : G in {Classes(alg`dedupe, c)| G : G in Gs}];
    end if;
    PGG_GlobalTimer_Pop();
    s`possible_groups := Gs;
  end if;
  s`choice_state := Start(alg`subgroup_choice, s);
  vprint PGG_GaloisGroup: "#possible groups =", #s`possible_groups;
  return s;
end intrinsic;

intrinsic Start(alg :: PGGAlg_ResGroups_RootsMaximal, emb :: PGGHomGrpPerm, upto :: PGGAlg_GaloisGroup_ARM_UpTo) -> PGGAlgState_ResGroups_RootsMaximal
  {"}
  s := New(PGGAlgState_ResGroups_RootsMaximal);
  s`algorithm := alg;
  s`overgroup_embedding := emb;
  s`upto := Start(upto, emb);
  s`galois_conjugacy := Domain(emb);
  s`resolvent_overgroup := Codomain(emb);
  s`conjugacy_classes := Classes(alg`dedupe, s`galois_conjugacy);
  s`cur_group := s`conjugacy_classes ! GalGrpOvergroup(s`galois_conjugacy);
  return s;
end intrinsic;

intrinsic Start(alg :: PGGAlg_ResGroups_Sequence, emb :: PGGHomGrpPerm, upto :: PGGAlg_GaloisGroup_ARM_UpTo) -> PGGAlgState_ResGroups_Sequence
  {"}
  s := New(PGGAlgState_ResGroups_Sequence);
  s`algorithm := alg;
  s`overgroup_embedding := emb;
  s`upto := Start(upto, emb);
  s`idx := 0;
  return s;
end intrinsic;

intrinsic Start(alg :: PGGAlg_ResGroups_Null, emb :: PGGHomGrpPerm, upto :: PGGAlg_GaloisGroup_ARM_UpTo) -> PGGAlgState_ResGroups_Null
  {"}
  s := New(PGGAlgState_ResGroups_Null);
  s`algorithm := alg;
  s`overgroup_embedding := emb;
  s`upto := Start(upto, emb);
  return s;
end intrinsic;

intrinsic Start(alg :: PGGAlg_ResGroups_Maximal, emb :: PGGHomGrpPerm, upto :: PGGAlg_GaloisGroup_ARM_UpTo) -> PGGAlgState_ResGroups_Maximal
  {"}
  s := New(PGGAlgState_ResGroups_Maximal);
  s`algorithm := alg;
  s`overgroup_embedding := emb;
  s`upto := Start(upto, emb);
  s`galois_conjugacy := Domain(emb);
  s`resolvent_overgroup := Codomain(emb);
  s`choice_state := Start(alg`subgroup_choice, s);
  s`graph := PGGAlgState_ResGroups_Maximal_Graph_Make(s);
  s`pool := { Root(s`graph) };
  s`resolvents := [**];
  s`useful_is_special := false;
  s`actually_useful_is_special := false;
  return s;
end intrinsic;

intrinsic Start(alg :: PGGAlg_ResGroups_Maximal2, emb :: PGGHomGrpPerm, upto :: PGGAlg_GaloisGroup_ARM_UpTo) -> PGGAlgState_ResGroups_Maximal2
  {"}
  s := New(PGGAlgState_ResGroups_Maximal2);
  s`algorithm := alg;
  s`overgroup_embedding := emb;
  s`upto := Start(upto, emb);
  s`galois_conjugacy := Domain(emb);
  s`resolvent_overgroup := Codomain(emb);
  s`choice_state := Start(alg`subgroup_choice, s);
  s`graph := PGGAlgState_ResGroups_Maximal_Graph_Make(s);
  s`pool := { Root(s`graph) };
  // s`subgroup_classes := Classes(alg`dedupe);
  return s;
end intrinsic;

intrinsic RestartFromScratch(alg :: PGGAlg_ResGroups, s :: PGGAlgState_ResGroups) -> PGGAlgState_ResGroups
  {Same as ``Start(alg, s`overgroup_embedding)`` but with a warning that it's happening.}
  vprint PGG_GaloisGroup: "*** starting ResGroups algorithm from scratch:", alg;
  return Start(alg, s`overgroup_embedding, s`upto`algorithm);
end intrinsic;

intrinsic AllowRestartFromScratch(a :: PGGAlg_ResGroups, b :: PGGAlg_ResGroups) -> BoolElt
  {True if we allow restarting an instance of b from a.}
  // by default, we allow it if the algorithms are of different type
  return Type(a) ne Type(b);
end intrinsic;

intrinsic Restart(alg :: PGGAlg_ResGroups, s :: PGGAlgState_ResGroups) -> PGGAlgState_ResGroups
  {Like ``Start(alg, s`overgroup_embedding)`` but tries to preserve information.}
  if AllowRestartFromScratch(alg, s`algorithm) then
    return RestartFromScratch(alg, s);
  else
    error "not implemented: Restart:", Type(alg);
  end if;
end intrinsic;

intrinsic Restart(alg :: PGGAlg_ResGroups_Null, s :: PGGAlgState_ResGroups) -> PGGAlgState_ResGroups_Null
  {"}
  s2 := Start(alg, s`overgroup_embedding, s`upto`algorithm);
  s2`previous := s;
  return s2;
end intrinsic;

intrinsic Restart(alg :: PGGAlg_ResGroups_Null, s :: PGGAlgState_ResGroups_Null) -> PGGAlgState_ResGroups_Null
  {"}
  return s;
end intrinsic;

intrinsic Restart(alg :: PGGAlg_ResGroups, s :: PGGAlgState_ResGroups_Null) -> PGGAlgState_ResGroups
  {"}
  if assigned s`previous then
    return Restart(alg, s`previous);
  else
    return RestartFromScratch(alg, s);
  end if;
end intrinsic;

intrinsic Restart(alg :: PGGAlg_ResGroups_All, s :: PGGAlgState_ResGroups_All) -> PGGAlgState_ResGroups_All
  {"}
  s`algorithm := alg;
  s`choice_state := Restart(alg`subgroup_choice, s`choice_state);
  return s;
end intrinsic;

intrinsic Restart(alg :: PGGAlg_ResGroups_Maximal2, s :: PGGAlgState_ResGroups_Maximal2) -> PGGAlgState_ResGroups_Maximal2
  {"}
  error if s`algorithm`dedupe ne alg`dedupe, "not implemented: Restart: PGGAlg_ResGroups_Maximal2: change dedupe";
  s`algorithm := alg;
  s`choice_state := Restart(alg`subgroup_choice, s`choice_state);
  return s;
end intrinsic;

intrinsic Restart(alg :: PGGAlg_ResGroups_Sequence, s :: PGGAlgState_ResGroups_Sequence) -> PGGAlgState_ResGroups_Sequence
  {"}
  s`algorithm := alg;
  s`idx := 0;
  return s;
end intrinsic;

intrinsic Init(X :: PGG_ResGroups_Graph, c :: PGGConj, Dedupe :: PGGAlg_SubgrpDedupe)
  {Initializes X.}
  X`galois_conjugacy := c;
  X`dedupe := Dedupe;
  X`conjugacy_classes := Classes(Dedupe, c);
  X`id_to_node := [];
  X`class_to_id := AssociativeArray();
end intrinsic;

intrinsic Init(X :: PGGAlgState_ResGroups_Maximal_Graph, s :: PGGAlgState_ResGroups)
  {Initializes X.}
  Init(X, s`galois_conjugacy, s`algorithm`dedupe);
  X`state := s;
end intrinsic;

intrinsic PGGAlgState_ResGroups_Maximal_Graph_Make(s :: PGGAlgState_ResGroups_Maximal) -> PGGAlgState_ResGroups_Maximal_Graph
  {The graph for s.}
  X := New(PGGAlgState_ResGroups_Maximal_Graph);
  Init(X, s);
  return X;
end intrinsic;

intrinsic PGGAlgState_ResGroups_Maximal_Graph_Make(s :: PGGAlgState_ResGroups_Maximal2) -> PGGAlgState_ResGroups_Maximal_Graph
  {The graph for s.}
  X := New(PGGAlgState_ResGroups_Maximal_Graph);
  Init(X, s);
  return X;
end intrinsic;

intrinsic Init(n :: PGG_ResGroups_GraphNode)
  {Initializes n.}
  return;
end intrinsic;

intrinsic Init(n :: PGGAlgState_ResGroups_Maximal_GraphNode)
  {"}
  n`pooled_children := {};
end intrinsic;

intrinsic Node(X :: PGG_ResGroups_Graph, G :: _PGGSubgrpcls) -> PGG_ResGroups_GraphNode
  {The node of X at G.}
  ok, i := IsDefined(X`class_to_id, G);
  if ok then
    return X`id_to_node[i];
  end if;
  n := New(ElementType(X));
  n`id := #X`id_to_node + 1;
  n`parent := X;
  n`class := G;
  n`possible_parents := {};
  n`maybe_equal := true;
  n`maybe_subgroup := true;
  X`class_to_id[G] := n`id;
  X`id_to_node[n`id] := n;
  Init(n);
  return n;
end intrinsic;

intrinsic Node(X :: PGG_ResGroups_Graph, G :: GrpPerm) -> PGG_ResGroups_GraphNode
  {"}
  return Node(X, X`conjugacy_classes ! G);
end intrinsic;

intrinsic Node(X :: PGG_ResGroups_Graph, G :: PGG_ResGroups_GraphNode) -> PGG_ResGroups_GraphNode
  {"}
  return G;
end intrinsic;

intrinsic IsCoercible(X :: PGG_ResGroups_Graph, G) -> BoolElt, PGG_ResGroups_Graph
  {True if G is coercible to a node of X.}
  return true, Node(X, G);
end intrinsic;

intrinsic Root(X :: PGG_ResGroups_Graph) -> PGG_ResGroups_GraphNode
  {The root node of X.}
  return Node(X, GalGrpOvergroup(X`galois_conjugacy));
end intrinsic;

intrinsic Parent(n :: PGG_ResGroups_GraphNode) -> PGG_ResGroups_Graph
  {Parent.}
  return n`parent;
end intrinsic;

intrinsic Hash(n :: PGG_ResGroups_GraphNode) -> RngIntElt
  {Hash.}
  return n`id;
end intrinsic;

intrinsic 'eq'(n1 :: PGG_ResGroups_GraphNode, n2 :: PGG_ResGroups_GraphNode) -> BoolElt
  {Equality.}
  return n1`id eq n2`id;
end intrinsic;

intrinsic Children(n :: PGG_ResGroups_GraphNode) -> {}
  {The child nodes of n.}
  if not assigned n`children then
    c := n`parent`galois_conjugacy;
    PGG_GlobalTimer_Push("subgroups");
    subgrps := [x`subgroup : x in MaximalSubgroups(Rep(n`class) : IsSolvable, IsTransitive:=GalGrpIsTransitive(c))];
    PGG_GlobalTimer_Swap("consistent");
    subgrps := [G : G in subgrps | GalGrpIsConsistentWithSubgroup(c, G)];
    PGG_GlobalTimer_Swap("nodes");
    n`children := { Node(n`parent, G) : G in subgrps };
    for m in n`children do
      m2 := m;
      // Include(~m2`parents, n);
      Include(~m2`possible_parents, n);
    end for;
    PGG_GlobalTimer_Pop();
  end if;
  return n`children;
end intrinsic;

intrinsic PossibleChildren(n :: PGG_ResGroups_GraphNode) -> {}
  {The possible children of n.}
  if not assigned n`possible_children then
    if n`maybe_subgroup then
      PGG_GlobalTimer_Push("children");
      n`possible_children := {m : m in Children(n) | m`maybe_subgroup};
      PGG_GlobalTimer_Pop();
      if #n`possible_children eq 0 and not n`maybe_equal then
        NotSubgroup(n);
      end if;
    else
      n`possible_children := {};
    end if;
    PossibleChildrenHook(n);
  end if;
  return n`possible_children;
end intrinsic;

intrinsic PossibleChildrenHook(n :: PGG_ResGroups_GraphNode)
  {Overload this to do something when `PossibleChildren(n)` is called for the first time.}
  return;
end intrinsic;

intrinsic PossibleChildrenHook(n :: PGGAlgState_ResGroups_Maximal_GraphNode)
  {"}
  if (n in n`parent`state`pool) and (not n`maybe_equal) and #PoolableChildren(n) eq 0 then
    RemoveFromPool(n);
  end if;
end intrinsic;

intrinsic NotEqual(n :: PGG_ResGroups_GraphNode)
  {Record that n is not equal to the Galois group.}
  if n`maybe_equal then
    n`maybe_equal := false;
    if assigned n`possible_children and #n`possible_children eq 0 then
      NotSubgroup(n);
    end if;
    NotEqualHook(n);
  end if;
end intrinsic;

intrinsic NotEqualHook(n :: PGG_ResGroups_GraphNode)
  {Overload this to do something when `NotEqual(n)` is called.}
  return;
end intrinsic;

intrinsic NotEqualHook(n :: PGGAlgState_ResGroups_Maximal_GraphNode)
  {"}
  if (n in n`parent`state`pool) and assigned n`possible_children and #PoolableChildren(n) eq 0 then
    RemoveFromPool(n);
  end if;
end intrinsic;

intrinsic NotSubgroup(n :: PGG_ResGroups_GraphNode)
  {Record that the Galois group is not a subgroup of n.}
  if n`maybe_subgroup then
    children := assigned n`possible_children select n`possible_children else {};
    parents := n`possible_parents;
    n`maybe_equal := false;
    n`maybe_subgroup := false;
    n`possible_children := {};
    n`possible_parents := {};
    for m in children do
      mm := m;
      Exclude(~mm`possible_parents, n);
    end for;
    for p in parents do
      pp := p;
      Exclude(~pp`possible_children, n);
    end for;
    for m in children do
      NotSubgroup(m);
    end for;
    for p in parents do
      if p`maybe_subgroup and (not p`maybe_equal) and #p`possible_children eq 0 then
        NotSubgroup(p);
      end if;
    end for;
    NotSubgroupHook(n);
  end if;
end intrinsic;

intrinsic NotSubgroupHook(n :: PGG_ResGroups_GraphNode)
  {Overload this to do something when `NotSubgroup(n)` is called.}
  return;
end intrinsic;

intrinsic NotSubgroupHook(n :: PGGAlgState_ResGroups_Maximal_GraphNode)
  {"}
  RemoveFromPool(n);
end intrinsic;

intrinsic PGG_Image(G :: GrpPerm, h :: Map) -> GrpPerm
  {The image of G under h.}
  return G @ h;
end intrinsic;

intrinsic PGG_Image(G :: GrpPerm, hs :: [Map]) -> GrpPerm, [], Map
  {The image of G under hs, the corresponding GSets, and the combined map.}
  ds := [Z| Degree(Codomain(h)) : h in hs];
  cds := [&+ds[1..i] : i in [0..#hs]];
  d := &+ds;
  Sd := SymmetricGroup(d);
  gens := [Sd| &cat[[cds[j] + k^gj : k in [1..ds[j]]] where gj:=hs[j](g) : j in [1..#hs]] where g:=G.i : i in [1..Ngens(G)]];
  G2 := sub<Sd | gens>;
  h := hom<G -> G2 | gens>;
  os := [GSet(G2, {cds[j]+1..cds[j+1]}) : j in [1..#hs]];
  return G2, os, h;
end intrinsic;

intrinsic GroupStat(alg :: PGGStat, G :: GrpPerm, h :: Map) -> .
  {"}
  return GroupStat(alg, G @ h);
end intrinsic;

intrinsic GroupStat(alg :: PGGStat, G :: GrpPerm, hs :: [Map]) -> .
  {"}
  G2, os := PGG_Image(G, hs);
  return GroupStat(alg, G2, os);
end intrinsic;

intrinsic GroupStat(n :: PGGAlgState_ResGroups_Maximal_GraphNode, h) -> .
  {The statistic of n under the action h.}
  return GroupStat(n`parent`state`algorithm`statistic, Rep(n`class) @ n`parent`state`overgroup_embedding, h);
end intrinsic;

intrinsic CurrentState(s :: PGGAlgState_ResGroups_Sequence) -> PGGAlgState_ResGroups
  {The current state.}
  if s`idx eq 0 then
    assert not assigned s`cur_state;
    PGG_GlobalTimer_Push("start groups");
    s2 := NextState(s);
    PGG_GlobalTimer_Pop();
    return s2;
  else
    return s`cur_state;
  end if;
end intrinsic;

intrinsic NextState(s :: PGGAlgState_ResGroups_Sequence) -> PGGAlgState_ResGroups
  {The next state. At the end of the sequence, returns Null.}
  assert s`idx ge 0;
  s`idx +:= 1;
  alg := AtEnd(s) select PGGAlg_ResGroups_Null_Make() else s`algorithm`seq[s`idx];
  if assigned s`cur_state then
    vprint PGG_GaloisGroup: "next groups algorithm =", alg;
    s`prev_state := s`cur_state;
    s`cur_state := Restart(alg, s`prev_state);
  else
    vprint PGG_GaloisGroup: "initial groups algorithm =", alg;
    s`cur_state := Start(alg, s`overgroup_embedding, s`upto`algorithm);
  end if;
  return s`cur_state;
end intrinsic;

intrinsic AtEnd(s :: PGGAlgState_ResGroups_Sequence) -> BoolElt
  {True if s has reached the end of the sequence.}
  return s`idx gt #s`algorithm`seq;
end intrinsic;

intrinsic IsDone(s :: PGGAlgState_ResGroups_All) -> BoolElt, GrpPerm
  {True if we have deduced the Galois group.}
  assert #s`possible_groups ge 1;
  if forall{G : G in s`possible_groups[2..#s`possible_groups] | IsEquivalent(s`upto, s`possible_groups[1], G)} then
    return true, s`possible_groups[1];
  else
    return false, _;
  end if;
end intrinsic;

intrinsic IsDone(s :: PGGAlgState_ResGroups_RootsMaximal) -> BoolElt, GrpPerm
  {"}
  if #PossibleSubgroups(s) eq 0 then
    return true, Rep(s`cur_group);
  else
    return false, _;
  end if;
end intrinsic;

intrinsic IsDone(s :: PGGAlgState_ResGroups_Sequence) -> BoolElt, GrpPerm
  {"}
  return IsDone(CurrentState(s));
end intrinsic;

intrinsic IsDone(s :: PGGAlgState_ResGroups_Null) -> BoolElt, GrpPerm
  {"}
  return false, _;
end intrinsic;

function pool_has_one_class(s)
  if #s`pool eq 0 then
    return false;
  else
    // OPTIONAL: do a linear pass through the data to eliminate conjugates
    s`pool := {pool_seq[i] : i in [1..#pool_seq] | i eq 1 or not (n eq n0 or IsConjugate(Group(s`galois_conjugacy), Rep(n`class), Rep(n0`class)) where n:=pool_seq[i] where n0:=pool_seq[i-1])} where pool_seq := SetToSequence(s`pool);
    // now see if there is one group up to UpTo
    return forall{i : i in [2..pool_seq] | IsEquivalent(s`upto, Rep(pool_seq[i]`class), Rep(n0`class))} where n0:=pool_seq[1] where pool_seq := SetToSequence(s`pool);
  end if;
end function;

intrinsic IsDone(s :: PGGAlgState_ResGroups_Maximal) -> BoolElt, GrpPerm
  {"}
  if forall{n : n in s`pool | #PossibleChildren(n) eq 0} then
    assert forall{n : n in s`pool | n`maybe_equal};
    assert #s`pool gt 0;
    if pool_has_one_class(s) then
      return true, Rep(Rep(s`pool)`class);
    end if;
  end if;
  return false;
end intrinsic;

intrinsic IsDone(s :: PGGAlgState_ResGroups_Maximal2) -> BoolElt, GrpPerm
  {"}
  if forall{n : n in s`pool | #PossibleChildren(n) eq 0} then
    assert forall{n : n in s`pool | n`maybe_equal};
    assert #s`pool gt 0;
    if pool_has_one_class(s) then
      return true, Rep(Rep(s`pool)`class);
    end if;
  end if;
  return false;
end intrinsic;

intrinsic TheGroup(s :: PGGAlgState_ResGroups) -> GrpPerm
  {The Galois group, assuming it is already known.}
  ok, G := IsDone(s);
  assert ok;
  return G;
end intrinsic;

intrinsic PossibleSubgroups(s :: PGGAlgState_ResGroups_RootsMaximal) -> []
  {The subgroups of s`cur_group which have not yet been ruled out.}
  if not assigned s`possible_subgroups then
    PGG_GlobalTimer_Push("maximal subgroups");
    Gs := [x`subgroup : x in MaximalSubgroups(Rep(s`cur_group))];
    PGG_GlobalTimer_Swap("consistent");
    Gs := [G : G in Gs | GalGrpIsConsistentWithSubgroup(s`galois_conjugacy, G)];
    PGG_GlobalTimer_Swap("dedupe");
    s`possible_subgroups := SetToSequence({s`conjugacy_classes!G : G in Gs});
    PGG_GlobalTimer_Pop();
    vprint PGG_GaloisGroup: "#possible_subgroups =", #s`possible_subgroups;
  end if;
  return s`possible_subgroups;
end intrinsic;

intrinsic ProcessResolvent(s :: PGGAlgState_ResGroups_All, resolvent, subgroup)
  {Use the resolvent to get information about the Galois group.}
  PGG_GlobalTimer_Push("resolvent stat");
  rstat := ResolventStat(s`algorithm`statistic, resolvent);
  vprint PGG_GaloisGroup: "resolvent statistic =", rstat;
  PGG_GlobalTimer_Swap("filter");
  h := PGG_CosetAction(Group(s`resolvent_overgroup), subgroup);
  stats := [<G, GroupStat(s`algorithm`statistic, G @ s`overgroup_embedding, h)> : G in s`possible_groups];
  vprint PGG_GaloisGroup: "group statistics =", {* x[2] : x in stats *};
  s`possible_groups := [x[1] : x in stats | x[2] eq rstat];
  vprint PGG_GaloisGroup: "#possible groups =", #s`possible_groups;
  PGG_GlobalTimer_Pop();
end intrinsic;

intrinsic ProcessResolvent(s :: PGGAlgState_ResGroups_RootsMaximal, resolvent :: PGGPol, subgroup :: GrpPerm)
  {"}
  PGG_GlobalTimer_Push("assert");
  poss := PossibleSubgroups(s);
  assert #poss gt 0;
  assert IsConjugate(Group(s`resolvent_overgroup), subgroup, Rep(poss[1]));
  PGG_GlobalTimer_Swap("has root");
  ok := HasRoot(resolvent);
  if ok then
    PGG_GlobalTimer_Swap("descend");
    vprint PGG_GaloisGroup: "descending";
    s`cur_group := poss[1];
    if assigned s`possible_subgroups then
      delete s`possible_subgroups;
    end if;
    vprint PGG_GaloisGroup: "index =", Index(Group(s`resolvent_overgroup), Rep(s`cur_group));
  else
    PGG_GlobalTimer_Swap("eliminate");
    vprint PGG_GaloisGroup: "eliminating";
    s`possible_subgroups := poss[2..#poss];
    vprint PGG_GaloisGroup: "#possible_subgroups =", #s`possible_subgroups;
  end if;
  PGG_GlobalTimer_Pop();
end intrinsic;

intrinsic ProcessResolvent(s ::PGGAlgState_ResGroups_Sequence, resolvent, subgroup)
  {"}
  ProcessResolvent(CurrentState(s), resolvent, subgroup);
end intrinsic;

intrinsic ProcessResolvent(s :: PGGAlgState_ResGroups_Maximal, resolvent, subgroup)
  {"}
  PGG_GlobalTimer_Push("resolvent stat");
  rstat := ResolventStat(s`algorithm`statistic, resolvent);
  vprint PGG_GaloisGroup: "resolvent statistic =", rstat;
  h := PGG_CosetAction(Group(s`resolvent_overgroup), subgroup);
  // vprint PGG_GaloisGroup: "pool statistics =", {* <GroupStat(n,h), n`maybe_equal, {*GroupStat(m,h) : m in PossibleChildren(n)*}> : n in s`pool *};
  ridx := #s`resolvents + 1;
  s`resolvents[ridx] := <resolvent, subgroup, rstat, h>;
  PGG_GlobalTimer_Swap("process");
  if s`actually_useful_is_special then
    old_pool := s`pool;
    old_children := &join[PossibleChildren(n) : n in old_pool];
    ProcessResolvent(s, ridx);
    if old_pool ne s`pool or old_children ne &join[PossibleChildren(n) : n in s`pool] then
      MarkSpecial(s`choice_state);
    end if;
  else
    ProcessResolvent(s, ridx);
  end if;
  vprint PGG_GaloisGroup: "#pool =", #s`pool, #{n : n in s`pool | n`maybe_equal}, #{m : m in PossibleChildren(n), n in s`pool}, #{m : m in PossibleChildren(n), n in s`pool | m`maybe_equal};
  while true do
    PGG_GlobalTimer_Swap("is done");
    if IsDone(s) then
      break;
    end if;
    PGG_GlobalTimer_Swap("should descend");
    if ShouldDescend(s) then
      PGG_GlobalTimer_Swap("descend");
      vprint PGG_GaloisGroup: "descending...";
      Descend(s);
      vprint PGG_GaloisGroup: "#pool =", #s`pool, #{n : n in s`pool | n`maybe_equal}, #{m : m in PossibleChildren(n), n in s`pool}, #{m : m in PossibleChildren(n), n in s`pool | m`maybe_equal};
      if s`algorithm`reprocess then
        PGG_GlobalTimer_Swap("reprocess");
        for ridx in [1..#s`resolvents] do
          if IsDone(s) then
            break;
          end if;
          ProcessResolvent(s, ridx);
        end for;
        vprint PGG_GaloisGroup: "#pool =", #s`pool, #{n : n in s`pool | n`maybe_equal}, #{m : m in PossibleChildren(n), n in s`pool}, #{m : m in PossibleChildren(n), n in s`pool | m`maybe_equal};
      end if;
      if s`algorithm`reset then
        PGG_GlobalTimer_Swap("reset");
        s`useful_is_special := false;
        s`actually_useful_is_special := false;
        Reset(s`choice_state);
      end if;
    else
      break;
    end if;
  end while;
  PGG_GlobalTimer_Pop();
end intrinsic;

intrinsic Stabilizers(G :: GrpPerm, h :: Map, v :: PGGStatVal) -> []
  {The preimages under h of the stabilizers of v in G@h}
  return [G meet (S@@h) : S in Stabilizers(G@h, v)];
end intrinsic;

intrinsic Stabilizers(G :: GrpPerm, h :: [Map], v :: PGGStatVal) -> []
  {"}
  G2, os, h2 := PGG_Image(G, h);
  return [G meet (S@@h2) : S in Stabilizers(G2, os, v)];
end intrinsic;

intrinsic ProcessResolvent(s :: PGGAlgState_ResGroups_Maximal2, resolvent, subgroup)
  {"}
  PGG_GlobalTimer_Push("resolvent stat");
  rstat := ResolventStat(s`algorithm`statistic, resolvent);
  vprint PGG_GaloisGroup: "resolvent statistic =", rstat;
  h := PGG_CosetAction(Group(s`resolvent_overgroup), subgroup);
  vprint PGG_GaloisGroup, 1: "pool stats =", {* GroupStat(n, h) : n in s`pool *};
  vprint PGG_GaloisGroup, 1: "child stats =", {* GroupStat(m, h) : m in PossibleChildren(n), n in s`pool *};
  // eliminate any children not containing the galois group
  PGG_GlobalTimer_Swap("filter children");
  ok_children := false;
  for n in s`pool join &join[PossibleChildren(n) : n in s`pool] do
    nstat := GroupStat(n, h);
    if not (rstat subset nstat) then
      NotSubgroup(n);
      ok_children := true;
    end if;
  end for;
  // mark pool nodes as not equal
  PGG_GlobalTimer_Swap("filter pool");
  ok_pool := false;
  for n in s`pool do
    nstat := GroupStat(n, h);
    if not (rstat eq nstat) then
      NotEqual(n);
      ok_pool := true;
    end if;
  end for;
  // eliminate the not equal pool nodes
  ns := {n : n in s`pool | not n`maybe_equal};
  if #ns gt 0 then
    case s`algorithm`descend:
    when "Strict":
      // for each P in the pool, we find it the maximal subgroups with the observed statistic (called "stabilizers")
      emb := s`overgroup_embedding;
      PGG_GlobalTimer_Swap("stabilizers");
      Hs := &cat[[(S @@ emb) meet P : S in Stabilizers(P @ emb, h, rstat)] where P:=Rep(n`class) : n in ns];
      PGG_GlobalTimer_Swap("valid");
      Hs := [H : H in Hs | GalGrpIsConsistentWithSubgroup(Domain(emb), H)];
      PGG_GlobalTimer_Swap("nodes");
      ms := {m : H in Hs | m`maybe_subgroup where m:=Node(s`graph, H)};
    when "StrictOld":
      // this version takes stabilizers once in Codomain(h), and then needs to intersect conjugates with groups in the pool
      emb := s`overgroup_embedding;
      W1 := Group(Domain(emb));
      W2 := Group(Codomain(emb));
      I2 := Group(Image(emb));
      PGG_GlobalTimer_Swap("stabilizers");
      S3s := Stabilizers(Codomain(h), rstat);
      S2s := [S3 @@ h : S3 in S3s];
      PGG_GlobalTimer_Swap("double coset representatives");
      S1s := [(S2^x) @@ emb : x in DoubleCosetRepresentatives(W2, Normalizer(W2, S2), I2), S2 in S2s];
      PGG_GlobalTimer_Swap("conjugates");
      Ss := [S : S in Conjugates(W1,S1), S1 in S1s];
      PGG_GlobalTimer_Swap("intersect");
      Hs := [S meet Rep(n`class) : S in Ss, n in ns];
      PGG_GlobalTimer_Swap("valid");
      Hs := [H : H in Hs | GalGrpIsConsistentWithSubgroup(Domain(emb), H)];
      PGG_GlobalTimer_Swap("nodes");
      ms := { m : H in Hs | m`maybe_equal and m`maybe_subgroup where m := Node(s`graph, H)};
    when "SimpleMagic":
      ms := {m : m in PossibleChildren(n), n in ns | GroupStat(m,h) eq rstat};
    when "Magic":
      ms := {};
      for n in ns do
        todo := [n];
        nms := {};
        while #todo ne 0 do
          nn := todo[1];
          todo := todo[2..#todo];
          if #PossibleChildren(nn) ne 0 then
            nnms := {m : m in PossibleChildren(nn) | GroupStat(m,h) eq rstat};
            if #nnms ne 0 then
              nms join:= nnms;
            else
              for m in PossibleChildren(nn) do
                if not (rstat subset GroupStat(m,h)) then
                  NotSubgroup(m);
                end if;
              end for;
              todo cat:= SetToSequence(PossibleChildren(nn));
            end if;
          end if;
        end while;
        ms join:= nms;
      end for;
    else
      assert false;
    end case;
    if s`algorithm`reset and #(ms diff s`pool) ne 0 then
      PGG_GlobalTimer_Swap("reset");
      Reset(s`choice_state);
    end if;
    s`pool diff:= ns;
    s`pool join:= ms;
  end if;
  assert ok_children or ok_pool;
  PGG_GlobalTimer_Swap("possible children");
  vprint PGG_GaloisGroup: "#pool =", #s`pool, #&join[PossibleChildren(n) : n in s`pool];
  vprint PGG_GaloisGroup: "pool indices =", {* Index(Group(s`galois_conjugacy), Rep(n`class)) : n in s`pool *};
  PGG_GlobalTimer_Pop();
end intrinsic;

intrinsic ProcessResolvent(s :: PGGAlgState_ResGroups_Maximal, ridx :: RngIntElt)
  {Processes the ridxth resolvent of s.}
  _, _, rstat, h := Explode(s`resolvents[ridx]);
  for n in s`pool join &join[PossibleChildren(n) : n in s`pool] do
    if n`maybe_equal then
      nstat := GroupStat(n, h);
      if rstat ne nstat then
        NotEqual(n);
      end if;
      if n`maybe_subgroup and not (rstat subset nstat) then
        NotSubgroup(n);
      end if;
    elif n`maybe_subgroup then
      nstat := GroupStat(n, h);
      if not (rstat subset nstat) then
        NotSubgroup(n);
      end if;
    end if;
  end for;
end intrinsic;

intrinsic MustDescend(s :: PGGAlgState_ResGroups_Maximal) -> BoolElt
  {True if we must descend.}
  assert #s`pool gt 0;
  if pool_has_one_class(s) then
    ok := forall{n : n in s`pool | n`maybe_equal};
    assert ok eq exists{n : n in s`pool | n`maybe_equal};
    return not ok;
  elif forall{n : n in s`pool | not n`maybe_equal} and #&join[PossibleChildren(n) : n in s`pool] eq 1 then
    return true;
  else
    return false;
  end if;
end intrinsic;

intrinsic CanDescend(s :: PGGAlgState_ResGroups_Maximal) -> BoolElt
  {True if we can descend.}
  return exists{n : n in s`pool | (not n`maybe_equal) and (#PoolableChildren(n) ne 0)};
end intrinsic;

intrinsic ShouldDescend(s :: PGGAlgState_ResGroups_Maximal) -> BoolElt
  {True if we should descend.}
  if MustDescend(s) then
    return true;
  elif not CanDescend(s) then
    return false;
  end if;
  case s`algorithm`descend_when:
  when "Always":
    return exists{n : n in s`pool | not n`maybe_equal};
  when "AllUnequal":
    return forall{n : n in s`pool | not n`maybe_equal};
  when "Sufficient":
    return not(#[n : n in s`pool | n`maybe_equal] ge 2 or exists{n : n in s`pool | n`maybe_equal and #PossibleChildren(n) ge 1});
  when "Necessary":
    error "not implemented";
  when "NoSubgroup":
    s`useful_is_special := true;
    return not HasSubgroup(s : NoForget);
  when "AllUnequalAndNoSubgroup":
    if forall{n : n in s`pool | not n`maybe_equal} then
      s`useful_is_special := true;
      return not HasSubgroup(s : NoForget);
    else
      return false;
    end if;
  when "Ask":
    while true do
      read ans, "Descend? (y/n) ";
      case ans:
      when "y":
        return true;
      when "n":
        return false;
      end case;
    end while;
  else
    assert false;
  end case;
end intrinsic;

intrinsic PooledChildren(n :: PGGAlgState_ResGroups_Maximal_GraphNode) -> {}
  {The pooled children of n.}
  return n`pooled_children;
end intrinsic;

intrinsic PoolableChildren(n :: PGGAlgState_ResGroups_Maximal_GraphNode) -> {}
  {The poolable children of n.}
  return PossibleChildren(n) diff PooledChildren(n);
end intrinsic;

// intrinsic PoolAllChildren(ns :: [PGGAlgState_ResGroups_Maximal_GraphNode])
//   {Pools all children from the given nodes.}
//   if #ns gt 0 then
//     s := ns[1]`parent`state;
//     for n in ns do
//       n2 := n;
//       n2`pooled_children join:= PossibleChildren(n) meet s`pool;
//       newpool := PoolableChildren(n);
//       assert forall{m : m in newpool | m`maybe_subgroup};
//       s`pool join:= newpool;
//       n2`pooled_children join:= newpool;
//     end for;
//     s`pool diff:= SequenceToSet(ns);
//     NormalizePool(s);
//   end if;
// end intrinsic;

// intrinsic PoolAllChildren(n :: PGGAlgState_ResGroups_Maximal_GraphNode)
//   {Pools all children of n.}
//   PoolAllChildren([n]);
// end intrinsic;

// intrinsic PoolChild(m :: PGGAlgState_ResGroups_Maximal_GraphNode)
//   {Pools the child m.}
//   assert m`maybe_subgroup;
//   s := m`parent`state;
//   Include(~s`pool, m);
//   for n in m`possible_parents do
//     nn := n;
//     nn`pooled_children join:= PossibleChildren(n) meet s`pool;
//     Include(~nn`pooled_children, m);
//     if #PoolableChildren(n) eq 0 then
//       Exclude(~s`pool, n);
//     end if;
//   end for;
//   NormalizePool(s);
// end intrinsic;

intrinsic RemoveFromPool(n :: PGGAlgState_ResGroups_Maximal_GraphNode)
  {Removes n from the pool.}
  assert (not n`maybe_subgroup) or ((not n`maybe_equal) and #PoolableChildren(n) eq 0);
  Exclude(~n`parent`state`pool, n);
end intrinsic;

intrinsic AddToPool(ns :: {PGGAlgState_ResGroups_Maximal_GraphNode})
  {Adds the given nodes to the pool.}
  if #ns gt 0 then
    s := Rep(ns)`parent`state;
    // each node must be a direct child of an element of the pool
    assert forall{n : n in ns | exists{m : m in n`possible_parents | m in s`pool}};
    // add to the pool
    s`pool join:= ns;
    // record that each one was pooled
    for n in ns do
      for m in n`possible_parents do
        mm := m;
        Include(~mm`pooled_children, n);
      end for;
    end for;
    // remove from the pool anything which is unequal and has no poolable children
    for m in &join[n`possible_parents : n in ns] do
      if (m in s`pool) and (not m`maybe_equal) and (#PoolableChildren(m) eq 0) then
        RemoveFromPool(m);
      end if;
    end for;
  end if;
end intrinsic;

intrinsic Descend(s :: PGGAlgState_ResGroups_Maximal)
  {Descends.}
  case s`algorithm`descend:
  when "All":
    AddToPool(&join[PoolableChildren(n) : n in s`pool | not n`maybe_equal]);
    // PoolAllChildren([n : n in s`pool | not n`maybe_equal]);
  when "OneNode":
    cands := [n : n in s`pool | not n`maybe_equal];
    assert #cands gt 0;
    _, i := Max([#Rep(n`class) : n in cands]);
    AddToPool(PoolableChildren(cands[i]));
  when "OneChild":
    cands := &join[PoolableChildren(n) : n in s`pool | not n`maybe_equal];
    assert #cands gt 0;
    _, i := Max([#Rep(n`class) : n in cands]);
    AddToPool(cands[i]);
  else
    assert false;
  end case;
end intrinsic;

function layer_has_consistent_subgroup(s, i, G)
  // check the inputs
  assert i ge 1;
  if i eq 1 then
    assert GalGrpIsConsistentWithSubgroup(s`conjugacy, G);
  end if;
  // base case
  if i eq #s`current_layer then
    assert GalGrpIsConsistentWithSubgroup(s`conjugacy, G);
    return true;
  end if;
  // general case
  assert i lt #s`current_layer;
  for H in Conjugates(s`conjugacy`group, Rep(s`current_layer[i+1]`class)) do
    G2 := H meet G;
    if GalGrpIsConsistentWithSubgroup(s`conjugacy, G2) then
      if layer_has_consistent_subgroup(s, i+1, G2) then
        return true;
      end if;
    end if;
  end for;
  return false;
end function;

intrinsic UsefulHook(s :: PGGAlgState_ResGroups, i :: PGGAlgState_TrancheItem)
  {This is called when i is useful. Overload this.}
  return;
end intrinsic;

intrinsic UsefulHook(s :: PGGAlgState_ResGroups_Maximal, i :: PGGAlgState_TrancheItem)
  {"}
  if s`useful_is_special then
    MarkSpecialTranche(i);
  end if;
end intrinsic;

intrinsic IsUseful(s :: PGGAlgState_ResGroups, i :: PGGAlgState_TrancheItem) -> BoolElt
  {True if i represents a useful group.}
  ok := IsUseful(s, i`subgroup);
  if ok then
    UsefulHook(s, i);
  else
    Forget(i);
  end if;
  return ok;
end intrinsic;

intrinsic IsUseful(s :: PGGAlgState_ResGroups_Maximal, i :: PGGAlgState_TrancheItem) -> BoolElt
  {"}
  ok := IsUseful(s, i`subgroup);
  if not ok then
    Forget(i);
  elif s`useful_is_special then
    MarkSpecialTranche(i);
  end if;
end intrinsic;

intrinsic IsUseful(s :: PGGAlgState_ResGroups_All, G) -> BoolElt
  {True if G is useful.}
  Gs := s`possible_groups;
  h := CosetAction(Group(s`resolvent_overgroup), G);
  statalg := s`algorithm`statistic;
  stat1 := GroupStat(statalg, Gs[1] @ s`overgroup_embedding, h);
  is_useful := exists{G : G in Gs[2..#Gs] | GroupStat(statalg, G @ s`overgroup_embedding, h) ne stat1};
  return is_useful;
end intrinsic;

intrinsic IsUseful(s :: PGGAlgState_ResGroups_Maximal, G) -> BoolElt
  {"}
  case s`algorithm`useful:
  when "Sufficient":
    h := CosetAction(Group(s`resolvent_overgroup), G);
    nstats := [<n, GroupStat(n, h)> : n in s`pool | n`maybe_equal];
    if exists{nstat : nstat in nstats[2..#nstats] | nstat[2] ne nstats[1][2]} then
      is_useful := true;
    else
      mstats := AssociativeArray();
      for m in &join[PossibleChildren(n) : n in s`pool] do
        mstats[m] := GroupStat(m, h);
      end for;
      is_useful := exists{nstat : m in PossibleChildren(nstat[1]), nstat in nstats | (mstat ne nstat[2]) and (mstat subset nstat[2]) where mstat:=mstats[m]};
      // mstats := [GroupStat(m, h) : m in &join[PossibleChildren(n) : n in s`pool]];
      // is_useful := exists{nstat : nstat in nstats, mstat in mstats | (nstat ne mstat) and (mstat subset nstat)};
    end if;
  when "Necessary":
    h := CosetAction(Group(s`resolvent_overgroup), G);
    nstats := [GroupStat(n, h) : n in s`pool | n`maybe_equal];
    if exists{nstat : nstat in nstats[2..#nstats] | nstat ne nstats[1]} then
      is_useful := true;
    else
      mstats := [GroupStat(m, h) : m in &join[PossibleChildren(n) : n in s`pool]];
      is_useful := exists{nstat : nstat in nstats, mstat in mstats | not (nstat subset mstat)};
    end if;
  when "Generous":
    h := CosetAction(Group(s`resolvent_overgroup), G);
    nstats := [GroupStat(n, h) : n in {n : n in s`pool | n`maybe_equal} join &join[PossibleChildren(n) : n in s`pool]];
    assert #nstats ge 2;
    is_useful := exists{nstat : nstat in nstats[2..#nstats] | nstat ne nstats[1]};
  when "All":
    is_useful := Index(Group(s`resolvent_overgroup), G) gt 1;
  else
    assert false;
  end case;
  return is_useful;
end intrinsic;

intrinsic IsUseful(s :: PGGAlgState_ResGroups_Maximal2, G) -> BoolElt
  {"}
  assert #s`pool ne 0;
  assert forall{n : n in s`pool | n`maybe_equal and n`maybe_subgroup};
  assert forall{m : m in &join[PossibleChildren(n) : n in s`pool] | m`maybe_subgroup};
  assert (not pool_has_one_class(s)) or exists{n : n in s`pool | #PossibleChildren(n) gt 0};
  h := CosetAction(Group(s`resolvent_overgroup), G);
  nstats := [<n, GroupStat(n, h)> : n in s`pool];
  is_useful := exists{nstat : nstat in nstats[2..#nstats] | nstat[2] ne nstats[1][2]}
    or exists{nstat : m in PossibleChildren(nstat[1]), nstat in nstats | GroupStat(m, h) ne nstat[2]};
  return is_useful;
end intrinsic;

intrinsic Information(i :: PGGAlgState_TrancheItem, s :: PGGAlgState_ResGroups_All) -> FldReElt
  {The information provided by this item.}
  Gs := s`possible_groups;
  h := CosetAction(s`overgroup, i`subgroup);
  statalg := s`algorithm`statistic;
  stats := {* GroupStat(statalg, h(G)) : G in Gs *};
  return multiplicities_to_information(Multiplicities(stats));
end intrinsic;

intrinsic Diversity(i :: PGGAlgState_TrancheItem, s :: PGGAlgState_ResGroups_All) -> RngIntElt
  {The information provided by this item.}
  Gs := s`possible_groups;
  h := CosetAction(s`overgroup, i`subgroup);
  statalg := s`algorithm`statistic;
  stats := { GroupStat(statalg, h(G)) : G in Gs };
  return #stats;
end intrinsic;

intrinsic HasSubgroup(s :: PGGAlgState_ResGroups_All) -> BoolElt, GrpPerm
  {True if s has a subgroup to form a resolvent from.}
  return HasSubgroup(s`choice_state);
end intrinsic;

intrinsic HasSubgroup(s :: PGGAlgState_ResGroups_RootsMaximal) -> BoolElt, GrpPerm
  {"}
  ps := PossibleSubgroups(s);
  if #ps eq 0 then
    return false, _;
  else
    return true, Rep(ps[1]);
  end if;
end intrinsic;

intrinsic HasSubgroup(s :: PGGAlgState_ResGroups_Sequence) -> BoolElt, GrpPerm
  {"}
  s2 := CurrentState(s);
  while true do
    ok, G := HasSubgroup(s2);
    if ok then
      return true, G;
    elif AtEnd(s) then
      return false, _;
    else
      PGG_GlobalTimer_Push("start groups");
      s2 := NextState(s);
      PGG_GlobalTimer_Pop();
    end if;
  end while;
end intrinsic;

intrinsic HasSubgroup(s :: PGGAlgState_ResGroups_Null) -> BoolElt, GrpPerm
  {"}
  return false, _;
end intrinsic;

intrinsic HasSubgroup(s :: PGGAlgState_ResGroups_Maximal : NoForget:=false) -> BoolElt, GrpPerm
  {"}
  return HasSubgroup(s`choice_state : NoForget:=NoForget);
end intrinsic;

intrinsic HasSubgroup(s :: PGGAlgState_ResGroups_Maximal2) -> BoolElt, GrpPerm
  {"}
  return HasSubgroup(s`choice_state);
end intrinsic;

intrinsic Subgroup(s :: PGGAlgState_ResGroups) -> GrpPerm
  {A subgroup to form a resolvent from.}
  ok, G := HasSubgroup(s);
  error if not ok, "no next subgroup";
  return G;
end intrinsic;