Field.mag

import "Utils.mag": is_extension_of, xdiv, ramification_tower, precision_error, trim_apr, not_implemented, Q, Z;

// Base type for any representation of a p-adic field
declare type PGGFld[PGGFldElt];
declare attributes PGGFld: polynomial_ring;
declare attributes PGGFldElt: parent;

// Base type for any p-adic field wrapping another representation
declare type PGGFldWrap[PGGFldWrapElt]: PGGFld;
declare attributes PGGFldWrap: actual, base_field, defining_polynomial, residue_class_field, generator;
declare type PGGFldWrapElt: PGGFldElt;
declare attributes PGGFldWrapElt: actual, eltseq;

// Wraps Magma's standard p-adics (FldPad)
declare type PGGFldStd[PGGFldStdElt]: PGGFldWrap;
declare type PGGFldStdElt: PGGFldWrapElt;

// Represents fields as the fixed field of some subgroup of some universal galois group
declare type PGGStrFldGrp[PGGFldGrp];
declare attributes PGGStrFldGrp: group, prime, splitting_polynomial, ramification_groups, top_field, base_field, ramification_fields, lower_breaks;
declare type PGGFldGrp: PGGFld;
declare attributes PGGFldGrp: universe, group, degree, normal_closure, defining_polynomial, ramification_tower, ramification_degree, inertia_degree;

intrinsic PGGFldStd_Make(K :: FldPad) -> PGGFldStd
  {The p-adic field K.}
  F := New(PGGFldStd);
  F`actual := K;
  return F;
end intrinsic;

intrinsic PGGFldStd_Make(p :: RngIntElt) -> PGGFldStd
  {The p-adic field.}
  require p gt 0 and IsPrime(p): "p must be prime";
  return PGGFldStd_Make(pAdicField(p));
end intrinsic;

intrinsic PGGFldStd_Make(K :: FldPad, F :: PGGFldStd) -> PGGFldStd
  {The p-adic field K as an extension of F.}
  t := [];
  L := K;
  while true do
    if L eq Actual(F) then
      break;
    elif IsPrimeField(L) then
      error "K is not an extension of F";
    else
      Append(~t, L);
      L := BaseField(L);
    end if;
  end while;
  E := F;
  for L in Reverse(t) do
    E0 := E;
    E := PGGFldStd_Make(L);
    E`base_field := E0;
  end for;
  return E;
end intrinsic;

intrinsic Print(F :: PGGFldWrap, lvl :: MonStgElt)
  {Print.}
  printf "%O", Actual(F), lvl;
end intrinsic;

intrinsic Actual(F :: PGGFldWrap) -> FldPad
  {The actual field.}
  return F`actual;
end intrinsic;

intrinsic DefaultPrecision(F :: PGGFldStd) -> RngIntElt
  {The default precision of F.}
  return F`actual`DefaultPrecision;
end intrinsic;

intrinsic SetDefaultPrecision(F :: PGGFldStd, k :: RngIntElt)
  {Sets the default precision of F to k.}
  F`actual`DefaultPrecision := k;
end intrinsic;

intrinsic 'eq'(F1 :: PGGFld, F2 :: PGGFld) -> BoolElt
  {Equality.}
  return IsIdentical(F1, F2);
end intrinsic;

intrinsic 'eq'(F1 :: PGGFldStd, F2 :: PGGFldStd) -> BoolElt
  {"}
  return Actual(F1) eq Actual(F2);
end intrinsic;

intrinsic IsPrimeField(F :: PGGFldStd) -> BoolElt
  {True if F is a completion of Q.}
  return IsPrimeField(Actual(F));
end intrinsic;

intrinsic Prime(F :: PGGFldStd) -> RngIntElt
  {The p in p-adic.}
  return Prime(Actual(F));
end intrinsic;

intrinsic Prime(F :: PGGFldGrp) -> RngIntElt
  {"}
  require assigned Universe(F)`prime: "prime not known";
  return Universe(F)`prime;
end intrinsic;

intrinsic Degree(F :: PGGFldStd) -> RngIntElt
  {The degree of F over its base field.}
  return Degree(Actual(F));
end intrinsic;

intrinsic InertiaDegree(F :: PGGFldStd) -> RngIntElt
  {The inertia degree of F over its base field.}
  return InertiaDegree(Actual(F));
end intrinsic;

intrinsic RamificationDegree(F :: PGGFldStd) -> RngIntElt
  {The ramification degree of F over its base field.}
  return RamificationDegree(Actual(F));
end intrinsic;

intrinsic Degree(E :: PGGFldStd, F :: PGGFldStd) -> RngIntElt
  {The degree of E over F.}
  return Degree(Actual(E), Actual(F));
end intrinsic;

intrinsic InertiaDegree(E :: PGGFldStd, F :: PGGFldStd) -> RngIntElt
  {The inertia degree of E over F.}
  return InertiaDegree(Actual(E), Actual(F));
end intrinsic;

intrinsic RamificationDegree(E :: PGGFldStd, F :: PGGFldStd) -> RngIntElt
  {The ramification degree of E over F.}
  return RamificationDegree(Actual(E), Actual(F));
end intrinsic;

intrinsic AbsoluteDegree(F :: PGGFldStd) -> RngIntElt
  {The absolute degree of F.}
  return AbsoluteDegree(Actual(F));
end intrinsic;

intrinsic AbsoluteInertiaDegree(F :: PGGFldStd) -> RngIntElt
  {The absolute inertia degree of F.}
  return AbsoluteInertiaDegree(Actual(F));
end intrinsic;

intrinsic AbsoluteRamificationDegree(F :: PGGFldStd) -> RngIntElt
  {The absolute ramification degree of F.}
  return AbsoluteRamificationDegree(Actual(F));
end intrinsic;

intrinsic Parent(x :: PGGFldElt) -> PGGFld
  {The parent field of x.}
  return x`parent;
