/* Quadratic Hensel lifting for Fp[x,y] */

Z := IntegerRing();
Zx<x> := PolynomialRing(Z);

height := function(f)
    m := 0;
    for c in Coefficients(f) do m := Max(Abs(c),m); end for;
    return m;
end function;

FactorBound := function(f)
    h := height(f);
    d := Degree(f);
    return 2^(d-1)*h*Ceiling(Sqrt(d));
end function;

getpoly := function(d,m,p)
    n := p^m;
    /*C := [ Random(-n,n) : i in [0..d-1] ];*/
    C := [ Random(1,n) : i in [0..d-1] ];
    g := Zx!C;
    return x^d+g;
end function;

DivideModm := function(f,g,m)
   Zm := ResidueClassRing(m);
   Zmx<x> := PolynomialRing(Zm);
   F := Zmx!f;
   G := Zmx!g;
   Q,R := Quotrem(F,G);
   q := Zx!Q;
   r := Zx!R;
   return q,r;
end function;

reduce := function( f, m ) 
    return Zx![ c mod m : c in Coefficients(f) ]; 
end function;
mods := function(a,p,p2) 
   if a ge p2 then return a-p; else return a; end if; 
end function;
mapmods := function(f,p,p2) 
   return Zx![ mods(c,p,p2) : c in Coefficients(f) ];
end function;


HenselLift := function( f, g0, h0, p ) 

  Fp := GaloisField(p);
  Fpy<y> := PolynomialRing(Fp);

  g0 := mapmods(g0,p,p/2);
  h0 := mapmods(h0,p,p/2);
  G0,S0,T0 := XGCD( Fpy!g0, Fpy!h0 );
  s := Zx!S0; t := Zx!T0;

  B := FactorBound(f); // height bound on factors of a

  g := g0;
  h := h0;
  m := p;
  e := f-g*h;
  k := 0;
  while m lt 2*B do

    k := k+1; /* print "Step", k; */
    m := m^2; m2 := (m-1)/2;

    e := reduce(e,m);
    q,r := DivideModm(s*e,h,m);
    u := t*e+q*g; u := reduce(u,m);
    g := g + u; g := mapmods( g, m, m2 );
    h := h + r; h := mapmods( h, m, m2 );
    e := f-g*h;

    if e eq 0 then break; end if;
    if m gt 2*B then break; end if;

    // Lift diophantine equation solutions s and t
    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;

  end while;
  return  e, g, h  ;
end function;


d := 800;
m := 800;
n := 10;
p := 2^31-1;
p := 2^50-27;
g := getpoly(d,m,p);
h := getpoly(d,m,p);
time for i in [1..n] do f := g*h; end for;
g0 := reduce(g,p);
h0 := reduce(h,p);
printf "Start Hensel lift:  d=%m  m=%m  p=%m\n", d,m,p;
time e, G, H := HenselLift(f,g0,h0,p);
e; G-g; H-h; // should all be 0