// Magma code for Quadratic Hensel lifting
// Michael Monagan and Garrett Paluck, January 2022
//
// Reference: Modern Computer Algebra, von zur Gathen and Gerhard
//     Algorithm 15.10 Hensel step
//     Algorithm 15.17 Multifactor Hensel lifting
// Note: we have modified it to terminate early if the error is zero

p := 2^31-1;

Fp := GaloisField(p);
Zpy<y> := PolynomialRing(Fp);
Zpxy<x> := PolynomialRing(Zpy);

function getexpons(dy,n)
   L := [];
   while #L lt Minimum(dy+1,n) do
      e := Random(dy);
      if not( e in L ) then L := Append(L,e); end if;
   end while;
   return L;
end function;

function getcoeff(dy,p)
   E := getexpons(dy,dy); 
   C := 0;
   for e in E do C := C + Random(1,p-1)*y^e; end for;
   return C;
end function;

function getpoly(dx,dy,p)
   C := [ getcoeff(dy,p) : i in [0..dx-1] ];
   f := Zpxy!C;
   return x^dx+f;
end function;

function reduce( f, m )
   return Zpxy![ c mod m : c in Coefficients(f) ];
end function;

function convertytox( f )
g := Zpxy!Coefficients(f);
return g;
end function;

function DivideModm(f,g,m)
   I := ideal<Zpy|m>;
   Q := quo<Zpy|I>;
   Rx := PolynomialRing(Q);
   F := Rx!f;
   G := Rx!g;
   Q,R := Quotrem(F,G);
   q := Zpxy!Q;
   r := Zpxy!R;
   return q,r ;
end function;

function BivariateHensel(f,F,n,m,L)
   
   if n eq 1 then
      return [reduce(f,m^L)];
   end if;

   k := Floor(n/2);

   g := 1; 
   F1 := [];
   for i := 1 to k do
      g := g*F[i];
      F1 := Append(F1,F[i]);
   end for;
   h := 1; 
   F2 := [];
   for i := k+1 to n do
      h := h*F[i];
      F2 := Append(F2,F[i]);
   end for;

   G0 := Zpy!g; H0 := Zpy!h;
   gcd,S0,T0 := XGCD(G0,H0);
   s := convertytox(S0); t := convertytox(T0);

   e := f-g*h;
   K := 1;
   M := m;
   while Degree(M) le L do
      //print "Step", K, "Degree", Degree(m);
      M := M^2;

      e := reduce(e,M);
      q,r := DivideModm(s*e,h,M);
      u := t*e+q*g; u := reduce(u,M);
      g := g + u;
      h := h + r;

      e := f-g*h;
      if e eq 0 then break; end if;

      // Lift diophantine equation solutions
      b := s*g + t*h - 1; b := reduce(b,M);
      c,d := DivideModm(s*b,h,M);
      u := t*b+c*g; u := reduce(u,M);
      s := s - d;
      t := t - u;

      K := K+1;
   end while;

   F1 := BivariateHensel(g,F1,k,m,L);
   F2 := BivariateHensel(h,F2,n-k,m,L);

   H := [];
   for i := 1 to k do
      H := Append(H,F1[i]);
   end for;
   for i := 1 to (n-k) do
      H := Append(H,F2[i]);
   end for;

   return H;

end function;


n := 4;
print "number of factors is",n;
dx := 20;
print "degree in x", dx;
dy := 20;
print "degree in y", dy;

alpha := 3;
print "evaluating at y =",alpha;
m := Zpy!(y-alpha);

f1 := getpoly(dx,dy,p);
f2 := getpoly(dx,dy,p);
f3 := getpoly(dx,dy,p);
f4 := getpoly(dx,dy,p);

f := f1*f2*f3*f4;
F := [reduce(f1,m),reduce(f2,m),reduce(f3,m),reduce(f4,m)];

time G := BivariateHensel(f,F,n,m,4*dy);

G[1] - f1; G[2] - f2; G[3] - f3; G[4] - f4;