end intrinsic;

intrinsic IsCoercible(F :: PGGFld, X) -> BoolElt, .
  {True if X is coercible into F.}
  return false, "wrong type";
end intrinsic;

intrinsic IsCoercible(F :: PGGFld, X :: PGGFldElt) -> BoolElt, .
  {"}
  if Parent(X) eq F then
    return true, X;
  else
    return false, "wrong parent";
  end if;
end intrinsic;

intrinsic IsCoercible(F :: PGGFldStd, X) -> BoolElt, .
  {"}
  ok, xx := IsCoercible(Actual(F), X);
  if ok then
    return true, Element(F, xx);
  elif assigned xx then
    return false, xx;
  else
    return false, "cannot coerce to actual field";
  end if;
end intrinsic;

intrinsic IsCoercible(F :: PGGFldStd, X :: PGGFldStdElt) -> BoolElt, .
  {"}
  if Parent(X) eq F then
    return true, X;
  end if;
  return IsCoercible(F, Actual(X));
end intrinsic;

intrinsic IsCoercible(F :: PGGFldStd, X :: FldPadElt) -> BoolElt, .
  {"}
  if Parent(X) eq Actual(F) then
    return true, Element(F, X);
  end if;
  ok, xx := IsCoercible(Actual(F), X);
  if ok then
    idx := xdiv(AbsoluteRamificationDegree(F), AbsoluteRamificationDegree(Parent(X)));
    if (not IsWeaklyZero(X)) /*and Precision(Actual(F)) eq Infinity() and Precision(xx) lt Precision(X)*idx*/ then
      // SetDefaultPrecision(F, Precision(X)*idx);
      // xx := Actual(F) ! X;
      assert Precision(xx) ge Precision(X)*idx;
    end if;
    return true, Element(F, xx);
  elif assigned xx then
    return false, xx;
  else
    return false, "cannot coerce to actual field";
  end if;
end intrinsic;

intrinsic IsCoercible(F :: PGGFldStd, X :: []) -> BoolElt, .
  {"}
  if IsPrimeField(F) then
    return false, "not an extension";
  end if;
  ok, Y := CanChangeUniverse(X, BaseField(F));
  if ok then
    return true, Element(F, Actual(F) ! [Actual(x) : x in Y]);
  elif assigned Y then
    return false, Y;
  else
    return false, "cannot coerce to base field";
  end if;
end intrinsic;

intrinsic Element(F :: PGGFldStd, actual :: FldPadElt) -> PGGFldStdElt
  {An element of F.}
  require Parent(actual) eq Actual(F): "actual must be an element of Actual(F)";
  x := New(PGGFldStdElt);
  x`parent := F;
  x`actual := actual;
  return x;
end intrinsic;

intrinsic Print(x :: PGGFldWrapElt, lvl :: MonStgElt)
  {Print.}
  printf "%o", Actual(x), lvl;
end intrinsic;

intrinsic RationalApproximation(x :: PGGFldStdElt) -> FldRatElt
  {A rational number approximating x.}
  require IsPrimeField(Parent(x)): "x must be in a prime field";
  return RationalField() ! Actual(x);
end intrinsic;

intrinsic Actual(x :: PGGFldWrapElt) -> FldPadElt
  {The actual value of x.}
  return x`actual;
end intrinsic;

intrinsic AutomorphismGroup(E :: PGGFldStd, F :: PGGFldStd) -> GrpPerm
  {The automorphism group of E/F.}
  return AutomorphismGroup(Actual(E), Actual(F));
end intrinsic;

intrinsic AutomorphismGroup(E :: PGGFldGrp, F :: PGGFldGrp) -> GrpPerm
  {"}
  require IsExtensionOf(E, F): "E must be an extension of F";
  return CosetImage(Normalizer(F`group, E`group), E`group);
end intrinsic;

intrinsic RamificationTower(E :: PGGFldStd, F :: PGGFldStd) -> []
  {The ramification tower of E/F.}
  t := ramification_tower(Actual(E), Actual(F));
  return [i eq 1 select F else PGGFldStd_Make(t[i], Self(i-1)) : i in [1..#t]];
end intrinsic;

intrinsic IsWeaklyZero(x :: PGGFldStdElt) -> BoolElt
  {True if x is weakly zero.}
  return IsWeaklyZero(Actual(x));
end intrinsic;

intrinsic WeakValuation(x :: PGGFldStdElt) -> RngIntElt
  {A lower bound on the valuation of x. It is correct if x is not weakly zero.}
  return Valuation(Actual(x));
end intrinsic;

intrinsic Valuation(x :: PGGFldStdElt) -> RngIntElt
  {The valuation of x. Returns an error if x is weakly but not precisely zero.}
  v := WeakValuation(x);
  if IsWeaklyZero(x) and v lt Infinity() then
    precision_error();
  end if;
  return v;
end intrinsic;

intrinsic AbsolutePrecision(x :: PGGFldStdElt) -> RngIntElt
  {Absolute precision.}
  return AbsolutePrecision(Actual(x));
end intrinsic;

intrinsic ValuationEq(x :: PGGFldStdElt, n :: RngIntElt) -> BoolElt
  {True if the valuation of x is n.}
  if WeakValuation(x) gt n then
    return false;
  elif IsWeaklyZero(x) then
    precision_error();
  else
    return WeakValuation(x) eq n;
  end if;
end intrinsic;

intrinsic ValuationGe(x :: PGGFldStdElt, n :: RngIntElt) -> BoolElt
  {True if the valuation of x is at least n.}
  if WeakValuation(x) ge n then
    return true;
  elif IsWeaklyZero(x) then
    precision_error();
  else
    return false;
  end if;
end intrinsic;

intrinsic ResidueClassField(F :: PGGFldWrap) -> FldFin, Map
  {The residue class field of F.}
  if not assigned F`residue_class_field then
    FF, m := _ResidueClassField(F);
    assert Type(FF) eq FldFin;
    assert Characteristic(FF) eq Prime(F);
    assert Type(m) eq Map;
    assert Domain(m) eq F;
    assert Codomain(m) eq FF;
    F`residue_class_field := [* FF, m *];
  end if;
  return Explode(F`residue_class_field);
end intrinsic;

intrinsic _ResidueClassField(F :: PGGFldStd) -> FldFin, Map
  {"}
  FF, m := ResidueClassField(Integers(Actual(F)));
  return FF, map<F -> FF | x :-> Actual(x) @ m, y :-> Actual(F) ! (y @@ m)>;
end intrinsic;

intrinsic ChangeAbsolutePrecision(x :: PGGFldStdElt, n :: RngIntElt) -> PGGFldStdElt
  {Changes the absolute precision of x to n.}
  if IsWeaklyZero(x) or ValuationGe(x, n) then
    return Zero(Parent(x), n);
  else
    return Parent(x) ! ChangePrecision(Actual(x), n - Valuation(x));
  end if;
end intrinsic;

intrinsic Zero(F :: PGGFldStd, n :: RngIntElt) -> PGGFldStdElt
  {The zero of F to absolute precision n.}
  z := (Actual(F)!1) - (Actual(F)!1);
  return F ! ShiftValuation(z, n - AbsolutePrecision(z));
end intrinsic;

intrinsic HasIsomorphism(L1 :: PGGFldStd, L2 :: PGGFldStd, K :: PGGFldStd : MaximizeAPr:=true) -> BoolElt, Map
  {True if there is a K-isomorphism L1 to L2.}

  // check inputs
  ok, t1 := IsExtensionOf(L1, K);
  assert ok;
  ok, t2 := IsExtensionOf(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);
    roots1 := Roots(ChangeRing(f1, L2));
    if #roots1 eq 0 then
      error if HasRoot(ChangeRing(f2, L1)), "precision error";
      return false, _;
    end if;
    root1 := roots1[1];
    roots2 := Roots(ChangeRing(f2, L1));
    error if #roots2 eq 0, "precision error";
    idxs := [i : i in [1..#roots2] | IsWeaklyEqual(trim_apr(Actual(L1).1,FUDGE), trim_apr(&+[Actual(L1)| Actual(cs[i]) * r^(i-1) : i in [1..#cs]],FUDGE) where cs:=Eltseq(root1)) where r:=Actual(roots2[i])];
    error if #idxs ne 1, "precision error";
    root2 := roots2[idxs[1]];
    if MaximizeAPr then
      root1 := MaximizeAbsolutePrecision(root1);
      root2 := MaximizeAbsolutePrecision(root2);
    end if;
    return true, map<L1 -> L2 |
      x :-> &+[Actual(L2)| Actual(cs[i]) * Actual(root1)^(i-1) : i in [1..#cs]] where cs:=Eltseq(x),
      y :-> &+[Actual(L1)| Actual(cs[i]) * Actual(root2)^(i-1) : i in [1..#cs]] where cs:=Eltseq(y)>;
  end if;

  // general case
  not_implemented("HasIsomorphism: general towers of extensions");
end intrinsic;

intrinsic IsExtensionOf(L :: PGGFldStd, K :: PGGFldStd) -> BoolElt, []
  {True if L is an extension of K. If so, returns the sequence of fields from K to L.}
  if L eq K then
    return true, [K];
  elif IsPrimeField(L) then
    return false, _;
  else
    ok, twr := IsExtensionOf(BaseField(L), K);
    if ok then
      return true, Append(twr, L);
    else
      return false, _;
    end if;
  end if;
end intrinsic;

intrinsic BaseField(F :: PGGFldWrap) -> PGGFldWrap
  {The base field of F.}
  if not assigned F`base_field then
    require not IsPrimeField(F): "F must be an extension";
    F`base_field := _BaseField(F);
  end if;
  return F`base_field;
end intrinsic;

intrinsic _BaseField(F :: PGGFldStd) -> PGGFldStd
  {"}
  return PGGFldStd_Make(BaseField(Actual(F)));
end intrinsic;

intrinsic MaximizeAbsolutePrecision(x :: PGGFldStdElt) -> PGGFldStdElt
  {Maximizes the absolute precision of x.}
  if IsWeaklyZero(x) then
    return Parent(x) ! 0;
  else
    return Parent(x) ! ChangePrecision(Actual(x), Precision(Parent(Actual(x))));
  end if;
end intrinsic;

intrinsic DefiningPolynomial(E :: PGGFldWrap) -> PGGPol
  {The defining polynomial of E over its base field.}
  if not assigned E`defining_polynomial then
    E`defining_polynomial := _DefiningPolynomial(E);
  end if;
  return E`defining_polynomial;
end intrinsic;

intrinsic _DefiningPolynomial(E :: PGGFldStd) -> PGGPolStd
  {"}
  return PolynomialRing(BaseField(E)) ! DefiningPolynomial(Actual(E));
end intrinsic;

intrinsic DefiningPolynomial(E :: PGGFldStd, F :: PGGFldStd) -> PGGPolStd
  {The defining polynomial of E over F.}
  return PolynomialRing(F) ! DefiningPolynomial(Actual(E), Actual(F));
end intrinsic;

intrinsic ExtConstructor(F :: PGGFldStd, t :: Tup) -> PGGFldStd
  {An extension of F defined by t.}
  if #t eq 1 and Type(t[1]) eq PGGPolStd then
    return PGGFldStd_Make(ext<Actual(F) | Actual(t[1])>, F);
  elif #t eq 1 and Type(t[1]) eq FldPad then
    PGGFldStd_Make(t[1], F);
  else
    error "bad constructor";
  end if;
end intrinsic;

intrinsic '&+'(xs :: [PGGFldElt]) -> PGGFldElt
  {Sum.}
  not_implemented("&+:", ElementType(Universe(xs)));
end intrinsic;

intrinsic '&+'(xs :: [PGGFldStdElt]) -> PGGFldStdElt
  {Sum.}
  F := Universe(xs);
  return F ! &+[Actual(F)| Actual(x) : x in xs];
end intrinsic;

intrinsic '+'(x :: PGGFldElt, y :: PGGFldElt) -> PGGFldElt
  {Plus.}
  ok, F := ExistsCoveringStructure(Parent(x), Parent(y));
  require ok: "different fields";
  return &+[F|x,y];
end intrinsic;

intrinsic '&*'(xs :: [PGGFldElt]) -> PGGFldElt
  {Product.}
  not_implemented("&*:", ElementType(Universe(xs)));
end intrinsic;

intrinsic '&*'(xs :: [PGGFldStdElt]) -> PGGFldStdElt
  {Product.}
  F := Universe(xs);
  return F ! &*[Actual(F)| Actual(x) : x in xs];
end intrinsic;

intrinsic '*'(x :: PGGFldElt, y :: PGGFldElt) -> PGGFldElt
  {Multiply.}
  ok, F := ExistsCoveringStructure(Parent(x), Parent(y));
  require ok: "different fields";
  return &*[F|x,y];
end intrinsic;

intrinsic '^'(x :: PGGFldStdElt, n :: RngIntElt) -> PGGFldElt
  {Power.}
  return Parent(x) ! (Actual(x)^n);
end intrinsic;

intrinsic '-'(x :: PGGFldStdElt) -> PGGFldStdElt
  {Negation.}
  return Parent(x) ! (-Actual(x));
end intrinsic;

intrinsic '/'(x :: PGGFldStdElt, y :: PGGFldStdElt) -> PGGFldStdElt
  {Divide.}
  ok, F := ExistsCoveringStructure(Parent(x), Parent(y));
  require ok: "different fields";
  return F ! (Actual(F ! x) / Actual(F ! y));
end intrinsic;

intrinsic '-'(x :: PGGFldStdElt, y :: PGGFldStdElt) -> PGGFldStdElt
  {Minus.}
  ok, F := ExistsCoveringStructure(Parent(x), Parent(y));
  require ok: "different fields";
  return F ! (Actual(F!x) - Actual(F!y));
end intrinsic;

intrinsic ExistsCoveringStructure(E :: PGGFldStd, F :: PGGFldStd) -> BoolElt, .
  {True if there is a structure containing both E and F.}
  ok, A := ExistsCoveringStructure(Actual(E), Actual(F));
  if ok then
    if A eq Actual(E) then
      return true, E;
    elif A eq Actual(F) then
      return true, F;
    else
      assert false;
    end if;
  else
    return false, _;
  end if;
end intrinsic;

intrinsic '.'(F :: PGGFldStd, n :: RngIntElt) -> BoolElt, .
  {The nth generator of F.}
  require n eq 1: "n must be 1";
  return F![0,1];
end intrinsic;

intrinsic PGGFldGrp_Make(G :: Grp : Prime:=false, RamificationGroups:=false, LowerBreaks:=false) -> PGGFldGrp
  {"}
  U := New(PGGStrFldGrp);
  U`group := G;
  if Prime cmpne false then
    p := Prime;
    require Type(p) eq RngIntElt and p gt 0 and IsPrime(p): "Prime must be a prime integer";
    U`prime := p;
  end if;
  if RamificationGroups cmpne false then
    R := RamificationGroups;
    require Type(R) eq SeqEnum: "RamificationGroups must be a sequence";
    require #R ge 1: "RamificationGroups must have length at least 1";
    require forall{N : N in R | Type(N) eq Type(G) and N subset G and IsNormal(G,N)}: "RamificationGroups must be normal subgroups of G";
    require R[1] eq G: "First RamificationGroup must be G";
    require R[#R] eq sub<G|Id(G)>: "Last RamificationGroup must be trivial";
    require forall{i : i in [2..#R] | R[i] ne R[i-1] and R[i] subset R[i-1]}: "RamificationGroups must be properly descending";
    U`ramification_groups := R;
  end if;
  if LowerBreaks cmpne false then
    vs := LowerBreaks;
    require Type(vs) eq SeqEnum and #vs ge 1 and Universe(vs) cmpeq Integers() and vs[1] eq -1 and forall{i : i in [2..#vs] | vs[i-1] lt vs[i]}: "LowerBreaks must be an increasing sequence of integers starting at -1";
    if assigned U`ramification_groups then
      require #vs eq #U`ramification_groups: "LowerBreaks must have the same length as RamificationGroups";
    end if;
    U`lower_breaks := vs;
  end if;
  U`top_field := Field(U, sub<U`group | Id(U`group)>);
  U`base_field := Field(U, U`group);
  return U`top_field;
end intrinsic;

intrinsic Parent(F :: PGGFldGrp) -> PGGStrFldGrp
  {The universe of F.}
  return F`universe;
end intrinsic;

intrinsic Universe(F :: PGGFldGrp) -> PGGStrFldGrp
  {"}
  return Parent(F);
end intrinsic;

intrinsic Print(U :: PGGStrFldGrp)
  {Print.}
  printf "Set of subfields of a field with Galois group %o", U`group;
end intrinsic;

intrinsic Field(U :: PGGStrFldGrp, G :: Grp) -> PGGFldGrp
  {The fixed field of G.}
  require G subset U`group: "G must be a subgroup of the universe group";
  F := New(PGGFldGrp);
  F`group := G;
  F`universe := U;
  return F;
end intrinsic;

intrinsic BaseField(F :: PGGFldGrp) -> PGGFldGrp
  {The base field of F.}
  return Universe(F)`base_field;
end intrinsic;

intrinsic UniverseField(F :: PGGFldGrp) -> PGGFldGrp
  {The universe field of F, inside which everything lives.}
  return Universe(F)`top_field;
end intrinsic;

intrinsic UniverseGroup(F :: PGGFldGrp) -> Grp
  {The universe group, the Galois group of the universe field.}
  return Universe(F)`group;
end intrinsic;

intrinsic UniverseSplittingPolynomial(F :: PGGFldGrp) -> PGGPolGrp
  {A polynomial with Galois group equal to the universe group.}
  if not assigned F`universe`splitting_polynomial then
    U := UniverseField(F);
    G := UniverseGroup(F);
    f := &*[PolynomialRing(BaseField(U))| DefiningPolynomial(Subfield(U, Stabilizer(G,Rep(o)))) : o in Orbits(G)];
    f`galois_group_quo := hom<G -> G | [G.i : i in [1..Ngens(G)]]>;
    F`universe`splitting_polynomial := f;
  end if;
  return F`universe`splitting_polynomial;
end intrinsic;

intrinsic Subfield(F :: PGGFldGrp, G :: Grp) -> PGGFldGrp
  {The subfield of F generated by the group G.}
  require F`group subset G: "G must contain the group defining F";
  return Field(Universe(F), G);
end intrinsic;

intrinsic Extension(F :: PGGFldGrp, G :: Grp) -> PGGFldGrp
  {The extension of F generated by G.}
  require G subset F`group: "The group defining F must contain G";
  return Field(Universe(F), G);
end intrinsic;

intrinsic Print(F :: PGGFldGrp, lvl :: MonStgElt)
  {Print.}
  case lvl:
  when "Magma":
    if F eq UniverseField(F) then
      printf "PGGFldGrp_Make(%m", F`group;
      params := [];
      if assigned F`prime then
        Append(~params, Sprintf("Prime:=%m", F`prime));
      end if;
      if assigned F`ramification_groups then
        Append(~params, Sprintf("RamificationGroups:=%m", F`ramification_groups));
      end if;
      if #params ne 0 then
        printf " : %o", Join(params, ", ");
      end if;
      printf ")";
    else
      printf "Subfield(%m, %m)", UniverseField(F), F`group;
    end if;
  else
    printf "Field of degree %o", Degree(F);
    if F eq UniverseField(F) then
      printf " (the universe)";
    end if;
  end case;
end intrinsic;

intrinsic 'eq'(U :: PGGStrFldGrp, V :: PGGStrFldGrp) -> BoolElt
  {Equality.}
  return IsIdentical(U, V);
end intrinsic;

intrinsic 'eq'(E :: PGGFldGrp, F :: PGGFldGrp) -> BoolElt
  {Equality.}
  require Universe(E) eq Universe(F): "different universes";
  return E`group eq F`group;
end intrinsic;

intrinsic Degree(F :: PGGFldGrp) -> RngIntElt
  {Degree of F over its base field.}
  if not assigned F`degree then
    F`degree := Degree(F, BaseField(F));
  end if;
  return F`degree;
end intrinsic;

intrinsic Degree(E :: PGGFldGrp, F :: PGGFldGrp) -> RngIntElt
  {Degree of E/F.}
  require IsExtensionOf(E, F): "E must be an extension of F";
  return Index(F`group, E`group);
end intrinsic;

intrinsic IsExtensionOf(E :: PGGFldGrp, F :: PGGFldGrp) -> BoolElt
  {True if E is an extension of F.}
  require Universe(E) eq Universe(F): "different universes";
  return E`group subset F`group;
end intrinsic;

intrinsic 'subset'(F :: PGGFldGrp, E :: PGGFldGrp) -> BoolElt
  {True if F is a subfield of E.}
  return IsExtensionOf(E, F);
end intrinsic;

intrinsic 'join'(F :: PGGFldGrp, E :: PGGFldGrp) -> PGGFldGrp
  {Compositum.}
  require Universe(E) eq Universe(F): "different universes";
  return Field(Universe(F), E`group meet F`group);
end intrinsic;

intrinsic 'meet'(F :: PGGFldGrp, E :: PGGFldGrp) -> PGGFldGrp
  {Intersection.}
  require Universe(E) eq Universe(F): "different universes";
  return Field(Universe(F), sub<UniverseGroup(F) | E`group, F`group>);
end intrinsic;

intrinsic NormalClosure(E :: PGGFldGrp, F :: PGGFldGrp) -> PGGFldGrp
  {Normal closure of E/F.}
  require IsExtensionOf(E, F): "E must be an extension of F";
  return Field(Universe(F), Core(F`group, E`group));
end intrinsic;

intrinsic NormalClosure(F :: PGGFldGrp) -> PGGFldGrp
  {Normal closure of F over its base field.}
  if not assigned F`normal_closure then
    F`normal_closure := NormalClosure(F, BaseField(F));
  end if;
  return F`normal_closure;
end intrinsic;

intrinsic IsIsomorphic(E1 :: PGGFldGrp, E2 :: PGGFldGrp, F :: PGGFldGrp) -> BoolElt, GrpElt
  {True if E1 is isomorphic to E2 over F. If so, also returns an element of the group defining F conjugating E1 to E2.}
  require IsExtensionOf(E1,F) and IsExtensionOf(E2,F): "E1 and E2 must be extensions of F";
  ok, g := IsConjugate(F`group, E1`group, E2`group);
  if ok then
    assert g in F`group;
    assert Conjugate(E1,g) eq E2;
    return true, g;
  else
    return false, _;
  end if;
end intrinsic;

intrinsic IsIsomorphic(E1 :: PGGFldGrp, E2 :: PGGFldGrp) -> BoolElt, GrpElt
  {True if E1 is isomorphic to E2 over the base field. If so, also returns an element of the group defining F conjugating E1 to E2.}
  require Universe(E1) eq Universe(E2): "different universes";
  return IsIsomorphic(E1, E2, Universe(E1)`base_field);
end intrinsic;

intrinsic Conjugate(F :: PGGFldGrp, g :: GrpElt) -> PGGFldGrp
  {Conjugation.}
  require g in UniverseGroup(F): "g must be an element of the universe group";
  return Field(Universe(F), F`group^g);
end intrinsic;

intrinsic IsCoercible(U :: PGGStrFldGrp, F) -> BoolElt, .
  {True if F is coercible into U.}
  return false, "wrong type";
end intrinsic;

intrinsic IsCoercible(U :: PGGStrFldGrp, F :: PGGFldGrp) -> BoolElt, .
  {"}
  if Universe(F) eq U then
    return true, F;
  else
    return false, _;
  end if;
end intrinsic;

intrinsic UniverseRamificationGroups(F :: PGGFldGrp) -> []
  {The ramification groups of the universe.}
  require assigned Universe(F)`ramification_groups: "not known";
  return Universe(F)`ramification_groups;
end intrinsic;

intrinsic UniverseRamificationFields(F :: PGGFldGrp) -> []
  {The ramification fields of the universe.}
  U := Universe(F);
  if not assigned U`ramification_fields then
    U`ramification_fields := [U| Field(U, N) : N in UniverseRamificationGroups(F)];
  end if;
  return U`ramification_fields;
end intrinsic;

intrinsic UniverseFindLowerBreak(F :: PGGFldGrp, v :: FldRatElt) -> RngIntElt
  {The index of the vth ramification group in the sequence of distinct ramification groups.}
  U := F`universe;
  require assigned U`ramification_groups and assigned U`lower_breaks: "universe must have ramification groups and lower breaks";
  if v le -1 then
    return 1;
  end if;
  i := Max([i : i in [1..#U`lower_breaks] | U`lower_breaks[i] le v]);
  return i;
end intrinsic;

intrinsic UniverseFindLowerBreak(F :: PGGFldGrp, v :: RngIntElt) -> RngIntElt
  {"}
  return UniverseFindLowerBreak(F, RationalField()!v);
end intrinsic;

intrinsic UniverseRamificationGroup(F :: PGGFldGrp, v) -> Grp
  {The vth ramification group.}
  return F`universe`ramification_groups[UniverseFindLowerBreak(F,v)];
end intrinsic;

intrinsic UniverseRamificationField(F :: PGGFldGrp, v) -> Grp
  {The vth ramification group.}
  return UniverseRamificationFields()[UniverseFindLowerBreak(F,v)];
end intrinsic;

intrinsic RamificationTower(E :: PGGFldGrp, F :: PGGFldGrp) -> [], []
  {The tower of ramification fields of E/F.}
  require IsExtensionOf(E,F): "E must be an extension of F";
  Ns := UniverseRamificationFields(F);
  T1 := [(F join N) meet E : N in Ns];
  T2 := [F join (N meet E) : N in Ns];
  assert T1 eq T2;
  idxs := [i : i in [1..#T1] | i eq 1 or T1[i] ne T1[i-1]];
  T := [T1[i] : i in idxs];
  if assigned Universe(F)`lower_breaks then
    vs := Universe(F)`lower_breaks;
    return T, [vs[i] : i in idxs];
  else
    return T, _;
  end if;
end intrinsic;

intrinsic RamificationTower(F :: PGGFldGrp) -> [], []
  {The ramification fields of F over its base field.}
  if not assigned F`ramification_tower then
    Fs, vs := RamificationTower(F, BaseField(F));
    F`ramification_tower := assigned vs select [*Fs,vs*] else [*Fs*];
  end if;
  if #F`ramification_tower eq 1 then
    return F`ramification_tower[1], _;
  else
    return F`ramification_tower[1], F`ramification_tower[2];
  end if;
end intrinsic;

intrinsic RamificationField(E :: PGGFldGrp, F :: PGGFldGrp, v) -> PGGFldGrp
  {The vth ramification field between E and F.}
  return RamificationTower(E,F)[UniverseFindLowerBreak(F,v)];
end intrinsic;

intrinsic RamificationField(F :: PGGFldGrp, v) -> PGGFldGrp
  {The vth ramification field below F.}
  return RamificationTower(F)[UniverseFindLowerBreak(F,v)];
end intrinsic;

intrinsic RamificationDegree(E :: PGGFldGrp, F :: PGGFldGrp, v) -> RngIntElt
  {The v-ramification degree of E/F.}
  return Degree(E, RamificationField(E,F,v));
end intrinsic;

intrinsic RamificationDegree(F :: PGGFldGrp, v) -> RngIntElt
  {The v-ramification degree of F over its base field.}
  return Degree(F, RamificationField(F,v));
end intrinsic;

intrinsic RamificationDegree(E :: PGGFldGrp, F :: PGGFldGrp) -> RngIntElt
  {The ramification degree of E/F.}
  return RamificationDegree(E, F, 0);
end intrinsic;

intrinsic RamificationDegree(F :: PGGFldGrp) -> RngIntElt
  {The ramification degree of F over its base field.}
  if not assigned F`ramification_degree then
    F`ramification_degree := RamificationDegree(F, 0);
  end if;
  return F`ramification_degree;
end intrinsic;

intrinsic InertiaDegree(E :: PGGFldGrp, F :: PGGFldGrp) -> RngIntElt
  {The inertia degree of E/F.}
  return xdiv(Degree(E,F), RamificationDegree(E,F));
end intrinsic;

intrinsic InertiaDegree(F :: PGGFldGrp) -> RngIntElt
  {The inertia degree of F over its base field.}
  if not assigned F`inertia_degree then
    F`inertia_degree := xdiv(Degree(F), RamificationDegree(F));
  end if;
  return F`inertia_degree;
end intrinsic;

intrinsic RationalApproximation(x :: PGGFldStdElt, pr :: RngIntElt) -> FldRatElt
  {A rational approximation of x, to absolute precision pr.}
  ok, xx := IsCoercible(Q, Actual(x));
  require ok: xx;
  require pr le AbsolutePrecision(x): "pr must be at most the absolute precision of x";
  return xx;
end intrinsic;

intrinsic Eltseq(x :: PGGFldWrapElt) -> []
  {Represents x as a sequence of elements of its base field.}
  if not assigned x`eltseq then
    require not IsPrimeField(Parent(x)): "x must lie in an extension";
    cs := _Eltseq(x);
    assert Universe(cs) eq BaseField(Parent(x));
    x`eltseq := cs;
  end if;
  return x`eltseq;
end intrinsic;

intrinsic _Eltseq(x :: PGGFldStdElt) -> []
  {The coefficients of x as a vector over the base field.}
  F := Parent(x);
  if AbsolutePrecision(x) eq Infinity() then
    assert IsWeaklyZero(x);
    return [BaseField(F)| ];
  else
    return [BaseField(F)| c : c in Eltseq(Actual(x))];
  end if;
end intrinsic;

intrinsic Generator(E :: PGGFldWrap) -> PGGFldWrapElt
  {The generator of E over its base field.}
  if not assigned E`generator then
    require not IsPrimeField(E): "E must be an extension";
    gen := _Generator(E);
    assert Parent(gen) eq E;
    E`generator := gen;
  end if;
  return E`generator;
end intrinsic;

intrinsic _Generator(E :: PGGFldWrap) -> PGGFldWrapElt
  {"}
  not_implemented("_Generator:", Type(E));
end intrinsic;

intrinsic _Generator(E :: PGGFldStd) -> PGGFldStdElt
  {"}
  return E ! (Actual(E).1);
end intrinsic;

intrinsic TameRParameter(E :: PGGFld, F :: PGGFld) -> RngIntElt, RngIntElt, RngIntElt
  {Given tame extension `E/F`, returns r such that `E=K(zeta,(pi*zeta^r)^(1/e))`.}
  not_implemented("TameRParameter:", Type(E));
end intrinsic;

intrinsic TameRParameter(E :: PGGFldWrap, F :: PGGFldWrap : Alg:="Roots") -> RngIntElt
  {"}
  f := InertiaDegree(E, F);
  e := RamificationDegree(E, F);
  p := Prime(F);
  require not IsDivisibleBy(e, p): "E/F must be tame";
  if e eq 1 then
    return 0;
  end if;
  assert e gt 1;
  q := p^AbsoluteInertiaDegree(F);
  piF := UniformizingElement(F);
  piE := UniformizingElement(E);
  zr := piE^e / piF;
  FF := ResidueClassField(F);
  assert q eq #FF;
  FE, modq := ResidueClassField(E);
  assert Degree(FE, FF) eq f;
  zetar := zr @ modq;
  if zetar eq 1 then
    return 0;
  end if;
  case Alg:
  when "Log":
    // This method selects a primitive root of unity `zeta` and then finds `r = log_zeta(zetar)`, but this logarithm gets slow as the residue field gets large
    zeta := FE.1;
    assert IsPrimitive(zeta);
    r := Log(zeta, zetar);
  when "Roots":
    // This method tries successive values of `r` until some `r`th root of `zetar` is primitive (this takes advantage of the fact that we don't care about the choice of the primitive root `zeta`)
    // TODO: Are there `r` that we can skip over?
    // TODO: Do we really need to check all roots? Does one suffice?
    r := 1;
    while true do
      if exists{zeta : zeta in AllRoots(zetar, r) | IsPrimitive(zeta)} then
        break;
      else
        r +:= 1;
      end if;
    end while;
  else
    not_implemented("TameRParameter: Alg:", Alg);
  end case;
  return r mod GCD(e, q^f-1);
end intrinsic;

intrinsic UniformizingElement(F :: PGGFldStd) -> PGGFldStdElt
  {A uniformizing element.}
  return F ! UniformizingElement(Actual(F));
end intrinsic;