PolRing.mag

import "Utils.mag": not_implemented;

// base type
declare type PGGRngPol[PGGPol];
declare attributes PGGRngPol: base_ring;
declare attributes PGGPol: parent, coefficients, roots, num_roots, has_root, factorization, factor_degrees, extension, degree, splitting_field;

// base wrapping type
declare type PGGRngPolWrap[PGGPolWrap]: PGGRngPol;
declare attributes PGGRngPolWrap: actual;
declare type PGGPolWrap: PGGPol;
declare attributes PGGPolWrap: actual;

// standard p-adics (FldPad)
declare type PGGRngPolStd[PGGPolStd]: PGGRngPolWrap;
declare type PGGPolStd: PGGPolWrap;

// group representation
declare type PGGRngPolGrp[PGGPolGrp]: PGGRngPol;
declare type PGGPolGrp: PGGPol;
declare attributes PGGPolGrp: factor_fields, galois_group_quo;

intrinsic PolynomialRing(F :: PGGFld) -> PGGRngPol
  {The univariate polynomial ring of F.}
  if not assigned F`polynomial_ring then
    F`polynomial_ring := _PolynomialRing(F);
  end if;
  return F`polynomial_ring;
end intrinsic;

intrinsic _PolynomialRing(F :: PGGFld) -> PGGRngPol
  {"}
  not_implemented("_PolynomialRing:", Type(F));
end intrinsic;

intrinsic _PolynomialRing(F :: PGGFldStd) -> PGGRngPolStd
  {"}
  R := New(PGGRngPolStd);
  R`base_ring := F;
  return R;
end intrinsic;

