diophantsolve := proc(a,b,s,t,c)
# Assume as + bt = 1. Solve a*sigma + b*tau = c with deg(sigma) < deg(b) and deg(tau) < deg(a).
local sigma,tau,q;
    sigma := mulrpoly(s,c);
    sigma := remrpoly(sigma, b, 'q');
    tau := addrpoly(mulrpoly(t,c), mulrpoly(q,a));
    return sigma,tau;
end proc:


hensellift_univar := proc(u0,v0,f,p)
# Assume u0*v0 = f mod p.
#        gcd(u0,v0) = 1 mod p.
#        f is over Q.
#        u0,v0 are over Z/p.
# Lift to u1*v1 = f mod p^2, etc.
#
# This version assumes f, u0, and v0 are monic like in the paper.
# But in RECDEN the polynomials in the triangular set T are primitive over ZZ.
# So to use this you've either got to make f, u0 and v0 monic on input
# or scale the output factorization f = u v as needed.
local u,v,c,c_p,S,Sigma,k,pk,pk_mod,uk,vk,e,u_rr,v_rr,i,T,s;

    printf("  Hensel: zero-divisor?  %a mod %d\n", rpoly(u0),p);
    printf("  Hensel: f = %a\n", rpoly(f) );
    T := getalgexts(f); # extracted from f
    # First iteration of Hensel lifting.
    u := retcharpoly(modsrpoly(u0)); # Lift from modulo a prime.
    v := retcharpoly(modsrpoly(v0));
    for i from 1 to nops(T) do
        u := liftrpoly(u);
        v := liftrpoly(v);
    od;
    u := reduce_by_tset(u,T);
    v := reduce_by_tset(v,T);
    e := subrpoly(f, mulrpoly(u,v));
    pk := p;
    pk_mod := p;
    #print(1,u,v,c);

    # Initial solve of extended Euclidean algorithm.
    S := traperror(gcdexrpoly(u0,v0));
    if S[1] = "inverse does not exist" then
        error(S);
    fi;

    # All subsequent iterations of Hensel lifting.
    for k from 2 to 10 do
        c_p := phirpoly(scarpoly(1/pk, e),p);
        Sigma := diophantsolve(u0,v0,S[2],S[3],c_p);
        uk := retcharpoly(Sigma[2]);
        vk := retcharpoly(Sigma[1]);
        for i from 1 to nops(T) do
            uk := liftrpoly(uk);
            vk := liftrpoly(vk);
        od;
        uk := reduce_by_tset(uk,T);
        vk := reduce_by_tset(vk,T);
        u := addrpoly(u, scarpoly(pk,uk));
        v := addrpoly(v, scarpoly(pk,vk));
        pk_mod := pk_mod*p;
        u_rr := irrrpoly(phirpoly(u,pk_mod));
        #print(k,u_rr);
        if u_rr <> FAIL then
            u_rr := retextsrpoly(u_rr);
            u_rr := reduce_by_tset(u_rr, T);
            # RECDEN requires the primitive associate over ZZ
	    # u_rr := ipprpoly(u_rr); # -21 z + 6 ==> 7 z - 2
            if divrpoly(f, u_rr, v_rr) then
                printf("  Hensel: zero-divisor over Q = %a\n",rpoly(u_rr));
                return u_rr, v_rr;
            fi;
        fi;

        # Set up for next iteration.
        pk := pk_mod;
        u := retcharpoly(u);
        v := retcharpoly(v);
        for i from 1 to nops(T) do
            u := liftrpoly(u);
            v := liftrpoly(v);
        od;
        u := reduce_by_tset(u,T);
        v := reduce_by_tset(v,T);
        e := subrpoly(f, mulrpoly(u,v));
        #print(k,u,v,e);
    od;
    printf("  Hensel: lifting failed\n");
    return FAIL;
end proc: