Factorization.mag

// This file is part of ExactpAdics
// Copyright (C) 2018 Christopher Doris
//
// ExactpAdics is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// ExactpAdics is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with ExactpAdics.  If not, see <http://www.gnu.org/licenses/>.


// Suppose we have f(x), g(x), Pi(x), Rho(x) such that g(x) approximates an irreducible factor of f(x) uniquely, and Pi(x) is a uniformizer of and Rho(x) generates the residue class field of the extension defined by the factor. Can we prove that g(x) is close to an irreducible factor? Can we (Hensel) lift it to be closer to the factor? Consider f mod g in the basis {Rho^i * Pi^j mod g : i<F, j<E}, we should be able to Hensel-lift this directly.

///# Univariate polynomials
///## Root-finding and factorization

import "ExactpAdics.mag": CAP_APR, CHANGE_APR, WZERO, WVAL, VAL, APR;
OO := Infinity();
Z := Integers();
declare verbose ExactpAdics_Factorization, 2;
declare verbose ExactpAdics_Roots, 2;

function reslope(f, slope : pivot:=Degree(f))
  deg := Degree(f);
  return Parent(f) ! [ShiftValuation(Coefficient(f, i), slope*(i-pivot)) : i in [0..deg]];
end function;

function reslope0(f, slope)
  return reslope(f, slope : pivot:=0);
end function;

function deslope(f)
  deg := Degree(f);
  if deg le 0 then
    return f, 0;
  end if;
  np := NewtonPolygon([<i, v> : i in [0..deg] | v lt Infinity() where v:=Valuation(Coefficient(f,i))] : Faces:="Lower");
  slope := Ceiling(Max(Slopes(np)));
  return reslope(f, -slope), slope;
end function;

function reshift(f, shift)
  return Parent(f) ! [ShiftValuation(c, shift) : c in Coefficients(f)];
end function;

function integralize(f)
  deg := Degree(f);
  if deg lt 0 then
    return f, 0;
  end if;
  val := Min([Valuation(c) : c in Coefficients(f)]);
  return reshift(f, -val), val;
end function;

function nicen(f)
  f2, slope := deslope(f);
  f3, val := integralize(f2);
  return f3, slope, val;
end function;

