Stat.mag

import "Utils.mag": has_root, roots, factorization, extension, not_implemented, is_subpartition_of, all_binnings, mset_apply, xdiv, Z;

declare type PGGStat;
declare type PGGStat_EqSimple: PGGStat;
declare type PGGStat_HasRoot: PGGStat_EqSimple;
declare type PGGStat_NumRoots: PGGStat_EqSimple;
declare type PGGStat_FactorDegrees: PGGStat_EqSimple;
declare type PGGStat_Degree: PGGStat_EqSimple;
declare type PGGStat_Order: PGGStat_EqSimple;
declare type PGGStat_AutGroup: PGGStat_EqSimple;
declare attributes PGGStat_AutGroup: conjugacy_classes;
declare type PGGStat_NumAuts: PGGStat_EqSimple;
declare type PGGStat_Tup: PGGStat_EqSimple;
declare attributes PGGStat_Tup: stats;
declare type PGGStat_Factors: PGGStat_EqSimple;
declare attributes PGGStat_Factors: stat;
declare type PGGStat_Factors2: PGGStat;
declare attributes PGGStat_Factors2: strict, stat1, stat2;
declare type PGGStat_Stab: PGGStat_EqSimple;
declare attributes PGGStat_Stab: stat;
declare type PGGStat_GaloisGroup: PGGStat;
declare attributes PGGStat_GaloisGroup: alg, degmax;
declare type PGGStat_SubfieldDegrees: PGGStat_EqSimple; // pair <(L:K), {*(M:K) : L/M/K*}>
declare type PGGStat_FactorDegreesSeq: PGGStat_EqSimple;

declare type PGGStatVal;
declare attributes PGGStatVal: stat, value, norbits, hash, hashable;

intrinsic MakeValue(S :: PGGStat, value, norbits :: RngIntElt) -> PGGStatVal
  {A value for the statistic S.}
  x := New(PGGStatVal);
  x`stat := S;
  x`value := value;
  x`norbits := norbits;
  return x;
end intrinsic;

intrinsic MakeValue(S :: PGGStat, value) -> PGGStatVal
  {"}
  return MakeValue(S, value, 1);
end intrinsic;

intrinsic Hash(v :: PGGStatVal) -> .
  {Hash of v.}
  if not assigned v`hash then
    v`hash := Hash(Hashable(v));
  end if;
  return v`hash;
end intrinsic;

