Utils.mag

Z := IntegerRing();
Q := RationalField();

function not_implemented(msg, ...)
  error "not implemented: " cat Join([Sprintf("%o", x) : x in msg], " ");  
end function;

function precision_error()
  error "precision error";
end function;

function user_error(msg, ...)
  error "error: " cat Join([Sprintf("%o", x) : x in msg], " ");
end function;

procedure isort_by(~xs, key)
  keys := [key(xs[i]) : i in [1..#xs]];
  Sort(~keys, ~permutation);
  xs := [xs[i^permutation] : i in [1..#xs]];
end procedure;

function sort_by(xs, key)
  isort_by(~xs, key);
  return xs;
end function;

function group_by(xs, key)
  A := AssociativeArray();
  for x in xs do
    k := key(x);
    if IsDefined(A, k) then
      Append(~A[k], x);
    else
      A[k] := [x];
    end if;
  end for;
end function;

function enumerate(xs)
  return [<xs[i], i> : i in [1..#xs] | IsDefined(xs, i)];
end function;

function fldpad_eltseq(x)
  if AbsolutePrecision(x) eq Infinity() then
    assert IsWeaklyZero(x);
    return [BaseField(Parent(x))| 0 : i in [1..Degree(Parent(x))]];
  else
    return Eltseq(x);
  end if;
end function;

function is_any_coercible(coercers, x)
  for i in [1..#coercers] do
    ok, y := coercers[i](x);
    if ok then
      return true, i, y;
    end if;
  end for;
  return false, _, _;
end function;

function all_coercible(dflts, coercers, xs)
  assert #dflts eq #coercers;
  dflts := [* d : d in dflts *];
  for x in xs do
    ok, i, y := is_any_coercible(coercers, x);
    if ok then
      dflts[i] := y;
    else
      return false, _;
    end if;
  end for;
  return true, dflts;
end function;

function the(xs)
  assert #xs eq 1;
  return [x : x in xs][1];
end function;

// the roots of f as a sequence (no multiplicities)
// f must be squarefree
function roots(f : Lift:=true)
  return PGG_Roots(f : Lift:=Lift);
end function;

// true if f has a root
// f must be squarefree
// also returns a root
function has_root(f)
  rs := roots(f);
  if #rs eq 0 then
    return false, _;
  else
    return true, rs[1];
  end if;
end function;

// the factorization of f as a sequence (no multiplicities)
// f must be squarefree
// also the certificates
function factorization(f : Extensions:=false, Lift:=true)
  return PGG_Factorization(f : Extensions:=Extensions, Lift:=Lift);
end function;

function extension(f)
  facs, certs := factorization(f : Extensions);
  error if #facs gt 1, "expecting at most 1 factor";
  return certs[1]`Extension;
end function;

function pol_invmod(g, f)
  d := Degree(f);
  // rows are coefficients of g^i mod f
  M := Matrix([[Coefficient(gg, j) : j in [0..d-1]] where gg := Polynomial([0 : j in [1..i]] cat Coefficients(g)) mod f : i in [0..d-1]]);
  // rows are linear combinations of g^i summing to x^j
  ok, Minv := IsInvertible(M);
  assert ok;
  // first row is h(x) so that h(x)g(x)=1 mod f(x)
  return Polynomial(Eltseq(Rows(Minv)[1]));
end function;

// the factorization of f into pieces corresponding to the slopes of its newton polygon as a sequence
function newton_polygon_factorization(f : alg:="Hensel")
  // special cases
  d := Degree(f);
  error if d lt 0, "f must be non-zero";
  error if IsWeaklyZero(Coefficient(f, d)), "precision error: leading coefficient is weakly zero";
  if d eq 0 then
    return [];
  end if;
  case alg:
  when "Cheat":
    // we cheat and do a full factorization and then aggregate the results
    // this will be very inefficient for some polynomials!
    facs := factorization(f);
    a := AssociativeArray();
    for fac in facs do
      s := the(Slopes(NewtonPolygon(fac)));
      if IsDefined(a, s) then
        Append(~a[s], fac);
      else
        a[s] := [fac];
      end if;
    end for;
    return [&*a[s] : s in Sort([s : s in Keys(a)])];
  when "Hensel":
    // use Hensel's lemma
    R := Parent(f);
    K := BaseRing(R);
    pr := Precision(K);
    assert pr lt Infinity();
    np := NewtonPolygon(f);
    vs := ChangeUniverse(Vertices(np), car<Z,Z>);
    assert #vs ge 2;
    assert vs[1][1] eq 0;
    assert vs[#vs][1] eq Degree(f);
    error if exists{v : v in vs | IsWeaklyZero(Coefficient(f, v[1]))}, "precision error: coefficient weakly zero";
    factors := [];
    for n in [2..#vs] do
      i0, v0 := Explode(vs[n-1]);
      i1, v1 := Explode(vs[n]);
      width := i1 - i0;
      slope := (v1 - v0) / width;
      h := -Numerator(slope);
      e := Denominator(slope);
      g0 := R ! [Coefficient(f, i) / Coefficient(f, i1) : i in [i0..i1]];
      g := R ! [IsWeaklyZero(c) select 0 else ChangePrecision(c, pr) : c in Coefficients(g0)];
      assert Degree(g) eq width;
      assert Valuation(Coefficient(g, 0)) eq v0-v1;
      assert Valuation(Coefficient(g, width)) eq 0;
      vprint PGG_GaloisGroup, 2: "hensel lifting slope =", slope;
      niters := 0;
      while true do
        niters +:= 1;
        if false then
          // use arithmetic in R mod g instead of reducing things mod g explicitly
          Q := quo<R | g>;
          PGG_GlobalTimer_Push("quotrem g");
          h, fmodg := Quotrem(f, g);
          PGG_GlobalTimer_Swap("gnew");
          gnew := g + R!(Q!fmodg/Q!h);
        else
          PGG_GlobalTimer_Push("quotrem g");
          h, fmodg := Quotrem(f, g);
          if true then
            // use Newton lifting to find hinv = h^-1 mod g after the first iteration
            // TODO: even better would be to only do pol_invmod() to limited precision at first, and then do many rounds of Newton lifting
            if niters eq 1 then
              PGG_GlobalTimer_Swap("invmod g");
              hinv := pol_invmod(h, g);
            else
              PGG_GlobalTimer_Swap("newton lift");
              hinv := (hinv * (2 - ((h * hinv) mod g))) mod g;
            end if;
          else
            PGG_GlobalTimer_Swap("invmod g");
            hinv := pol_invmod(h, g);
          end if;
          PGG_GlobalTimer_Swap("mult");
          tmp := fmodg * hinv;
          PGG_GlobalTimer_Swap("mod g");
          tmp := tmp mod g;
          PGG_GlobalTimer_Swap("add");
          gnew := g + tmp;
        end if;
        PGG_GlobalTimer_Swap("check");
        assert Degree(gnew) eq width;
        assert Valuation(Coefficient(g, 0)) eq v0-v1;
        assert Valuation(Coefficient(g, width)) eq 0;
        delta := g - gnew;
        PGG_GlobalTimer_Pop();
        if forall{i : i in [0..width] | IsWeaklyZero(Coefficient(delta, i))} then
          vprint PGG_GaloisGroup, 2: "f mod g =", fmodg;
          Append(~factors, g);
          continue n;
        else
          g := R ! [IsWeaklyZero(c) select 0 else ChangePrecision(c, pr) : c in Coefficients(gnew)];
        end if;
      end while;
    end for;
    return factors;
  else
    error "invalid algorithm:", alg;
  end case;
end function;

function multiplicities_to_information(mults)
  mults := [m : m in mults | m ne 0];
  assert forall{m : m in mults | m gt 0};
  assert #mults gt 0;
  if #mults eq 1 then
    return 0.0;
  end if;
  total := &+mults;
  ans := &+[-p*Log(2,p) where p:=m/total : m in mults];
  assert ans gt 0;
  return ans;
end function;

function largest_coefficient(x)
  if IsPrimeField(Parent(x)) then
    return Abs(x);
  else
    return Max([largest_coefficient(c) : c in Eltseq(x)]);
  end if;
end function;

function reduce_coefficients(x, M)
  if IsPrimeField(Parent(x)) then
    return Parent(x) ! ((Z!x) mod M);
  else
    return Parent(x) ! [reduce_coefficients(c, M) : c in Eltseq(x)];
  end if;
end function;

function startswith(x, y)
  if #x ge #y and x[1..#y] eq y then
    return true, x[#y+1..#x];
  else
    return false, _;
  end if;
end function;

procedure seq_idiff(~xs, ys)
  ys := {y : y in ys};
  xs := [xs[i] : i in [1..#xs] | xs[i] notin ys];
end procedure;

// a bound on |f(a_1,...,a_n)| where |a_i| le b.
// f is a RngSLPolElt
function slpol_bound(f, b)
  op := Operator(f);
  f1, f2 := Operands(f);
  case op:
  when "var":
    return b;
  when "const":
    return Abs(f1);
  when "+", "-":
    return slpol_bound(f1, b) + slpol_bound(f2, b);
  when "*":
    return slpol_bound(f1, b) * slpol_bound(f2, b);
  when "^":
    return slpol_bound(f1, b) ^ f2;
  else
    assert false;
  end case;
end function;

function permute_seq(g, xs)
  d := Degree(Parent(g));
  assert d eq #xs;
  return [xs[i^g] : i in [1..d]];
end function;

procedure pop_start(~x, ~xs)
  assert #xs ne 0;
  x := xs[1];
  xs := xs[2..#xs];
end procedure;

function polynomial_with_roots(rs)
  U := Universe(rs);
  R := PolynomialRing(U);
  return &*[R| [-r, 1] : r in rs];
end function;

function left_coset_representatives(G, H)
  assert H subset G;
  cosets, indices := DoubleCosetRepresentatives(G, H, sub<G | Id(G)>);
  return cosets;
end function;

function right_coset_representatives(G, H)
  assert H subset G;
  cosets, indices := DoubleCosetRepresentatives(G, sub<G | Id(G)>, H);
  return cosets;
end function;

function is_extension_of(L, K)
  if L eq K then
    return true, [K];
  elif IsPrimeField(L) then
    return false, _;
  else
    ok, twr := is_extension_of(BaseField(L), K);
    if ok then
      return true, Append(twr, L);
    else
      return false, _;
    end if;
  end if;
end function;

function tower(L, K)
  ok, twr := is_extension_of(L, K);
  assert ok;
  return twr;
end function;

function is_in_standard_form(L, K)
  t := tower(L, K);
  if #t eq 1 then
    return true, K;
  elif #t eq 2 then
    if RamificationDegree(t[2]) gt 1 then
      return true, K;
    elif InertiaDegree(t[2]) gt 1 then
      return true, L;
    end if;
  elif #t eq 3 then
    U := t[2];
    if RamificationDegree(L) gt 1 and InertiaDegree(U) gt 1 then
      return true, U;
    end if;
  end if;
  return false, _;
end function;

function valuation_eq(c, n)
  if Valuation(c) gt n then
    return false;
  elif IsWeaklyZero(c) then
    precision_error();
  else
    return Valuation(c) eq n;
  end if;
end function;

function valuation_ge(c, n)
  if Valuation(c) ge n then
    return true;
  elif IsWeaklyZero(c) then
    precision_error();
  else
    return false;
  end if;
end function;

function is_eisenstein(f)
  d := Degree(f);
  return (d ge 0) and valuation_eq(Coefficient(f, d), 0) and valuation_eq(Coefficient(f, 0), 1) and forall{i : i in [1..d-1] | valuation_ge(Coefficient(f, i), 1)};
end function;

function is_inertial(f)
  d := Degree(f);
  if exists{c : c in Coefficients(f) | AbsolutePrecision(c) lt 1} then
    precision_error();
  end if;
  return (d ge 1) and valuation_eq(Coefficient(f, d), 0) and valuation_eq(Coefficient(f, 0), 0) and forall{c : c in Coefficients(f) | valuation_ge(c, 0)} and IsIrreducible(Polynomial([c@m : c in Coefficients(f)]) where _,m:=ResidueClassField(Integers(BaseRing(f))));
end function;

function xdiv(x, y)
  ok, z := IsDivisibleBy(x, y);
  assert ok;
  return z;
end function;

// solves S*M=V for S
function solve(M, V)
  return V * M^-1;
  ok, S := IsSolvable(M, V);
  if ok then
    return S;
  else
    return V * M^-1;
  end if;
end function;

// given L/K and a sequence pis of elements of L which are uniformizing elements of subextensions K(pis[#pis])/.../K(pis[1]), return this tower from the bottom up
function tower_from_uniformizers(L, K, pis : alg:="LA")
  assert Universe(pis) eq L;
  assert not IsPrimeField(L);
  assert BaseField(L) eq K;
  t := [];
  K2 := K;
  L2 := L;
  LtoL2 := map<L -> L2 | x :-> x>;
  f := DefiningPolynomial(L);
  Lpi := L.1;
  for i in [1..#pis] do
    PGG_GlobalTimer_Push("minimal polynomial");
    pi2 := LtoL2(pis[i]);
    m := MinimalPolynomial(pi2, K2);
    PGG_GlobalTimer_Swap("factorization");
    assert is_eisenstein(m);
    d := Degree(m);
    n := xdiv(Degree(L2, K2), d);
    if d gt 1 then
      K2new := ext<K2 | ChangePrecision(m, Precision(K2))>;
      case alg:
      when "Factorization":
        f2 := factorization(ChangeRing(f, K2new))[1];
      when "LA":
        V := VectorSpace(K2, Degree(L2, K2));
        vmap := map<L2 -> V | x :-> V!Eltseq(x), y :-> L2!Eltseq(y)>;
        L2pi := LtoL2(Lpi);
        vec := vmap(L2pi^n);
        mat := Matrix([vmap(L2pi^i * pi2^j) : j in [0..d-1], i in [0..n-1]]);
        coeffs := Eltseq(solve(mat, vec));
        f2 := Polynomial([K2new| i eq n select 1 else -&+[coeffs[i*d+j+1] * K2new.1^j : j in [0..d-1]] : i in [0..n]]);
      else
        assert false;
      end case;
      assert Degree(f2) eq n;
      assert is_eisenstein(f2);
      L2new := ext<K2new | ChangePrecision(f2, Precision(K2new))>;
      LtoL2new := map<L -> L2new | x :-> &+[L2new| L2new.1^(i-1) * cs[i] : i in [1..#cs]] where cs:=Eltseq(x)>;
      // rename for iterating
      K2 := K2new;
      L2 := L2new;
      LtoL2 := LtoL2new;
    end if;
    Append(~t, K2);
    PGG_GlobalTimer_Pop();
  end for;
  return t;
end function;

function ramification_tower(L, K)
  ok, U := is_in_standard_form(L, K);
  assert ok;
  if L eq K then
    return [K];
  elif L eq U then
    return [K, U];
  elif U eq K then
    f := DefiningPolynomial(L);
    R := Parent(f);
    assert is_eisenstein(f);
    pi := L.1;
    assert IsWeaklyZero(Evaluate(f, pi));
    r := Evaluate(f, PolynomialRing(L) ! [pi, 1]);
    assert IsWeaklyZero(Coefficient(r, 0));
    PGG_GlobalTimer_Push("newton polygon factorization");
    rfacs := [R.1] cat newton_polygon_factorization(r div R.1);
    PGG_GlobalTimer_Swap("get uniformizers");
    assert Degree(rfacs[1]) eq 1;
    assert forall{i : i in [1..#rfacs] | IsDivisibleBy(Degree(f), &+[Degree(rfacs[j]) : j in [1..i]])};
    crfacs := [&*rfacs[1..i] : i in [1..#rfacs]];
    cffacs := [Evaluate(crfac, PolynomialRing(L) ! [-pi, 1]) : crfac in crfacs];
    pis := [Coefficient(cffac, 0) : cffac in cffacs];
    PGG_GlobalTimer_Swap("tower from uniformizers");
    twr := tower_from_uniformizers(L, K, Reverse(pis));
    PGG_GlobalTimer_Pop();
    return twr;
  else
    return [K] cat ramification_tower(L, U);
  end if;
end function;

function random_element(E, B)
  if IsPrimeField(E) then
    return E ! Random(B);
  else
    return E ! [random_element(BaseField(E), B) : i in [1..Degree(E)]];
  end if;
end function;

function random_primitive_element(E, F)
  d := Degree(E, F);
  if d eq 1 then
    return E!0, PolynomialRing(F)![0,1];
  end if;
  B := 1;
  while true do
    x := random_element(E, B);
    minpol := MinimalPolynomial(x, F);
    if Degree(minpol) eq d then
      return x, minpol;
    else
      B +:= 1;
    end if;
  end while;
end function;

function zero(K, apr)
  z := (K!1) - (K!1);
  return ShiftValuation(z, apr - AbsolutePrecision(z));
end function;

function change_apr(x, apr)
  if IsWeaklyZero(x) or apr le Valuation(x) then
    return zero(Parent(x), apr);
  else
    return ChangePrecision(x, apr - Valuation(x));
  end if;
end function;

function trim_apr(x, trim)
  return change_apr(x, AbsolutePrecision(x) - trim);
end function;

function maximize_apr(x)
  if IsWeaklyZero(x) then
    return Parent(x) ! 0;
  else
    return ChangePrecision(x, Precision(Parent(x)));
  end if;
end function;

function residually_primitive_element(K)
  F, m := ResidueClassField(Integers(K));
  return PrimitiveElement(F) @@ m;
end function;

function primitive_element(L, K)
  if Degree(L, K) eq 1 then
    return L!1;
  elif RamificationDegree(L, K) eq 1 then
    return residually_primitive_element(L);
  elif InertiaDegree(L, K) eq 1 then
    return UniformizingElement(L);
  else
    return UniformizingElement(L) + residually_primitive_element(L);
  end if;
end function;

// the following works, but looks like it was fixed anyway in V2.23-6
// function automorphism_group(L, K)
//   ok, twr := is_extension_of(L, K);
//   assert ok;
//   if #twr eq 1 or Degree(L, K) eq 1 then
//     return SymmetricGroup(1);
//   end if;
//   hs := Automorphisms(L, K);
//   us := [h(u) : h in hs] where u:=primitive_element(L, K);
//   huss := [[h(u) : u in us] : h in hs];
//   perms := [[the([i : i in [1..#us] | IsWeaklyEqual(hu, us[i])]) : hu in hus] : hus in huss];
//   return sub<SymmetricGroup(#us) | perms>;
// end function;

DEFAULT_TSCHIRNHAUS_TRANSFORMATIONS := [PolynomialRing(Z) | [0,1], [0,0,1], [0,-1,1], [0,1,1], [0,0,0,1], [0,1,0,1], [0,-1,0,1], [0,0,1,1], [0,0,-1,1], [0,1,1,1], [0,-1,-1,1]];

function tschirnhaus_transformation(ntries, degree)
  assert ntries gt 0;
  if ntries le #DEFAULT_TSCHIRNHAUS_TRANSFORMATIONS then
    return DEFAULT_TSCHIRNHAUS_TRANSFORMATIONS[ntries];
  else
    return PolynomialRing(Z) ! [Random(1,ntries) : i in [1..ntries-#DEFAULT_TSCHIRNHAUS_TRANSFORMATIONS+2]];
  end if;
end function;

function dflt(what, d)
  return (what cmpne false) select what else d;
end function;

function dedupe_conjugage_subgroups(G, Hs)
  CG := PGG_SubgroupClasses(G);
  return [Rep(c) : c in {CG ! H : H in Hs}];
end function;

function dedupe_conjugage_subgroups_simple(G, Hs)
  newHs := [];
  for H in Hs do
    if not exists{H2 : H2 in newHs | IsConjugate(G, H, H2)} then
      Append(~newHs, H);
    end if;
  end for;
  return newHs;
end function;

// non-checking version
// TODO: consider _is_subpartition_of({*1^^9*},{*2,3,4*}); currently we try assigning 1 to each of 2, 3, 4, and then assign 1 again, and keep going; but since we are assigning the same thing over and over, the order that these occur in is unimportant, so we are repeating work; hence, we should do a strategy which assigns all of the same element in one go
function _is_subpartition_of(P, Q)
  if #P eq 0 then
    assert #Q eq 0;
    return true;
  end if;
  assert #P gt 0;
  assert #Q gt 0;
  // pick an element of p
  // hopefully dealing with the largest ones first is optimal
  p := Max(P);
  P2 := Exclude(P, p);
  // pick an element of Q it could come from
  Rs := [];
  for q in MultisetToSet(Q) do
    if p le q then
      Q2 := Exclude(Q, q);
      if p lt q then
        Include(~Q2, q-p);
      end if;
      ok := _is_subpartition_of(P2, Q2);
      if ok then
        return true;
      end if;
    end if;
  end for;
  return false;
end function;

// sequence of different ks such that k = sum_i ks[i] ms[i]
function partitions_of_shape(k, ms)
  assert k ge 0;
  if k eq 0 then
    return [[0 : m in ms]];
  elif #ms eq 0 then
    return [];
  elif #ms eq 1 then
    ok, k1 := IsDivisibleBy(k, ms[1]);
    if ok then
      return [[k1]];
    else
      return [];
    end if;
  else
    m1 := ms[1];
    ret := [];
    for k1 in [0..k div m1] do
      k2 := k - k1 * m1;
      assert k2 ge 0;
      ret cat:= [[k1] cat p : p in partitions_of_shape(k2, ms[2..#ms])];
    end for;
    return ret;
  end if;
end function;

// given two multisets of positive integers, determine if P may be obtained by partitioning each element of Q
function is_subpartition_of(P, Q)
  assert Type(P) eq SetMulti;
  assert Type(Q) eq SetMulti;
  assert Universe(P) cmpeq Z;
  assert Universe(Q) cmpeq Z;
  assert forall{x : x in MultisetToSet(P) | x gt 0};
  assert forall{x : x in MultisetToSet(Q) | x gt 0};
  assert ((#P eq 0) and (#Q eq 0)) or ((#P gt 0) and (#Q gt 0) and (&+P eq &+Q));
  return _is_subpartition_of(P, Q);
end function;

function _all_partition_groupings(P, Q)
  ptodo := [<Q, [{**} : q in Q]>];
  for p in Reverse(Sort(SetToSequence(MultisetToSet(P)))) do
    n := Multiplicity(P, p);
    ptodo_new := [];
    for x in ptodo do
      Qrem, Ps := Explode(x);
      todo := [<n, []>];
      for q in Qrem do
        todo := [<x[1]-k, Append(x[2],k)> : k in [0..Min(x[1], q div p)], x in todo];
      end for;
      nss := [x[2] : x in todo | x[1] eq 0];
      for ns in nss do
        Append(~ptodo_new, <[Qrem[i]-ns[i]*p : i in [1..#Q]], [Ps[i] join {*p^^ns[i]*} : i in [1..#Q]]>);
      end for;
    end for;
    ptodo := ptodo_new;
  end for;
  return [x[2] : x in ptodo | x[1] eq [0 : q in Q]];
end function;

// given a multiset P of integers, and a sequence Q of integers, returns all sequences Ps of length #Q of multisets of integers so that &join Ps eq P and &+Ps[i] eq Qs[i]; that is, assigns the elements of P to bins of size Q
function all_partition_groupings(P, Q)
  assert Type(P) eq SetMulti;
  assert Type(Q) eq SeqEnum;
  assert Universe(P) eq Z;
  assert Universe(Q) eq Z;
  assert forall{x : x in MultisetToSet(P) | x gt 0};
  assert forall{x : x in Q | x gt 0};
  assert ((#P eq 0) and (#Q eq 0)) or ((#P gt 0) and (#Q gt 0) and (&+P eq &+Q));
  return _all_partition_groupings(P, Q);
end function;

intrinsic PGG_all_partition_groupings(P :: {* RngIntElt *}, Q :: [RngIntElt]) -> []
  {All sequences Ps of length #Q of multisets of integers such that &join Ps eq P and &+Ps[i] eq Qs[i].}
  return all_partition_groupings(P, Q);
end intrinsic;

function extend_binning(binning, items : limit:=Infinity(), is_valid:=func<i,b | true>, is_semivalid:=func<i,b | true>)
  // assert #binning eq #bins;
  // assert forall{i : i in [1..#bins] | #binning[i] eq bins[i]};
  assert limit ge 0;
  if limit eq 0 then
    return [];
  end if;
  // base case: no items to bin
  n := #items;
  if n eq 0 then
    if forall{i : b in binning[i], i in [1..#binning] | is_semivalid(i,b) and is_valid(i,b)} then
      return [binning];
    else
      return [];
    end if;
  end if;
  // base case: multiplicity zero
  if items[n] eq 0 then
    return extend_binning(binning, items[1..n-1] : limit:=limit, is_valid:=is_valid, is_semivalid:=is_semivalid);
  end if;
  // main case: add an item into a bin
  ret := [];
  items[n] -:= 1;
  for i in [1..#binning] do
    b := binning[i];
    // group the bin by excluding n, for each group record the highest multiplicity m, and whether m-1 was seen
    A := AssociativeArray();
    for x in MultisetToSet(b) do
      m := Multiplicity(x,n);
      x0 := Exclude(x, n^^m);
      assert n notin x0;
      if IsDefined(A, x0) then
        if m eq A[x0][1]-1 then
          A[x0][2] := true;
        elif m eq A[x0][1]+1 then
          A[x0] := <m, true>;
        elif m gt A[x0][1] then
          A[x0] := <m, false>;
        end if;
      else
        A[x0] := <m, false>;
      end if;
    end for;
    // add an item into the largest or largest but one bin
    for x0 in Keys(A) do
      for m in [A[x0][1]-(A[x0][2] select 1 else 0) .. A[x0][1]] do
        x := Include(x0, n^^m);
        assert x in b;
        x2 := Include(x0, n^^(m+1));
        b2 := Include(Exclude(b,x),x2);
        if is_semivalid(i,x2) then
          binning2 := binning;
          binning2[i] := b2;
          ret cat:= extend_binning(binning2, items : limit:=limit-#ret, is_valid:=is_valid, is_semivalid:=is_semivalid);
          if #ret ge limit then
            return ret[1..limit];
          end if;
        end if;
      end for;
    end for;
  end for;
  return ret;
end function;

function all_binnings(items, bins : limit:=Infinity(), is_valid:=func<i,b | true>, is_semivalid:=func<i,b | true>)
  return extend_binning([PowerMultiset(PowerMultiset(Integers()))| {*{**}^^bins[i]*} : i in [1..#bins]], items : limit:=limit, is_valid:=is_valid, is_semivalid:=is_semivalid);
end function;

intrinsic PGG_all_binnings(items :: [RngIntElt], bins :: [RngIntElt] : limit:=Infinity(), is_valid:=func<i,b | true>, is_semivalid:=func<i,b | true>) -> []
  {If items is a sequence of multiplicities of different items, and bins is a sequence of multiplicities of different bins, returns a sequence of possible binnings, where a binning is a sequence (of distinct bins) of multisets (of bins) of multisets (of items) of valid binnings of the items. A multiset b of items into the ith bin is valid if is_valid(i,b) is true. If b is a partial binning extendable to a valid binning, then is_semivalid(i,b) must be true; this is used to terminate branches of the algorithm early.}
  return all_binnings(items, bins : limit:=limit, is_valid:=is_valid, is_semivalid:=is_semivalid);
end intrinsic;

function mset_apply(xs, f)
  return {*f(x)^^Multiplicity(xs,x) : x in MultisetToSet(xs)*};
end function;

function sprint_with_parens(x : Level:="Default", Chars:=" \t\n")
  str := Sprintf("%O", x, Level);
  if exists{i : i in [1..#Chars] | Chars[i] in str} then
    return "(" cat str cat ")";
  else
    return str;
  end if;
end function;

function fldpad_has_isomorphism(L1, L2, K : MaximizeAPr:=true)
  // check inputs
  ok, t1 := is_extension_of(L1, K);
  assert ok;
  ok, t2 := is_extension_of(L2, K);
  assert ok;

  // case degrees unequal
  if Degree(L1, K) ne Degree(L2, K) then
    return false, _;
  elif RamificationDegree(L1, K) ne RamificationDegree(L2, K) then
    return false, _;
  end if;
  d := Degree(L1, K);

  // case d=1
  assert d eq Degree(L2, K);
  if d eq 1 then
    return true, map<L1 -> L2 | x :-> L2!K!x, y :-> L1!K!y>;
  end if;

  // case K is the direct base field
  if #t1 eq 2 and #t2 eq 2 then
    assert BaseField(L1) eq K;
    assert BaseField(L2) eq K;
    f1 := DefiningPolynomial(L1);
    f2 := DefiningPolynomial(L2);
    FUDGE := (e eq 1 select 0 else 2*e) where e:=AbsoluteRamificationDegree(L1);
    ok, root1 := has_root(ChangeRing(f1, L2));
    if not ok then
      error if has_root(ChangeRing(f2, L1)), "precision error";
      return false, _;
    end if;
    roots2 := roots(ChangeRing(f2, L1));
    error if #roots2 eq 0, "precision error";
    idxs := [i : i in [1..#roots2] | IsWeaklyEqual(trim_apr(L1.1,FUDGE), trim_apr(&+[L1| cs[i] * r^(i-1) : i in [1..#cs]],FUDGE) where cs:=fldpad_eltseq(root1)) where r:=roots2[i]];
    error if #idxs ne 1, "precision error";
    root2 := roots2[idxs[1]];
    if MaximizeAPr then
      root1 := maximize_apr(root1);
      root2 := maximize_apr(root2);
    end if;
    return true, map<L1 -> L2 |
      x :-> &+[L2| cs[i] * root1^(i-1) : i in [1..#cs]] where cs:=fldpad_eltseq(x),
      y :-> &+[L1| cs[i] * root2^(i-1) : i in [1..#cs]] where cs:=fldpad_eltseq(y)>;
  end if;

  // general case
  not_implemented("fldpad_has_isom: general towers of extensions");
end function;

procedure print_recursive(x)
  if ISA(Type(x), MonStgElt) then
    printf "%o", x;
  elif ISA(Type(x), {List,SeqEnum,Tup}) then
    for y in x do
      print_recursive(y);
    end for;
  else
    Print(x);
  end if;
end procedure;

procedure print_header_then_indent(hdr, xs)
  print_recursive(hdr);
  if #xs gt 0 then
    print "";
    IndentPush();
    for i in [1..#xs] do
      print_recursive(xs[i]);
      if i ne #xs then
        print "";
      end if;
    end for;
    IndentPop();
  end if;
end procedure;

function seq_shuffle(xs)
  return [xs[i^g] : i in [1..#xs]] where g:=Random(SymmetricGroup(#xs));
end function;

intrinsic PGG_has_random_subgroup(G :: Grp, is_valid, is_semivalid : MaxTries:=Infinity()) -> BoolElt, Grp
  {Tries to find a subgroup H of G such that `is_valid(H)` is true. For all subgroups K of all such H, we must have `is_semivalid(K)` true. On success, returns true and a random such H. Otherwise, returns false. Can fail if MaxTries is exceeded.}
  ntries := 0;
  while ntries lt MaxTries do
    ntries +:= 1;
    H := sub<G | Id(G)>;
    while true do
      if is_valid(H) then
        return true, H;
      elif H eq G then
        break;
      elif is_semivalid(H) then
        H := sub<G | H, Random(G)>;
      else
        break;
      end if;
    end while;
  end while;
  return false, _;
end intrinsic;

intrinsic PGG_has_random_subgroup_of_index(G :: Grp, n :: RngIntElt : MaxTries:=Infinity()) -> BoolElt, Grp
  {"}
  if not IsDivisibleBy(#G, n) then
    return false, _;
  end if;
  return PGG_has_random_subgroup(G, func<H | Index(G,H) eq n>, func<H | IsDivisibleBy(Index(G,H), n)> : MaxTries:=MaxTries);
end intrinsic;

intrinsic PGG_linear_divisions(length :: RngIntElt, lengths :: {* RngIntElt *} : Bound:=false) -> []
  {Given an interval of given length, find all possible divisions of it from among the given lengths. A division is a sequence of integers (sorted largest first) summing to the given length, whose elements are a subset of lengths. If Bound is given, the divisions must be lexicographically no more than it (i.e. the first element to disagree must be smaller).}
  require length ge 0: "length must be at least 0";
  require forall{len : len in MultisetToSet(lengths) | len gt 0}: "lengths must be at least 0";
  if &+lengths lt length then
    return [];
  elif length eq 0 then
    return [[]];
  end if;
  // parse the Bound
  if Bound cmpeq false then
    Bound := [length];
  else
    require forall{b : b in Bound | b gt 0}: "Bound must have positive entries";
    require &+Bound eq length: "Bound must sum to length";
    require forall{i : i in [2..#Bound] | Bound[i] le Bound[i-1]}: "Bound must be decreasing";
  end if;
  // loop over the first len
  ret := [];
  for len in Reverse(Sort(SetToSequence(MultisetToSet(lengths)))) do
    if len le length and len le Bound[1] then
      for lens in PGG_linear_divisions(length - len, lengths diff {* len *} : Bound := len eq Bound[1] select Bound[2..#Bound] else [Min(len,length-len-i) : i in [0..length-len-1 by len]]) do
        Append(~ret, [len] cat lens);
      end for;
    end if;
  end for;
  return ret;
end intrinsic;

intrinsic PGG_rectangle_divisions(width :: RngIntElt, height :: RngIntElt, areas :: {* RngIntElt *} : Bound:=false) -> []
  {Given a rectangle of given width and height, returns all possible integer divisions of the rectangle with the given areas. A division is a series of vertical cuts, through the whole rectangle, followed by a series of horizontal cuts within each resulting rectangle. It is returned in the form `[<w1,[h11,h12,...]>,<w2,[h21,...]>,...]` where the `wi` sum to width, and for each `i` the `hij` sum to height, and areas is the multiset of `wi*hij`.}
  require width ge 0: "width must be at least 0";
  require height ge 0: "height must be at least 0";
  require forall{a : a in MultisetToSet(areas) | a gt 0}: "each area must be positive";
  if &+areas lt width * height then
    return [];
  elif width eq 0 or height eq 0 then
    return [[]];
  end if;
  // parse the Bound
  if Bound cmpne false then
    max_w, max_hs := Explode(Bound);
  else
    max_w := width;
    max_hs := [height];
  end if;
  // make a w->hs lookup table
  poss_hss := AssociativeArray();
  for a in MultisetToSet(areas) do
    m := Multiplicity(areas, a);
    for h in Divisors(a) do
      w := xdiv(a, h);
      if h le height and w le width and w le max_w then
        if not IsDefined(poss_hss, w) then
          poss_hss[w] := {* h^^m *};
        else
          poss_hss[w] join:= {* h^^m *};
        end if;
      end if;
    end for;
  end for;
  // loop over possible w
  ret := [];
  for w in Reverse(Sort(SetToSequence(Keys(poss_hss)))) do
    poss_hs := poss_hss[w];
    // now find possible hs
    for hs in PGG_linear_divisions(height, poss_hs : Bound:=max_hs) do
      for rdiv in PGG_rectangle_divisions(width-w, height, areas diff {* w*h : h in hs*} : Bound:=<w, w eq max_w select max_hs else [height]>) do
        Append(~ret, [<w, hs>] cat rdiv);
      end for;
    end for;
  end for;
  return ret;
end intrinsic;

function realcomplexpairs(xs)
  rs := [Z|];
  cs := [car<Z,Z>|];
  todo := [<true,x> : x in xs];
  for i in [1..#todo] do
    if todo[i][1] then
      x := todo[i][2];
      sz,j := Min([todo[j][1] select Abs(x - Conjugate(todo[j][2])) else Infinity() : j in [i..#todo]]);
      j +:= i-1;
      assert todo[j][1];
      if j eq i then
        Append(~rs, i);
      else
        todo[j][1] := false;
        Append(~cs, <i,j>);
      end if;
    end if;
  end for;
  assert #rs + 2*#cs eq #xs;
  assert ({r:r in rs} join {r:r in c, c in cs}) eq {1..#xs};
  return rs, cs;
end function;

intrinsic PGG_realcomplexpairs(xs :: [FldComElt]) -> [], []
  {Returns the sequence of `i` such that `xs[i]` are real, and the sequence of pairs `<i,j>` such that `xs[i]` and `xs[j]` are a complex conjugate pair. Assumes that xs are all the roots of a squarefree real polynomial in some order.}
  return realcomplexpairs(xs);
end intrinsic;

intrinsic PGG_CloseVector(L :: Lat, vec : Close:=0) -> LatElt
  {A vector of L close to vec.}
  V := AmbientSpace(L);
  vec := V ! vec;
  n := Dimension(V);
  m := #Basis(L);
  M := Matrix([Eltseq(V!b) cat [0] : b in Basis(L)] cat [Eltseq(vec) cat [0]]);
  assert Nrows(M) eq m+1;
  assert Ncols(M) eq n+1;
  weight := BaseRing(L) ! 1.0;
  bestnorm := false;
  bestvec := false;
  normbound := Close^2;
  while true do
    M[m+1,n+1] := weight;
    M2, T := LLL(M);
    goodrows := [i : i in [1..m+1] | T[i,m+1] in [1,-1]];
    for i in goodrows do
      coeffs := [-T[i,j]*T[i,m+1] : j in [1..m]];
      lvec := &+[coeffs[j]*Basis(L)[j] : j in [1..m]];
      norm := Norm(vec - V!lvec);
      if norm le normbound then
        return lvec;
      elif bestnorm cmpeq false or norm lt bestnorm then
        bestnorm := norm;
        bestvec := lvec;
      end if;
    end for;
    if #goodrows gt 0 then
      weight *:= 0.5;
    elif bestnorm cmpeq false then
      weight *:= 2;
    else
      return bestvec;
    end if;
  end while;
end intrinsic;