Chapter 6, Exercise 15, page 276, GCL
-------------------------------------

For this question, to do multivariate Hensel lifting we do not want to
use Algorithm 6.4 at the end of Chapter 6.
(I am not discussing all of the details of Algorithm 6.4.)
We want to use a more direct generalization of Algorithm 6.1, generalized
to the multivariate case.

What you should understand is:

  (i) The basic idea of univariate Hensel lifting (Algorithm 6.1).

 (ii) The idea of how we generalize Newton's iteration from the p-adic
      case (as in univariate Hensel lifting) to the I-adic case where
      I is an ideal such as I = <y-1, z> .

I present below the Multivariate Hensel Lifting Algorithm which we want to
use in Exercise 15. It is based on Theorem 6.6 in the text.


The multivariate Hensel construction of Theorem 6.6
---------------------------------------------------

Condition: a(x,y,z) must be primitive w.r.t. the variable x.

Note the values mentioned in the statement of Theorem 6.6:
  We are using  p=13, l=1 .
  We will lift with respect to the polynomial a(x,y,z).
  We are using the evaluation homomorphism with kernel I = <y-1,z> .
  Note that p does not divide the leading coefficient of a (mod I).
  In part (a) you determine u1(x) = GCD(a,b) (mod I);  then define
    the cofactor w1(x) = a/u1(x) (mod I).

Keep in mind that Theorem 6.6 is talking about a Hensel lifting process
for lifting from the solution mod I, to solutions mod I^2, mod I^3, ... .
The lifting is with respect to I.
All computation will be done mod p at all times.
(We have determined that, for the given problem, we can equate Zp with Z.)

Theorem 6.6 states that we can compute polynomials

        u[k], w[k]  in  Zp[x,y,z] / I^k

such that
        a(x,y,z) == u[k]*w[k]  mod <I^k,p>
and
        u[k] == u1(x) mod <I,p>;   w[k] == w1(x) mod <I,p> .


Multivariate Hensel Lifting Algorithm
=====================================

Note: The textbook presents a Hensel lifting process in the proof of
Theorem 6.6 but then goes on to develop more efficient variations of it
to be presented in the Algorithms at the end of the chapter.
We are ignoring the developments beyond Theorem 6.6.
We just want to use the basic Hensel algorithm suggested in Theorem 6.6.

This multivariate algorithm follows the same outline as Algorithm 6.1.

Note: In Maple you should assign `mod` := mods;  so that the symmetric
representation will be used throughout the session.
Also, the algorithm below assumes that the ideal I = <y-1,z> is represented
by a Maple list Id := [y=1,z=0].


# Step 1: Define new polynomial and its mod <I,p> factors
# ======  (i.e. as in Alg. 6.1, multiply through by the leading coefficient)

alpha  <--  lcoeff(a(x,y,z), x)
gamma  <--  alpha
a(x,y,z)  <--  expand( gamma*a(x,y,z) ) mod p
u1(x)  <--  eval(gamma, Id) * N(u1(x)) mod p
w1(x)  <--  eval(alpha, Id) * N(w1(x)) mod p
         # where the operation N(..) simply means the normalization
	 # "divide by lcoeff to make monic"


# Step 2: Apply extended Euclidean algorithm to u1(x), w1(x) in Zp[x]
# ======

s(x), t(x)  <--  polynomials in Zp[x] computed by Algorithm 2.2 such that
                       s(x)*u1(x) + t(x)*w1(x) == 1  (mod p)


# Step 3: Initialization for the iteration
# ======

u(x,y,z)  <--  replace_lc(u1(x), gamma)
w(x,y,z)  <--  replace_lc(w1(x), alpha)
e(x,y,z)  <--  expand( a(x,y,z) - u(x,y,z)*w(x,y,z) ) mod p


# Step 4: Iterate until either the factorization in Zp[x,y,z] is obtained or
# ======  else the degree bound is reached.

#  (NOTE: This main iteration loop of the algorithm is presented in the box
#   at the top of page 263 of the text.  We are ignoring where it mentions
#   "Substitute" since that is a detail which Maple's "mtaylor" command takes
#   care of for you.)

for k from 1 to degree(a(x,y,z), {y,z}) while e(x,y,z) <> 0 do

   For each term of total degree k in  mtaylor(e(x,y,z), Id, k+1) mod p  do

      Pick off the coefficient c(x)

      Solve in the domain Zp[x] the polynomial equation
                sigma(x)*u1(x) + tau(x)*w1(x) == c(x)  (mod p)
         # (solving this equation is done just as in Algorithm 6.1 -- use
         #  Step 4.1 of Algorithm 6.1, ignoring its definition of c(x))

      u(x,y,z)  <--  u(x,y,z)  +  tau(x) * (y-1)^k1 * z^k2
      w(x,y,z)  <--  w(x,y,z)  +  sigma(x) * (y-1)^k1 * z^k2
         # where k1, k2 are the powers of (y-1), z in the term from which c(x)
	 # was picked off
   end do

   e(x,y,z)  <--  expand( a(x,y,z) - u(x,y,z)*w(x,y,z) ) mod p
end do


# Step 5: Check termination status
# ======

if e(x,y,z) = 0 then
   # Factorization obtained - remove contents
   u(x,y,z)  <--  expand( u(x,y,z) ) mod p
   w(x,y,z)  <--  expand( w(x,y,z) ) mod p
   delta  <--  content(u(x,y,z), x)
   u(x,y,z)  <--  u(x,y,z) / delta          # polynomial division mod p
   w(x,y,z)  <--  w(x,y,z) / (gamma/delta)  # polynomial divisions mod p
   return u(x,y,z), w(x,y,z)
else
   return "no such factorization exists"
fi