intrinsic Hashable(v :: PGGStatVal) -> .
  {Hashable value for v.}
  if not assigned v`hashable then
    v`hashable := _Hashable(v`stat, v);
  end if;
  return v`hashable;
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat, S2 :: PGGStat) -> BoolElt
  {Equality.}
  return IsIdentical(S1, S2);
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat_HasRoot, S2 :: PGGStat_HasRoot) -> BoolElt
  {"}
  return true;
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat_NumRoots, S2 :: PGGStat_NumRoots) -> BoolElt
  {"}
  return true;
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat_FactorDegrees, S2 :: PGGStat_FactorDegrees) -> BoolElt
  {"}
  return true;
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat_Degree, S2 :: PGGStat_Degree) -> BoolElt
  {"}
  return true;
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat_AutGroup, S2 :: PGGStat_AutGroup) -> BoolElt
  {"}
  return true;
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat_NumAuts, S2 :: PGGStat_NumAuts) -> BoolElt
  {"}
  return true;
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat_Order, S2 :: PGGStat_Order) -> BoolElt
  {"}
  return true;
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat_Factors, S2 :: PGGStat_Factors) -> BoolElt
  {"}
  return S1`stat eq S2`stat;
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat_Stab, S2 :: PGGStat_Stab) -> BoolElt
  {"}
  return S1`stat eq S2`stat;
end intrinsic;

intrinsic 'eq'(S1 :: PGGStat_Tup, S2 :: PGGStat_Tup) -> BoolElt
  {"}
  return #S1`stats eq #S2`stats and forall{i : i in [1..#S1`stats] | S1`stats[i] eq S2`stats[i]};
end intrinsic;

intrinsic 'eq'(v1 :: PGGStatVal, v2 :: PGGStatVal) -> BoolElt
  {Equality.}
  require v1`stat eq v2`stat: "must be values of the same statistic";
  require v1`norbits eq v2`norbits: "must have the same number of orbits";
  return 'eq'(v1`stat, v1, v2);
end intrinsic;

intrinsic 'subset'(v1 :: PGGStatVal, v2 :: PGGStatVal) -> BoolElt
  {True if v1 is possibly a statistic for a subgroup of the group v2 is a subset for.}
  require v1`stat eq v2`stat: "must be values of the same statistic";
  require v1`norbits eq v2`norbits: "must have the same number of orbits";
  return 'subset'(v1`stat, v1, v2);
end intrinsic;

intrinsic Print(v :: PGGStatVal)
  {Prints.}
  _Print(v`stat, v);
end intrinsic;

intrinsic PGGStat_HasRoot_Make() -> PGGStat_HasRoot
  {Makes a new PGGStat.}
  return New(PGGStat_HasRoot);
end intrinsic;

intrinsic PGGStat_NumRoots_Make() -> PGGStat_NumRoots
  {Makes a new PGGStat.}
  return New(PGGStat_NumRoots);
end intrinsic;

intrinsic PGGStat_FactorDegrees_Make() -> PGGStat_FactorDegrees
  {Makes a new PGGStat.}
  return New(PGGStat_FactorDegrees);
end intrinsic;

intrinsic PGGStat_FactorDegreesSeq_Make() -> PGGStat_FactorDegreesSeq
  {"}
  return New(PGGStat_FactorDegreesSeq);
end intrinsic;

intrinsic PGGStat_AutGroup_Make() -> PGGStat_AutGroup
  {"}
  return New(PGGStat_AutGroup);
end intrinsic;

intrinsic PGGStat_NumAuts_Make() -> PGGStat_NumAuts
  {"}
  return New(PGGStat_NumAuts);
end intrinsic;

intrinsic PGGStat_Degree_Make() -> PGGStat_Degree
  {Makes a new PGGStat.}
  return New(PGGStat_Degree);
end intrinsic;

intrinsic PGGStat_Order_Make() -> PGGStat_Order
  {"}
  return New(PGGStat_Order);
end intrinsic;

intrinsic PGGStat_Factors_Make(:Stat:=false) -> PGGStat_Factors
  {"}
  S := New(PGGStat_Factors);
  S`stat := Stat cmpne false select Stat else PGGStat_Degree_Make();
  return S;
end intrinsic;

intrinsic PGGStat_Factors2_Make(:Stat1:=false,Stat2:=false,Strict:=false) -> PGGStat_Factors2
  {"}
  S := New(PGGStat_Factors2);
  S`stat2 := Stat2 cmpne false select Stat2 else PGGStat_Degree_Make();
  S`stat1 := Stat1 cmpne false select Stat1 else PGGStat_Degree_Make();
  S`strict := Strict;
  return S;
end intrinsic;

intrinsic PGGStat_Tup_Make(args) -> PGGStat_Tup
  {Makes a new PGGStat.}
  S := New(PGGStat_Tup);
  S`stats := args;
  return S;
end intrinsic;

intrinsic PGGStat_Stab_Make(:Stat:=false) -> PGGStat_Stab
  {"}
  S := New(PGGStat_Stab);
  S`stat := Stat cmpne false select Stat else PGGStat_Degree_Make();
  return S;
end intrinsic;

intrinsic PGGStat_GaloisGroup_Make(:Alg:=false,DegreeLimit:=false) -> PGGStat_GaloisGroup
  {"}
  S := New(PGGStat_GaloisGroup);
  S`degmax := DegreeLimit cmpne false select DegreeLimit else Infinity();
  if Alg cmpne false then
    S`alg := Alg;
  end if;
  return S;
end intrinsic;

intrinsic PGGStat_SubfieldDegrees_Make() -> PGGStat_SubfieldDegrees
  {"}
  S := New(PGGStat_SubfieldDegrees);
  return S;
end intrinsic;

intrinsic Print(S :: PGGStat_HasRoot)
  {Print.}
  printf "has root";
end intrinsic;

intrinsic Print(S :: PGGStat_NumRoots)
  {"}
  printf "number of roots";
end intrinsic;

intrinsic Print(S :: PGGStat_FactorDegrees)
  {"}
  printf "degrees of factors";
end intrinsic;

intrinsic Print(S :: PGGStat_FactorDegreesSeq)
  {"}
  printf "degrees of factors (sequence)";
end intrinsic;

intrinsic Print(S :: PGGStat_AutGroup)
  {"}
  printf "automorphism group";
end intrinsic;

intrinsic Print(S :: PGGStat_NumAuts)
  {"}
  printf "number of automorphisms";
end intrinsic;

intrinsic Print(S :: PGGStat_Degree)
  {"}
  printf "degree";
end intrinsic;

intrinsic Print(S :: PGGStat_Factors)
  {"}
  printf "factors -> ";
  Print(S`stat);
end intrinsic;

intrinsic Print(S :: PGGStat_Factors2)
  {"}
  print "factors of factors";
  IndentPush();
  printf "stat2 = "; Print(S`stat2);
  print "";
  printf "stat1 = "; Print(S`stat1);
  print "";
  printf "strict = %o", S`strict;
  IndentPop();
end intrinsic;

intrinsic Print(S :: PGGStat_Tup)
  {"}
  if #S`stats eq 0 then
    printf "tuple (empty)";
  else
    print "tuple";
    IndentPush();
    for i in [1..#S`stats] do
      Print(S`stats[i]);
      if i lt #S`stats then
        print "";
      end if;
    end for;
    IndentPop();
  end if;
end intrinsic;

intrinsic Print(S :: PGGStat_Stab)
  {"}
  printf "stabilizer -> "; Print(S`stat);
end intrinsic;

intrinsic Print(S :: PGGStat_GaloisGroup)
  {"}
  printf "galois group";
  if assigned S`alg then
    print "";
    IndentPush();
    printf "alg = ";
    Print(S`alg);
    IndentPop();
  end if;
  if S`degmax lt Infinity() then
    print "";
    IndentPush();
    printf "max degree = %o", S`degmax;
    IndentPop();
  end if;
end intrinsic;

intrinsic Print(S :: PGGStat_Order)
  {"}
  printf "galois group order";
end intrinsic;

intrinsic Print(S :: PGGStat_SubfieldDegrees)
  {"}
  printf "degrees of subfields";
end intrinsic;

intrinsic ConjugacyClasses(S :: PGGStat_AutGroup, d :: RngIntElt) -> PGGSetSubgrpcls
  {The conjugacy classes of Sym(d).}
  if not assigned S`conjugacy_classes then
    S`conjugacy_classes := AssociativeArray();
  end if;
  if not IsDefined(S`conjugacy_classes, d) then
    S`conjugacy_classes[d] := PGG_SubgroupClasses(SymmetricGroup(d));
  end if;
  return S`conjugacy_classes[d];
end intrinsic;

intrinsic GroupStat(S :: PGGStat, G :: GrpPerm) -> PGGStatVal
  {The value of the statistic on the group.}
  return GroupStat(S, G, [GSet(G)]);
end intrinsic;

intrinsic GroupStat(S :: PGGStat, G :: GrpPerm, os :: [GSetIndx[RngIntElt]]) -> PGGStatVal
  {The value of the statistic on the group, with the distinguished orbits os.}
  not_implemented("GroupStat:", Type(S));
end intrinsic;

intrinsic GroupStat(S :: PGGStat_HasRoot, G :: GrpPerm) -> PGGStatVal
  {"}
  return MakeValue(S, exists{o : o in Orbits(G) | #o eq 1});
end intrinsic;

intrinsic GroupStat(S :: PGGStat_NumRoots, G :: GrpPerm) -> PGGStatVal
  {"}
  return MakeValue(S, #[o : o in Orbits(G) | #o eq 1]);
end intrinsic;

intrinsic GroupStat(S :: PGGStat_FactorDegrees, G :: GrpPerm) -> PGGStatVal
  {"}
  return MakeValue(S, {* #o : o in Orbits(G) *});
end intrinsic;

intrinsic GroupStat(S :: PGGStat_FactorDegreesSeq, G :: GrpPerm, os :: [GSetIndx[RngIntElt]]) -> PGGStatVal
  {"}
  assert #os ge 1;
  if #os eq 1 then
    return MakeValue(S, {* #o : o in Orbits(OrbitImage(G, os[1])) *});
  elif #os eq 2 then
    return MakeValue(S, {* <#o, GroupStat(S, Go, os[2..#os])> where Go:=Stabilizer(G,Rep(o)) : o in Orbits(G,os[1]) *}, #os);
  else
    not_implemented("GroupStat: PGGStat_FactorDegreesSeq: more than 2 orbits");
  end if;
end intrinsic;

intrinsic GroupStat(S :: PGGStat_Degree, G :: GrpPerm) -> PGGStatVal
  {"}
  return MakeValue(S, Degree(G));
end intrinsic;

intrinsic GroupStat(S :: PGGStat_Order, G :: GrpPerm) -> PGGStatVal
  {"}
  return MakeValue(S, #G);
end intrinsic;

intrinsic GroupStat(S :: PGGStat_AutGroup, G :: GrpPerm) -> PGGStatVal
  {"}
  assert IsTransitive(G);
  stab := Stabilizer(G, 1);
  norm := Normalizer(G, stab);
  aut := CosetImage(norm, stab);
  return MakeValue(S, ConjugacyClasses(S, Degree(aut)) ! aut);
end intrinsic;

intrinsic GroupStat(S :: PGGStat_NumAuts, G :: GrpPerm) -> PGGStatVal
  {"}
  assert IsTransitive(G);
  stab := Stabilizer(G, 1);
  norm := Normalizer(G, stab);
  return MakeValue(S, Index(norm, stab));
end intrinsic;

intrinsic GroupStat(S :: PGGStat_Tup, G :: GrpPerm) -> PGGStatVal
  {"}
  return MakeValue(S, <GroupStat(S2, G) : S2 in S`stats>);
end intrinsic;

intrinsic GroupStat(S :: PGGStat_Factors, G :: GrpPerm) -> PGGStatVal
  {"}
  return MakeValue(S, {* GroupStat(S`stat, OrbitImage(G, o)) : o in Orbits(G) *});
end intrinsic;

intrinsic GroupStat(S :: PGGStat_Factors2, G :: GrpPerm) -> PGGStatVal
  {"}
  return MakeValue(S,
    < [GroupStat(S`stat1, OrbitImage(G,o)) : o in os]
    , [ [ {* GroupStat(S`stat2, OrbitImage(Go,o2)) : o2 in Orbits(Go) *} 
          where Go := OrbitImage(So,o) 
          : o in os
        ]
        where So := Stabilizer(G,Rep(o))
        : o in os 
        ]
    >
    where os := [o : o in Orbits(G)]
  );
end intrinsic;

intrinsic GroupStat(S :: PGGStat_Stab, G :: GrpPerm) -> PGGStatVal
  {"}
  error if not IsTransitive(G), "only applies to transitive groups";
  return MakeValue(S, GroupStat(S`stat, Stabilizer(G,1)));
end intrinsic;

intrinsic GroupStat(S :: PGGStat_GaloisGroup, G :: GrpPerm) -> PGGStatVal
  {"}
  return MakeValue(S, Degree(G) le S`degmax select G else sub<G | []>);
end intrinsic;

intrinsic GroupStat(S :: PGGStat_SubfieldDegrees, G :: GrpPerm) -> PGGStatVal
  {"}
  error if not IsTransitive(G), "only applies to transitive groups";
  d := Degree(G);
  return MakeValue(S, <d, {* xdiv(d, #X) : X in AllPartitions(G) join {{1},{1..d}}*}>);
end intrinsic;

intrinsic ResolventStat(S :: PGGStat, R :: PGGPol) -> PGGStatVal
  {The value of the statistic on the resolvent.}
  return ResolventStat(S, [R]);
end intrinsic;

intrinsic ResolventStat(S :: PGGStat, Rs :: [PGGPol]) -> PGGStatVal
  {The value of the statistic on the sequence of resolvents.}
  not_implemented("ResolventStat:", Type(S));
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_HasRoot, R :: PGGPol) -> PGGStatVal
  {"}
  return MakeValue(S, HasRoot(R));
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_NumRoots, R :: PGGPol) -> PGGStatVal
  {"}
  return MakeValue(S, NumRoots(R));
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_FactorDegrees, R :: PGGPol) -> PGGStatVal
  {"}
  return MakeValue(S, {* d : d in FactorDegrees(R) *});
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_FactorDegreesSeq, R :: [PGGPol]) -> PGGStatVal
  {"}
  assert #R ge 1;
  if #R eq 1 then
    return MakeValue(S, {* d : d in FactorDegrees(R[1]) *});
  else
    return MakeValue(S, {* <Degree(fac), ResolventStat(S, [PolynomialRing(Extension(fac))| f : f in R[2..#R]])> : fac in Factorization(R[1]) *}, #R);
  end if;
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_Degree, R :: PGGPol) -> PGGStatVal
  {"}
  return MakeValue(S, Degree(R));
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_AutGroup, R :: PGGPol) -> PGGStatVal
  {"}
  K := BaseRing(R);
  L := Extension(R);
  aut := AutomorphismGroup(L, K);
  return MakeValue(S, ConjugacyClasses(S, Degree(aut)) ! aut);
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_NumAuts, R :: PGGPol) -> PGGStatVal
  {"}
  K := BaseRing(R);
  L := Extension(R);
  f := DefiningPolynomial(L,K);
  naut := #Roots(ChangeRing(f, L));
  return MakeValue(S, naut);
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_Tup, R :: PGGPol) -> PGGStatVal
  {"}
  return MakeValue(S, <ResolventStat(S2, R) : S2 in S`stats>);
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_Factors, R :: PGGPol) -> PGGStatVal
  {"}
  return MakeValue(S, {* ResolventStat(S`stat, fac) : fac in Factorization(R) *});
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_Factors2, R :: PGGPol) -> PGGStatVal
  {"}
  return MakeValue(S,
    < [ ResolventStat(S`stat1, fac) : fac in facs ]
    , [ [ {* ResolventStat(S`stat2, fac2) : fac2 in Factorization(ChangeRing(fac, L)) *}
          : fac in facs
        ]
        where L := certs[i]`Extension
        : i in [1..#facs]
      ]
    >
    where facs, certs := Factorization(R : Extensions)
  );
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_Stab, R :: PGGPol) -> PGGStatVal
  {"}
  L := Extension(R);
  // facs, certs := Factorization(R : Extensions);
  // error if #facs ne 1, "only applies to irreducible polynomials";
  return MakeValue(S, ResolventStat(S`stat, ChangeRing(R, L)));
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_GaloisGroup, R :: PGGPol) -> PGGStatVal
  {"}
  if Degree(R) gt S`degmax then
    return MakeValue(S, sub<SymmetricGroup(Degree(R)) | []>);
  elif assigned S`alg then
    PGG_GlobalTimer_Push("galois group");
    IndentPush();
    G := PGG_GaloisGroup(R : Alg:=S`alg);
    IndentPop();
    PGG_GlobalTimer_Pop();
    return MakeValue(S, G);
  else
    error "ResolventStat: Galois group algorithm not specified";
  end if;
end intrinsic;

intrinsic ResolventStat(S :: PGGStat_SubfieldDegrees, R :: PGGPol) -> PGGStatVal
  {"}
  K := BaseRing(R);
  L := Extension(R);
  return MakeValue(S, <Degree(L,K), {*1, Degree(L,K)*}>); // TODO!!!
end intrinsic;

intrinsic 'eq'(S :: PGGStat, v1, v2) -> BoolElt
  {Equality.}
  not_implemented("equality:", Type(S));
end intrinsic;

intrinsic 'eq'(S :: PGGStat_EqSimple, v1, v2) -> BoolElt
  {Equality.}
  return v1`value eq v2`value;
end intrinsic;

intrinsic 'eq'(S :: PGGStat_AutGroup, v1, v2) -> BoolElt
  {"}
  return Degree(Rep(v1`value)) eq Degree(Rep(v2`value)) and v1`value eq v2`value;
end intrinsic;

intrinsic 'eq'(S :: PGGStat_GaloisGroup, v1, v2) -> BoolElt
  {"}
  return Degree(v1`value) eq Degree(v2`value) and IsConjugate(SymmetricGroup(Degree(v1`value)), v1`value, v2`value);
end intrinsic;

procedure assoc_append(~A, k, x)
  if IsDefined(A, k) then
    Append(~A[k], x);
  else
    A[k] := [x];
  end if;
end procedure;

function is_conjugate2(G, X, xs, ys)
  H := G;
  for i in [1..#xs] do
    ok, g := IsConjugate(H, X, xs[i], ys[i]);
    if ok then
      for j in [i..#xs] do
        xs[j] := xs[j]^g;
      end for;
      assert xs[i] eq ys[i];
      H := Stabilizer(H, X, xs[i]);
    else
      return false;
    end if;
  end for;
  return true;
end function;

intrinsic 'eq'(S :: PGGStat_Factors2, v1, v2) -> BoolElt
  {"}
  if Hash(v1) ne Hash(v2) then
    return false;
  elif Hashable(v1) ne Hashable(v2) then
    return false;
  elif S`strict then
    xs1, yss1 := Explode(v1`value);
    xs2, yss2 := Explode(v2`value);
    n := #xs1;
    assert n eq #yss1;
    assert n eq #xs2;
    assert n eq #yss2;
    A1 := AssociativeArray();
    B1 := AssociativeArray();
    A2 := AssociativeArray();
    B2 := AssociativeArray();
    for i in [1..n] do
      assoc_append(~A1, xs1[i], [i,n+1]);
      assoc_append(~A2, xs2[i], [i,n+1]);
      for j in [1..n] do
        assoc_append(~B1, yss1[i][j], [i,j]);
        assoc_append(~B2, yss2[i][j], [i,j]);
      end for;
    end for;
    S := DirectProduct(SymmetricGroup(n), SymmetricGroup(1));
    as := SetToSequence(Keys(A1));
    bs := SetToSequence(Keys(B1));
    X := {[i,j] : i,j in [1..n+1]};
    GX := GSet(S, X, map<car<X, S> -> X | p :-> [p[1][1]^p[2], p[1][2]^p[2]] >);
    ok := is_conjugate2(S, GX, [SequenceToSet(B1[b]) : b in bs] cat [SequenceToSet(A1[a]) : a in as], [SequenceToSet(B2[b]) : b in bs] cat [SequenceToSet(A2[a]) : a in as]);
    if not ok then
      vprint PGG_GaloisGroup: "***** Factors2", v1, v2;
    end if;
    return ok;
  else
    return true;
  end if;
end intrinsic;

intrinsic 'subset'(S :: PGGStat, v1, v2) -> BoolElt
  {True if v1 is a refinement of v2.}
  not_implemented("subset:", Type(S));
end intrinsic;

intrinsic 'subset'(S :: PGGStat_HasRoot, v1, v2) -> BoolElt
  {"}
  return (not v2`value) or (v1`value);
end intrinsic;

intrinsic 'subset'(S :: PGGStat_NumRoots, v1, v2) -> BoolElt
  {"}
  return v1`value ge v2`value;
end intrinsic;

intrinsic 'subset'(S :: PGGStat_FactorDegrees, v1, v2) -> BoolElt
  {"}
  return is_subpartition_of(v1`value, v2`value);
end intrinsic;

intrinsic 'subset'(S :: PGGStat_FactorDegreesSeq, v1, v2 : Partial:=false) -> BoolElt
  {"}
  // put an ordering on the unique values
  uvs1 := SetToSequence(MultisetToSet(v1`value));
  uvs2 := SetToSequence(MultisetToSet(v2`value));
  // try to bin them
  bs := all_binnings(
    [Multiplicity(v1`value, x) : x in uvs1],
    [Multiplicity(v2`value, x) : x in uvs2]:
    limit := 1,
    is_semivalid := v1`norbits eq 1
    select func<i,b |
      // check the degrees
      &+mset_apply(b,func<j | uvs1[j]>) le uvs2[i]
    >
    else func<i,b |
      // check the first degrees
      (&+mset_apply(b,func<j | uvs1[j][1]>) le uvs2[i][1]) and
      // check the inner statistics, after joining and multiplying by first degrees
      'subset'(S,
        '&join'(S, mset_apply(b, func<j | '*'(S, uvs1[j][2], uvs1[j][1])>)),
        '*'(S, uvs2[i][2], uvs2[i][1])
      : Partial)
    >,
    is_valid := Partial
    select func<i,b | true>
    else v1`norbits eq 1
    select func<i,b |
      // check the degrees
      &+mset_apply(b,func<j | uvs1[j]>) eq uvs2[i]
    >
    else func<i,b |
      // check the first degrees
      (&+mset_apply(b,func<j | uvs1[j][1]>) eq uvs2[i][1]) and
      // check the inner statistics
      'subset'(S,
        '&join'(S, mset_apply(b, func<j | '*'(S, uvs1[j][2], uvs1[j][1])>)),
        '*'(S, uvs2[i][2], uvs2[i][1])
      )
    >
  );
  return #bs ne 0;
end intrinsic;

intrinsic 'subset'(S :: PGGStat_Degree, v1, v2) -> BoolElt
  {"}
  return v1`value eq v2`value;
end intrinsic;

intrinsic 'subset'(S :: PGGStat_Tup, v1, v2) -> BoolElt
  {"}
  return forall{i : i in [1..#S`stats] | 'subset'(S`stats[i], v1`value[i], v2`value[i])};
end intrinsic;

intrinsic 'subset'(S :: PGGStat_Stab, v1, v2) -> BoolElt
  {"}
  return 'subset'(S`stat, v1`value, v2`value);
end intrinsic;

// intrinsic 'subset'(S :: PGGStat_Factors, v1, v2) -> BoolElt
//   {"}
//   w1s := SetToSequence(MultisetToSet(v1`value));
//   w2s := SetToSequence(MultisetToSet(v2`value));
//   bs := all_binnings(
//     [Multiplicity(v1`value, w1) : w1 in w1s],
//     [Multiplicity(v2`value, w2) : w2 in w2s]:
//     limit:=1,
//     is_semivalid:=func<i,b | PartialLe(ProductStat(S`stat, mset_apply(b, func<j | w1s[j]>)), w2s[i])>,
//     is_valid:=func<i,b | ProductStat(S`stat, mset_apply(b,func<j | w1s[j]>)) subset w2s[i]>
//   );
//   return #bs gt 0;
// end intrinsic;

intrinsic 'subset'(S :: PGGStat_Factors2, v1, v2) -> BoolElt
  {"}
  // v`value[1][i] corresponds to factor i
  // v`value[2][i][i] corresponds to the stabilizer of factor i
  // v`value[2][i][j] for i ne j do not behave well in the ordering, so are ignored 
  S2 := PGGStat_Stab_Make(S`stat2);
  S3 := PGGStat_Factors_Make(S2);
  S4 := PGGStat_Tup_Make(<S`stat1, S3>);
  S5 := PGGStat_Factors_Make(:Stat:=S4);
  m := func<v | MakeValue(S5, {* MakeValue(S4, <v`value[1][i], MakeValue(S3, {*MakeValue(S2,x) : x in v`value[2][i][i]*})>) : i in [1..n]*}) where n:=#v`value[1]>;
  return m(v1) subset m(v2);
end intrinsic;

intrinsic 'subset'(S :: PGGStat_GaloisGroup, v1, v2) -> BoolElt
  {"}
  return Degree(v1`value) eq Degree(v2`value) and IsConjugateSubgroup(SymmetricGroup(Degree(v1`value)), v2`value, v1`value);
end intrinsic;

intrinsic 'subset'(S :: PGGStat_SubfieldDegrees, v1, v2) -> BoolElt
  {"}
  return v1`value[1] eq v2`value[1] and v2`value[2] subset v1`value[2];
end intrinsic;

intrinsic '&join'(S :: PGGStat_FactorDegreesSeq, vs :: SetMulti) -> PGGStatVal
  {The values of this statistic are multisets; this joins them together.}
  norbits := Rep(vs)`norbits;
  return MakeValue(S, &join mset_apply(vs, func<v | v`value>), norbits);
end intrinsic;

intrinsic '*'(S :: PGGStat_FactorDegreesSeq, v :: PGGStatVal, d :: RngIntElt) -> PGGStatVal
  {Multiplies the degrees by d.}
  if v`norbits eq 1 then
    return MakeValue(S, mset_apply(v`value, func<x | x*d>));
  else
    return MakeValue(S, mset_apply(v`value, func<x | <x[1]*d, x[2]>>), v`norbits);
  end if;
end intrinsic;

// intrinsic ProductStat(S :: PGGStat_Degree, vs :: {*PGGStatVal*}) -> PGGStatVal
//   {Statistict of the direct product of groups or product of polynomials.}
//   return MakeValue(S, &+ChangeUniverse(mset_apply(vs, func<v | v`value>), Integers()));
// end intrinsic;

// intrinsic ProductStat(S :: PGGStat_HasRoot, vs :: {*PGGStatVal*}) -> PGGStatVal
//   {"}
//   return MakeValue(S, exists{v : v in MultisetToSet(vs) | v`value});
// end intrinsic;

// intrinsic ProductStat(S :: PGGStat_NumRoots, vs :: {*PGGStatVal*}) -> PGGStatVal
//   {"}
//   return MakeValue(S, &+ChangeUniverse(mset_apply(vs, func<v | v`value>), Integers()));
// end intrinsic;

// intrinsic ProductStat(S :: PGGStat_Tup, vs :: {*PGGStatVal*}) -> PGGStatVal
//   {"}
//   return MakeValue(S, <ProductStat(S`stats[i], mset_apply(vs, func<v | v`value[i]>)) : i in [1..#S`stats]>);
// end intrinsic;

// intrinsic ProductStat(S :: PGGStat_FactorDegrees, vs :: {*PGGStatVal*}) -> PGGStatVal
//   {"}
//   return MakeValue(S, &join mset_apply(vs, func<v | v`value>));
// end intrinsic;

// intrinsic ProductStat(S :: PGGStat_Factors, vs :: {*PGGStatVal*}) -> PGGStatVal
//   {"}
//   return MakeValue(S, &join mset_apply(vs, func<v | v`value>));
// end intrinsic;

// intrinsic ProductStat(S :: PGGStat_AutGroup, vs :: {*PGGStatVal*}) -> PGGStatVal
//   {"}
//   aut := DirectProduct([Rep(v`value) : v in vs]);
//   return MakeValue(S, ConjugacyClasses(S, Degree(aut)) ! aut);
// end intrinsic;

// intrinsic ProductStat(S :: PGGStat_Stab, vs :: {*PGGStatVal*}) -> PGGStatVal
//   {"}
//   return MakeValue(S, ProductStat(S`stat, mset_apply(vs, func<v | v`value>)));
// end intrinsic;

// intrinsic ProductStat(S :: PGGStat_GaloisGroup, vs :: {*PGGStatVal*}) -> PGGStatVal
//   {"}
//   return MakeValue(S, DirectProduct([v`value : v in vs]));
// end intrinsic;

// intrinsic PartialLe(v1 :: PGGStatVal, v2 :: PGGStatVal) -> BoolElt
//   {True if G1 is a subgroup the action of G2 on some of its orbits, and v1 and v2 are corresponding statistics.}
//   require v1`stat eq v2`stat: "must be values of the same statistic";
//   return PartialLe(v1`stat, v1, v2);
// end intrinsic;

// intrinsic PartialLe(S :: PGGStat, v1, v2) -> BoolElt
//   {"}
//   error "not implemented: PartialLe:", Type(S);
// end intrinsic;

// intrinsic PartialLe(S :: PGGStat_HasRoot, v1, v2) -> BoolElt
//   {"}
//   return true;
// end intrinsic;

// intrinsic PartialLe(S :: PGGStat_NumRoots, v1, v2) -> BoolElt
//   {"}
//   return true;
// end intrinsic;

// intrinsic PartialLe(S :: PGGStat_Degree, v1, v2) -> BoolElt
//   {"}
//   return v1`value le v2`value;
// end intrinsic;

// intrinsic PartialLe(S :: PGGStat_Tup, v1, v2) -> BoolElt
//   {"}
//   return forall{i : i in [1..#S`stats] | PartialLe(v1`value[i], v2`value[i])};
// end intrinsic;

// intrinsic PartialLe(S :: PGGStat_FactorDegrees, v1, v2) -> BoolElt
//   {"}
//   d1s := SetToSequence(MultisetToSet(v1`value));
//   d2s := SetToSequence(MultisetToSet(v2`value));
//   bs := all_binnings(
//     [Multiplicity(v1`value, d1) : d1 in d1s],
//     [Multiplicity(v2`value, d2) : d2 in d2s]:
//     limit:=1,
//     is_semivalid:=func<i,b | &+mset_apply(b,func<j | d1s[j]>) le d2s[i]>
//   );
//   return #bs gt 0;
// end intrinsic;

// intrinsic PartialLe(S :: PGGStat_Stab, v1, v2) -> BoolElt
//   {"}
//   return PartialLe(v1`value, v2`value);
// end intrinsic;

intrinsic PossibleIntermediateStatistics(S1 :: PGGStat, S2 :: PGGStat) -> []
  {A sequence of statistics to try deducing to from S1, and to S2 from.}
  return [];
end intrinsic;

intrinsic PossibleIntermediateStatistics(S1 :: PGGStat, S2 :: PGGStat_HasRoot) -> []
  {"}
  return [PGGStat_NumRoots_Make()];
end intrinsic;

intrinsic PossibleIntermediateStatistics(S1 :: PGGStat, S2 :: PGGStat_NumRoots) -> []
  {"}
  return [PGGStat_FactorDegrees_Make()];
end intrinsic;

intrinsic PossibleIntermediateStatistics(S1 :: PGGStat, S2 :: PGGStat_FactorDegrees) -> []
  {"}
  return [PGGStat_Factors_Make(:Stat:=PGGStat_Degree_Make())];
end intrinsic;

intrinsic PossibleIntermediateStatistics(S1 :: PGGStat, S2 :: PGGStat_Degree) -> []
  {"}
  return [PGGStat_FactorDegrees_Make()];
end intrinsic;

intrinsic Implies(S1 :: PGGStat, S2 :: PGGStat) -> BoolElt, .
  {True if values of S2 can be deduced from values of S1. If so, also returns the corresponding map.}
  if S1 eq S2 then
    return true, func<v | v>;
  end if;
  for S in PossibleIntermediateStatistics(S1,S2) do
    ok, m1 := Implies(S1, S);
    if ok then
      ok, m2 := Implies(S, S2);
      if ok then
        return true, func<v | m2(m1(v))>;
      end if;
    end if;
  end for;
  return false, _;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_NumRoots, S2 :: PGGStat_HasRoot) -> BoolElt, .
  {"}
  return true, func<v | MakeValue(S2, v`value gt 0)>;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_FactorDegrees, S2 :: PGGStat_NumRoots) -> BoolElt, .
  {"}
  return true, func<v | MakeValue(S2, Multiplicity(v`value, 1))>;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_Factors, S2 :: PGGStat_FactorDegrees) -> BoolElt, .
  {"}
  ok, m := Implies(S1`stat, PGGStat_Degree_Make());
  if ok then
    return true, func<v | MakeValue(S2, mset_apply(v`value, func<w | m(w)`value>))>;
  end if;
  return false, _;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_FactorDegrees, S2 :: PGGStat_Degree) -> BoolElt, .
  {"}
  return true, func<v | MakeValue(S2, &+v`value)>;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_SubfieldDegrees, S2 :: PGGStat_Degree) -> BoolElt, .
  {"}
  return true, func<v | MakeValue(S2, &+[Z| x[1] : x in v`value])>;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_GaloisGroup, S2 :: PGGStat_Degree) -> BoolElt, .
  {"}
  return true, func<v | MakeValue(S2, Degree(v`value))>;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_FactorDegrees, S2 :: PGGStat_Factors) -> BoolElt, .
  {"}
  SD := PGGStat_Degree_Make();
  ok, m := Implies(SD, S2`stat);
  if ok then
    return true, func<v | MakeValue(S2, mset_apply(v`value, func<d | m(MakeValue(SD,d))>))>;
  end if;
  return false, _;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_Factors, S2 :: PGGStat_Factors) -> BoolElt, .
  {"}
  if S1 eq S2 then
    return true, func<v | v>;
  end if;
  ok, m := Implies(S1`stat, S2`stat);
  if ok then
    return true, func<v | MakeValue(S2, mset_apply(v`value, m))>;
  end if;
  return false, _;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_Tup, S2 :: PGGStat_Tup) -> BoolElt, .
  {"}
  if S1 eq S2 then
    return true, func<v | v>;
  end if;
  ms := [];
  for S in S2`stats do
    ok, m := Implies(S1, S);
    if ok then
      Append(~ms, m);
    else
      return false, _;
    end if;
  end for;
  return true, func<v | MakeValue(S2, <m(v) : m in ms>)>;
end intrinsic;

intrinsic Implies(S1 :: PGGStat, S2 :: PGGStat_Tup) -> BoolElt, .
  {"}
  ms := [];
  for S in S2`stats do
    ok, m := Implies(S1, S);
    if ok then
      Append(~ms, m);
    else
      return false, _;
    end if;
  end for;
  return true, func<v | MakeValue(S2, <m(v) : m in ms>)>;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_Tup, S2 :: PGGStat) -> BoolElt, .
  {"}
  for i in [1..#S1`stats] do
    ok, m := Implies(S1`stats[i], S2);
    if ok then
      return true, func<v | m(v`value[i])>;
    end if;
  end for;
  return false, _;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_Stab, S2 :: PGGStat_Stab) -> BoolElt, .
  {"}
  if S1 eq S2 then
    return true, func<v | v>;
  end if;
  ok, m := Implies(S1`stat, S2`stat);
  if ok then
    return true, func<v | MakeValue(S2, m(v`value))>;
  end if;
  return false, _;
end intrinsic;

intrinsic Implies(S1 :: PGGStat_Factors2, S2 :: PGGStat_Factors) -> BoolElt, .
  {"}
  ok, m := Implies(S1`stat1, S2`stat);
  if ok then
    return true, func<v | MakeValue(S2, {*m(x) : x in v`value[1]*})>;
  end if;
  return false, _;
end intrinsic;

// Is it true that if H<G then N_G(Stab_G(1))/Stab_G(1) < N_H(Stab_H(1))/Stab_H(1)?
// AG < AH <=> NG/SG < NH/SH = NH/(SG meet NH) ~ (NH SG)/SG <=> NG < NH SG ?
// It is true for all transitive groups of degree up to 15.
// 
// SG = {s in G : s(1)=1}
// SH = {s in H : s(1)=1} = SG meet H
// NG = {n in G : s in SG => s^n in SG} = {n in G : s in G, s(1)=1 => (s^n)(1)=1 i.e. s(n(1))=n(1)}
// 
// intrinsic 'subset'(S :: PGGStat_AutGroup, v1, v2) -> BoolElt
//   {"}
// end intrinsic;

intrinsic _Hashable(S :: PGGStat, v) -> .
  {A hashable value derived from v.}
  not_implemented("_Hashable:", Type(S));
end intrinsic;

intrinsic _Hashable(S :: PGGStat_EqSimple, v) -> .
  {"}
  return v`value;
end intrinsic;

intrinsic _Hashable(S :: PGGStat_AutGroup, v) -> .
  {"}
  return <Degree(Rep(v`value)), v`value>;
end intrinsic;

intrinsic _Hashable(S :: PGGStat_Factors2, v) -> .
  {"}
  xs, yss := Explode(v`value);
  n := #xs;
  return 
    < {* x : x in xs *}
    , {* <xs[i], {* y : y in yss[i] *}> : i in [1..n] *}
    , {* <xs[i], {* ys[i] : ys in yss *}> : i in [1..n] *}
    >;
end intrinsic;

intrinsic _Hashable(S :: PGGStat_GaloisGroup, v) -> .
  {"}
  return <Degree(v`value), #v`value, Hash(v`value)>;
end intrinsic;

intrinsic _Print(S :: PGGStat, v)
  {"}
  printf "%o", v`value;
end intrinsic;

intrinsic _Print(S :: PGGStat_AutGroup, v)
  {"}
  printf "%o", GroupName(Rep(v`value));
end intrinsic;

function stabilizers(G, S, v)
  Gv := GroupStat(S, G);
  if v eq Gv then
    return [G];
  elif v subset Gv then
    Hs := [x`subgroup : x in MaximalSubgroups(G)];
    return #Hs eq 0 select [] else &cat[stabilizers(H, S, v) : H in Hs];
  else
    return [];
  end if;
end function;

function stabilizers_seq(G, os, S, v)
  Gv := GroupStat(S, G, os);
  if v eq Gv then
    return [G];
  elif v subset Gv then
    Hs := [x`subgroup : x in MaximalSubgroups(G)];
    return #Hs eq 0 select [] else &cat[stabilizers_seq(H, os, S, v) : H in Hs];
  else
    return [];
  end if;
end function;

intrinsic Stabilizers(G :: GrpPerm, v :: PGGStatVal) -> []
  {Minimal subgroups of G whose statistic is v. One per conjugacy class.}
  return _Stabilizers(G, v`stat, v);
end intrinsic;

intrinsic Stabilizers(G :: GrpPerm, os :: [GSetIndx[RngIntElt]], v :: PGGStatVal) -> []
  {"}
  assert #os eq v`norbits;
  return _Stabilizers(G, os, v`stat, v);
end intrinsic;

intrinsic _Stabilizers(G :: GrpPerm, os :: [], S :: PGGStat, v) -> []
  {"}
  not_implemented("_Stabilizers:", Type(S));
end intrinsic;

intrinsic _Stabilizers(G :: GrpPerm, S :: PGGStat, v) -> []
  {"}
  S2s := [* <S2,m> : S2 in [* PGGStat_HasRoot_Make(), PGGStat_FactorDegrees_Make() *] | ok where ok,m:=Implies(S,S2) *];
  vprint PGG_GaloisGroup: "WARNING: using generic stabilizers algorithm";
  if #S2s gt 0 then
    vprint PGG_GaloisGroup: "...after using specialist algorithms for:";
    IndentPush();
    for x in S2s do
      vprint PGG_GaloisGroup: x[1];
    end for;
    IndentPop();
  end if;
  Gs := [G];
  for x in S2s do
    if #Gs eq 0 then
      break;
    end if;
    S2, m := Explode(x);
    v2 := m(v);
    Gs := &cat[_Stabilizers(G, S2, v2) : G in Gs];
  end for;
  if #Gs eq 0 then
    return [];
  else
    return &cat[stabilizers(G, S, v) : G in Gs];
  end if;
end intrinsic;

intrinsic _Stabilizers(G :: GrpPerm, S :: PGGStat_HasRoot, v) -> []
  {"}
  return v`value select [Stabilizer(G,Rep(o)) : o in Orbits(G)] else [G];
end intrinsic;

function partition_system(G)
  parts := AllPartitions(G) join {{1}, {1..Degree(G)}};
  sys := [{1}];
  while true do
    last := sys[#sys];
    parts := [P : P in parts | last ne P and last subset P];
    if #parts eq 0 then
      assert last eq {1..Degree(G)};
      return sys;
    else
      _, i := Min([#P : P in parts]);
      Append(~sys, parts[i]);
    end if;
  end while;
end function;

intrinsic PGG_DirectProductEmbedding(G :: GrpPerm) -> []
  {A sequence of transitive groups Gs and a permutation g such that G^g subset DirectProduct(Gs).}
  d := Degree(G);
  Sd := SymmetricGroup(d);
  if IsTransitive(G) then
    Gs := [G];
    g := Id(Sd);
  else
    os := [[x:x in o]:o in Orbits(G)];
    g := (Sd ! [x : x in o, o in os])^-1;
    Gs := [OrbitImage(G^g, {@ x^g : x in o @}) : o in os];
  end if;
  for i in [1..#Gs] do
    ReduceGenerators(~Gs[i]);
  end for;
  assert G^g subset DirectProduct(Gs);
  return Gs, g;
end intrinsic;

intrinsic PGG_WreathProductEmbedding(G :: GrpPerm) -> [], GrpPermElt
  {A sequence of groups Gs and a permutation g such that G^g subset WreathProduct(Gs).}
  d := Degree(G);
  Sd := SymmetricGroup(d);
  if IsPrimitive(G) then
    // terminal case: no non-trivial partitions
    ReduceGenerators(~G);
    Gs := [G];
    g := Id(Sd);
  else
    // recursive case
    // select a minimal partition (there is a choice here, maybe some choices are better than others?)
    P := MinimalPartition(G);
    n := #P;
    Sn := SymmetricGroup(n);
    ok, m := IsDivisibleBy(d, n);
    assert ok;
    Sm := SymmetricGroup(m);
    assert forall{B : B in P | #B eq m};
    // choose an ordering on each block; the choice of ordering determines how small the resulting wreath product will be; taking an arbitrary order is bad, because even if the factor should be cyclic, it could end up huge; on the other hand, we don't bother to find the smallest possible wreath product, which could be infeasible to compute; instead we order the first block arbitrarily, and then order the remaining blocks 
    Pseq := [#Self() eq 0 select [x : x in B] else ([x^g : x in B1] where ok,g:=IsConjugate(G,B1[1],Rep(B)) where B1:=Self(1)) : B in P];
    // find g1 such that G^g1 has the partition [[1..m],[m+1..2m],..,[..d]]
    g1 := (Sd ! [Pseq[i][j] : j in [1..m], i in [1..n]])^-1;
    G1 := G^g1;
    P1seq := [[(i-1)*m+1..i*m] : i in [1..n]];
    P1 := [SequenceToSet(B) : B in P1seq];
    assert forall{B : B in P1 | IsBlock(G1, B)};
    b := BlocksAction(G1, P1);
    Gquo := G1 @ b;
    Glift := WreathProduct(sub<Sm | Id(Sm)>, Gquo);
    lift := map<Sn -> Sd | g :-> [j+(i^g-1)*m : j in [1..m], i in [1..n]]>;
    topgens := [Sm|];
    for i in [1..Ngens(G1)] do
      gen := G1.i;
      gen0 := gen * (gen @ b @ lift)^-1;
      topgens cat:= [Sm| [(j+(i-1)*m)^gen0 - (i-1)*m : j in [1..m]] : i in [1..n]];
    end for;
    Gtop := sub<Sm | topgens>;
    assert G1 subset WreathProduct(Gtop, Gquo);
    // recurse on the quotient
    Gquos, gquo := PGG_WreathProductEmbedding(Gquo);
    glift := gquo @ lift;
    ReduceGenerators(~Gtop);
    Gs := [Gtop] cat Gquos;
    g := g1 * glift;
  end if;
  assert G^g subset WreathProduct(Gs);
  return Gs, g;
end intrinsic;

function is_partition_stabilizer(G, S)
  return Stabilizer(G, [{x:x in o} : o in Orbits(S)]) eq S;
end function;

function is_partition_stabilizer_of_shape(G, S, V)
  return ({*#o : o in Orbits(S)*} eq V) and is_partition_stabilizer(G, S);
end function;

function partitions_of_shape(G, V : simple:=false)
  os := Orbits(G);
  sh := {*#o:o in os*};
  if sh eq V then
    return [{{x:x in o} : o in os}];
  elif is_subpartition_of(V, sh) then
    if simple then
      return SetToSequence({{{x:x in o}:o in os} : x in Subgroups(G) | {*#o:o in os*} eq V where os:=Orbits(x`subgroup)});
    else
      return [P : P in partitions_of_shape(x`subgroup, V), x in MaximalSubgroups(G)];
    end if;
  else
    return [];
  end if;
end function;

function partitions_above_shape(G, V, d : avoid:=false)
  os := Orbits(G);
  sh := {*d*#o:o in os*};
  if sh eq V then
    return [{{x:x in o}:o in os}];
  elif is_subpartition_of(V, sh) then
    return (avoid cmpeq false or avoid ne os select [{{x:x in o}:o in os}] else []) cat [P : P in partitions_above_shape(x`subgroup, V, d : avoid:=os), x in MaximalSubgroups(G)];
  else
    return [];
  end if;
end function;

function partitions_of_wr_above_shape(Gs, V, d)
  assert d*&*[Degree(G):G in Gs] eq &+V;
  if #Gs eq 1 then
    return partitions_above_shape(Gs[1], V, d);
  end if;
  // partitions on the bottom of the wreath product
  P0s := partitions_of_wr_above_shape(Gs[2..#Gs], V, d*Degree(Gs[1]));
  // partitions on the top
  P1s := [{{x:x in o}:o in Orbits(S)} : x in Subgroups(Gs[1]) | is_partition_stabilizer(Gs[1],S) where S:=x`subgroup];
  Ps := [];
  D := Degree(Gs[1]);
  for P0 in P0s do
    P0seq := [B : B in P0];
    for x in CartesianPower(P1s, #P0seq) do
      P := {{y+(x-1)*D : y in P1, x in P0} where P0:=P0seq[i] : P1 in x[i], i in [1..#P0seq]};
      if is_subpartition_of(V, {*d*#o : o in P*}) then
        Append(~Ps, P);
      end if;
    end for;
  end for;
  return Ps;
end function;

function partitions_of_wr_of_shape(Gs, V)
  assert #Gs gt 0;
  if #Gs eq 1 then
    return partitions_of_shape(Gs[1], V);
  end if;
  return [P : P in partitions_of_wr_above_shape(Gs, V, 1) | {*#o : o in P*} eq V];
end function;

intrinsic _Stabilizers(G :: GrpPerm, S :: PGGStat_FactorDegrees, v) -> []
  {"}
  // We embed G into a wreath product, solve the problem there, and then intersect conjugates of those stabilizers with G
  // TODO: this only uses one system of partitions for G to find a wreath product; the information from the other partitions could probably be used, though this would make the algorithm more complicated, since the partitions used now form a DAG not a single path; might be worth doing if this routine becomes the bottleneck
  V := v`value;
  d := Degree(G);
  assert &+V eq d;
  os := [{x:x in o} : o in Orbits(G)];
  assert is_subpartition_of(V, {* #o : o in os *});
  if V eq {* #o : o in os *} then
    return [G];
  end if;
  // embed G into a direct product of transitive groups
  PGG_GlobalTimer_Push("direct product embedding");
  Ds, gD := PGG_DirectProductEmbedding(G);
  gDinv := gD^-1;
  // find the ways in which sub-partitions of V correspond to the factors
  PGG_GlobalTimer_Swap("matching");
  vprint PGG_GaloisGroup: "matching", V, "to", [Degree(D):D in Ds];
  Vss := PGG_all_partition_groupings(V, [Degree(D) : D in Ds]);
  vprint PGG_GaloisGroup: "#Vss =", #Vss;
  Sss := [];
  // for each direct factor, find all its stabilizers for each relevant partition
  PGG_GlobalTimer_Swap("solve each factor");
  for i in [1..#Ds] do
    D := Ds[i];
    PGG_GlobalTimer_Push("wreath product embedding");
    Gs, g := PGG_WreathProductEmbedding(D);
    vprint PGG_GaloisGroup: "#Gs =", #Gs;
    vprint PGG_GaloisGroup: "#W =", #WreathProduct(Gs);
    vprint PGG_GaloisGroup: "#D =", #D;
    PGG_GlobalTimer_Pop();
    ginv := g^-1;
    Sss[i] := AssociativeArray();
    for VD in {Vs[i] : Vs in Vss} do
      if (#D)^3 lt #WreathProduct(Gs) then
        PGG_GlobalTimer_Push("partitions of small group");
        Ss := [Stabilizer(D, SetToSequence(P)) : P in partitions_of_shape(D, VD : simple)];
        PGG_GlobalTimer_Pop();
      else
        PGG_GlobalTimer_Push("partitions of wr");
        P0s := [{{x^ginv : x in B} : B in P} : P in partitions_of_wr_of_shape(Gs, VD)];
        PGG_GlobalTimer_Pop();
        assert forall{P0 : P0 in P0s | {*#o : o in P0*} eq VD};
        W := WreathProduct(Gs)^ginv;
        assert D subset W;
        S0s := [Stabilizer(W, SetToSequence(P0)) : P0 in P0s];
        assert forall{i : i in [1..#P0s] | P0s[i] eq {{x : x in o} : o in Orbits(S0s[i])}};
        Ss := [];
        PGG_GlobalTimer_Push("images");
        for S0 in S0s do
          PGG_GlobalTimer_Push("double coset representatives");
          ws := DoubleCosetRepresentatives(W, Normalizer(W,S0), D);
          vprint PGG_GaloisGroup: "#ws =", #ws;
          PGG_GlobalTimer_Pop();
          for w in ws do
            PGG_GlobalTimer_Push("meet");
            S := (S0^w) meet D;
            PGG_GlobalTimer_Swap("orbits");
            VS := {* #o : o in Orbits(S) *};
            PGG_GlobalTimer_Pop();
            if VS eq VD then
              PGG_GlobalTimer_Push("dedupe conjugates");
              for S2 in Ss do
                if IsConjugate(D, S2, S) then
                  PGG_GlobalTimer_Pop();
                  continue w;
                end if;
              end for;
              PGG_GlobalTimer_Pop();
              Append(~Ss, S);
            else
              ; // the following assert can be expensive
              // PGG_GlobalTimer_Push("assert is subpartition");
              // assert is_subpartition_of(VS, VD);
              // PGG_GlobalTimer_Pop();
            end if;
          end for;
        end for;      
        PGG_GlobalTimer_Pop();
      end if;
      Sss[i][VD] := Ss;
    end for;
  end for;
  // take combinations to find the overall factors
  PGG_GlobalTimer_Swap("combinations");
  Ss := [];
  W := DirectProduct(Ds)^gDinv;
  assert G subset W;
  for Vs in Vss do
    for x in CartesianProduct(<Sss[i][Vs[i]] : i in [1..#Ds]>) do
      S0 := DirectProduct([S : S in x])^gDinv;
      assert S0 subset W;
      for w in DoubleCosetRepresentatives(W, Normalizer(W,S0), G) do
        S := (S0^w) meet G;
        VS := {* #o : o in Orbits(S) *};
        if VS eq V then
          PGG_GlobalTimer_Push("dedupe conjugates");
          for S2 in Ss do
            if IsConjugate(G, S2, S) then
              PGG_GlobalTimer_Pop();
              continue w;
            end if;
          end for;
          PGG_GlobalTimer_Pop();
          Append(~Ss, S);
        else
          assert is_subpartition_of(VS, V);
        end if;
      end for;
    end for;
  end for;
  PGG_GlobalTimer_Pop();
  return Ss;
end intrinsic;

intrinsic _Stabilizers(G :: GrpPerm, S :: PGGStat_FactorDegreesSeq, v) -> []
  {"}
  return _Stabilizers(G, [GSet(G)], S, v);
end intrinsic;

intrinsic _Stabilizers(G :: GrpPerm, os :: [], S :: PGGStat_FactorDegreesSeq, v) -> []
  {"}
  // first coerce the statistic to FactorDegrees, and get stabilizers there
  o1 := os[1];
  q1 := OrbitAction(G, o1);
  G1 := G @ q1;
  S1 := PGGStat_FactorDegrees_Make();
  v1 := MakeValue(S1, #os eq 1 select v`value else mset_apply(v`value, func<x | x[1]>));
  Hs1 := [H @@ q1 : H in Stabilizers(G1, v1)];
  // now use the generic algorithm
  if (#os eq 1) or (#Hs1 eq 0) then
    return Hs1;
  else
    vprint PGG_GaloisGroup: "WARNING: PGGStat_FactorDegreesSeq: using generic stabilizers algorithm";
    return &cat[stabilizers_seq(H, os, S, v) : H in Hs1];
  end if;
end intrinsic;