GloMod.mag

import "Utils.mag": Q, not_implemented, largest_coefficient, reduce_coefficients, change_apr, maximize_apr, has_root, factorization, ramification_tower, random_primitive_element, fldpad_has_isomorphism, seq_shuffle, roots, fldpad_eltseq, is_eisenstein, is_inertial, precision_error, Z;

declare type PGGGloMod;

declare type PGGGloMod_Fld: PGGGloMod;
declare attributes PGGGloMod_Fld: local_field, global_field, embeddings, approximations, primes;

declare type PGGGloMod_UPol: PGGGloMod;
declare attributes PGGGloMod_UPol: local_pol, global_pol;

declare type PGGGloMod_Ext: PGGGloMod;
declare attributes PGGGloMod_Ext: base_model;

declare type PGGGloMod_FldExt: PGGGloMod_Fld, PGGGloMod_Ext;

declare type PGGGloMod_Rational: PGGGloMod_Fld;

declare type PGGGloMod_Symmetric: PGGGloMod_FldExt, PGGGloMod_UPol;

declare type PGGGloMod_Factors: PGGGloMod_Ext;
declare attributes PGGGloMod_Factors: factors;

declare type PGGGloMod_Tower: PGGGloMod_FldExt;
declare attributes PGGGloMod_Tower: tower, global_pol, global_pol_root;

declare type PGGGloMod_RootOfUnity: PGGGloMod_FldExt, PGGGloMod_UPol;
declare attributes PGGGloMod_RootOfUnity: primitive_root, aut_powers, sum_powers, zeta_pol;

declare type PGGGloMod_RootOfUniformizer: PGGGloMod_FldExt, PGGGloMod_UPol;
declare attributes PGGGloMod_RootOfUniformizer: degree, global_uniformizer, global_zeta_pol;

declare type PGGGloMod_PthRoots: PGGGloMod_FldExt, PGGGloMod_UPol;
declare attributes PGGGloMod_PthRoots: p, n, local_gens, global_gens, global_root_coeffs, global_zetap;

declare type PGGGloMod_UPol_Cheat: PGGGloMod_UPol;
declare attributes PGGGloMod_UPol_Cheat: model, overgroup_embedding;

