Factorization_new.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/>.


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

import "common/Factorization_new.mag": CERT, STATE, BRANCH, APR, factorize, irred_done_branches, add_flatten_errs, absdeg, mb_has_value, mb_seq, mb_approximate_factor, br_lift, approximate_factor, base_flatten, elt_with_valuation, elt_with_residue_class, absramdeg, absresdeg, mb_eisenstein_poly, resfield, not_implemented, eisensteinPolyDefinesUniqueExtension;

function state_make_exact(f)
  assert Type(f) eq RngUPolElt_FldPadExact;
  return rec<STATE |
    xpol_state := 1,
    xpol_get := function (xst)
      return true, EpochApproximation(f, xst), xst+1;
    end function,
    xpols := [**],
    brs_todo := [rec<BRANCH |
      on_pop := [* "next_xpol" *],
      xpol_index := 0,
      basis := [],
      offset := 0,
      slopes_done := {},
      resfield := ResidueClassField(BaseRing(f)),
      width := WeakDegree(f)
    >],
    brs_done := []
  >;
end function;

function exact_remove_zeros(f)
  n := 0;
  while n lt WeakDegree(f) and WeakValuation(Coefficient(f, n)) eq Infinity() do
    n +:= 1;
  end while;
  if n eq 0 then
    return f, n;
  else
    return Parent(f) ! Coefficients(f)[n+1..WeakDegree(f)+1], n;
  end if;
end function;