intrinsic _PolynomialRing(F :: PGGFldGrp) -> PGGRngPolGrp
  {"}
  R := New(PGGRngPolGrp);
  R`base_ring := F;
  return R;
end intrinsic;

intrinsic BaseRing(R :: PGGRngPol) -> PGGFld
  {The base ring of R.}
  return R`base_ring;
end intrinsic;

intrinsic Parent(f :: PGGPol) -> PGGRngPol
  {The parent ring of f.}
  return f`parent;
end intrinsic;

intrinsic BaseRing(f :: PGGPol) -> PGGFld
  {The base ring of f.}
  return BaseRing(Parent(f));
end intrinsic;

intrinsic 'eq'(R1 :: PGGRngPol, R2 :: PGGRngPol) -> BoolElt
  {Equality.}
  return IsIdentical(R1, R2);
end intrinsic;

intrinsic 'eq'(R1 :: PGGRngPolStd, R2 :: PGGRngPolStd) -> BoolElt
  {"}
  return BaseRing(R1) eq BaseRing(R2);
end intrinsic;

intrinsic 'eq'(R1 :: PGGRngPolGrp, R2 :: PGGRngPolGrp) -> BoolElt
  {"}
  return BaseRing(R1) eq BaseRing(R2);
end intrinsic;

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

intrinsic Print(f :: PGGPolStd)
  {"}
  printf "%o", f`actual;
end intrinsic;

intrinsic Actual(R :: PGGRngPolStd) -> RngUPol
  {The actual ring.}
  if not assigned R`actual then
    R`actual := PolynomialRing(Actual(BaseRing(R)));
  end if;
  return R`actual;
end intrinsic;

intrinsic Element(R :: PGGRngPolStd, xf :: RngUPolElt[FldPad]) -> PGGPolStd
  {An element of R.}
  require Parent(xf) eq Actual(R): "xf must be an element of Actual(R)";
  f := New(PGGPolStd);
  f`actual := xf;
  f`parent := R;
  return f;
end intrinsic;

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

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

intrinsic IsCoercible(R :: PGGRngPolStd, X) -> BoolElt, .
  {"}
  ok, Y := IsCoercible(Actual(R), X);
  if ok then
    return true, Element(R, Y);
  else
    return false, "not coercible to the actual ring" cat (assigned Y select ": " cat Y else "");
  end if;
end intrinsic;

intrinsic IsCoercible(R :: PGGRngPolStd, X :: PGGPolStd) -> BoolElt, .
  {"}
  if Parent(X) eq R then
    return true, X;
  else
    return IsCoercible(R, Actual(X));
  end if;
end intrinsic;

intrinsic IsCoercible(R :: PGGRngPolStd, X :: []) -> BoolElt, .
  {"}
  ok, Y := CanChangeUniverse(X, BaseRing(R));
  if ok then
    return true, Element(R, Actual(R)![Actual(x) : x in Y]);
  else
    return false, "not coercible to base ring" cat (assigned Y select ": " cat Y else "");
  end if;
end intrinsic;

intrinsic IsCoercible(R :: PGGRngPolGrp, X :: PGGPolGrp) -> BoolElt, .
  {"}
  F := BaseRing(X);
  E := BaseRing(R);
  if Universe(F) eq Universe(E) then
    if Parent(X) eq R then
      return true, X;
    elif F subset E then
      _, q := GaloisGroup(X);
      G := E`group @ q;
      return true, &*[R| DefiningPolynomial(Field(Universe(F), Stabilizer(G, Rep(o)) @@ q), E) : o in Orbits(G)];
    end if;
  end if;
  return false, _;
end intrinsic;

intrinsic PGG_Polynomial(f) -> PGGPol
  {Converts f to a PGGPol.}
  error "not coercible to a PGGPol:", Type(f);
end intrinsic;

intrinsic PGG_Polynomial(f :: PGGPol) -> PGGPol
  {"}
  return f;
end intrinsic;

intrinsic PGG_Polynomial(f :: RngUPolElt[FldPad]) -> PGGPol
  {"}
  return PolynomialRing(PGGFldStd_Make(BaseRing(f))) ! f;
end intrinsic;

intrinsic Actual(f :: PGGPolStd : FixPr:=false) -> RngUPolElt
  {The current actual value of f.}
  xf := f`actual;
  if FixPr and Precision(BaseRing(xf)) eq Infinity() then
    pr := Max([1] cat [IsWeaklyZero(c) select 0 else Precision(c) : c in Coefficients(xf)]);
    xf := PolynomialRing(ChangePrecision(BaseRing(xf), pr)) ! xf;
  end if;
  return xf;
end intrinsic;

intrinsic Degree(f :: PGGPol) -> RngIntElt
  {Degree.}
  if not assigned f`degree then
    f`degree := _Degree(f);
  end if;
  return f`degree;
end intrinsic;

intrinsic _Degree(f :: PGGPol) -> RngIntElt
  {"}
  not_implemented("_Degree:", Type(f));
end intrinsic;

intrinsic _Degree(f :: PGGPolStd) -> RngIntElt
  {Degree.}
  return Degree(Actual(f));
end intrinsic;

intrinsic _Degree(f :: PGGPolGrp) -> RngIntElt
  {"}
  return &+[Integers()| Degree(E, BaseRing(f)) : E in f`factor_fields];
end intrinsic;

intrinsic Coefficients(f :: PGGPol) -> []
  {The coefficients of f.}
  if not assigned f`coefficients then
    f`coefficients := _Coefficients(f);
  end if;
  return f`coefficients;
end intrinsic;

intrinsic _Coefficients(f :: PGGPol) -> []
  {"}
  not_implemented("_Coefficients:", Type(f));
end intrinsic;

intrinsic _Coefficients(f :: PGGPolStd) -> []
  {"}
  return [BaseRing(f)| c : c in Coefficients(Actual(f))];
end intrinsic;

intrinsic Roots(f :: PGGPol) -> []
  {The roots of f, as a sequence. Assumes f is squarefree, and so does not return multiplicities.}
  if not assigned f`roots then
    f`roots := _Roots(f);
  end if;
  return f`roots;
end intrinsic;

intrinsic _Roots(f :: PGGPol) -> []
  {"}
  not_implemented("_Roots:", Type(f));
end intrinsic;

intrinsic _Roots(f :: PGGPolStd) -> []
  {"}
  return [BaseRing(f)| root : root in PGG_Roots(Actual(f))];
end intrinsic;

intrinsic NumRoots(f :: PGGPol) -> RngIntElt
  {The number of roots of f.}
  if not assigned f`num_roots then
    f`num_roots := _NumRoots(f);
  end if;
  return f`num_roots;
end intrinsic;

intrinsic _NumRoots(f :: PGGPol) -> RngIntElt
  {"}
  return #Roots(f);
end intrinsic;

intrinsic _NumRoots(f :: PGGPolStd) -> RngIntElt
  {"}
  if assigned f`roots then
    return #f`roots;
  else
    return #PGG_Roots(Actual(f) : Lift:=false);
  end if;
end intrinsic;

intrinsic _NumRoots(f :: PGGPolGrp) -> RngIntElt
  {"}
  return Multiplicity(FactorDegrees(f), 1);
end intrinsic;

intrinsic HasRoot(f :: PGGPol) -> BoolElt
  {True if f has a root.}
  if not assigned f`has_root then
    f`has_root := _HasRoot(f);
  end if;
  return f`has_root;
end intrinsic;

intrinsic _HasRoot(f :: PGGPol) -> BoolElt
  {"}
  return NumRoots(f) ne 0;
end intrinsic;

CERT := recformat<F, E, Pi, Rho, Extension>;

intrinsic Factorization(f :: PGGPol : Extensions:=false) -> [], []
  {The factorization of f, as a sequence of irreducible factors. Assumes f is squarefree, and so does not return multiplicities.}
  if (not assigned f`factorization) or (Extensions and not f`factorization[1]) then
    f`factorization := [*Extensions, facs, certs*] where facs, certs := _Factorization(f : Extensions:=Extensions);
  end if;
  return f`factorization[2], f`factorization[3];
end intrinsic;

intrinsic _Factorization(f :: PGGPol : Extensions:=false) -> [], []
  {"}
  not_implemented("_Factorization:", Type(f));
end intrinsic;

intrinsic _Factorization(f :: PGGPolStd : Extensions:=false) -> [], []
  {"}
  facs0, certs0 := PGG_Factorization(Actual(f) : Extensions:=Extensions);
  facs := [Parent(f)| fac : fac in facs0];
  certs := [];
  for c0 in certs0 do
    c := rec<CERT | F:=c0`F, E:=c0`E, Pi:=Parent(f)!c0`Pi, Rho:=Parent(f)!c0`Rho>;
    if Extensions then
      c`Extension := PGGFldStd_Make(c0`Extension, BaseRing(f));
    end if;
    Append(~certs, c);
  end for;
  if Extensions then
    for i in [1..#facs] do
      facs[i]`extension := certs[i]`Extension;
    end for;
  end if;
  return facs, certs;
end intrinsic;

intrinsic FactorDegrees(f :: PGGPol) -> {**}
  {The degrees of factors of f.}
  if not assigned f`factor_degrees then
    f`factor_degrees := _FactorDegrees(f);
  end if;
  return f`factor_degrees;
end intrinsic;

intrinsic _FactorDegrees(f :: PGGPol) -> {**}
  {"}
  return {*Integers()| Degree(fac) : fac in Factorization(f)*};
end intrinsic;

intrinsic _FactorDegrees(f :: PGGPolStd) -> {**}
  {"}
  if assigned f`factorization then
    return {*Degree(fac) : fac in f`factorization*};
  else
    return {*Degree(fac) : fac in PGG_Factorization(Actual(f) : Lift:=false)*};
  end if;
end intrinsic;

intrinsic NumFactors(f :: PGGPol) -> RngIntElt
  {The number of factors of f.}
  return #FactorDegrees(f);
end intrinsic;

intrinsic IsIrreducible(f :: PGGPol) -> RngIntElt
  {True if f is irreducible.}
  return NumFactors(f) eq 1;
end intrinsic;

intrinsic Extension(f :: PGGPol) -> PGGFld
  {The extension defined by f, which must be irreducible.}
  if not assigned f`extension then
    f`extension := _Extension(f);
  end if;
  return f`extension;
end intrinsic;

intrinsic _Extension(f :: PGGPol) -> PGGFld
  {"}
  facs, certs := Factorization(f : Extensions);
  require #facs eq 1: "f must be irreducible";
  return certs[1]`Extension;  
end intrinsic;

intrinsic ChangeRing(f :: PGGPol, F :: PGGFld) -> PGGPol
  {Changes the base ring of f to F.}
  return PolynomialRing(F) ! f;
end intrinsic;

intrinsic IsEisenstein(f :: PGGPol) -> BoolElt
  {True if f is Eisenstein.}
  d := Degree(f);
  return d ge 1 and ValuationEq(Coefficient(f,0),1) and ValuationEq(Coefficient(f,d),0) and forall{i : i in [1..d-1] | ValuationGe(Coefficient(f,i),1)};
end intrinsic;

intrinsic IsInertial(f :: PGGPol) -> BoolElt
  {True if f is inertial.}
  d := Degree(f);
  return d ge 1 and ValuationEq(Coefficient(f,0),0) and ValuationEq(Coefficient(f,d),0) and forall{i : i in [1..d-1] | ValuationGe(Coefficient(f,i),0)} and IsIrreducible(Polynomial([c@m : c in Coefficients(f)] where _,m:=ResidueClassField(BaseRing(f))));
end intrinsic;

intrinsic Coefficient(f :: PGGPol, i :: RngIntElt) -> PGGFldElt
  {The ith coefficient of f.}
  require i ge 0: "i must be at least 0";
  if i le Degree(f) then
    return Coefficients(f)[i+1];
  else
    return BaseRing(f) ! 0;
  end if;
end intrinsic;

intrinsic IsWeaklyEqual(f :: PGGPolStd, g :: PGGPolStd) -> BoolElt
  {True if f and g are weakly equal.}
  return IsWeaklyEqual(Actual(f), Actual(g));
end intrinsic;

intrinsic Polynomial(cs :: [PGGFldElt]) -> PGGPol
  {The polynomial with the given coefficients.}
  return PolynomialRing(Universe(cs)) ! cs;
end intrinsic;

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

intrinsic Print(f :: PGGPolGrp)
  {Print.}
  printf "Polynomial of degree %o", Degree(f);
end intrinsic;

intrinsic SplittingField(f :: PGGPolGrp) -> PGGFldGrp
  {The splitting field of f.}
  if not assigned f`splitting_field then
    F := BaseRing(f);
    f`splitting_field := &join[Universe(F)| NormalClosure(E,F) : E in f`factor_fields];
  end if;
  return f`splitting_field;
end intrinsic;

intrinsic GaloisGroup(f :: PGGPolGrp) -> GrpPerm, Map
  {The Galois group of f, and the quotient map from the defining group of the base field to it.}
  if not assigned f`galois_group_quo then
    F := BaseRing(f);
    G := F`group;
    require Degree(f) gt 0: "must be a non-constant polynomial";
    if Degree(f) eq 1 then
      A := SymmetricGroup(1);
      q := hom<G -> A | [Id(A) : i in [1..Ngens(G)]]>;
    else
      Hs := [E`group : E in f`factor_fields];
      qs := [CosetAction(G,H) : H in Hs];
      As := [Codomain(q) : q in qs];
      A0, incls, projs := DirectProduct(As);
      genimgs := [&*[incls[j](qs[j](g)) : j in [1..#Hs]] where g:=G.i : i in [1..Ngens(G)]];
      A := sub<A0 | genimgs>;
      q := hom<G -> A | genimgs>;
    end if;
    f`galois_group_quo := q;
  end if;
  return Codomain(f`galois_group_quo), f`galois_group_quo;
end intrinsic;

intrinsic '&*'(fs :: [PGGPolGrp]) -> PGGPolGrp
  {Product.}
  g := New(PGGPolGrp);
  g`parent := Universe(fs);
  g`factor_fields := &cat[f`factor_fields : f in fs];
  g`factorization := [*true, &cat[f`factorization[2] : f in fs], &cat[f`factorization[3] : f in fs]*];
  return g;
end intrinsic;

intrinsic '*'(f :: PGGPolGrp, g :: PGGPolGrp) -> PGGPolGrp
  {Multiply.}
  require Parent(f) eq Parent(g): "different parents";
  return &*[f,g];
end intrinsic;

intrinsic DefiningPolynomial(E :: PGGFldGrp, F :: PGGFldGrp) -> PGGPolGrp
  {The defining polynomial of E over F.}
  f := New(PGGPolGrp);
  f`parent := PolynomialRing(F);
  f`factor_fields := [E];
  f`factorization := [*true,[f],[rec<CERT | Extension:=E>]*];
  return f;
end intrinsic;

intrinsic DefiningPolynomial(F :: PGGFldGrp) -> PGGPolGrp
  {The defining polynomial of F over the base field.}
  if not assigned F`defining_polynomial then
    F`defining_polynomial := DefiningPolynomial(F, BaseField(F));
  end if;
  return F`defining_polynomial;
end intrinsic;

intrinsic AssignNames(~R :: PGGRngPolStd, names :: [MonStgElt])
  {Assigns names to generators of R.}
  AssignNames(~R`actual, names);
end intrinsic;