// Quadratic Hensel lifting for Zp[x,y] for p prime.
dx :=  80; dx;
dy :=  80; dy;
n := 1;
p := 1125899906842597;
p := 2^31-1; p;
alpha := 3;

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

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

getcoeff := function(dy,p)
    E := getexpons(dy,3);
    C := 0;
    for e in E do C := C + Random(p)*y^e; end for;
    return C; // C has 3 distinct terms
end function;

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

G := getpoly(dx,dy,p);
H := getpoly(dx,dy,p);
time f := G*H;

print "Multiplication done";
    
reduce := function( f, m ) return Zpxy![ c mod m : c in Coefficients(f) ]; end function;

g0 := reduce(G,y-alpha); G0 := Zpy!g0;
h0 := reduce(H,y-alpha); H0 := Zpy!h0;
gcd,S0,T0 := XGCD( G0, H0 );

print "gcd", gcd ;

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

DivideModm := function(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;

time for i := 1 to n do

  s := convertytox(S0);
  t := convertytox(T0);

  g := g0;
  h := h0;
  m := Zpy!(y-alpha);
  e := f-g*h;

  k := 1;
  while Degree(m) le 2*dy 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;
end for;

e; G-g; H-h; // should be 0