/*

  Absolute factorization of bivariate polynomials
  for zero or "high" characteristic
 
  (c) G. Chèze, G. Lecerf, March 2005.

  Guillaume.Cheze (at) math.unice.fr
  Gregoire.Lecerf (at) math.uvsq.fr

  Require Magma >= V2.10-4 (untested with older version)

Manual.

F is a polynomial in K[x,y] of total degree d, where K is a field
of characteristic 0 or at least d(d-1)+1.

intrinsic 
AbsoluteFactorization(F::RngMPolElt) -> [* <RngMPolElt, RngIntElt> *], RngElt

Returns the absolutely irreducible factors together with their
multiplicities and the corresponding scalar.

The internal algorithms are randomized but the output is always
correct. 

The sparse case is handled by Gao's algorithm.

Verbose flag: AbsoluteFactorization in [0..5]

Reference.

G. Chèze and G. Lecerf, "Lifting and Recombination Techniques
  for Absolute Factorization", manuscript, 2005.
*/

declare verbose AbsoluteFactorization,5;

/***********************************************************************
  PART I.

  Multi-factor Hensel lifting
 
  We follow Shoup's lifting algorith as described in the first
  edition of "Modern Computer Algebra" by von zur Gather and Gerhard.

**********************************************************************/ 

hensel_step:= procedure(f,~g,~h,~s,~t)
/*
  precision is the one of f;
  m=(maximal ideal)^(precision/2 +1);
  
  f = g*h mod m, s*g+t*h=1 mod m
  h is monic, deg(f)=n=deg(g)+deg(h), deg(s) < deg(h) and deg(t)<deg(g)
  
  output:
  f = G*H mod m^2, S*G+T*H=1 mod m^2
  H is monic, G=g mod m, H=h mod m, S=s mod m, T=t mod m,
  deg(G)=deg(g), deg(H)=deg(h), deg(S) < deg(H) and deg(T)<deg(G)
*/

  trim:= function(z,d)
    return Parent(z)![Coefficient(z,i) : i in [0..d]];
  end function;
  
  d:=Degree(f);
  dg:=Degree(g);
  dh:=Degree(h);
  e:=trim(f-g*h,d);

  q,r:= Quotrem(s*e,h);
  r:=trim(r,dh);
  G:=trim(g+t*e+q*g,dg);
  H:=trim(h+r,dh);
  
  b:=trim(s*G+t*H-1,d);
  c,d:= Quotrem(s*b,H);
  S:=trim(s-d,dh-1);
  T:=trim(t-t*b-c*G,dg-1);

  g:=G;
  h:=H;
  s:=S;
  t:=T;
end procedure;

/*
  The hensel step is to be applied to the subproduct tree enriched
  with the Bézout identities.
*/

node:=recformat< f, g, h, s, t >;

is_leaf:= function(n)
  return not assigned n`g;
end function;

get_leaves:=function(tree)
  if is_leaf(tree) then return [tree`f]; end if;
  return $$(tree`g) cat $$(tree`h);
end function;

build_tree:= function(f)
/*
  f is a list of factors.
*/
  n:=#f;
  assert n ge 1;
  if n eq 1 then
    return rec<node | f:=f[1]>;
  end if;
  n1:= n div 2;
  t1:=$$(f[1..n1]);
  t2:=$$(f[n1+1..n]);
  
  d,c1,c2:=XGCD(t1`f,t2`f);
  assert d eq 1;
  return rec<node|
    f:= t1`f * t2`f,
    g:= t1,
    h:= t2,
    s:= c1,
    t:= c2
  >;
end function;

change_base_ring_tree:= procedure(~node,R)
  node`f:=R!node`f;
  if is_leaf(node) then
    return;
  end if;
  node`s:=R!node`s;  
  node`t:=R!node`t;
  $$(~node`g,R);
  $$(~node`h,R);
end procedure;

map_tree:= procedure(~node,h)
  node`f:=h(node`f);
  if is_leaf(node) then
    return;
  end if;
  node`s:=h(node`s);
  node`t:=h(node`t);
  $$(~node`g,h);
  $$(~node`h,h);
end procedure;