function wrap_factorization(xf : Certificates:=false, Extensions:=false, Ideals:=false)
  // it seems to be faster to compute extensions separately, so we do this
  Certificates or:= Extensions or Ideals;
  facs, lc, certs := Factorization(xf : IsSquarefree:=true, Certificates:=Certificates);
  if Extensions or Ideals then
    for i in [1..#facs] do
      facs2, lc2, certs2 := Factorization(facs[i][1] : IsSquarefree:=true, Certificates:=Certificates, Extensions:=Extensions, Ideals:=Ideals);
      assert #facs2 eq 1;
      assert facs2[1][2] eq 1;
      if Extensions then
        certs[i]`Extension := certs2[1]`Extension;
      end if;
      if Ideals then
        certs[i]`IdealGen1 := certs2[1]`IdealGen1;
        certs[i]`IdealGen2 := certs2[1]`IdealGen2;
      end if;
    end for;
  end if;
  if Certificates then
    return facs, lc, certs;
  else
    return facs, lc, _;
  end if;
end function;

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

/// The factorization of f as a sequence of `<factor, multiplicity>` pairs.
/// 
/// It is only possible to factorize squarefree polynomials, so `multiplicity` is always 1. The factors are monic; also returns the leading coefficient of `f`.
/// 
///param Strategy The precision strategy.
///param Alg:="OM" The algorithm to use, either "Builtin" to use Magma's bulitin algorithm or "OM" to use our implementation.
///param Certificates When `true`, also returns a corresponding sequence of certificates including fields `F`, `E`, `Rho` and `Pi`.
///param Extensions When `true` implies `Certificates` and also includes `Extension` in each certificate.
intrinsic Factorization(f :: RngUPolElt_FldPadExact : Strategy:="default", Alg:="OM", Certificates:=false, Extensions:=false) -> [], FldPadElt, []
  {The factorization of f.}
  case Alg:
  when "Builtin":
    return Builtin_Factorization(f : Strategy:=Strategy, Certificates:=Certificates, Extensions:=Extensions);
  when "OM":
    return ExactpAdics_Factorization(f : Strategy:=Strategy, Certificates:=Certificates, Extensions:=Extensions);
  else
    require false: "Alg must be Builtin or OM";
  end case;
end intrinsic;

intrinsic Roots(f :: RngUPolElt_FldPadExact : Strategy:="default", Alg:="OM") -> []
  {The roots of f.}
  case Alg:
  when "Builtin":
    return Builtin_Roots(f : Strategy:=Strategy);
  when "OM":
    return ExactpAdics_Roots(f : Strategy:=Strategy);
  else
    require false: "Alg must be Builtin or OM";
  end case;
end intrinsic;

/// The factorization of `f` according to its Newton polygon, as a sequence of factors.
///
///param Residual:=false When `true` then further factorizes according to the factorization of each residual polynomial.
intrinsic NewtonPolygonFactorization(f :: RngUPolElt_FldPadExact : Strategy:="default", Residual:=false) -> []
  {The factorization of `f` along its Newton polygon, that is, one factor per face of the Newton polygon. If `Residual` is true, then further factorizes according to the coprime factorization of each residual polynomial.}
  if not assigned f`npfactorization then
    f`npfactorization := AssociativeArray();
  end if;
  if not IsDefined(f`npfactorization, Residual) then
    R := Parent(f);
    K := BaseRing(f);
    facs := [R|];
    np := NewtonPolygon(f : Strategy:=Strategy, Support:=<1,WeakDegree(f)>);
    vs := ChangeUniverse(Vertices(np), car<Z,Z>);
    for i in [1..#vs] do
      if i eq 1 then
        if vs[1][1] eq 0 then
          continue;
        else
          assert vs[1][1] eq 1;
          v1 := vs[1][2];
          v0 := Valuation(Coefficient(f`approximation, 0));
          if v0 eq OO then
            fac := R![0,1];
          else
            ok, fac := IsHenselLiftable(f, R![0,1] : Strategy:=Strategy, Slope:=v1-v0, fShift:=v0, gShift:=v0-v1);
            assert ok;
          end if;
          Append(~facs, fac);
        end if;
      else
        x0,v0 := Explode(vs[i-1]);
        x1,v1 := Explode(vs[i]);
        w := x1 - x0;
        slope := (v1-v0) / w;
        h := Numerator(-slope);
        e := Denominator(-slope);
        ok, d := IsDivisibleBy(w, e);
        assert ok;
        resfld, resmap := ResidueClassField(K);
        respol := Polynomial([ShiftValuation(Coefficient(f, x0+i*e), i*h-v0)@resmap : i in [0..d]]);
        assert Degree(respol) eq d;
        assert Coefficient(respol, d) ne 0;
        assert Coefficient(respol, 0) ne 0;
        if Residual then
          resfacs := [x[1]^x[2] : x in Factorization(respol)];
        else
          resfacs := [respol / Coefficient(respol, d)];
        end if;
        for resfac in resfacs do
          d2 := Degree(resfac);
          xfac := R![K| ok select ShiftValuation(WeakApproximation(Coefficient(resfac, j)@@resmap), d2*h-j*h) else 0 where ok,j:=IsDivisibleBy(i,e) : i in [0..e*d2]];
          ok, fac := IsHenselLiftable(f, xfac : Strategy:=Strategy, Slope:=slope, fShift:=v0-slope*x0, gShift:=d2*h);
          assert ok;
          Append(~facs, fac);
        end for;
      end if;
    end for;
    f`npfactorization[Residual] := facs;
  end if;
  return f`npfactorization[Residual];
end intrinsic;