// populates the primes attribute
procedure make_primes(~m)
  m`primes := [ideal<Integers(m`global_field) | Prime(m`local_field), a(UniformizingElement(m`local_field), 2)> : a in m`approximations];
end procedure;

// the number field defined by pol
// works around the bug that if you just do NumberField(pol) then MaximalOrder sometimes doesn't work, so we make the maximal order up front
function number_field(pol : base:=BaseRing(pol))
  return NumberField(ext<(Type(base) in [FldRat,RngInt]) select Integers() else EquationOrder(base) | pol>);
end function;

intrinsic Approximation(m :: PGGGloMod_Fld, x :: PGGFldElt, pr :: RngIntElt) -> FldElt
  {Takes a global approximation to x across all embeddings.}
  if #m`approximations eq 1 then
    return m`approximations[1](x, pr);
  else
    return m`global_field ! CRT([Integers(m`global_field)| a(x,pr) : a in m`approximations], [p^pr : p in m`primes]);
  end if;
end intrinsic;

intrinsic Approximation(m :: PGGGloMod_Fld, x :: PGGFldElt, pr :: Infty) -> FldElt
  {Takes a global approximation to x across all embeddings.}
  require pr eq Infinity(): "expecting positive infinity";
  require AbsolutePrecision(x) eq pr and IsWeaklyZero(x): "expecting x to be precisely zero";
  return m`global_field ! 0;
end intrinsic;

intrinsic GlobalModel(F :: PGGFldWrap) -> PGGGloMod_Fld
  {A global model for F.}
  if IsPrimeField(F) then
    m := New(PGGGloMod_Rational);
    m`local_field := F;
    m`global_field := Q;
    m`embeddings := [map<Q -> F | x :-> x>];
    m`approximations := [func<y, pr | RationalApproximation(F!y, pr)>];
    m`primes := [ideal<Z | Prime(F)>];
    return m;
  else
    bm := GlobalModel(BaseField(F));
    alg := PGGAlg_ResEval_Global_Model_Symmetric_Make();
    m := GlobalModel(alg, DefiningPolynomial(F), bm : TopField:=F, Extendible);
    return m;
  end if;
end intrinsic;

intrinsic GlobalModel(F :: PGGFldGrp) -> PGGGloMod_Fld
  {"}
  m := New(PGGGloMod_Fld);
  m`local_field := F;
  return m;
end intrinsic;

intrinsic GlobalModel(alg :: PGGAlg_ResEval_Global_Model, f :: PGGPol, bm :: PGGGloMod_Fld : TopField:=false) -> PGGGloMod, PGGHomGrpPerm
  {Global model with the given base model, and the corresponding hom W -> W2 where Gal(f)<W and Gal(model)<W2.}
  error "not implemented: GlobalModel:", Type(alg), Type(f), Type(bm);
end intrinsic;

intrinsic GlobalModel(alg :: PGGAlg_ResEval_Global_Model_Cheat, f :: PGGPol, bm :: PGGGloMod_Fld : TopField:=false) -> PGGGloMod, PGGHomGrpPerm
  {"}
  m := New(PGGGloMod_UPol_Cheat);
  m`local_pol := f;
  m`model, m`overgroup_embedding := GlobalModel(alg`next, f, bm : TopField:=TopField);
  return m, m`overgroup_embedding;
end intrinsic;

intrinsic GlobalModel(alg :: PGGAlg_ResEval_Global_Model_Symmetric, f :: PGGPolWrap, bm :: PGGGloMod_Fld : TopField:=false, Extendible:=false) -> PGGGloMod, PGGHomGrpPerm
  {"}
  locW := PGGConj_Symmetric_Make(f : GaloisGroupAlg:=alg`galois_group_alg, TopField:=TopField);
  gloW := PGGGrpPerm_Symmetric_Make(Degree(f));
  Whom := IdentityHomomorphism(locW, gloW);
  m := New(PGGGloMod_Symmetric);
  m`base_model := bm;
  PGG_GlobalTimer_Push("global polynomial");
  if assigned locW`top_field then
    locfld := BaseRing(f);
    locext := locW`top_field;
    if locext eq locfld then
      locpol := Parent(f) ! [0,1];
      glopol := Polynomial([bm`global_field|0,1]);
      gloext := bm`global_field;
      m`local_pol := locpol;
      m`local_field := locext;
      m`global_pol := glopol;
      m`global_field := bm`global_field;
      m`embeddings := bm`embeddings;
      m`approximations := bm`approximations;
      m`primes := bm`primes;
    else
      if BaseField(locext) eq locfld then
        locpol := DefiningPolynomial(locext);
        assert BaseRing(locpol) eq locfld;
        if InertiaDegree(locext, locfld) eq 1 then
          assert IsEisenstein(locpol);
          ty := "E";
        elif RamificationDegree(locext, locfld) eq 1 then
          assert IsInertial(locpol);
          ty := "I";
        else
          assert false;
        end if;
      else
        if Extendible then
          not_implemented("extendible symmetric global models for ramified but not totally ramified extensions");
        end if;
        ty := "M";
        locpol := DefiningPolynomial(locext, locfld);
        assert BaseRing(locpol) eq locfld;
        assert IsIrreducible(locpol);
      end if;
      nembs := #bm`embeddings;
      // minimize the coefficients of glopol
      K := BaseRing(locpol);
      p := Prime(K);
      ramdeg := AbsoluteRamificationDegree(K);
      max_pr := Ceiling(Max([AbsolutePrecision(c) : c in Coefficients(locpol)]) / ramdeg);
      // max_pr := Ceiling(Log(Prime(K), Max([largest_coefficient(c) : c in Coefficients(glopol0)])));
      pr := 0;
      while true do
        if pr lt 4 then
          pr +:= 1;
        elif pr lt max_pr then
          pr := Min(max_pr, 2*pr);
        else
          precision_error();
        end if;
        vprint PGG_GaloisGroup, 1: "pr =", pr;
        PGG_GlobalTimer_Push("approximation");
        glopol := Polynomial([reduce_coefficients(Approximation(bm, c, pr*ramdeg), p^pr) : c in Coefficients(locpol)]);
        // glopol := Polynomial([reduce_coefficients(c, Prime(K)^pr) : c in Coefficients(glopol0)]);
        vprint PGG_GaloisGroup, 2: "glopol =", glopol;
        PGG_GlobalTimer_Swap("embed");
        locpol2s := [Polynomial([MaximizeAbsolutePrecision(c@e) : c in Coefficients(glopol)]) : e in bm`embeddings];
        if exists{locpol2 : locpol2 in locpol2s | not case<ty | "E": IsEisenstein(locpol2), "I": IsInertial(locpol2), "M": IsIrreducible(locpol2), default: not_implemented()>} then
          PGG_GlobalTimer_Pop();
          continue;
        end if;
        locext2s := [case<ty | "E": ext<K | pol>, "I": ext<K | pol>, "M": Extension(pol), default: not_implemented()> : pol in locpol2s];
        PGG_GlobalTimer_Swap("isomorphic");
        isoms := [**];
        ok := true;
        for locext2 in locext2s do
          ok, isom := HasIsomorphism(locext, locext2, K);
          if ok then
            Append(~isoms, isom);
          else
            break;
          end if;
        end for;
        PGG_GlobalTimer_Pop();
        if ok then
          break;
        end if;
      end while;
      // make the embeddings
      gloext := number_field(glopol);
      m`local_pol := locpol;
      m`local_field := locext;
      m`global_pol := glopol;
      m`global_field := gloext;
      if Extendible then
        assert ty in ["E","I"];
        m`embeddings := [
          map<gloext -> locext | x :-> (locext2![e(c) : c in Eltseq(x)]) @@ isom>
            where locext2:=locext2s[i]
            where isom:=isoms[i]
            where e:=bm`embeddings[i]
          : i in [1..nembs]
        ];
        m`approximations := [function (y, pr)
          y := locext ! y;
          pr0 := Ceiling(pr / ramdeg);
          return gloext ! [a(c, pr0) : c in Eltseq(y @ isom)];
        end function where ramdeg:=RamificationDegree(locext) where isom:=isoms[i] where a:=bm`approximations[i] : i in [1..nembs]];
        PGG_GlobalTimer_Push("primes");
        make_primes(~m);
        PGG_GlobalTimer_Pop();
      end if;
      // TODO: can we construct the primes directly, instead of by factoring?
      //       does it suffice to take the ideal generated by p and by a uniformizer?
      // m`primes := [Factorization(Integers(m`global_field) !! p)[1][1] : p in bm`primes];
    end if;
  else
    vprint PGG_GaloisGroup: "WARNING: not minimizing global coefficients (if the complex precision explodes, or complex root finding fails, this might be why) (it might even hang right here)";
    vprint PGG_GaloisGroup: "WARNING: assuming coefficients to current precision are sufficient";
    m`local_pol := f;
    m`global_pol := Polynomial([Approximation(bm, c, AbsolutePrecision(c)) : c in Coefficients(f)]);
  end if;
  PGG_GlobalTimer_Pop();
  return m, Whom;
end intrinsic;

intrinsic GlobalModel(alg :: PGGAlg_ResEval_Global_Model_Symmetric, f :: PGGPolGrp, bm :: PGGGloMod_Fld : TopField:=false) -> PGGGloMod, PGGHomGrpPerm
  {"}
  locW := PGGConj_Symmetric_Make(f : GaloisGroupAlg:=alg`galois_group_alg, TopField:=TopField);
  gloW := PGGGrpPerm_Symmetric_Make(Degree(f));
  Whom := IdentityHomomorphism(locW, gloW);
  m := New(PGGGloMod_Symmetric);
  m`local_pol := f;
  return m, Whom;
end intrinsic;

intrinsic GlobalModel(alg :: PGGAlg_ResEval_Global_Model_Factors, f :: PGGPol, bm :: PGGGloMod_Fld : TopField:=false) -> PGGGloMod, PGGHomGrpPerm
  {"}
  if TopField cmpeq false then
    PGG_GlobalTimer_Push("factorize");
    facs, certs := Factorization(f : Extensions);    
    PGG_GlobalTimer_Pop();
  else
    facs := [f];
    certs := [rec<recformat<Extension> | Extension:=TopField>];
  end if;
  PGG_GlobalTimer_Push("recurse");
  data := [* <m, h> where m,h:=GlobalModel(alg`next, facs[i], bm : TopField:=certs[i]`Extension) : i in [1..#facs] *];
  PGG_GlobalTimer_Pop();
  m := New(PGGGloMod_Factors);
  m`factors := [x[1] : x in data];
  m`base_model := bm;
  dom := PGGConj_Factors_Make(f, <Domain(x[2]) : x in data>);
  hom := PGGHomGrpPerm_DirProd_Make(<x[2] : x in data> : Domain:=dom);
  return m, hom;
end intrinsic;

intrinsic Tower(alg :: PGGAlg_ResEval_Global_Model_Tower, L :: PGGFld, K :: PGGFld) -> []
  {The tower from K up to L.}
  not_implemented("Tower:", Type(alg), Type(L));
end intrinsic;

intrinsic Tower(alg :: PGGAlg_ResEval_Global_Model_RamTower, L :: PGGFld, K :: PGGFld) -> []
  {"}
  return RamificationTower(L, K);
end intrinsic;

intrinsic GlobalModel(alg :: PGGAlg_ResEval_Global_Model_Tower, f :: PGGPol, bm :: PGGGloMod_Fld : TopField:=false, Extendible:=false) -> PGGGloMod, PGGHomGrpPerm
  {"}
  if TopField cmpeq false then
    PGG_GlobalTimer_Push("factorize");
    facs, certs := Factorization(f : Extensions);
    PGG_GlobalTimer_Pop();
    assert #facs eq 1;
    L := certs[1]`Extension;
  else
    L := TopField;
  end if;
  PGG_GlobalTimer_Push("tower");
  t := Tower(alg, L, BaseRing(f));
  PGG_GlobalTimer_Swap("recurse");
  m := New(PGGGloMod_Tower);
  m`tower := [**];
  m`base_model := bm;
  m2 := bm;
  hs := [**];
  for i in [2..#t] do
    m2, h2 := GlobalModel(alg`next, DefiningPolynomial(t[i], t[i-1]), m2 : TopField:=t[i], Extendible:=Extendible or i lt #t);
    Append(~m`tower, m2);
    Append(~hs, h2);
  end for;
  if assigned m2`global_field then
    PGG_GlobalTimer_Swap("primitive element");
    m`global_field := m2`global_field;
    gloroot, glopol := random_primitive_element(m2`global_field, bm`global_field);
    m`global_pol := glopol;
    m`global_pol_root := gloroot;
  end if;
  PGG_GlobalTimer_Pop();
  if assigned m2`local_field then
    m`local_field := m2`local_field;
  end if;
  if assigned m2`embeddings then
    m`embeddings := m2`embeddings;
  end if;
  if assigned m2`approximations then
    m`approximations := m2`approximations;
  end if;
  if assigned m2`primes then
    m`primes := m2`primes;
  end if;
  dom := PGGConj_Tower_Make(f, <Domain(h) : h in hs>);
  hom := PGGHomGrpPerm_WrProd_Make(<h : h in hs> : Domain:=dom);
  return m, hom;
end intrinsic;

intrinsic GlobalModel(alg :: PGGAlg_ResEval_Global_Model_Select, f :: PGGPol, bm :: PGGGloMod_Fld : TopField:=false, Extendible:=false) -> PGGGloMod, PGGHomGrpPerm
  {"}
  PGG_GlobalTimer_Push("predicates");
  K := BaseRing(f);
  L := TopField;
  for i in [1..#alg`models] do
    if (i gt #alg`predicates)
    or EvaluateLazy(
        alg`predicates[i],
        func<n, v | case<n
          | "p": Prime(BaseRing(f))
          , "faccerts": certs where _,certs := Factorization(f : Extensions)
          , "irr": L cmpne false or #v("faccerts") eq 1
          , "deg": Degree(f)
          , "unram": L cmpne false select RamificationDegree(L,K) eq 1 else forall{c : c in v("faccerts") | c`E eq 1}
          , "tame": (L cmpne false select not IsDivisibleBy(RamificationDegree(L,K), p) else forall{c : c in v("faccerts") | not IsDivisibleBy(c`E, p)}) where p:=Prime(BaseRing(f))
          , "ram": not v("unram")
          , "wild": not v("tame")
          , "totram": L cmpne false select InertiaDegree(L,K) eq 1 else forall{c : c in v("faccerts") | c`F eq 1}
          , "totwild": (L cmpne false select InertiaDegree(L,K) eq 1 and IsPowerOf(RamificationDegree(L,K),p) else forall{c : c in v("faccerts") | c`F eq 1 and IsPowerOf(c`E, p)}) where p:=Prime(BaseRing(f))
          , default: not_implemented()
          >
        > : Recursive)
    then
      PGG_GlobalTimer_Swap("recurse");
      m, h := GlobalModel(alg`models[i], f, bm : TopField:=TopField, Extendible:=Extendible);
      PGG_GlobalTimer_Pop();
      return m, h;
    end if;
  end for;
  PGG_GlobalTimer_Pop();
  error "cases were not exhaustive";
end intrinsic;

intrinsic GlobalModel(alg :: PGGAlg_ResEval_Global_Model_RootOfUnity, f :: PGGPolWrap, bm :: PGGGloMod_Fld : TopField:=false, Extendible:=false) -> PGGGloMod, PGGHomGrpPerm
  {"}
  // check the inputs
  K := BaseRing(f);
  L := TopField;
  if L cmpeq false then
    not_implemented("expecting f to be irreducible");
  end if;
  if RamificationDegree(L, K) ne 1 then
    not_implemented("ramified extensions");
  end if;

  // choose n such that we are adjoining an nth root of unity
  d := Degree(L, K);
  assert InertiaDegree(L, K) eq d;
  p := Prime(K);
  q := p^AbsoluteInertiaDegree(K);
  n0 := q^d-1;
  if alg`minimize then
    n := Min([n : n in Divisors(n0) | forall{i : i in [1..d-1] | not IsDivisibleBy(q^i-1, n)}]);
  else
    n := n0;
  end if;
  // choose a factor of the nth cyclotomic polynomial over the global field
  Cn := CyclotomicPolynomial(n);
  gK := Domain(bm`embeddings[1]);
  assert forall{e : e in bm`embeddings | Domain(e) eq gK and Codomain(e) eq K};
  C := Factorization(ChangeRing(Cn, gK))[1][1];
  dN := Degree(C);
  ok, sN := IsDivisibleBy(dN, d);
  assert ok;
  // extend by this factor, so z is a primitive nth root of unity
  gN<zeta> := number_field(C : base:=gK);//ext<gK | C>;
  assert Evaluate(C, zeta) eq 0;
  assert Degree(gN, gK) eq dN;
  // find the i such that zeta->zeta^i are automorphisms
  auts := {@ i : i in [1..n-1] | Evaluate(C,zeta^i) eq 0 @};
  assert #auts eq dN;
  assert forall{i : i in [0..d-1] | ((q^i) mod n) in auts};
  // put these into orbits under z->z^q, corresponding to the local Galois group (each orbit corresponds to an embedding of gM into K)
  auts_todo := {i : i in auts};
  orbits := [];
  for j in [1..sN] do
    orbit := [];
    i := Min(auts_todo);
    for k in [1..d] do
      assert i in auts_todo;
      Exclude(~auts_todo, i);
      Append(~orbit, i);
      i := (q*i) mod n;
    end for;
    Append(~orbits, orbit);
  end for;
  // make the global field actually forming the model, and the intermediate field embedding into the local base field
  if sN eq 1 then
    dL := dN;
    sL := 1;
    gL<z> := gN;
    gM := gK;
    sum_powers := [1];
    Lorbits := orbits;
  elif alg`complement then
    M, Mmap := MultiplicativeGroup(Integers(n));
    G := sub<M | [i @@ Mmap : i in auts]>;
    assert #G eq #auts;
    H := sub<G | q @@ Mmap>;
    assert #H eq d;
    ok, Hc := HasComplement(G, H);
    if not ok then
      not_implemented("RootOfUnity: Complement: imperfect complement");
    end if;
    dL := Index(G, Hc);
    ok, sL := IsDivisibleBy(dL, d);
    assert ok;
    sum_powers := [Z| g @ Mmap : g in Hc];
    gL<z> := sub<gN | &+[zeta^i : i in sum_powers]>;
    assert Degree(gL, gK) eq dL;
    if sL eq 1 then
      gM := gK;
    else
      not_implemented("RootOfUnity: Complement: imperfect complements");
    end if;
    Lorbits := [o : o in orbits | not exists{o2 : o2 in Self() | exists{i : i in o2 | (o[1]@@Mmap)-(i@@Mmap) in Hc}}];
  else
    gL<z> := gN;
    gM<w> := sub<gL | &+[z^(q^i) : i in [0..d-1]]>;
    dL := dN;
    sL := sN;
    assert Degree(gM, gK) eq sL;
    sum_powers := [1];
    Lorbits := orbits;
  end if;
  vprint PGG_GaloisGroup, 2: "orbits =", Lorbits;
  // choose a factor of C over M
  zpol := DefiningPolynomial(gL);
  assert Evaluate(zpol, z) eq 0;
  zpolM := Factorization(ChangeRing(zpol, gM))[1][1];
  assert Degree(zpolM) eq d;
  // choose primes of L above the primes of K
  p0s := [Factorization(Integers(gL) !! p)[1][1] : p in bm`primes];
  // embed gM into K
  if sL eq 1 then
    Membs := bm`embeddings;
  else
    wpol := DefiningPolynomial(gM);
    assert Evaluate(wpol, w) eq 0;
    Membs := [];
    for e in bm`embeddings do
      pol := Polynomial([K| c@e : c in Coefficients(wpol)]);
      roots := Roots(pol);
      assert #roots ne 0;
      root := MaximizeAbsolutePrecision(roots[1]);
      Append(~Membs, map<gM -> K
        | x :-> &+[K| (cs[i] @ e) * root^(i-1) : i in [1..#cs]] where cs:=Eltseq(x)
      >);
    end for;
  end if;
  // make automorphisms of gL corresponding to the orbits
  gLauts := [hom<gL -> gL | &+[zeta^((o[1]*j) mod n) : j in sum_powers]> : o in Lorbits];
  // make primes of gL corresponding to embeddings
  ps := [ideal<Integers(gL) | [gen @ aut : gen in Generators(p0)]> : aut in gLauts, p0 in p0s];
  // make representative embeddings of gL to L
  Lembs0 := [];
  Lapprs0 := [];
  for i in [1..#bm`embeddings] do
    e := bm`embeddings[i];
    a := bm`approximations[i];
    e0 := Membs[i];
    loczpol := Polynomial([K| MaximizeAbsolutePrecision(c @ e0) : c in Coefficients(zpolM)]);
    L2 := ext<K | loczpol>;
    ok, isom := HasIsomorphism(L, L2, K);
    assert ok;
    Append(~Lembs0, map<gL -> L
      | x :-> (&+[L2| Generator(L2)^(i-1) * (L2!(cs[i]@e)) : i in [1..#cs]])@@isom where cs:=Eltseq(x)
    >);
    Append(~Lapprs0, function (y, pr)
      y := L!y;
      cs := Eltseq(y @ isom);
      return &+[gL| (gL.1)^(i-1) * a(cs[i], pr) : i in [1..#cs]];
    end function);
  end for;
  // make all embeddings
  Lembs := [map<gL -> L | x :-> x @ aut @ e> : aut in gLauts, e in Lembs0];
  Lapprs := [func<y, pr | a(y,pr) @@ aut> : aut in gLauts, a in Lapprs0];

  // create the model
  m := New(PGGGloMod_RootOfUnity);
  m`base_model := bm;
  m`embeddings := Lembs;
  m`approximations := Lapprs;
  m`primes := ps;
  m`local_field := L;
  m`global_field := gL;
  // m`local_pol := C2;
  m`global_pol := zpol;
  m`primitive_root := n;
  m`zeta_pol := C;
  m`aut_powers := [i : i in o, o in Lorbits];
  m`sum_powers := sum_powers;

  // create the corresponding groups
  // W is cyclic of order d
  W := PGGConj_Unramified_Make(f : TopField:=L);
  // W2wr = S_d wr S_s
  W2wr := PGGGrpPerm_WrProd_Make(<PGGGrpPerm_Symmetric_Make(sL), PGGGrpPerm_Symmetric_Make(d)>);
  // W2 is the subgroup of W2wr corresponding to the global Galois group (identifying (i->z^i) with the index of i in auts2)
  W2 := PGGGrpPerm_RawSubgroup_Make(W2wr, sub<Group(W2wr) | [[dN eq dL select Index(m`aut_powers, (i*j) mod n) else [k : k in [1..#m`aut_powers] | (m`aut_powers[k] @@ Mmap) - ((i*j) @@ Mmap) in Hc][1] : j in m`aut_powers] : i in m`aut_powers]>);
  // the diagonal embedding of W into W2
  Whom := PGGHomGrpPerm_WrDiagonal_Make(W, W2);
  assert #Group(W2) eq dL;
  assert IsAbelian(Group(W2));
  vprint PGG_GaloisGroup: "W2 =", GroupName(Group(W2));
  return m, Whom;
end intrinsic;

intrinsic GlobalModel(alg :: PGGAlg_ResEval_Global_Model_RootOfUniformizer, f :: PGGPolWrap, bm :: PGGGloMod_Fld : TopField:=false, Extendible:=false) -> PGGGloMod, PGGHomGrpPerm
  {"}

  // check the inputs
  if TopField cmpeq false then
    not_implemented("f must be known irreducible");
  end if;
  L := TopField;
  K := BaseRing(f);
  d := Degree(f);
  p := Prime(K);
  assert Degree(L, K) eq d;
  if (InertiaDegree(L, K) ne 1) or IsDivisibleBy(d, p) then
    not_implemented("expecting a totally tame extension");
  end if;

  // get the global base field
  gK := Domain(bm`embeddings[1]);
  assert forall{e : e in bm`embeddings | Domain(e) eq gK};
  assert forall{e : e in bm`embeddings | Codomain(e) eq K};

  // find Aut(gK(zeta_d)/gK) in terms of i s.t. zeta -> zeta^i
  C := Factorization(ChangeRing(CyclotomicPolynomial(d), gK))[1][1];
  if Degree(C) eq 1 then
    gM := gK;
    z := Roots(C)[1][1];
  else
    gM<z> := number_field(C : base:=gK);//ext<gK | C>;
  end if;
  zauts := {@ i : i in [1..d-1] | Evaluate(C, z^i) eq 0 @};
  // zauts := {@ i : i in [1..d-1] | GCD(i,d) eq 1 @};
  vprint PGG_GaloisGroup, 1: "zauts =", zauts;
  assert 1 in zauts;
  assert IsDivisibleBy(EulerPhi(d), #zauts);

  // the global field
  piK := -(Generator(L)^d);
  while Parent(piK) ne K do
    piK := Eltseq(piK)[1];
  end while;
  // piK := Coefficient(DefiningPolynomial(L, K), 0);
  assert ValuationEq(piK, 1);
  gpiK0 := Approximation(bm, piK, 2);
  gpiK := reduce_coefficients(gpiK0, p^2);
  gL := number_field(Polynomial([case<i | 0:gpiK, d:1, default: 0> : i in [0..d]]) : base:=gK);//ext<gK | Polynomial([case<i | 0:gpiK, d:1, default: 0> : i in [0..d]])>;

  // extend each embedding
  Lembs := [];
  Lapprs := [];
  for i in [1..#bm`embeddings] do
    e := bm`embeddings[i];
    a := bm`approximations[i];
    L2 := ext<K | Polynomial([K| case<i | 0:e(gpiK), d:1, default: 0> : i in [0..d]])>;
    ok, isom := HasIsomorphism(L, L2, K);
    assert ok;
    Append(~Lembs, map<gL -> L
      | x :-> (&+[L2| (L2!(cs[i]@e)) * Generator(L2)^(i-1) : i in [1..#cs]])@@isom where cs:=Eltseq(x)
    >);
    Append(~Lapprs, function (y, pr)
      y := L ! y;
      pr0 := Ceiling(pr / d);
      cs := Eltseq(y @ isom);
      return &+[gL| a(cs[i], pr0) * (gL.1)^(i-1) : i in [1..#cs]];
    end function);
  end for;

  // make the model
  m := New(PGGGloMod_RootOfUniformizer);
  m`base_model := bm;
  m`embeddings := Lembs;
  m`approximations := Lapprs;
  m`local_field := L;
  m`global_field := gL;
  m`global_pol := DefiningPolynomial(gL);
  m`global_uniformizer := -gpiK;
  m`global_zeta_pol := C;
  m`degree := d;
  make_primes(~m);

  // the Galois groups
  Sd := PGGGrpPerm_Symmetric_Make(d);
  gW := PGGGrpPerm_RawSubgroup_Make(Sd, sub<Group(Sd) | [[(((k-1)*j+i) mod d)+1 : k in [1..d]] : j in zauts, i in [0..d-1]]>);
  assert #Group(gW) eq d*#zauts;
  W := PGGConj_TotRam_Make(f : TopField:=L);
  Whom := PGGHomGrpPerm_Id_Make(W, gW);
  // error "DEBUG";

  return m, Whom;
end intrinsic;

intrinsic ComplexEmbeddings(m :: PGGGloMod_Rational, embs :: [Map]) -> []
  {The complex embeddings of m extending embs.}
  return embs;
end intrinsic;

intrinsic ComplexEmbeddings(m :: PGGGloMod_Symmetric, embs :: [Map]) -> []
  {"}
  rss := ComplexRoots(m, embs);
  return [map<m`global_field -> C | x :-> &+[(cs[j] @ e) * r^(j-1) : j in [1..#cs]] where cs:=Eltseq(x)> where C:=Codomain(e) where e:=embs[i] : r in rss[i], i in [1..#embs]];
end intrinsic;

intrinsic ComplexEmbeddings(m :: PGGGloMod_RootOfUnity, embs :: [Map]) -> []
  {"}
  rss := ComplexRoots(m, embs);
  return [map<m`global_field -> C | x :-> &+[(cs[j] @ e) * r^(j-1) : j in [1..#cs]] where cs:=Eltseq(x)> where C:=Codomain(e) where e:=embs[i] : r in rss[i], i in [1..#embs]];
end intrinsic;

intrinsic ComplexEmbeddings(m :: PGGGloMod_RootOfUniformizer, embs :: [Map]) -> []
  {"}
  rss := ComplexRoots(m, embs);
  return [map<m`global_field -> C | x :-> &+[(cs[j] @ e) * r^(j-1) : j in [1..#cs]] where cs:=Eltseq(x)> where C:=Codomain(e) where e:=embs[i] : r in rss[i], i in [1..#embs]];
end intrinsic;

intrinsic ComplexEmbeddings(m :: PGGGloMod_PthRoots, embs :: [Map]) -> []
  {"}
  p := m`p;
  n := m`n;
  ret := [];
  for e in embs do
    C := Codomain(e);
    zetap := m`global_zetap @ e;
    for J in CartesianPower([0..p-1], n) do
      e2 := e;
      for i in [1..n] do
        xroot := Root(m`global_gens[i] @ e, p) * zetap^J[i];
        e2 := func<y | &+[e0(cs[i]) * xroots[i] : i in [1..p]] where cs:=Eltseq(y)> where e0:=e2 where xroots:=[xroot^i : i in [0..p-1]];
      end for;
      Append(~ret, map<m`global_field -> C | x :-> e2(x)>);      
    end for;
  end for;
  return ret;
end intrinsic;

intrinsic ComplexEmbeddings(m :: PGGGloMod_Tower, embs :: [Map]) -> []
  {"}
  es := embs;
  for m2 in m`tower do
    es := ComplexEmbeddings(m2, es);
  end for;
  return es;
end intrinsic;

intrinsic ComplexRoots(m :: PGGGloMod_Symmetric, embs :: [Map]) -> []
  {The complex roots of the polynomial model.}
  return [[r[1] : r in Roots(Polynomial([c@e : c in Coefficients(m`global_pol)]))] : e in embs];
end intrinsic;

intrinsic ComplexRoots(m :: PGGGloMod_RootOfUnity, embs :: [Map]) -> []
  {"}
  n := m`primitive_root;
  zetas := [r where ok, r := HasRoot(Polynomial([c@e : c in Coefficients(m`zeta_pol)])) : e in embs];
  assert forall{zeta : zeta in zetas | Abs(Abs(zeta)-1) lt 1e-10};
  assert forall{zeta : zeta in zetas | Abs(zeta^n-1) lt 1e-10};
  return [[&+[zeta^((i*j) mod n) : j in m`sum_powers] : i in m`aut_powers] : zeta in zetas];
end intrinsic;

intrinsic ComplexRoots(m :: PGGGloMod_RootOfUniformizer, embs :: [Map]) -> []
  {"}
  // rs := [Roots(Polynomial([c@e : c in Coefficients(m`global_pol)]))[1][1] : e in embs];
  rs := [(m`global_uniformizer @ e)^(1/m`degree) : e in embs];
  if m`degree eq 1 then
    return [[r] : r in rs];
  else
    // zs := [((C!-1)^(1/m`degree))^2 where C:=Codomain(e) : e in embs];
    // THE ABOVE LINE WILL CAUSE ERRORS!!
    // We don't have free choice over the root of unity, it needs to be consistent between different models, i.e. a root of a fixed polynomial over the base field
    zs := [Roots(Polynomial([c@e : c in Coefficients(m`global_zeta_pol)]))[1][1] : e in embs];
    return [[j eq 0 select r else r * z^j : j in [0..m`degree-1]] where z:=zs[i] where r:=rs[i] : i in [1..#embs]];
  end if;
end intrinsic;

intrinsic ComplexRoots(m :: PGGGloMod_PthRoots, embs :: [Map]) -> []
  {"}
  p := m`p;
  n := m`n;
  return [[&+[zetapows[J[i]+1] * roots[i] : i in [1..n]] : J in CartesianPower([0..p-1], n)]
    where zetapows := [zetap^i : i in [0..p-1]]
    where zetap := m`global_zetap @ e
    where roots := [Root(m`global_gens[i] @ e, p) * m`global_root_coeffs[i] : i in [1..n]]
    : e in embs
  ];
end intrinsic;

intrinsic ComplexRoots(m :: PGGGloMod_Factors, embs :: [Map]) -> []
  {"}
  rsss := [ComplexRoots(fac, embs) : fac in m`factors];
  return [&cat[rss[i] : rss in rsss] : i in [1..#embs]];
end intrinsic;

intrinsic ComplexRoots(m :: PGGGloMod_Tower, embs :: [Map]) -> []
  {"}
  es := ComplexEmbeddings(m, embs);
  d := #m`tower eq 0 select 1 else Degree(m`tower[#m`tower]`global_field, BaseRing(m`tower[1]`global_pol));
  assert #es eq #embs * d;
  // todo: also find complex roots of the global_pol, which should be more accurate but unordered, and pair them up
  return [[e(m`global_pol_root) : e in es[1+d*(i-1)..d*i]] : i in [1..#embs]];
end intrinsic;

intrinsic AllComplexEmbeddings(m :: PGGGloMod_Rational, C :: FldCom) -> []
  {A complex embedding into C.}
  return [map<m`global_field -> C | x :-> C!x>];
end intrinsic;

intrinsic AllComplexEmbeddings(m :: PGGGloMod_Ext, C :: FldCom) -> []
  {"}
  return ComplexEmbeddings(m, AllComplexEmbeddings(m`base_model, C));
end intrinsic;