lift_tree:=procedure(~node)
  if is_leaf(node) then return; end if;
  hensel_step(node`f,~node`g`f,~node`h`f,~node`s,~node`t);
  $$(~node`g);
  $$(~node`h);
end procedure;

multi_factor_hensel_lifting:= function(lf,F,eta)
/*
  Lift the factorization F= &*lf at precision eta.
  F and the elements of lf must be monic.
*/
  UPy:= Parent(F);
  UPx:= BaseRing(UPy);
  K:= BaseRing(UPx);

  time_begin := Cputime();
  
  tree:=build_tree(lf);

  LS<t>:=LaurentSeriesRing(K,1);
  UP<T>:=PolynomialRing(LS);

  change_base_ring_tree(~tree,UP);
  tree`f:= UP![ LS!z : z in Coefficients(F)];

  // compute the sequence of successive precisions
  k:= eta;
  precisions:=[];
  while k gt 1 do
    Append(~precisions,k);
    if k mod 2 eq 0 then k:= (k div 2); else k:= (k div 2) +1; end if;
  end while;

  // main loop
  while not IsEmpty(precisions) do
    precision:=precisions[#precisions];
    precisions:=precisions[1..#precisions-1];
    vprintf AbsoluteFactorization, 2: "multi_factor_hensel_lifting, entering precision %o\n",precision;
    LS<t>:= LaurentSeriesRing(K,precision);
    UP<T>:= PolynomialRing(LS);

    change_precision:= function(z)
      ans:=[LS|];
      for i:=0 to Degree(z) do
        tmp:=Coefficient(z,i);
        ChangePrecision(~tmp,precision);
        Append(~ans,tmp);
      end for;
      return UP!ans;
    end function;

    map_tree(~tree,change_precision);
    tree`f:= UP![ LS!z : z in Coefficients(F)];
    lift_tree(~tree);
  end while;

  return [ UPy![ UPx![K|Coefficient(z,i): i in [0..eta-1]] : z in Coefficients(f)] : f in get_leaves(tree)];
end function;


/***********************************************************************
  Part II.

  In this part, F \in K[x,y] is assumed to satisfy:
  H_1) monic in y.
  H_2) Resultant_y(F, Derivative(F,y))(0) <> 0.

  Notation:
    d denotes the total degree of F.
    A<a> := K[y]/F(0,y)
    phi in A[[x]] is defined by phi(0)=a and F(x,phi(x))=0.
***********************************************************************/

/*
 * Computation of phi(x)
 */

evaluate_ps_phi:= function(P, psi)
/*
  Paterson & Stockmeyer's evaluation scheme, optimized for the lifting of phi.
  psi is an approximation of phi.
*/
  n:= Degree(P);
  k:= Iroot(n+1,2);
  kp:= Ceiling((n+1)/k);

  S:= Parent(psi);
  A:= BaseRing(S);
  d:= Degree(Modulus(A));
  sigma:= S`Precision;
  K:= BaseRing(A);
  R:= PowerSeriesRing(K, sigma);
  
  tmp:= [S| i eq 0 select S!1 else psi*Self(i) : i in [0..k-1]];
  M:= Matrix([ [ R![Coefficient(Coefficient(tmp[j+1],z),i)
      : z in [0..sigma-1]]
         : j in [0..k-1]]
	    : i in [0..d-1]]);
  N:= Matrix([ [ R!Coefficients(Coefficient(P,k*j +i))
      : j in [0..kp-1]]
         : i in [0..k-1]]);
  C:= M * N;
  delete M;
  delete N;
  D:= [ &+[S| S![A| Coefficient(C[j+1,i],z)
      : z in [0..sigma-1] ] * (A.1)^j
         : j in [0..d-1]]
	    : i in [1..kp]];
  tmp:= psi*tmp[k];
  tmp:= [S| i eq 0 select S!1 else tmp*Self(i) : i in [0..kp-1]];
  return &+[S| D[i] * tmp[i] : i in [1..kp]];
end function;

lifting_initialization_phi:= function(F)
/*
  Returns phi and D_y F(x,phi)^{-1} at precision (x).
*/
  Fp:= Derivative(F);
  UPy:= Parent(F);
  UPx:= BaseRing(F);
  K:= BaseRing(UPx);
  d:= Degree(F);
  UP:= PolynomialRing(K);
  f:= UP![K| Coefficient(z,0): z in Coefficients(F)];
  fp:= Derivative(f);
  r, ifp:= XGCD(fp,f);
  ifp:= ifp;
  A<a>:= quo<UP|f>;
  E:= PowerSeriesRing(A,1);
  phi:= E!a;
  ifp:= E!A!ifp;
  return phi, ifp;
end function;

lifting_step_phi:= procedure(~phi, ~ifp, F, Fp, sigma)
/*
  Updates phi and ifp = D_y F(x,phi)^{-1} up to precision sigma.
  The input precision is assumed (and not checked) to be sufficient.
*/
  vprintf AbsoluteFactorization, 2: "lifting_step_phi: going to precision %o\n", sigma;
  begin_time:= Cputime();
  E:= Parent(phi);
  A:= BaseRing(E);
  ifp:= E!Coefficients(ifp);
  t1:= evaluate_ps_phi(Fp, phi);
  ifp := ifp + ifp * (1 - ifp * t1);
  E:= PowerSeriesRing(A,sigma);
  phi:=E!Coefficients(phi);
  t2:=evaluate_ps_phi(F, phi);
  vprintf AbsoluteFactorization, 3: "lifting_step_phi: valuation = %o\n", Valuation(t2);
  phi:= phi - t2 * E!Coefficients(ifp);
  vprintf AbsoluteFactorization, 4: "lifting_step_phi: time = %o\n", Cputime()-begin_time;
end procedure;

/*
 * Linear system solving
 */

probabilistic_solve_L:= function(F, Fp, phi, a, b);
  vprintf AbsoluteFactorization, 2: "probabilistic_solve_L, entering probabilistic resolution of the linear system, a = %o, b = %o\n",a,b;
  begin_time:= Cputime();
  d:= Degree(F);
  UPy:= Parent(F);
  UPx:= BaseRing(UPy);
  K:= BaseRing(UPx);

  LS:= Parent(phi);
  A:= BaseRing(LS);
  tau:= LS`Precision;
  M:=[];
  for u in {a,b} do 
    if Coefficient(Coefficient(F,0),0) ne 0 then
      N:= Evaluate(F,u*UPx.1);
      N:= LS!ChangeUniverse(Coefficients(N), A);
      D:= (A!u)*LS.1 - phi;
      T:= N * D^(-1);
    else
	K:= BaseRing(UPx);
        LS1:= LaurentSeriesRing(K, tau);
	phi1:= LS1![K| Coefficient(z,0) : z in Coefficients(phi)];
	UP1:= PolynomialRing(LS1);
	N1:= UP1!ChangeUniverse(Coefficients(F), LS1);
	D1:= UP1.1 - phi1;
	T1:= N1 div D1;
	T1:= Evaluate(T1, u*LS.1);

	f:= UPx![K| Coefficient(z,0) : z in Coefficients(F)];
	fh:= f div UPx.1;
	A2:= quo<UPx |fh>;
	LS2:= LaurentSeriesRing(A2, tau);

        N2:= Evaluate(F,u*UPx.1);
	N2:= LS2!ChangeUniverse(Coefficients(N2), A2);
	D2:= u * LS2.1 - LS2![A2| A2!((UPx!z) mod fh) : z in Coefficients(phi)];
	T2:= N2 * D2^(-1);
	iy:= (A2.1)^(-1);
	crt:= function(b1, b2)
            p1:= UPx!b1;
	    p2:= UPx!b2;
	    p3:= UPx!(A2!(p2 - p1)*iy);
	    return p1 + UPx.1 * p3;
	end function;
	T:= LS![A| A!crt(Coefficient(T1,i),Coefficient(T2,i)): i in [0..tau-1]];
    end if;
    Tx:= T * Derivative(phi);
    for j:=d to tau-2 do
      Append(~M,[ Coefficient(Coefficient(T,j),l): l in [0..d-1]]);
      Append(~M,[ Coefficient(Coefficient(Tx,j),l): l in [0..d-1]]);
    end for;
  end for;
  vprintf AbsoluteFactorization, 2: "probabilistic_solve_L, size of the matrix: %o x %o\n", #M, d;
  M:= Matrix(M);
  L:= Kernel(Transpose(M));
  vprintf AbsoluteFactorization, 5: "probabilistic_solve_L, L = %o\n", L;
  vprintf AbsoluteFactorization, 4: "probabilistic_solve_L, total time: %o\n", Cputime() - begin_time;
  return L;
end function;

/*
 * Partial fraction decomposition
 */

characteristic_polynomial_of_the_residues:=function(g,f)
  UP:= Parent(f);
  K:= BaseRing(UP);
  d:= Degree(f);
  fp:= Derivative(f);
  floc:= f/ LeadingCoefficient(f);
  return Interpolation([K|j : j in [0..d]], [ Resultant(floc, j*fp - g): j in [0..d]]);
end function;

univariate_partial_fraction_decomposition:= function(g,f)
  UP:= Parent(f);
  fp:= Derivative(f);
  Q:= characteristic_polynomial_of_the_residues(g,f);
  tmp, unit:= Factorization(Q);
  ans:=[* *];
  for z in tmp do
    q:= z[1];
    FF<a>:= quo<UP | q>;
    UPa:= PolynomialRing(FF);
    fA:= UP!Coefficients(f);
    fpA:= UPa!Coefficients(fp);
    g:= UPa!Coefficients(g);
    ft:= GCD(fA, a*fpA - g);
    Append(~ans,ft);
  end for;
  return ans;
end function;

/*
 * Lifting of Ft from ft: ft is a factor of f = F(0,y) 
 * with coefficients in an algebraic extension FF of K of degree r.
 */

evaluate_ps_Ft:= function(P, psi)
/*
  Paterson & Stockmeyer's evaluation scheme
*/
  n:= Degree(P);
  k:= Iroot(n+1,2);
  kp:= Ceiling((n+1)/k);
  A:= Parent(psi);
  UP:= PreimageRing(A);
  S:= BaseRing(UP);
  R:= BaseRing(S);
  d:= Degree(Modulus(A));
  tmp:= [A| i eq 0 select A!1 else psi*Self(i) : i in [0..k-1]];
  M:= Matrix([ [ S!Coefficient(tmp[j+1],i)
       : j in [0..k-1]]
          : i in [0..d-1]]);
  N:= Matrix([ [ S!ChangeUniverse(Coefficients(Coefficient(P,k*j +i)), R)
       : j in [0..kp-1]]
          : i in [0..k-1]]);
  C:= M * N;
  D:= [A| &+[A| C[j+1,i] * (A.1)^j: j in [0..d-1]] : i in [1..kp]];
  tmp:= psi*tmp[k];
  tmp:= [A| i eq 0 select A!1 else tmp*Self(i) : i in [0..kp-1]];
  return &+[A| D[i] * tmp[i] : i in [1..kp]];
end function;

lifting_initialization_Ft:= function(F,ft)
  Fp:= Derivative(F);
  UPy:= Parent(F);
  UPx:= BaseRing(F);
  UP:= Parent(ft);
  FF:= BaseRing(UP);
  d:= Degree(F);
  f:= UP![FF| FF!Coefficient(z,0): z in Coefficients(F)];
  fp:= Derivative(f);
  r, ifp:= XGCD(fp,ft);
  ifp:= ifp;
  S:= LaurentSeriesRing(FF, 1);
  E:= PolynomialRing(S);
  d:= Degree(ft);
  ftloc:= E![S| S!z : z in Coefficients(ft)];
  ifploc:= E![S| S!z : z in Coefficients(ifp)];
  B:= quo<E|ftloc>;
  return B.1, B!ifploc, Fp;
end function;

change_precision:= function(z,S)
  w:=z;
  ChangePrecision(~w,S`Precision);
  return S!w;
end function;

lifting_step_Ft:= procedure(~psi, ~iFp, F, Fp, sigma)
  vprintf AbsoluteFactorization, 2: "lifting_step_Ft: going to precision %o\n", sigma;
  begin_time := Cputime();
  E:= Parent(psi);
  UP:= PreimageRing(E);
  S:= BaseRing(UP);
  A:= BaseRing(S);

  iFp:= E!UP![S| change_precision(z,S) : z in Coefficients(iFp)];
  t1:= evaluate_ps_Ft(Fp, psi);
  iFp := iFp + iFp * (1 - iFp * t1);
  R:= LaurentSeriesRing(A,sigma);
  UP:= PolynomialRing(R);
  q:= Modulus(E);
  q:= UP![R| change_precision(z,R) : z in Coefficients(q)];

  E:= quo<UP|q>;
  psi:= E!UP![R| change_precision(z,R) : z in Coefficients(psi)];
  t2:= evaluate_ps_Ft(F, psi);
  vprintf AbsoluteFactorization, 3: "lifting_step_Ft: valuation = %o\n", Min([Valuation(z): z in Coefficients(t2)]);
  tmp:= t2 * E![R| change_precision(z,R) : z in Coefficients(iFp)];
  tmp:= E!Derivative(q) * tmp;
  q:= q + UP!tmp;
  E:= quo<UP|q>;
  psi:= E.1;
  vprintf AbsoluteFactorization, 4: "lifting_step_Ft: time = %o\n", Cputime()-begin_time;
end procedure;

lift_Ft:= function(F,ft)
  begin_time:= Cputime();
  sigma:= Degree(ft)+1;
  psi, iFp, Fp:= lifting_initialization_Ft(F,ft);
  precision:= [ sigma ];
  while precision[#precision] gt 1 do
    Append(~precision, Ceiling(precision[#precision]/2));
  end while;
  for i:= #precision-1 to 1 by -1 do
    lifting_step_Ft(~psi,~iFp, F, Fp, precision[i]);
  end for;
  FF:= BaseRing(Parent(ft));
  UPx:= PolynomialRing(FF);
  UPy:= PolynomialRing(UPx);
  AssignNames(~UPx, [Sprint(BaseRing(Parent(F)).1)]);
  AssignNames(~UPy, [Sprint(Parent(F).1)]);
  Ft:= Modulus(Parent(psi));
  ans:= UPy![ UPx![FF| Coefficient(z,i) : i in [0..sigma-1]] : z in Coefficients(Ft)];
  vprintf AbsoluteFactorization, 4: "lift_Ft: total time = %o\n", Cputime()-begin_time;
  return ans;
end function;  

norm:= function(ft)
  UP:= Parent(ft);
  E:= BaseRing(ft);
  K:= BaseRing(E);
  q:= Modulus(E);
  d:= Degree(q) * Degree(ft);
  return Interpolation([K| i : i in[0..d]], [K| Norm(Evaluate(ft,i)) : i in [0..d]]);
end function;

check_Ft_div_F:= function(F,Ft)
/*
  Tells if Ft divides F.
*/
  UPy:= Parent(Ft);
  UPx:= BaseRing(UPy);
  FF:= BaseRing(UPx);
  d:= Degree(F);
  S:= LaurentSeriesRing(FF, d+1);
  UP:= PolynomialRing(S);
  tmp:= UP![S| S!Coefficients(z): z in Coefficients(Ft)];
  E:= quo<UP | tmp>;
  tmp:=evaluate_ps_Ft(F,E.1);
  return Min([ Valuation(z) : z in Coefficients(tmp)]) ge d+1;
end function;

/*
 * Absolute multi factor hensel lifting and verification
 */

absolute_multi_factor_hensel_lifting_and_check:=function(ft, F)
/*
  Lift ft and verify that the lifted factors are true factors.
*/
  p:= [norm(z) : z in ft];
  eta:= 1 + Max([ Degree(z) : z in p]);
  P:= multi_factor_hensel_lifting(p,F,eta);
  Ft:=[* *];
  for z in P do
    d:= Degree(z);
    for i:=0 to d do
      if i + Degree(Coefficient(z,i)) gt d then
        return P,Ft,false;
      end if;
    end for;
  end for;
  if &*P ne F then
    return P,Ft,false;
  end if;
  for i:=1 to #ft do
    z:=lift_Ft(P[i],ft[i]);
    if not check_Ft_div_F(P[i],z) then
      return P,Ft,false;
    end if;
    Append(~Ft, z);
  end for;
  return P,Ft,true;
end function;

/*
 * Main routines
 */

vectors_of_bounded_weight:=function(k,r)
  if r eq 1 then return [ [1] ]; end if;
  ans:= [ ];
  for i := 0 to k do
    ans cat:= [ Append(z,i) : z in $$(k-i, r-1)];
  end for;
  return ans;
end function;

probabilistic_recombination:= function(F,L,cs);
  d:= Degree(F);
  UPy:= Parent(F);
  UPx:= BaseRing(UPy);
  K:= BaseRing(UPx);
  r:=Rank(L);
  vprintf AbsoluteFactorization, 2: "probabilistic_recombination, expected number of irreducible absolute factors: %o\n", r;
  if r eq 1 then return [*F*],true; end if;

  UP:= UPx;
  f:= UP![Coefficient(z,0) : z in Coefficients(F)];
  fp:= Derivative(f);
  A:= quo<UP|f>;
  UPA:= PolynomialRing(A);
  T:= UPA!ChangeUniverse(Coefficients(f),A) div (UPA.1 - A.1);

  for c in cs do
    if &and[ IsZero(z): z in c] then continue; end if;
    vprintf AbsoluteFactorization, 2: "probabilistic_recombination, trying c = %o\n", c;
    nu:= [ &+[ (L.j)[i] * c[j] : j in [1..r] ] : i in [1..d]];
    g:= UP![ &+[ nu[l+1]*Coefficient(Coefficient(T,i),l) : l in [0..d-1]] : i in [0..d-1]];
    ft:= univariate_partial_fraction_decomposition(g,f);
    if &+[ Degree(Modulus(BaseRing(Parent(z)))) : z in ft] lt r then
      continue;
    end if;
    P, Ft, t:= absolute_multi_factor_hensel_lifting_and_check(ft, F);
    delete ft;
    if t then
      return Ft, true;
    end if;
  end for;
  return [* *], false;
end function;

absolute_factorization_normalized:= function(F)
/*
  call the previous function with different values of the parameters b and c
  until the correct factorization is found.
*/
  UPy:= Parent(F);
  UPx:= BaseRing(F);
  K:= BaseRing(UPx);
  d:= Degree(F);
  tau:= 2*d+1;
  
  // sequence of the successive precisions
  precision:= [ tau ];
  while precision[#precision] gt 1 do
    Append(~precision, Ceiling( precision[#precision]/2));
  end while;

  // lifting
  phi, ifp := lifting_initialization_phi(F);
  Fp:= Derivative(F);
  lifting_begin_time:= Cputime();
  for i:= #precision-1 to 1 by -1 do
    lifting_step_phi(~phi,~ifp, F, Fp, precision[i]);
  end for;

  // recombination
  vprintf AbsoluteFactorization,4: "probabilistic_absolute_factorization, total time for lifting phi up to precision %o: %o\n", tau, Cputime() - lifting_begin_time; 

  // Heuristic
  for b:=0 to d*(d-1) do
    L:= probabilistic_solve_L(F, Fp, phi, 0, b);
    r:= Rank(L);
    c:= [0 : i in [1..r]];
    for i:=1 to r do
      tmp:=c;
      tmp[i]:=1;
      af,t:= probabilistic_recombination(F,L, [ tmp ] );
      if t then return af; end if;
    end for;
    af,t:= probabilistic_recombination(F,L, [ [1..r] ] );
    if t then return af; end if;
  end for;

  // Worst case. TODO << improve
  for b:=0 to d*(d-1) do
    L:= probabilistic_solve_L(F, Fp, phi, 0, b);
    r:= Rank(L);
    af,t:= probabilistic_recombination(F,L, vectors_of_bounded_weight((r*(r-1)) div 2, r));
    if t then return af; end if;
  end for;
  error "absolute_factorization_normalized: unexpected case";
end function;

/*
 * Find good coordinates
 */

/*************************************************************************
From now on, F is not assumed to satisfy H_1 and H_2
*************************************************************************/

// Sparse case: here we implement Gao's algorithm

gao_matrix:=function(F)
/*
  Returns the matrix of D_x(G/F) = D_y(H/F) and
  the supports of G and H.
*/

  vprintf AbsoluteFactorization,2: "gao_matrix, entering\n";
  begin_time:= Cputime();

  MP:= Parent(F);
  K:= BaseRing(MP);

  d_x:= Degree(F,1);
  d_y:= Degree(F,2);

  S_F:= [ < Degree(m,1), Degree(m,2)> : m in Monomials(F)];
  C_F:= Coefficients(F);

  P:= NewtonPolygon(S_F);

  S_P:= {@ @}; // Support of P
  if IsHomogeneous(F) then
    // This case should be optimized
    // but I think that the Magma functions
    // should be fixed to behave well with degenerated
    // polygons
    for i:=0 to d_x do
      for j:=0 to d_y do
        if IsBoundary(P,<i,j>) then
          Include(~S_P, Monomial(MP,[i,j]));
        end if;
      end for;
    end for;
  else
    for i:=0 to d_x do
      for j:=0 to d_y do
        if IsPoint(P,<i,j>) then
          Include(~S_P, Monomial(MP,[i,j]));
        end if;
      end for;
    end for;
  end if;

  S_G:= {@ Exponents(Derivative(m,2)) : m in S_P | Derivative(m,2) ne 0 @};
  S_H:= {@ Exponents(Derivative(m,1)) : m in S_P | Derivative(m,1) ne 0 @};
  delete S_P;

  // Target support
  S_T:= {@ @};
  for m1 in S_G do
    for m2 in S_H do
      Include(~S_T, < m1[1] + m2[1], m1[2]+ m2[2]>);
    end for;
  end for;

  s_G:= #S_G;
  s_H:= #S_H;

  M:= SparseMatrix(K, s_G+s_H, #S_T); // Seems faster than dense here.

  time_10pc:= Cputime();
  for i:=1 to #S_F do
    if i eq (#S_F div 10)  then
      vprintf AbsoluteFactorization, 4: "10 percents of the matrix filled in %o seconds\n", Cputime()-time_10pc;
    end if;
    for j:=1 to s_G do
      l:=Index(S_T, < S_G[j][1]+S_F[i][1]-1, S_G[j][2]+S_F[i][2] >);
      if l gt 0 then
        M[j,l]:= S_G[j][1]*C_F[i]-S_F[i][1]*C_F[i];
      end if;
    end for;

    for j:=1 to s_H do
      l:=Index(S_T, < S_H[j][1]+S_F[i][1], S_H[j][2]+S_F[i][2]-1 >);
      if l gt 0 then
        M[s_G+j,l]:= S_F[i][2]*C_F[i]-S_H[j][2]*C_F[i];
      end if;
    end for;
  end for;

  vprintf AbsoluteFactorization,4: "gao_matrix, matrix built in %o seconds\n",Cputime()-begin_time;
  return M, S_G, S_H;
end function;

absolute_factorization_squarefree_sparse_gao:=function(F)
/*
  F is a non-constant squarefree polynomial in K(x)[y].
*/
  vprintf AbsoluteFactorization,2: "absolute_factorization_squarefree_sparse_gao, entering\n";
  MP:=Parent(F);
  d:= TotalDegree(F);
  K:= BaseRing(MP);
  M, S_G, _ := gao_matrix(F);
  vprintf AbsoluteFactorization,3: "absolute_factorization_squarefree_sparse_gao, solving linear system\n";
  ns_time:= Cputime();
  E:= Nullspace(M);
  ns_time:= Cputime() - ns_time;
  r:= Rank(E);
  vprintf AbsoluteFactorization,4: "absolute_factorization_squarefree_sparse_gao, system solved in %o seconds\n", ns_time;
  vprintf AbsoluteFactorization,4: "absolute_factorization_squarefree_sparse_gao, number of factors: %o\n", r;
  vprintf AbsoluteFactorization,3: "absolute_factorization_squarefree_sparse_gao, building G_1,...,G_r\n";
  UP:= PolynomialRing(K);
  G_basis:= [MP| &+[MP| (E.j)[i] * MP.1^S_G[i][1] * MP.2^S_G[i][2]  :  i in [1..#S_G]] : j in [1..r]];

  vprintf AbsoluteFactorization,3: "absolute_factorization_squarefree_sparse_gao, looking for a primitive element\n";

  // Heuristic
  for b:=0 to 3 do
    f:= UnivariatePolynomial(Evaluate(F,MP.1,K!b));
    to_test:= [ [j eq i select 1 else 0 : j in [1..r]] : i in [1..r]]
      cat [ [1 : i in [1..r]], [i : i in [1..r]]]; 
    for c in to_test do
      G:= &+[c[i]*G_basis[i]: i in [1..r]];
      g:= UnivariatePolynomial(Evaluate(G,MP.1,K!b));
      Q:= characteristic_polynomial_of_the_residues(g, f);
      q:= SquarefreePart(Q);
      if Degree(q) eq r then
        vprintf AbsoluteFactorization,3: "absolute_factorization_squarefree_sparse_gao, primitive element found: %o\n",c;
        break;
      end if;
    end for;
    if Degree(q) eq r then break; end if;
  end for;

  // Worst case. TODO << improve
  if Degree(q) ne r then 
    for b:=0 to d*(d-1) do
      f:= UnivariatePolynomial(Evaluate(F,MP.1,K!b));
      for c in vectors_of_bounded_weight((r*(r-1)) div 2, r) do
        G:= &+[c[i]*G_basis[i]: i in [1..r]];
        g:= UnivariatePolynomial(Evaluate(G,MP.1,K!b));
        Q:= characteristic_polynomial_of_the_residues(g, f);
        q:= SquarefreePart(Q);
        if Degree(q) eq r then
          vprintf AbsoluteFactorization,3: "absolute_factorization_squarefree_sparse_gao, primitive element found: %o\n",c;
          break;
        end if;
      end for;
    if Degree(q) eq r then break; end if;
    end for;
  end if;

  if Degree(q) ne r then
    error "absolute_factorization_normalized: unexpected case";
  end if;

  vprintf AbsoluteFactorization,3: "absolute_factorization_squarefree_sparse_gao, factoring q\n";
  facts_q:= Factorization(q);
  vprintf AbsoluteFactorization,3: "absolute_factorization_squarefree_sparse_gao, computing the factors\n";
  ans:=[* *];

  Fy:=Derivative(F,2);
  for q in facts_q do
    E<a>:= quo<UP|q[1]>;
    MP_E:= PolynomialRing(E,2);
    AssignNames(~MP_E, [Sprint(MP.1), Sprint(MP.2)]);
    Append(~ans, Gcd(MP_E!F, a*MP_E!Fy - MP_E!G));
  end for;
  return ans;
end function;

// Dense case: need to change the variables

total_degree:= function(F);
  return Max([Degree(Coefficient(F,i)) + i : i in [0..Degree(F)]]);
end function;

find_projective_noether_position:= function(F)
/*
  find a correct change of variables x <- x + u*y
*/
  vprintf AbsoluteFactorization, 2 : "find_projective_noether_position, entering\n";
  UPy:= Parent(F);
  UPx:= BaseRing(UPy);
  K:= BaseRing(UPx);
  d:= total_degree(F);
  if Degree(Coefficient(F,d)) eq 0 then
     u:=K!0;
     vprintf AbsoluteFactorization, 2 : "find_projective_noether_position: found u = %o\n", u;
     return u, F; end if;
  for u := 1 to d do
    G:= &+[UPy| Evaluate(Coefficient(F,i), UPx.1 + u*UPy.1)* UPy.1^i : i in [0..Degree(F)]];
    if Degree(Coefficient(G,d)) eq 0 then
       vprintf AbsoluteFactorization, 2 : "find_projective_noether_position: found u = %o\n", u;
       return K!u, G;
    end if;
  end for;
  error "find_strong_noether_position, BUG: unexpected case";
end function;

find_smooth_fiber:= function(F)
/*
  find a correct specialization x <- v
  F is assumed to be in projective Noether position.
*/
  vprintf AbsoluteFactorization, 2 : "find_smooth_fiber, entering\n";
  UPy:= Parent(F);
  UPx:= BaseRing(UPy);
  K:= BaseRing(UPx);
  d:= Degree(F);
  for v:=0 to d*(d-1) do
    tmp := UPx![ Evaluate(z,v): z in Coefficients(F)];
    if not IsZero(Discriminant(tmp)) then
      G := UPy![ Evaluate(z, v + UPx.1): z in Coefficients(F)];
      vprintf AbsoluteFactorization, 2 : "find_smooth_fiber, found v = %o\n", v;
      return v, G;
    end if;
  end for;
  error "find_smooth_fiber, BUG: unexpected case";
end function;

find_good_coordinates:= function(F)
  u, G:= find_projective_noether_position(F);
  v, G:= find_smooth_fiber(G);
  return u,v,G;
end function;

change_coordinates_back:= function(G,u,v)
/*
  Changes x <- x - uy -v
*/
  UPy:= Parent(G);
  UPx:= BaseRing(UPy);
  K:= BaseRing(UPx);
  if not IsZero(v) then
    F := UPy![ Evaluate(z, UPx.1 - K!v): z in Coefficients(G)];
  else
    F := G;
  end if;
  if not IsZero(u) then
    F:= &+[Evaluate(Coefficient(F,i), UPx.1 - (K!u)*UPy.1)* UPy.1^i : i in [0..Degree(F)]];
  end if;
  return F;
end function;

/*
 * Main sub-routines
 */

absolute_factorization_squarefree_dense:= function(G)
  vprintf AbsoluteFactorization,1: "absolute_factorization_squarefree_dense, entering\n";
  // coercion to K[x][y]
  tmp:= [ UnivariatePolynomial(z) : z in Coefficients(G,2)];
  UPx<x>:= Universe(tmp);
  UPy<y>:= PolynomialRing(UPx);
  F:= UPy!tmp;
  delete tmp;
  u,v,F:= find_good_coordinates(F);
  unit:= LeadingCoefficient(LeadingCoefficient(F));
  F/:= unit;
  ans:= [* *];
  for z in absolute_factorization_normalized(F) do
    tmp:= change_coordinates_back(z,u,v);
    MP:= PolynomialRing(BaseRing(BaseRing(tmp)),2);
    AssignNames(~MP, [Sprint(Parent(G).1), Sprint(Parent(G).2)]);
    tmp:= &+[MP| MultivariatePolynomial(MP, Coefficient(tmp,i), 1) * MP.2^i : i in [0..Degree(tmp)]];
    Append(~ans, tmp);
  end for;
  return ans, unit;
end function;

absolute_factorization_squarefree_sparse:= function(G)
  vprintf AbsoluteFactorization,1: "absolute_factorization_squarefree_sparse, entering\n";
  if not IsZero(Gcd(G,Derivative(G,2))) then
    ans:= absolute_factorization_squarefree_sparse_gao(G);
  else
    vprintf AbsoluteFactorization,1: "absolute_factorization_squarefree_sparse, swap variables\n";
    swap:= function(z)
      MP:= Parent(z);
      return hom<MP->MP | MP.2, MP.1>(z);
    end function;
    tmp:= absolute_factorization_squarefree_sparse_gao(swap(G));
    ans:=[* *];
    for z in tmp do Append(~ans, swap(z)); end for;
  end if;
  K:=BaseRing(Parent(G));
  unit:= K!1;
  for z in ans do unit*:= Norm(LeadingCoefficient(z)); end for;
  unit:= LeadingCoefficient(G)/unit;
  return ans, unit;
end function;

absolute_factorization_squarefree:= function(G)
  // special case for univariate polynomials
  b, q, i:= IsUnivariate(G);
  if b then
    unit:= LeadingCoefficient(q);
    UP:=Parent(q);
    E<a>:= quo<UP | q>;
    MP:= PolynomialRing(E,2);
    AssignNames(~MP, [Sprint(Parent(G).1), Sprint(Parent(G).2)]);
    return [* MP.i - a *], unit;
  end if;

  // select sparse or dense strategy
  d:=TotalDegree(G);
  if Length(G) le 2*d^(1.3) then
    return absolute_factorization_squarefree_sparse(G);
  else
    return absolute_factorization_squarefree_dense(G);
  end if;
end function;

absolute_factorization_two_variables:= function(G)
  MP:= Parent(G);
  K:= BaseRing(MP);
  vprintf AbsoluteFactorization,2: "absolute_factorization_two_variables, computing the squarefree factorization\n";
  sqf:= SquareFreeFactorization(G);
  ans:= [* *];
  unit:= K!1;
  for f in sqf do
    vprintf AbsoluteFactorization,1: "absolute_factorization_two_variables, dealing with a squarefree factor of degree %o\n", TotalDegree(f[1]);
    tmp1, tmp2:= absolute_factorization_squarefree(f[1]);
    for z in tmp1 do Append(~ans, < z, f[2] >); end for;
    unit *:= tmp2;
  end for;  
  return ans, unit;  
end function;

/*
 * Export
 */

absfact_random_example:= function(MP, d, r, H)
/*
  Returns a random irreducible polynomial in MP of degree d*r 
  with r absolutely irreducible factors.
  It is computed as Res_z(q(z), P(x,y,z)), where
  deg_{y,z} P = d, deg_z P = r-1, deg_z q = r,
  The coefficients of q and P are randomly taken in [0..H-1].
*/
  K:=BaseRing(MP);
  UP:= PolynomialRing(K);
  q:= UP![ K!Random(H) : i in [0..r-1]] + (UP.1)^r;
  while not IsIrreducible(q) do
    q:= UP![ K!Random(H) : i in [0..r-1]] + (UP.1)^r;
  end while;
  UPx<x>:= PolynomialRing(K);
  UPy<y>:= PolynomialRing(UPx);
  UPz<z>:= PolynomialRing(UPy);
  Ft:= UPz![ UPy![ UPx![ K!Random(H) : j in [0..d-i]] : i in [0..d-1]] : k in [0..r-1]] + UPz!((UPy.1)^d);
  F:= Resultant(UPz!ChangeUniverse(Coefficients(q),UPy), Ft);
  return &+[MP| MultivariatePolynomial(MP, Coefficient(F,i), 2) * MP.1^i : i in [0..Degree(F)]];
end function;

absfact_print_with_parent:= procedure(facts,u)
  "Absolute factorization";
  "Scalar : ",u;
  "Factors : ";
  for f in facts do
    f,BaseRing(Parent(f[1]));
  end for;
end procedure;

intrinsic AbsoluteFactorization(G::RngMPolElt) -> List, RngElt
{
  Factorization of the polynomial G into absolutely irreducible factors
  together with the appropriate scalar.
}
  vprintf AbsoluteFactorization: "--> AbsoluteFactorization\n";
begin_time:= Cputime();
  IndentPush();
  MP:= Parent(G);
  require Rank(MP) eq 2: "AbsoluteFactorization, the input polynomial must be bivariate";
  K:= BaseRing(MP);
  require IsField(K): "AbsoluteFactorization, the base ring must be a field";
  require HasPolynomialFactorization(K): "AbsoluteFactorization, the base field must support irreducible of univariate polynomials";
  d:= TotalDegree(G);
  p:= Characteristic(K);
  vprintf AbsoluteFactorization: "AbsoluteFactorization, total degree = %o\n", d;
  require (p eq 0 or p gt d*(d-1)): "AbsoluteFactorization, the characteristic is too small";
  ans,unit:= absolute_factorization_two_variables(G);
  IndentPop();
  vprintf AbsoluteFactorization: "<-- AbsoluteFactorization, elapsed time: %o\n", Cputime() - begin_time;  
  return ans,unit;
end intrinsic;