/// Called internally by `Factorization`.
/// 
///param Proof:=true When true, the output is proven.
///param Certificates:=false When true, implies `Proof:=true` and returns certificates.
///param Extensions:=false When true, implies `Certificates:=true` and certificates include extensions.
///param InternalData:=false When true, implies `Certificates:=true` and certificates include the state of the branch leading to this factor.
///param SimpleLift:=false When true, uses "single factor lifting" just enough to use simple Hensel lifting on the factors. Probably slower, but possibly more reliable.
///param DegreeMle:=0 Only returns factors whose degree divides this, i.e. the degree is "multiplicatively less than" this.
///param DegreeGe:=1 Only returns factors whose degree is at least this.
///param Limit A bound on the number of factors returned.
intrinsic ExactpAdics_Factorization(f :: RngUPolElt_FldPadExact : Proof:=true, Certificates:=false, Extensions:=false, InternalData:=false, SimpleLift:=false, DegreeMle:=0, DegreeGe:=1, Limit:=Infinity()) -> [], FldPadExactElt, []
  {Factorization of f.}
  require WeakDegree(f) ge 0: "f must be non zero";
  Certificates or:= Extensions;
  Certificates or:= InternalData;
  Proof or:= Certificates;
  R := Parent(f);
  K := BaseRing(R);
  f, nzeros := exact_remove_zeros(f);
  if nzeros gt 0 then
    zero_facs := [<R![0,1], nzeros>];
    zero_certs := [rec<CERT | F:=1, E:=1, Rho:=R!1, Pi:=R!UniformizingElement(K), Extension:=K>];
  else
    zero_facs := [];
    zero_certs := [];    
  end if;
  deg := Degree(f);
  st := state_make_exact(f);
  factorize(~st : DegreeMle:=DegreeMle, DegreeGe:=DegreeGe, Limit:=Limit);
  lc := Coefficient(f, deg);
  brs := irred_done_branches(st : TotalDegree:=deg);
  if Proof then
    for i in [1..#brs] do
      add_flatten_errs(~brs[i]);
    end for;
    if #brs eq 1 and absdeg(brs[1]) eq deg then
      facs := [<f / lc, 1>];
    elif SimpleLift then
      epoch := 0;
      for epoch in [1..99999] do
        // approximate f
        xf := EpochApproximation(f, epoch);
        // get approximate factors
        ok, xfacs := mb_has_value(mb_seq([mb_approximate_factor(br_lift(br, xf)) : br in brs]));
        if not ok then
          continue epoch;
        end if;
        // ensure they are pairwise distinct (this could be optimized by classifying factors by degree and slope and only checking pairs within classes)
        for i in [1..#xfacs] do
          for j in [i+1..#xfacs] do
            xfac1 := xfacs[i];
            xfac2 := xfacs[j];
            if Degree(xfac1) eq Degree(xfac2) and forall{i : i in [0..Degree(xfac1)] | IsWeaklyEqual(Coefficient(xfac1, i), Coefficient(xfac2, i))} then
              continue epoch;
            end if;
          end for;
        end for;
        // ensure they are Hensel liftable (again this could be optimized by remebering the lifts from earlier iterations of the strategy if all factors within a class are liftable)
        facs := [];
        for xfac in xfacs do
          d := Degree(xfac);
          slope := (Valuation(Coefficient(xfac, d)) - Valuation(Coefficient(xfac, 0)))/d;
          fShift := Min([Valuation(Coefficient(xf, i)) - i*slope : i in [0..Degree(xf)]]);
          gShift := Min([Valuation(Coefficient(xfac, i)) - i*slope : i in [0..Degree(xfac)]]);
          ok, fac := IsHenselLiftable(f, WeakApproximation(R!xfac) : Slope:=slope, fShift:=fShift, gShift:=gShift);
          if ok then
            Append(~facs, <fac, 1>);
          else
            continue epoch;
          end if;
        end for;
        break epoch;
      end for;
    else
      facs := [];
      for i in [1..#brs] do
        br := brs[i];
        xfac := approximate_factor(br);
        fac := New(RngUPolElt_FldPadExact);
        fac`parent := R;
        fac`dependencies := [* f *];
        fac`internal_data := br;
        fac`get_approximation := function (n, xds)
          locfac := fac;
          xf := xds[1];
          br := fac`internal_data;
          br2 := br_lift(br, xf);
          xfac := approximate_factor(br2);
          locfac`internal_data := br2;
          return xfac;
        end function;
        fac`min_epoch := br`xpol_index;
        Init(fac);
        Append(~facs, <fac, 1>);
      end for;
    end if;
  else
    not_implemented(:msg:="Proof:=false");
    facs := [];
  end if;
  if Certificates then
    certs := [rec<CERT
      | F:=absresdeg(br)
      , E:=absramdeg(br)
      , Rho:=WeakApproximation(R!base_flatten(br, elt_with_residue_class(br, NormalElement(resfield(br), br`resfield))))
      , Pi:=WeakApproximation(R!base_flatten(br, elt_with_valuation(br, 1/absramdeg(br))))
      , InternalData:=br
      >
      : br in brs
    ];
    if Extensions then
      // we compute extensions by factoring the polynomial again over the unramified extension of the right degree
      // TODO: this seems like the wrong thing to do, we should be able to alter the basis as if we had worked over the unramified extension from the beginning
      for i in [1..#certs] do
        if certs[i]`F eq 1 then
          if certs[i]`E eq 1 then
            certs[i]`Extension := K;
          else
            for epoch in [1..99999] do
              xf := EpochApproximation(f, epoch);
              br := br_lift(certs[i]`InternalData, xf);
              ok, epol := mb_has_value(mb_eisenstein_poly(br));
              if ok and forall{c : c in Coefficients(epol) | APR(c) gt 2} and eisensteinPolyDefinesUniqueExtension(epol) then
                L := TotallyRamifiedExtension(K, WeakApproximation(R ! epol));
                assert #Roots(ChangeRing(facs[i][1], L)) ne 0;
                break epoch;
              end if;
            end for;
            certs[i]`Extension := L;
          end if;
        else
          U := UnramifiedExtension(K, certs[i]`F);
          ufacs, _, ucerts := ExactpAdics_Factorization(ChangeRing(facs[i][1], U) : Extensions, Limit:=1);
          assert #ufacs eq 1 /*certs[i]`F*/;
          assert forall{c : c in ucerts | c`F eq 1 and c`E eq certs[i]`E};
          certs[i]`Extension := ucerts[1]`Extension;
        end if;
        assert Degree(certs[i]`Extension, K) eq Degree(facs[i][1]);
        assert InertiaDegree(certs[i]`Extension, K) eq certs[i]`F;
        assert RamificationDegree(certs[i]`Extension, K) eq certs[i]`E;
      end for;
    end if;
    if not InternalData then
      for i in [1..#certs] do
        delete certs[i]`InternalData;
      end for;
    end if;
    return zero_facs cat facs, lc, zero_certs cat certs;
  else
    return zero_facs cat facs, lc, _;
  end if;
end intrinsic;

/// Called internally by `Roots`.
intrinsic ExactpAdics_Roots(f :: RngUPolElt_FldPadExact : Limit:=Infinity()) -> []
  {The roots of f.}
  vprint ExactpAdics_Factorization: Cputime(), "start";
  R := Parent(f);
  K := BaseRing(R);
  f, nzeros := exact_remove_zeros(f);
  st := state_make_exact(f);
  factorize(~st : JustRoots, Limit:=Limit);
  vprint ExactpAdics_Factorization: Cputime(), "check";
  brs := irred_done_branches(st);
  assert forall{br : br in brs | (absdeg(br) eq 1) and (#br`basis eq 0)};
  vprint ExactpAdics_Factorization: Cputime(), "hensel lift";
  rs := (nzeros eq 0 select [] else [<K!0,nzeros>]) cat [<r, 1> where _,r:=IsHenselLiftable(f, WeakApproximation(K!br`offset)) : br in brs];
  vprint ExactpAdics_Factorization: Cputime(), "done";
  return rs;
end intrinsic;