/// Called internally by `Factorization` with `Alg:="Builtin"`.
/// 
///param Strategy The precision strategy.
///param Certificates When `true`, also returns a corresponding sequence of certificates including fields `F`, `E`, `Rho` and `Pi`.
///param Extensions When `true` implies `Certificates` and also includes `Extension` in each certificate.
///param Ideals When `true` implies `Certificates` and also includes `IdealGen1` and `IdealGen2` in each certificate.
///param UseNP:=true When `true`, factorizes `f` first by its Newton polygon. This can be a significant performance improvement for large degree polynomials with several faces.
intrinsic Builtin_Factorization(f :: RngUPolElt_FldPadExact : Strategy:="default", Certificates:=false, Extensions:=false, Ideals:=false, UseNP:=true) -> [], FldPadElt, []
  {The factorization of f.}
  Certificates or:= Extensions or Ideals;
  // factor by Newton polygons first
  if UseNP then
    npfacs := NewtonPolygonFactorization(f : Strategy:=Strategy, Residual);
    facs := [];
    certs := [];
    for npfac in npfacs do
      facs1, _, certs1 := Builtin_Factorization(npfac : Strategy:=Strategy, Certificates:=Certificates, Extensions:=Extensions, Ideals:=Ideals, UseNP:=false);
      facs cat:= facs1;
      if Certificates then
        certs cat:= certs1;
      end if;
    end for;
    lc := LeadingCoefficient(f : Strategy:=Strategy);
    if Certificates then
      return facs, lc, certs;
    else
      return facs, lc, _;
    end if;
  end if;
  // approximate factorization
  degf := WeakDegree(f);
  R := Parent(f);
  K := BaseRing(R);
  ok, _, data := ExactpAdics_ExecutePrecisionStrategy(function (pr)
    vprint ExactpAdics_Factorization: "precision =", pr;
    xf0 := Approximation(f, BaselineValuation(f) + pr : FixPr);
    xR := Parent(xf0);
    xK := BaseRing(xR);
    // ensure full degree
    if WZERO(Coefficient(xf0, degf)) then
      return false, _;
    end if;
    if degf eq 1 then
      xf := xf0;
      slope := 0;
      shift := 0;
      xfacs := [<xf / Coefficient(xf, degf), 1>];
      xcerts := [rec<CERT | F:=1, E:=1, Pi:=xR!UniformizingElement(xK), Rho:=xR!1>];
      if Extensions then
        xcerts[1]`Extension := xK;
      end if;
      if Ideals then
        error "not implemented";
      end if;
    else
      xf, slope, shift := nicen(xf0);
      // ensure monic
      xf /:= Coefficient(xf, degf);
      // try to factorize
      try
        xfacs, _, xcerts := wrap_factorization(xf : Certificates:=true, Extensions:=Extensions, Ideals:=Ideals);
      catch err
        if _ExactpAdics_IsPrecisionError(err) then
          vprint ExactpAdics_Factorization: "precision error";
          return false, _;
        else
          // run it again without catching the error, so we can enter the debugger
          // (it would be nice if instead magma allowed us to just catch some errors and let others propagate)
          xfacs, _, xcerts := wrap_factorization(xf : Certificates:=true, Extensions:=Extensions, Ideals:=Ideals);
        end if;
      end try;
    end if;
    // the factorization should be squarefree
    if exists{fac : fac in xfacs | fac[2] ne 1} then
      vprint ExactpAdics_Factorization: "not squarefree";
      return false, _;
    end if;
    assert &+[Degree(fac[1]) : fac in xfacs] eq Degree(xf);
    // check if each factor is hensel liftable
    facs := [];
    for i in [1..#xfacs] do
      xfac := xfacs[i][1];
      deg := Degree(xfac);
      xfac /:= Coefficient(xfac, deg);
      xcofac := xf div xfac;
      codeg := degf - deg;
      xfmodfac := xf mod xfac;
      s := Min([Valuation(c) : c in Coefficients(xfmodfac)]);
      M := Matrix([[j lt i select 0 else Coefficient(xcofac, j-i) : j in [0..degf-1]] : i in [0..deg-1]] cat [[j lt i select 0 else Coefficient(xfac, j-i) : j in [0..degf-1]] : i in [0..codeg-1]]);
      J := Determinant(M);
      if WZERO(J) then
        vprint ExactpAdics_Factorization: "resultant weakly zero";
        return false, _;
      end if;
      t := VAL(J);
      if s le 2*t then
        vprint ExactpAdics_Factorization: "not hensel liftable";
        return false, _;
      end if;
      // now do the lift
      init := reslope(Parent(xfac) ! [CAP_APR(c, s-t) : c in Coefficients(xfac)], slope);
      mkupdate := func<z | function (apr)
        xfac0 := GetData(z);
        pr := &join(apr - WVAL(z)) join (2*t+1);
        return ExactpAdics_GeneralGetter(pr,
          procedure (~pr, ~dep)
            dep := Approximation_Lazy(f, BaselineValuation(f) + pr) mod func<xf | reshift(reslope(xf, -slope), -shift)>;
          end procedure,
          procedure (xf, ~pr, ~value)
            vprint ExactpAdics_Factorization, 2: "apr =", apr;
            vprint ExactpAdics_Factorization, 2: "pr =", pr;
            xR := Parent(xf);
            xK := BaseRing(xR);
            xpr := Max([VAL(c) lt OO select APR(c)-VAL(c) else 0 : c in Coefficients(xf)]);
            vprint ExactpAdics_Factorization, 2: "xpr =", xpr;
            incpr := func<pol | xR ! [WZERO(c) select 0 else CHANGE_APR(xK!c, VAL(c)+xpr) : c in Coefficients(pol)]>;
            xfac := incpr(xfac0);
            xcofac := incpr(xf div xfac);
            xerr := xf - xfac * xcofac;
            s := Min([Valuation(c) : c in Coefficients(xerr)]);
            vprint ExactpAdics_Factorization, 2: "s =", s;
            vprint ExactpAdics_Factorization, 2: "t =", t;
            assert s gt 2*t;
            while true do
              M := Matrix([[j lt i select 0 else Coefficient(xcofac, j-i) : j in [0..degf-1]] : i in [0..deg-1]] cat [[j lt i select 0 else Coefficient(xfac, j-i) : j in [0..degf-1]] : i in [0..codeg-1]]);
              J := Determinant(M);
              assert not WZERO(J);
              assert VAL(J) eq t;
              Minv := M^-1;
              delta := Vector([Coefficient(xerr, i) : i in [0..degf-1]]) * M^-1;
              new_xfac := incpr(xfac + xR!Eltseq(delta)[1..deg]);
              new_xcofac := incpr(xcofac + xR!Eltseq(delta)[deg+1..deg+codeg]);
              new_xerr := xf - new_xfac * new_xcofac;
              new_s := Min([Valuation(c) : c in Coefficients(new_xerr)]);
              vprint ExactpAdics_Factorization, 2: "new_s =", new_s;
              if new_s gt s then
                xfac := new_xfac;
                xcofac := new_xcofac;
                xerr := new_xerr;
                s := new_s;
              else
                xfac2 := reslope(xR![CAP_APR(c, new_s-t) : c in Coefficients(new_xfac)], slope);
                d := IntegerValue(&join[apr(i) - APR(Coefficient(xfac2, i)) : i in [0..deg]]);
                if d gt 0 then
                  pr +:= d;
                else
                  Update(z, xfac2);
                  SetData(z, xfac);
                  value := true;
                  vprint ExactpAdics_Factorization, 2: "success";
                end if;
                return;
              end if;
            end while;
          end procedure
        );
      end function>;
      fac := R ! <init, mkupdate, xfac>;
      Append(~facs, <fac, 1>);
    end for;
    // now make the certificates
    if Certificates then
      certs := [];
      for xcert in xcerts do
        cert := rec<CERT |
          E := xcert`E,
          F := xcert`F,
          Pi := WeakApproximation(R!reslope0(xcert`Pi, slope)),
          Rho := WeakApproximation(R!reslope0(xcert`Rho, slope))
        >;
        if Ideals then
          cert`IdealGen1 := WeakApproximation(K!xcert`IdealGen1);
          cert`IdealGen2 := WeakApproximation(R!reslope0(xcert`IdealGen2, slope));
        end if;
        if Extensions then
          try
            cert`Extension := ExactpAdicField(xcert`Extension, xK, K);
          catch err
            if Type(err`Object) eq MonStgElt and "does not define a unique extension" in err`Object then
              return false, _;
            else
              cert`Extension := ExactpAdicField(xcert`Extension, xK, K);
            end if;
          end try;
        end if;
        Append(~certs, cert);
      end for;
    else
      certs := false;
    end if;
    return true, [*facs, certs*];
  end function, Strategy);
  if not ok then
    error "precision error";
  end if;
  facs, certs := Explode(data);
  lc := WeakLeadingCoefficient(f);
  if Certificates then
    return facs, lc, certs;
  else
    return facs, lc, _;
  end if;
end intrinsic;

/// Called internally by `Roots` with `Alg:="Builtin"`.
/// 
///param Strategy The precision strategy.
///param UseNP:=true When `true`, factorizes `f` first by its Newton polygon.
intrinsic Builtin_Roots(f :: RngUPolElt_FldPadExact : Strategy:="default", UseNP:=true) -> []
  {The roots of f.}
  if not assigned f`roots then
    if UseNP then
      facs := NewtonPolygonFactorization(f : Strategy:=Strategy);
      f`roots := &cat[PowerSequence(car<BaseRing(f), Z>) | Builtin_Roots(fac : Strategy:=Strategy, UseNP:=false) : fac in facs];
    else
      R := Parent(f);
      K := BaseRing(R);
      df := WeakDegree(f);
      ok, _, roots := ExactpAdics_ExecutePrecisionStrategy(function (pr)
        vprint ExactpAdics_Roots: "precision =", pr;
        // get an approximation
        xf0 := Approximation(f, BaselineValuation(f) + pr : FixPr);
        // require that we can see the full degree of xf
        if WZERO(Coefficient(xf0, df)) then
          vprint ExactpAdics_Roots: "not full degree";
          return false, _;
        end if;
        if df eq 1 then
          xf := xf0;
          slope := 0;
          shift := 0;
          xroots := [<-Coefficient(xf,0)/Coefficient(xf,1), 1>];
        else
          xf, slope, shift := nicen(xf0);
          vprint ExactpAdics_Roots: "xf =", xf;
          // call Magma's root-finding routine
          try
            xroots := Roots(xf : IsSquarefree);
          catch err
            if _ExactpAdics_IsPrecisionError(err) then
              vprint ExactpAdics_Roots: "precision error";
              return false, _;
            else
              error err;
            end if;
          end try;
        end if;
        vprint ExactpAdics_Roots: "xroots =", xroots;
        vprint ExactpAdics_Roots: "D =", Discriminant(xf);
        // the roots should be distinct
        if exists{r : r in xroots | r[2] ne 1} then
          vprint ExactpAdics_Roots: "repeated roots";
          return false, _;
        end if;
        // check they are hensel-liftable
        roots := [];
        for xr in xroots do
          ok, r := IsHenselLiftable(f, WeakApproximation(K ! ShiftValuation(xr[1], -slope)) : Strategy:=Strategy);
          if ok then
            Append(~roots, <r, 1>);
          else
            vprint ExactpAdics_Roots: "not Hensel liftable";
            return false, _;
          end if;
        end for;
        vprint ExactpAdics_Roots: "success";
        return true, roots;
      end function, Strategy);
      if ok then
        f`roots := roots;
      else
        error "precision error";
      end if;
    end if;
  end if;
  return f`roots;
end intrinsic;