Index: factory/facHensel.cc
===================================================================
--- factory/facHensel.cc	(revision ad3c3ff7e9858e8efdced4631a25aebfcf700f36)
+++ factory/facHensel.cc	(revision ad3c3ff7e9858e8efdced4631a25aebfcf700f36)
@@ -0,0 +1,1557 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file facHensel.cc
+ * 
+ * This file implements functions to lift factors via Hensel lifting and 
+ * functions for modular multiplication and division with remainder.
+ * 
+ * ABSTRACT: Hensel lifting is described in "Efficient Multivariate 
+ * Factorization over Finite Fields" by L. Bernardin & M. Monagon. Division with 
+ * remainder is described in "Fast Recursive Division" by C. Burnikel and 
+ * J. Ziegler. Karatsuba multiplication is described in "Modern Computer 
+ * Algebra" by J. von zur Gathen and J. Gerhard.   
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#include "assert.h"
+#include "debug.h"
+#include "timing.h"
+
+#include "facHensel.h"
+#include "cf_util.h"
+
+#ifdef HAVE_NTL
+#include <NTL/lzz_pEX.h>
+#include "NTLconvert.h"
+
+static inline
+CanonicalForm
+mulNTL (const CanonicalForm& F, const CanonicalForm& G)
+{
+  if (F.inCoeffDomain() || G.inCoeffDomain())
+    return F*G;
+  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
+  ASSERT (F.level() == G.level(), "expected polys of same level");
+  if (CFFactory::gettype() == GaloisFieldDomain)
+    return F*G;
+  zz_p::init (getCharacteristic());
+  Variable alpha;
+  CanonicalForm result;
+  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+  {
+    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+    zz_pE::init (NTLMipo);
+    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
+    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
+    mul (NTLF, NTLF, NTLG);
+    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
+  }
+  else
+  {
+    zz_pX NTLF= convertFacCF2NTLzzpX (F);
+    zz_pX NTLG= convertFacCF2NTLzzpX (G);
+    mul (NTLF, NTLF, NTLG);
+    result= convertNTLzzpX2CF(NTLF, F.mvar());
+  }
+  return result;
+}
+
+static inline
+CanonicalForm
+modNTL (const CanonicalForm& F, const CanonicalForm& G)
+{
+  if (F.inCoeffDomain() && G.isUnivariate())
+    return F;
+  else if (F.inCoeffDomain() && G.inCoeffDomain())
+    return mod (F, G);
+  else if (F.isUnivariate() && G.inCoeffDomain())
+    return mod (F,G);
+    
+  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
+  ASSERT (F.level() == G.level(), "expected polys of same level");
+  if (CFFactory::gettype() == GaloisFieldDomain)
+    return mod (F, G);
+  zz_p::init (getCharacteristic());
+  Variable alpha;
+  CanonicalForm result;
+  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+  {
+    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
+    zz_pE::init (NTLMipo);
+    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
+    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
+    rem (NTLF, NTLF, NTLG);
+    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
+  }
+  else
+  {
+    zz_pX NTLF= convertFacCF2NTLzzpX (F);
+    zz_pX NTLG= convertFacCF2NTLzzpX (G);
+    rem (NTLF, NTLF, NTLG);
+    result= convertNTLzzpX2CF(NTLF, F.mvar());
+  }
+  return result;
+}
+
+static inline
+CanonicalForm
+divNTL (const CanonicalForm& F, const CanonicalForm& G)
+{
+  if (F.inCoeffDomain() && G.isUnivariate())
+    return F;
+  else if (F.inCoeffDomain() && G.inCoeffDomain())
+    return div (F, G);
+  else if (F.isUnivariate() && G.inCoeffDomain())
+    return div (F,G);
+
+  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
+  ASSERT (F.level() == G.level(), "expected polys of same level");
+  if (CFFactory::gettype() == GaloisFieldDomain)
+    return div (F, G);
+  zz_p::init (getCharacteristic());
+  Variable alpha;
+  CanonicalForm result;
+  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+  {
+    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
+    zz_pE::init (NTLMipo);
+    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
+    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
+    div (NTLF, NTLF, NTLG);
+    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
+  }
+  else
+  {
+    zz_pX NTLF= convertFacCF2NTLzzpX (F);
+    zz_pX NTLG= convertFacCF2NTLzzpX (G);
+    div (NTLF, NTLF, NTLG);
+    result= convertNTLzzpX2CF(NTLF, F.mvar());
+  }
+  return result;
+}
+
+/*static inline
+void
+divremNTL (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q, 
+           CanonicalForm& R)
+{
+  if (F.inCoeffDomain() && G.isUnivariate())
+  {
+    R= F;
+    Q= 0;
+  }
+  else if (F.inCoeffDomain() && G.inCoeffDomain())
+  {
+    divrem (F, G, Q, R);
+    return;
+  }
+  else if (F.isUnivariate() && G.inCoeffDomain())
+  {
+    divrem (F, G, Q, R);
+    return;
+  }
+
+  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
+  ASSERT (F.level() == G.level(), "expected polys of same level");
+  if (CFFactory::gettype() == GaloisFieldDomain)
+  {
+    divrem (F, G, Q, R);
+    return;
+  }
+  zz_p::init (getCharacteristic());
+  Variable alpha;
+  CanonicalForm result;
+  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+  {
+    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
+    zz_pE::init (NTLMipo);
+    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
+    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
+    zz_pEX NTLQ; 
+    zz_pEX NTLR;
+    DivRem (NTLQ, NTLR, NTLF, NTLG);
+    Q= convertNTLzz_pEX2CF(NTLQ, F.mvar(), alpha);
+    R= convertNTLzz_pEX2CF(NTLR, F.mvar(), alpha);
+    return;
+  }
+  else
+  {
+    zz_pX NTLF= convertFacCF2NTLzzpX (F);
+    zz_pX NTLG= convertFacCF2NTLzzpX (G);
+    zz_pX NTLQ; 
+    zz_pX NTLR;
+    DivRem (NTLQ, NTLR, NTLF, NTLG);
+    Q= convertNTLzzpX2CF(NTLQ, F.mvar());
+    R= convertNTLzzpX2CF(NTLR, F.mvar());
+    return;
+  }
+}*/
+
+CanonicalForm mod (const CanonicalForm& F, const CFList& M) 
+{
+  CanonicalForm A= F;
+  for (CFListIterator i= M; i.hasItem(); i++) 
+    A= mod (A, i.getItem());
+  return A;
+}
+
+CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B, 
+                       const CanonicalForm& M)
+{
+  if (A.isZero() || B.isZero())
+    return 0;
+
+  ASSERT (M.isUnivariate(), "M must be univariate");
+
+  CanonicalForm F= mod (A, M);
+  CanonicalForm G= mod (B, M);
+  if (F.inCoeffDomain() || G.inCoeffDomain())
+    return F*G;
+  Variable y= M.mvar();
+  int degF= degree (F, y);
+  int degG= degree (G, y);
+
+  if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) && 
+      (F.level() == G.level()))
+  { 
+    CanonicalForm result= mulNTL (F, G);
+    return mod (result, M);
+  }
+  else if (degF <= 1 && degG <= 1)
+  { 
+    CanonicalForm result= F*G;
+    return mod (result, M);
+  }
+ 
+  int m= (int) ceil (degree (M)/2.0);
+  if (degF >= m || degG >= m)
+  {
+    CanonicalForm MLo= power (y, m);
+    CanonicalForm MHi= power (y, degree (M) - m);
+    CanonicalForm F0= mod (F, MLo);
+    CanonicalForm F1= div (F, MLo);
+    CanonicalForm G0= mod (G, MLo);
+    CanonicalForm G1= div (G, MLo);
+    CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
+    CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
+    CanonicalForm F0G0= mulMod2 (F0, G0, M);
+    return F0G0 + MLo*(F0G1 + F1G0);
+  }
+  else
+  {
+    m= (int) ceil (tmax (degF, degG)/2.0);
+    CanonicalForm yToM= power (y, m);
+    CanonicalForm F0= mod (F, yToM);
+    CanonicalForm F1= div (F, yToM);
+    CanonicalForm G0= mod (G, yToM);
+    CanonicalForm G1= div (G, yToM);
+    CanonicalForm H00= mulMod2 (F0, G0, M);
+    CanonicalForm H11= mulMod2 (F1, G1, M);
+    CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
+    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
+  } 
+  DEBOUTLN (cerr, "fatal end in mulMod2");
+}
+
+CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B, 
+                      const CFList& MOD)
+{
+  if (A.isZero() || B.isZero())
+    return 0;
+
+  if (MOD.length() == 1)
+    return mulMod2 (A, B, MOD.getLast());
+
+  CanonicalForm M= MOD.getLast();
+  CanonicalForm F= mod (A, M);
+  CanonicalForm G= mod (B, M);
+  if (F.inCoeffDomain() || G.inCoeffDomain())
+    return F*G;
+  Variable y= M.mvar();
+  int degF= degree (F, y);
+  int degG= degree (G, y);
+
+  if ((degF <= 1 && F.level() <= M.level()) && 
+      (degG <= 1 && G.level() <= M.level()))
+  {
+    CFList buf= MOD;
+    buf.removeLast();
+    if (degF == 1 && degG == 1)
+    {
+      CanonicalForm F0= mod (F, y);
+      CanonicalForm F1= div (F, y);
+      CanonicalForm G0= mod (G, y);
+      CanonicalForm G1= div (G, y);
+      if (degree (M) > 2)
+      {
+        CanonicalForm H00= mulMod (F0, G0, buf);
+        CanonicalForm H11= mulMod (F1, G1, buf);
+        CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
+        return H11*y*y + (H01 - H00 - H11)*y + H00;
+      }
+      else //here degree (M) == 2
+      {
+        buf.append (y);
+        CanonicalForm F0G1= mulMod (F0, G1, buf);
+        CanonicalForm F1G0= mulMod (F1, G0, buf);
+        CanonicalForm F0G0= mulMod (F0, G0, MOD);
+        CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
+        return result;
+      }
+    }
+    else if (degF == 1 && degG == 0)
+      return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
+    else if (degF == 0 && degG == 1)
+      return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
+    else
+      return mulMod (F, G, buf);
+  }
+  int m= (int) ceil (degree (M)/2.0);
+  if (degF >= m || degG >= m)
+  {
+    CanonicalForm MLo= power (y, m);
+    CanonicalForm MHi= power (y, degree (M) - m);
+    CanonicalForm F0= mod (F, MLo);
+    CanonicalForm F1= div (F, MLo);
+    CanonicalForm G0= mod (G, MLo);
+    CanonicalForm G1= div (G, MLo);
+    CFList buf= MOD;
+    buf.removeLast();
+    buf.append (MHi);
+    CanonicalForm F0G1= mulMod (F0, G1, buf);
+    CanonicalForm F1G0= mulMod (F1, G0, buf);
+    CanonicalForm F0G0= mulMod (F0, G0, MOD);
+    return F0G0 + MLo*(F0G1 + F1G0);
+  }
+  else
+  {
+    m= (int) ceil (tmax (degF, degG)/2.0);
+    CanonicalForm yToM= power (y, m);
+    CanonicalForm F0= mod (F, yToM);
+    CanonicalForm F1= div (F, yToM);
+    CanonicalForm G0= mod (G, yToM);
+    CanonicalForm G1= div (G, yToM);
+    CanonicalForm H00= mulMod (F0, G0, MOD);
+    CanonicalForm H11= mulMod (F1, G1, MOD);
+    CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
+    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
+  } 
+  DEBOUTLN (cerr, "fatal end in mulMod");
+}
+
+CanonicalForm prodMod (const CFList& L, const CanonicalForm& M) 
+{
+  if (L.isEmpty())
+    return 1;
+  int l= L.length();
+  if (l == 1)
+    return mod (L.getFirst(), M);
+  else if (l == 2) {
+    CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
+    return result;
+  }
+  else
+  {
+    l /= 2;
+    CFList tmp1, tmp2;
+    CFListIterator i= L;
+    CanonicalForm buf1, buf2;
+    for (int j= 1; j <= l; j++, i++) 
+      tmp1.append (i.getItem());
+    tmp2= Difference (L, tmp1);
+    buf1= prodMod (tmp1, M);
+    buf2= prodMod (tmp2, M);
+    CanonicalForm result= mulMod2 (buf1, buf2, M);
+    return result;
+  } 
+}
+
+CanonicalForm prodMod (const CFList& L, const CFList& M) 
+{
+  if (L.isEmpty())
+    return 1;
+  else if (L.length() == 1)
+    return L.getFirst();
+  else if (L.length() == 2)
+    return mulMod (L.getFirst(), L.getLast(), M);
+  else
+  {
+    int l= L.length()/2;
+    CFListIterator i= L;
+    CFList tmp1, tmp2;
+    CanonicalForm buf1, buf2;
+    for (int j= 1; j <= l; j++, i++) 
+      tmp1.append (i.getItem());
+    tmp2= Difference (L, tmp1);
+    buf1= prodMod (tmp1, M);
+    buf2= prodMod (tmp2, M);
+    return mulMod (buf1, buf2, M);
+  }
+}
+
+static inline
+CFList split (const CanonicalForm& F, const int m, const Variable& x) 
+{
+  CanonicalForm A= F;
+  CanonicalForm buf= 0;
+  bool swap= false;
+  if (degree (A, x) <= 0)
+    return CFList(A);
+  else if (x.level() != A.level())
+  {
+    swap= true;
+    A= swapvar (A, x, A.mvar());
+  }
+
+  int j= (int) floor ((double) degree (A)/ m);
+  CFList result; 
+  CFIterator i= A;
+  for (; j >= 0; j--)
+  {
+    while (i.hasTerms() && i.exp() - j*m >= 0) 
+    {
+      if (swap)
+        buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
+      else
+        buf += i.coeff()*power (x, i.exp() - j*m);
+      i++;
+    }
+    if (swap)
+      result.append (swapvar (buf, x, F.mvar()));
+    else
+      result.append (buf);  
+    buf= 0;
+  }
+  return result; 
+}
+
+static inline
+void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q, 
+               CanonicalForm& R, const CFList& M);
+
+static inline 
+void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q, 
+               CanonicalForm& R, const CFList& M)
+{
+  CanonicalForm A= mod (F, M);
+  CanonicalForm B= mod (G, M);
+  Variable x= Variable (1);
+  int degB= degree (B, x);
+  int degA= degree (A, x);
+  if (degA < degB)
+  {
+    Q= 0;
+    R= A;
+    return;
+  }
+  ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
+  if (degB < 1)
+  {
+    divrem (A, B, Q, R);
+    Q= mod (Q, M);
+    R= mod (R, M);
+    return;
+  }
+  
+  int m= (int) ceil ((double) (degB + 1)/2.0) + 1; 
+  CFList splitA= split (A, m, x);
+  CFList splitB= split (B, m, x);
+  if (splitA.length() == 3)
+    splitA.insert (0);
+  if (splitA.length() == 2)
+  {
+    splitA.insert (0);
+    splitA.insert (0);
+  } 
+  if (splitA.length() == 1)
+  {
+    splitA.insert (0);
+    splitA.insert (0);
+    splitA.insert (0);
+  }
+
+  CanonicalForm xToM= power (x, m);
+
+  CFListIterator i= splitA;
+  CanonicalForm H= i.getItem();
+  i++;
+  H *= xToM;
+  H += i.getItem();
+  i++;
+  H *= xToM;
+  H += i.getItem();
+  i++;
+
+  divrem32 (H, B, Q, R, M);
+
+  CFList splitR= split (R, m, x);
+  if (splitR.length() == 1)
+    splitR.insert (0);
+
+  H= splitR.getFirst();
+  H *= xToM;
+  H += splitR.getLast();
+  H *= xToM;
+  H += i.getItem();
+
+  CanonicalForm bufQ;
+  divrem32 (H, B, bufQ, R, M);
+
+  Q *= xToM;
+  Q += bufQ;
+  return;
+}
+
+static inline
+void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q, 
+               CanonicalForm& R, const CFList& M) 
+{
+  CanonicalForm A= mod (F, M);
+  CanonicalForm B= mod (G, M);
+  Variable x= Variable (1);
+  int degB= degree (B, x);
+  int degA= degree (A, x);
+  if (degA < degB) 
+  {
+    Q= 0;
+    R= A;
+    return;
+  }
+  ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
+  if (degB < 1)
+  {
+    divrem (A, B, Q, R);
+    Q= mod (Q, M);
+    R= mod (R, M);
+    return;
+  }
+  int m= (int) ceil ((double) (degB + 1)/ 2.0);
+
+  CFList splitA= split (A, m, x);
+  CFList splitB= split (B, m, x);
+
+  if (splitA.length() == 2)
+  {
+    splitA.insert (0);
+  }
+  if (splitA.length() == 1)
+  {
+    splitA.insert (0);
+    splitA.insert (0);
+  }
+  CanonicalForm xToM= power (x, m);
+
+  CanonicalForm H;
+  CFListIterator i= splitA;
+  i++;
+
+  if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
+  {
+    H= splitA.getFirst()*xToM + i.getItem();
+    divrem21 (H, splitB.getFirst(), Q, R, M);
+  }
+  else 
+  {
+    R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() - 
+       splitB.getFirst()*xToM;
+    Q= xToM - 1;
+  }
+
+  H= mulMod (Q, splitB.getLast(), M);
+
+  R= R*xToM + splitA.getLast() - H;
+
+  while (degree (R, x) >= degB)
+  {
+    xToM= power (x, degree (R, x) - degB);
+    Q += LC (R, x)*xToM;
+    R -= mulMod (LC (R, x), B, M)*xToM;
+    Q= mod (Q, M);
+    R= mod (R, M);
+  }
+
+  return;
+}
+
+void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+              CanonicalForm& R, const CanonicalForm& M) 
+{
+  CanonicalForm A= mod (F, M);
+  CanonicalForm B= mod (G, M);
+  Variable x= Variable (1);
+  int degB= degree (B, x);
+  if (degB > degree (A, x)) 
+  {
+    Q= 0;
+    R= A;
+    return;
+  }
+
+  CFList splitA= split (A, degB, x);
+
+  CanonicalForm xToDegB= power (x, degB);
+  CanonicalForm H, bufQ;
+  Q= 0;
+  CFListIterator i= splitA;
+  H= i.getItem()*xToDegB;
+  i++;
+  H += i.getItem();
+  CFList buf;
+  while (i.hasItem()) 
+  {
+    buf= CFList (M);
+    divrem21 (H, B, bufQ, R, buf);
+    i++;
+    if (i.hasItem())
+      H= R*xToDegB + i.getItem();
+    Q *= xToDegB;
+    Q += bufQ; 
+  }
+  return;
+}
+
+void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+             CanonicalForm& R, const CFList& MOD)
+{
+  CanonicalForm A= mod (F, MOD);
+  CanonicalForm B= mod (G, MOD);
+  Variable x= Variable (1);
+  int degB= degree (B, x);
+  if (degB > degree (A, x)) 
+  {
+    Q= 0;
+    R= A;
+    return;
+  }
+
+  if (degB == 0)
+  {
+    divrem (A, B, Q, R);
+    Q= mod (Q, MOD);
+    R= mod (R, MOD);
+    return;
+  }
+  CFList splitA= split (A, degB, x);
+
+  CanonicalForm xToDegB= power (x, degB);
+  CanonicalForm H, bufQ;
+  Q= 0;
+  CFListIterator i= splitA;
+  H= i.getItem()*xToDegB;
+  i++;
+  H += i.getItem();
+  while (i.hasItem()) 
+  {
+    divrem21 (H, B, bufQ, R, MOD);
+    i++;
+    if (i.hasItem())
+      H= R*xToDegB + i.getItem();
+    Q *= xToDegB;
+    Q += bufQ; 
+  }
+  return;
+}
+
+void sortList (CFList& list, const Variable& x) 
+{ 
+  int l= 1;
+  int k= 1;
+  CanonicalForm buf;
+  CFListIterator m;
+  for (CFListIterator i= list; l <= list.length(); i++, l++) 
+  {
+    for (CFListIterator j= list; k <= list.length() - l; k++) 
+    {
+      m= j;
+      m++;
+      if (degree (j.getItem(), x) > degree (m.getItem(), x)) 
+      {
+        buf= m.getItem();
+        m.getItem()= j.getItem();
+        j.getItem()= buf;
+        j++;
+        j.getItem()= m.getItem();
+      }  
+      else
+        j++;
+    }
+    k= 1;
+  }
+}
+
+static inline
+CFList diophantine (const CanonicalForm& F, const CFList& factors) 
+{
+  CanonicalForm buf1, buf2, buf3, S, T;
+  CFListIterator i= factors;
+  CFList result;
+  if (i.hasItem())
+    i++;
+  buf1= F/factors.getFirst(); 
+  buf2= divNTL (F, i.getItem());
+  buf3= extgcd (buf1, buf2, S, T);
+  result.append (S);
+  result.append (T);
+  if (i.hasItem())
+    i++;
+  for (; i.hasItem(); i++)
+  {
+    buf1= divNTL (F, i.getItem());
+    buf3= extgcd (buf3, buf1, S, T);
+    CFListIterator k= factors;
+    for (CFListIterator j= result; j.hasItem(); j++, k++)
+    {
+      j.getItem()= mulNTL (j.getItem(), S);
+      j.getItem()= modNTL (j.getItem(), k.getItem());
+    }
+    result.append (T);
+  }
+  return result;
+}
+
+void 
+henselStep12 (const CanonicalForm& F, const CFList& factors, 
+              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
+              CFArray& Pi, int j)
+{
+  CanonicalForm E;
+  CanonicalForm xToJ= power (F.mvar(), j);
+  Variable x= F.mvar();
+  // compute the error
+  if (j == 1)
+    E= F[j]; 
+  else
+  {
+    if (degree (Pi [factors.length() - 2], x) > 0) 
+      E= F[j] - Pi [factors.length() - 2] [j];
+    else
+      E= F[j]; 
+  }
+
+  CFArray buf= CFArray (diophant.length());
+  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
+  int k= 0;
+  CanonicalForm remainder;
+  // actual lifting
+  for (CFListIterator i= diophant; i.hasItem(); i++, k++) 
+  {
+    if (degree (bufFactors[k], x) > 0)
+    {
+      if (k > 0)
+        remainder= modNTL (E, bufFactors[k] [0]);
+      else
+        remainder= E;
+    }
+    else
+      remainder= modNTL (E, bufFactors[k]);
+
+    buf[k]= mulNTL (i.getItem(), remainder);
+    if (degree (bufFactors[k], x) > 0)
+      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
+    else 
+      buf[k]= modNTL (buf[k], bufFactors[k]); 
+  }
+  for (k= 1; k < factors.length(); k++)
+    bufFactors[k] += xToJ*buf[k];
+
+  // update Pi [0]
+  int degBuf0= degree (bufFactors[0], x);
+  int degBuf1= degree (bufFactors[1], x);
+  if (degBuf0 > 0 && degBuf1 > 0)
+    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
+  CanonicalForm uIZeroJ;
+  if (j == 1)
+  {
+    if (degBuf0 > 0 && degBuf1 > 0)
+      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
+                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
+    else if (degBuf0 > 0)
+      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
+    else if (degBuf1 > 0)
+      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
+    else
+      uIZeroJ= 0;
+    Pi [0] += xToJ*uIZeroJ; 
+  } 
+  else
+  {
+    if (degBuf0 > 0 && degBuf1 > 0)
+      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
+                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
+    else if (degBuf0 > 0)
+      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
+    else if (degBuf1 > 0)
+      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
+    else
+      uIZeroJ= 0;
+    Pi [0] += xToJ*uIZeroJ; 
+  }
+  CFArray tmp= CFArray (factors.length() - 1);
+  for (k= 0; k < factors.length() - 1; k++)
+    tmp[k]= 0;
+  CFIterator one, two;
+  one= bufFactors [0];
+  two= bufFactors [1];
+  if (degBuf0 > 0 && degBuf1 > 0)
+  {
+    for (k= 1; k <= (int) ceil (j/2.0); k++)
+    {  
+      if (k != j - k + 1)
+      {
+	if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
+	{
+	  tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), (bufFactors[1] [k] +
+                     two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
+	  one++;
+	  two++;
+	}
+	else if (one.exp() == j - k + 1)
+	{
+	  tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), bufFactors[1] [k]) - 
+                     M (k + 1, 1);
+	  one++;
+	}
+	else if (two.exp() == j - k + 1)
+	{
+	  tmp[0] += mulNTL (bufFactors[0] [k], (bufFactors[1] [k] + two.coeff())) - 
+                    M (k + 1, 1);
+	  two++;
+	}
+      }
+      else
+      {
+	tmp[0] += M (k + 1, 1);
+      }
+    }
+  }
+  Pi [0] += tmp[0]*xToJ*F.mvar();
+
+  // update Pi [l] 
+  int degPi, degBuf;
+  for (int l= 1; l < factors.length() - 1; l++) 
+  {
+    degPi= degree (Pi [l - 1], x);
+    degBuf= degree (bufFactors[l + 1], x);
+    if (degPi > 0 && degBuf > 0)
+      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);  
+    if (j == 1)
+    {
+      if (degPi > 0 && degBuf > 0)
+        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
+                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) - 
+                  M (1, l + 1));
+      else if (degPi > 0)
+        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
+      else if (degBuf > 0)
+        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
+    }
+    else
+    {
+      if (degPi > 0 && degBuf > 0)
+      {
+        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
+        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]); 
+      }
+      else if (degPi > 0)
+        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
+      else if (degBuf > 0)
+      {
+        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
+        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
+      }
+      Pi[l] += xToJ*uIZeroJ;
+    }
+    one= bufFactors [l + 1];
+    two= Pi [l - 1];
+    if (two.exp() == j + 1)
+    {
+      if (degBuf > 0 && degPi > 0)
+      {
+	  tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
+	  two++;
+      }
+      else if (degPi > 0)
+      {      
+	  tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
+	  two++;
+      } 
+    }
+    if (degBuf > 0 && degPi > 0) 
+    {
+      for (k= 1; k <= (int) ceil (j/2.0); k++)
+      {  
+	if (k != j - k + 1)
+	{
+	  if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
+	  {
+	    tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), (Pi[l - 1] [k] +
+                       two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
+	    one++;
+	    two++;
+	  }
+	  else if (one.exp() == j - k + 1)
+	  {
+	    tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), Pi[l - 1] [k]) - 
+                       M (k + 1, l + 1);     
+	    one++;
+	  }
+	  else if (two.exp() == j - k + 1)
+	  {
+	    tmp[l] += mulNTL (bufFactors[l + 1] [k], (Pi[l - 1] [k] + two.coeff())) - 
+                      M (k + 1, l + 1);
+	    two++;
+	  }
+	}
+	else
+	  tmp[l] += M (k + 1, l + 1);
+      }  
+    }  
+    Pi[l] += tmp[l]*xToJ*F.mvar();
+  }
+  return;
+}
+
+void 
+henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi, 
+              CFList& diophant, CFMatrix& M) 
+{
+  sortList (factors, Variable (1));
+  Pi= CFArray (factors.length() - 1);
+  CFListIterator j= factors;
+  diophant= diophantine (F[0], factors);
+  DEBOUTLN (cerr, "diophant= " << diophant);
+  j++;
+  Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()));
+  M (1, 1)= Pi [0];
+  int i= 1;
+  if (j.hasItem())
+    j++;
+  for (j; j.hasItem(); j++, i++) 
+  {
+    Pi [i]= mulNTL (Pi [i - 1], j.getItem());
+    M (1, i + 1)= Pi [i];
+  }
+  CFArray bufFactors= CFArray (factors.length());
+  i= 0;
+  for (CFListIterator k= factors; k.hasItem(); i++, k++)
+  {
+    if (i == 0)
+      bufFactors[i]= mod (k.getItem(), F.mvar());
+    else
+      bufFactors[i]= k.getItem();
+  }
+  for (i= 1; i < l; i++)
+    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
+
+  CFListIterator k= factors;
+  for (i= 0; i < factors.length (); i++, k++)
+    k.getItem()= bufFactors[i];
+  factors.removeFirst();
+  return;
+}
+
+void 
+henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
+                    end, CFArray& Pi, const CFList& diophant, CFMatrix& M)
+{
+  CFArray bufFactors= CFArray (factors.length());
+  int i= 0;
+  CanonicalForm xToStart= power (F.mvar(), start);
+  for (CFListIterator k= factors; k.hasItem(); k++, i++)
+  {
+    if (i == 0)
+      bufFactors[i]= mod (k.getItem(), xToStart);
+    else
+      bufFactors[i]= k.getItem();
+  }
+  for (i= start; i < end; i++)
+    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i); 
+  
+  CFListIterator k= factors;
+  for (i= 0; i < factors.length(); k++, i++)
+    k.getItem()= bufFactors [i];
+  factors.removeFirst();
+  return; 
+} 
+
+static inline
+CFList
+biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
+{
+  Variable y= F.mvar();
+  CFList result;
+  if (y.level() == 1)
+  {
+    result= diophantine (F, factors);
+    return result;
+  }
+  else
+  {
+    CFList buf= factors;
+    for (CFListIterator i= buf; i.hasItem(); i++)
+      i.getItem()= mod (i.getItem(), y);
+    CanonicalForm A= mod (F, y);
+    int bufD= 1;
+    CFList recResult= biDiophantine (A, buf, bufD);
+    CanonicalForm e= 1;
+    CFList p;
+    CFArray bufFactors= CFArray (factors.length());
+    CanonicalForm yToD= power (y, d);
+    int k= 0;
+    for (CFListIterator i= factors; i.hasItem(); i++, k++)
+    {
+      bufFactors [k]= i.getItem();
+    }
+    CanonicalForm b;
+    for (k= 0; k < factors.length(); k++) //TODO compute b's faster
+    {
+      b= 1;
+      if (fdivides (bufFactors[k], F))
+        b= F/bufFactors[k];
+      else 
+      {
+	for (int l= 0; l < factors.length(); l++)
+	{
+	  if (l == k)
+	    continue;
+	  else
+	  {
+	    b= mulMod2 (b, bufFactors[l], yToD);
+	  }
+	}
+      }
+      p.append (b);
+    }  
+
+    CFListIterator j= p;
+    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
+      e -= i.getItem()*j.getItem();
+
+    if (e.isZero())
+      return recResult; 
+    CanonicalForm coeffE;
+    CFList s;
+    result= recResult;
+    CanonicalForm g;
+    for (int i= 1; i < d; i++)
+    {
+      if (degree (e, y) > 0)
+        coeffE= e[i];
+      else
+        coeffE= 0;
+      if (!coeffE.isZero())
+      {
+	CFListIterator k= result;
+	CFListIterator l= p;
+        int ii= 0;
+	j= recResult;
+	for (; j.hasItem(); j++, k++, l++, ii++)
+	{
+	  g= coeffE*j.getItem();
+          if (degree (bufFactors[ii], y) <= 0)
+            g= mod (g, bufFactors[ii]);
+          else
+            g= mod (g, bufFactors[ii][0]);
+	  k.getItem() += g*power (y, i); 
+	  e -= mulMod2 (g*power(y, i), l.getItem(), yToD);
+          DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " << 
+                    mod (e, power (y, i + 1)));
+	} 
+      }  
+      if (e.isZero())
+        break;
+    }
+
+    DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y));
+
+#ifdef DEBUGOUTPUT
+    CanonicalForm test= 0;
+    j= p; 
+    for (CFListIterator i= result; i.hasItem(); i++, j++)
+      test += mod (i.getItem()*j.getItem(), power (y, d));
+    DEBOUTLN (cerr, "test= " << test);
+#endif
+    return result;
+  }
+}
+
+static inline
+CFList
+multiRecDiophantine (const CanonicalForm& F, const CFList& factors, 
+                     const CFList& recResult, const CFList& M, const int d)
+{
+  Variable y= F.mvar();
+  CFList result;
+  CFListIterator i;
+  CanonicalForm e= 1;
+  CFListIterator j= factors;
+  CFList p;
+  CFArray bufFactors= CFArray (factors.length());
+  CanonicalForm yToD= power (y, d);
+  int k= 0;
+  for (CFListIterator i= factors; i.hasItem(); i++, k++)
+    bufFactors [k]= i.getItem();
+  CanonicalForm b;
+  CFList buf= M;
+  buf.removeLast();
+  buf.append (yToD);
+  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
+  {
+    b= 1;
+    if (fdivides (bufFactors[k], F))
+      b= F/bufFactors[k];
+    else 
+    {
+      for (int l= 0; l < factors.length(); l++)
+      {
+	if (l == k)
+	  continue;
+	else
+	{
+	  b= mulMod (b, bufFactors[l], buf);
+	}
+      }
+    }
+    p.append (b);
+  }  
+  j= p;
+  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
+    e -= mulMod (i.getItem(), j.getItem(), M);
+
+  if (e.isZero())
+    return recResult; 
+  CanonicalForm coeffE;
+  CFList s;
+  result= recResult;
+  CanonicalForm g;
+  for (int i= 1; i < d; i++)
+  {
+    if (degree (e, y) > 0)
+      coeffE= e[i];
+    else
+      coeffE= 0;
+    if (!coeffE.isZero())
+    {
+      CFListIterator k= result;
+      CFListIterator l= p;
+      j= recResult;
+      int ii= 0;
+      CanonicalForm dummy;
+      for (; j.hasItem(); j++, k++, l++, ii++)
+      {
+	g= mulMod (coeffE, j.getItem(), M);
+	if (degree (bufFactors[ii], y) <= 0)
+	  divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy, 
+                  g, M);
+	else
+	  divrem (g, bufFactors[ii][0], dummy, g, M);
+	k.getItem() += g*power (y, i);
+	e -= mulMod (g*power (y, i), l.getItem(), M);
+      } 
+    }
+ 
+    if (e.isZero())
+      break;
+  }
+
+#ifdef DEBUGOUTPUT
+    CanonicalForm test= 0;
+    j= p; 
+    for (CFListIterator i= result; i.hasItem(); i++, j++)
+      test += mod (i.getItem()*j.getItem(), power (y, d));
+    DEBOUTLN (cerr, "test= " << test);
+#endif
+  return result;
+}
+
+static inline
+void 
+henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors, 
+            const CFList& diophant, CFMatrix& M, CFArray& Pi, int j, 
+            const CFList& MOD) 
+{
+  CanonicalForm E;
+  CanonicalForm xToJ= power (F.mvar(), j);
+  Variable x= F.mvar();
+  // compute the error
+  if (j == 1)
+  {
+    E= F[j];
+#ifdef DEBUGOUTPUT 
+    CanonicalForm test= 1;
+    for (int i= 0; i < factors.length(); i++)
+    {
+      if (i == 0)
+        test= mulMod (test, mod (bufFactors [i], xToJ), MOD);
+      else
+        test= mulMod (test, bufFactors[i], MOD);
+    }
+    CanonicalForm test2= mod (F-test, xToJ);
+
+    test2= mod (test2, MOD);
+    DEBOUTLN (cerr, "test= " << test2);
+#endif
+  } 
+  else
+  {
+#ifdef DEBUGOUTPUT 
+    CanonicalForm test= 1;
+    for (int i= 0; i < factors.length(); i++)
+    {
+      if (i == 0)
+        test *= mod (bufFactors [i], power (x, j));
+      else
+        test *= bufFactors[i];
+    }
+    test= mod (test, power (x, j));
+    test= mod (test, MOD);
+    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
+    DEBOUTLN (cerr, "test= " << test2);
+#endif
+
+    if (degree (Pi [factors.length() - 2], x) > 0)
+      E= F[j] - Pi [factors.length() - 2] [j];
+    else
+      E= F[j];
+  }
+
+  CFArray buf= CFArray (diophant.length());
+  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
+  int k= 0;
+  // actual lifting
+  CanonicalForm dummy, rest1;
+  for (CFListIterator i= diophant; i.hasItem(); i++, k++) 
+  {
+    if (degree (bufFactors[k], x) > 0)
+    {
+      if (k > 0)
+        divrem (E, bufFactors[k] [0], dummy, rest1, MOD);
+      else
+        rest1= E;
+    }
+    else
+      divrem (E, bufFactors[k], dummy, rest1, MOD);
+
+    buf[k]= mulMod (i.getItem(), rest1, MOD);
+
+    if (degree (bufFactors[k], x) > 0)
+      divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD);
+    else 
+      divrem (buf[k], bufFactors[k], dummy, buf[k], MOD);
+  }
+  for (k= 1; k < factors.length(); k++)
+    bufFactors[k] += xToJ*buf[k];
+
+  // update Pi [0]
+  int degBuf0= degree (bufFactors[0], x);
+  int degBuf1= degree (bufFactors[1], x);
+  if (degBuf0 > 0 && degBuf1 > 0)
+    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
+  CanonicalForm uIZeroJ;
+  if (j == 1)
+  {
+    if (degBuf0 > 0 && degBuf1 > 0)
+      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
+                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
+    else if (degBuf0 > 0)
+      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
+    else if (degBuf1 > 0)
+      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
+    else
+      uIZeroJ= 0;
+    Pi [0] += xToJ*uIZeroJ; 
+  } 
+  else
+  {
+    if (degBuf0 > 0 && degBuf1 > 0)
+      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]), 
+                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
+    else if (degBuf0 > 0)
+      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
+    else if (degBuf1 > 0)
+      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
+    else
+      uIZeroJ= 0;
+    Pi [0] += xToJ*uIZeroJ; 
+  }
+
+  CFArray tmp= CFArray (factors.length() - 1);
+  for (k= 0; k < factors.length() - 1; k++)
+    tmp[k]= 0;
+  CFIterator one, two;
+  one= bufFactors [0];
+  two= bufFactors [1];
+  if (degBuf0 > 0 && degBuf1 > 0)
+  {
+    for (k= 1; k <= (int) ceil (j/2.0); k++)
+    {  
+      if (k != j - k + 1)
+      {
+	if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
+	{
+	  tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()), 
+                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) - 
+                    M (j - k + 2, 1);
+	  one++;
+	  two++;
+	}
+	else if (one.exp() == j - k + 1)
+	{
+	  tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()), 
+                    bufFactors[1] [k], MOD) - M (k + 1, 1);
+	  one++;
+	}
+	else if (two.exp() == j - k + 1)
+	{
+	  tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] + 
+                    two.coeff()), MOD) - M (k + 1, 1);
+	  two++;
+	}
+      }
+      else
+      {
+	tmp[0] += M (k + 1, 1);
+      }
+    }
+  }
+  Pi [0] += tmp[0]*xToJ*F.mvar();
+
+  // update Pi [l] 
+  int degPi, degBuf;
+  for (int l= 1; l < factors.length() - 1; l++) 
+  {
+    degPi= degree (Pi [l - 1], x);
+    degBuf= degree (bufFactors[l + 1], x);
+    if (degPi > 0 && degBuf > 0)
+      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);  
+    if (j == 1)
+    {
+      if (degPi > 0 && degBuf > 0)
+        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]), 
+                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)- 
+                  M (1, l + 1));
+      else if (degPi > 0)
+        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
+      else if (degBuf > 0)
+        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
+    }
+    else
+    {
+      if (degPi > 0 && degBuf > 0)
+      {
+        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
+        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD); 
+      }
+      else if (degPi > 0)
+        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
+      else if (degBuf > 0)
+      {
+        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
+        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
+      }
+      Pi[l] += xToJ*uIZeroJ;
+    }
+    one= bufFactors [l + 1];
+    two= Pi [l - 1];
+    if (two.exp() == j + 1)
+    {
+      if (degBuf > 0 && degPi > 0)
+      {
+	  tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
+	  two++;
+      }
+      else if (degPi > 0)
+      {      
+	  tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
+	  two++;
+      } 
+    }
+    if (degBuf > 0 && degPi > 0) 
+    {
+      for (k= 1; k <= (int) ceil (j/2.0); k++)
+      {  
+	if (k != j - k + 1)
+	{
+	  if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
+	  {
+	    tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()), 
+                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) - 
+                      M (j - k + 2, l + 1);
+	    one++;
+	    two++;
+	  }
+	  else if (one.exp() == j - k + 1)
+	  {
+	    tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()), 
+                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);     
+	    one++;
+	  }
+	  else if (two.exp() == j - k + 1)
+	  {
+	    tmp[l] += mulMod (bufFactors[l + 1] [k], 
+                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
+	    two++;
+	  }
+	}
+	else
+	  tmp[l] += M (k + 1, l + 1);
+      }  
+    }  
+    Pi[l] += tmp[l]*xToJ*F.mvar();
+  }
+
+  return;
+}
+
+CFList 
+henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
+              diophant, CFArray& Pi, CFMatrix& M)
+{
+  CFList buf= factors;
+  int k= 0;
+  int liftBound;
+  int liftBoundBivar= l[k];
+  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
+  CFList MOD;
+  MOD.append (power (Variable (2), liftBoundBivar));
+  CFArray bufFactors= CFArray (factors.length());
+  k= 0;
+  CFListIterator j= eval;
+  j++;
+  buf.removeFirst();
+  buf.insert (LC (j.getItem(), 1));
+  for (CFListIterator i= buf; i.hasItem(); i++, k++)
+    bufFactors[k]= i.getItem();
+  Pi= CFArray (factors.length() - 1);
+  CFListIterator i= buf;
+  i++;
+  Variable y= j.getItem().mvar();
+  Pi [0]= mulMod (i.getItem(), mod (buf.getFirst(), y), MOD);
+  M (1, 1)= Pi [0];
+  k= 1;
+  if (i.hasItem())
+    i++;
+  for (i; i.hasItem(); i++, k++) 
+  {
+    Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
+    M (1, k + 1)= Pi [k];
+  }
+
+  for (int d= 1; d < l[1]; d++)
+    henselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, d, MOD);
+  CFList result;
+  for (k= 1; k < factors.length(); k++)
+    result.append (bufFactors[k]);
+  return result;
+} 
+
+void 
+henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end, 
+                  CFArray& Pi, const CFList& diophant, CFMatrix& M, 
+                  const CFList& MOD)
+{
+  CFArray bufFactors= CFArray (factors.length());
+  int i= 0;
+  CanonicalForm xToStart= power (F.mvar(), start);
+  for (CFListIterator k= factors; k.hasItem(); k++, i++)
+  {
+    if (i == 0)
+      bufFactors[i]= mod (k.getItem(), xToStart);
+    else
+      bufFactors[i]= k.getItem();
+  }
+  for (i= start; i < end; i++)
+    henselStep (F, factors, bufFactors, diophant, M, Pi, i, MOD); 
+  
+  CFListIterator k= factors;
+  for (i= 0; i < factors.length(); k++, i++)
+    k.getItem()= bufFactors [i];
+  factors.removeFirst();
+  return; 
+} 
+
+CFList
+henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
+            diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
+            lNew)
+{
+  diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
+  int k= 0;
+  CFArray bufFactors= CFArray (factors.length());
+  for (CFListIterator i= factors; i.hasItem(); i++, k++)
+  {
+    if (k == 0)
+      bufFactors[k]= LC (F.getLast(), 1);
+    else
+      bufFactors[k]= i.getItem();
+  }
+  CFList buf= factors;
+  buf.removeFirst();
+  buf.insert (LC (F.getLast(), 1));
+  CFListIterator i= buf;
+  i++;
+  Variable y= F.getLast().mvar();
+  Variable x= F.getFirst().mvar();
+  CanonicalForm xToLOld= power (x, lOld);
+  Pi [0]= mod (Pi[0], xToLOld);
+  M (1, 1)= Pi [0];
+  k= 1;
+  if (i.hasItem())
+    i++;
+  for (i; i.hasItem(); i++, k++) 
+  {
+    Pi [k]= mod (Pi [k], xToLOld);
+    M (1, k + 1)= Pi [k];
+  }
+
+  for (int d= 1; d < lNew; d++)
+    henselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, d, MOD);
+  CFList result;
+  for (k= 1; k < factors.length(); k++)
+    result.append (bufFactors[k]);
+  return result;
+}
+
+CFList
+henselLift (const CFList& eval, const CFList& factors, const int* l, const int
+            lLength)
+{
+  CFList diophant; 
+  CFList buf= factors; 
+  buf.insert (LC (eval.getFirst(), 1));
+  sortList (buf, Variable (1));
+  CFArray Pi;
+  CFMatrix M= CFMatrix (l[1], factors.length());
+  CFList result= henselLift23 (eval, buf, l, diophant, Pi, M);
+  if (eval.length() == 2)
+    return result;
+  CFList MOD;
+  for (int i= 0; i < 2; i++)
+    MOD.append (power (Variable (i + 2), l[i]));
+  CFListIterator j= eval;
+  j++;
+  CFList bufEval;
+  bufEval.append (j.getItem());
+  j++;
+ 
+  for (int i= 2; i <= lLength && j.hasItem(); i++, j++)
+  {
+    result.insert (LC (bufEval.getFirst(), 1));
+    bufEval.append (j.getItem());
+    M= CFMatrix (l[i], factors.length());
+    result= henselLift (bufEval, result, MOD, diophant, Pi, M, l[i - 1], l[i]);
+    MOD.append (power (Variable (i + 2), l[i]));
+    bufEval.removeFirst();
+  }
+  return result; 
+}
+
+#endif
+/* HAVE_NTL */
+
Index: factory/facHensel.h
===================================================================
--- factory/facHensel.h	(revision ad3c3ff7e9858e8efdced4631a25aebfcf700f36)
+++ factory/facHensel.h	(revision ad3c3ff7e9858e8efdced4631a25aebfcf700f36)
@@ -0,0 +1,368 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file facHensel.h
+ * 
+ * This file defines functions for Hensel lifting and modular multiplication. 
+ *
+ * ABSTRACT: function are used for bi- and multivariate factorization over 
+ * finite fields
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#ifndef FAC_HENSEL_H
+#define FAC_HENSEL_H
+
+#include <config.h>
+#include "assert.h"
+#include "debug.h"
+#include "timing.h"
+
+#include "canonicalform.h"
+#include "cf_iter.h"	
+#include "ftmpl_functions.h"
+#include "algext.h"
+
+CanonicalForm mul (const CanonicalForm& F, const CanonicalForm& G);
+CanonicalForm mulMod3 (const CanonicalForm& F, const CanonicalForm& G, const CFList& MOD);
+/// multiplication of univariate polys over a finite field using NTL, if we are 
+/// in GF factory's default multiplication is used.
+///
+/// @return @a mulNTL returns F*G
+static inline
+CanonicalForm
+mulNTL (const CanonicalForm& F, ///< [in] a univariate poly 
+        const CanonicalForm& G  ///< [in] a univariate poly
+       );
+
+/// mod of univariate polys over a finite field using NTL, if we are 
+/// in GF factory's default mod is used.
+///
+/// @return @a modNTL returns F mod G
+static inline
+CanonicalForm
+modNTL (const CanonicalForm& F, ///< [in] a univariate poly 
+        const CanonicalForm& G  ///< [in] a univariate poly
+       );
+
+/// division of univariate polys over a finite field using NTL, if we are 
+/// in GF factory's default division is used.
+///
+/// @return @a divNTL returns F/G
+static inline
+CanonicalForm
+divNTL (const CanonicalForm& F, ///< [in] a univariate poly 
+        const CanonicalForm& G  ///< [in] a univariate poly
+       );
+
+/*/// division with remainder of univariate polys over a finite field using NTL, 
+/// if we are in GF factory's default division with remainder is used.
+static inline
+void
+divremNTL (const CanonicalForm& F, ///< [in] a univariate poly 
+           const CanonicalForm& G, ///< [in] a univariate poly
+           CanonicalForm& Q,       ///< [in,out] dividend
+           CanonicalForm& R        ///< [in,out] remainder
+          );*/
+
+/// splits @a F into degree (F, x)/m polynomials each of degree less than @a m 
+/// in @a x.
+///
+/// @return @a split returns a list of polynomials of degree less than @a m in
+///         @a x. If degree (F, x) <= 0, F is returned.  
+/// @sa divrem32(), divrem21()
+static inline
+CFList split (const CanonicalForm& F, ///< [in] some poly
+              const int m,            ///< [in] some int
+              const Variable& x       ///< [in] some Variable
+             );      
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a M.
+/// 
+/// @sa divrem(), divrem21(), divrem2()
+static inline
+void divrem32 (const CanonicalForm& F, ///< [in] poly, s.t. 3*(degree (G, 1)/2)> 
+                                       ///< degree (F, 1)
+               const CanonicalForm& G, ///< [in] some poly
+               CanonicalForm& Q,       ///< [in,out] dividend
+               CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) < 
+                                       ///< degree (G, 1)
+               const CFList& M         ///< [in] only contains powers of
+                                       ///< Variables of level higher than 1
+              );
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a M.
+/// 
+/// @sa divrem(), divrem32(), divrem2()
+static inline
+void divrem21 (const CanonicalForm& F, ///< [in] poly, s.t. 2*degree (G, 1) >
+                                       ///< degree (F, 1)
+               const CanonicalForm& G, ///< [in] some poly
+               CanonicalForm& Q,       ///< [in,out] dividend
+               CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) < 
+                                       ///< degree (G, 1)
+               const CFList& M         ///< [in] only contains powers of
+                                       ///< Variables of level higher than 1
+              );
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a M. 
+/// 
+/// @return @a Q returns the dividend, @a R returns the remainder.
+/// @sa divrem()
+void divrem2 (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial 
+              const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
+              CanonicalForm& Q,       ///< [in,out] dividend
+              CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) < 
+                                      ///< degree (G, 1)
+              const CanonicalForm& M  ///< [in] power of Variable (2)
+             );
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a MOD.
+///
+/// @sa divrem2()
+void divrem (
+           const CanonicalForm& F, ///< [in] multivariate, compressed polynomial
+           const CanonicalForm& G, ///< [in] multivariate, compressed polynomial
+           CanonicalForm& Q,       ///< [in,out] dividend
+           CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) < 
+                                   ///< degree (G, 1)
+           const CFList& MOD       ///< [in] only contains powers of
+                                   ///< Variables of level higher than 1
+            );
+
+/// reduce @a F modulo elements in @a M.
+///
+/// @return @a mod returns @a F modulo @a M
+CanonicalForm mod (const CanonicalForm& F, ///< [in] compressed polynomial
+                   const CFList& M         ///< [in] list containing only 
+                                           ///< univariate polynomials
+                  );
+
+/// Karatsuba style modular multiplication for bivariate polynomials.
+///
+/// @return @a mulMod2 returns @a A * @a B mod @a M.
+CanonicalForm 
+mulMod2 (const CanonicalForm& A, ///< [in] bivariate, compressed polynomial
+         const CanonicalForm& B, ///< [in] bivariate, compressed polynomial
+         const CanonicalForm& M  ///< [in] power of Variable (2)
+        );
+
+/// Karatsuba style modular multiplication for multivariate polynomials.
+///
+/// @return @a mulMod2 returns @a A * @a B mod @a MOD.
+CanonicalForm 
+mulMod (const CanonicalForm& A, ///< [in] multivariate, compressed polynomial
+        const CanonicalForm& B, ///< [in] multivariate, compressed polynomial
+        const CFList& MOD       ///< [in] only contains powers of
+                                ///< Variables of level higher than 1
+       );
+
+/// product of all elements in @a L modulo @a M via divide-and-conquer.
+///
+/// @return @a prodMod returns product of all elements in @a L modulo @a M.
+CanonicalForm 
+prodMod (const CFList& L,       ///< [in] contains only bivariate, compressed 
+                                ///< polynomials
+         const CanonicalForm& M ///< [in] power of Variable (2)
+        );
+
+/// product of all elements in @a L modulo @a M via divide-and-conquer.
+///
+/// @return @a prodMod returns product of all elements in @a L modulo @a M.
+CanonicalForm 
+prodMod (const CFList& L, ///< [in] contains multivariate, compressed 
+                          ///< polynomials
+         const CFList& M  ///< [in] contains only powers of Variables
+        );
+ 
+/// sort a list of polynomials by their degree in @a x.
+///
+void sortList (CFList& list,     ///< [in, out] list of polys, sorted list 
+               const Variable& x ///< [in] some Variable 
+              );
+
+/// solve the univariate diophantine equation 
+/// \f$ 1\equiv \sum_{i= 1}^{r} {\delta_{i} F/g_{i}} \f$.
+/// Where \f$ F= \prod_{i= 1}^{r} {g_{i}} \f$ and \f$ F \f$ is squarefree 
+/// the \f$ \delta_{i} \f$ have degree less than the degree of \f$ g_{i} \f$.
+///
+/// @return @a diophantine returns a list of polynomials \f$ \delta_{i} \f$ as
+/// specified above
+/// @sa biDiophantine(), multiRecDiophantine()
+static inline
+CFList diophantine (const CanonicalForm& F, ///< [in] compressed, bivariate 
+                                            ///< polynomial
+                    const CFList& factors   ///< [in] a list of factors, as 
+                                            ///< specified above, including
+                                            ///< LC (F, Variable (1)) as first 
+                                            ///< element
+                   );
+
+/// Hensel lift from univariate to bivariate.
+/// 
+/// @sa henselLiftResume12(), henselLift23(), henselLiftResume(), henselLift()
+void 
+henselLift12 (const CanonicalForm& F, ///< [in] compressed, bivariate poly
+              CFList& factors,        ///< [in, out] monic univariate factors of  
+                                      ///< F, including leading coefficient as
+                                      ///< first element. Returns monic lifted 
+                                      ///< factors without the leading 
+                                      ///< coefficient
+              int l,                  ///< [in] lifting precision
+              CFArray& Pi,            ///< [in,out] stores intermediate results
+              CFList& diophant,       ///< [in,out] result of diophantine()
+              CFMatrix& M             ///< [in,out] stores intermediate results
+             );
+
+/// resume Hensel lift from univariate to bivariate. Assumes factors are lifted 
+/// to precision Variable (2)^start and lifts them to precision Variable (2)^end
+/// 
+/// @sa henselLift12(), henselLift23(), henselLiftResume(), henselLift()
+void 
+henselLiftResume12 (const CanonicalForm& F, ///< [in] compressed, bivariate poly
+                    CFList& factors,        ///< [in,out] monic factors of F, 
+                                            ///< lifted to precision start, 
+                                            ///< including leading coefficient 
+                                            ///< as first element. Returns monic  
+                                            ///< lifted factors without the  
+                                            ///< leading coefficient
+                    int start,              ///< [in] starting precision
+                    int end,                ///< [in] end precision
+                    CFArray& Pi,            ///< [in,out] stores intermediate
+                                            ///< results
+                    const CFList& diophant, ///< [in] result of diophantine
+                    CFMatrix& M             ///< [in,out] stores intermediate
+                                            ///< results
+                   );
+
+/// solves the bivariate polynomial diophantine equation
+/// \f$ 1\equiv \sum_{i= 1}^{r} {\delta_{i} F/g_{i}} \ mod\ y^{d} \f$,
+/// where \f$ F= \prod_{i= 1}^{r} {g_{i}} \ mod\ y^{d}\f$ and 
+/// \f$ F \in K[x][y]\f$ is squarefree, the \f$ \delta_{i} \f$ have degree less
+/// than the degree of \f$ g_{i} \f$ in x.
+///
+/// @return @a biDiophantine returns a list of polynomials \f$ \delta_{i} \f$ as
+///         specified above
+/// @sa diophantine(), multiRecDiophantine()
+static inline
+CFList
+biDiophantine (const CanonicalForm& F, ///< [in] compressed, bivariate poly
+               const CFList& factors,  ///< [in] list of monic bivariate factors
+                                       ///< including LC (F, Variable (1)) as  
+                                       ///< first element
+               const int d             ///< [in] precision
+              );
+
+/// solve the multivariate polynomial diophantine equation
+/// \f$ 1\equiv \sum_{i= 1}^{r} {\delta_{i} F/g_{i}} \ mod\ <M,F.mvar()^{d}>\f$,
+/// where \f$ F= \prod_{i= 1}^{r} {g_{i}} \ mod\ <M,F.mvar()^{d}>\f$ and 
+/// \f$ F \in K[x][x_1,\ldots , x_n]\f$ is squarefree, the \f$ \delta_{i} \f$ 
+/// have degree less than the degree of \f$ g_{i} \f$ in x. 
+/// 
+/// @return @a multiDiophantine returns a list of polynomials \f$ \delta_{i} \f$ 
+///         as specified above
+/// @sa diophantine(), biDiophantine()
+static inline
+CFList
+multiRecDiophantine (
+               const CanonicalForm& F,   ///< [in] compressed, 
+                                         ///< multivariate polynomial
+               const CFList& factors,    ///< [in] list of monic factors
+                                         ///< including LC (F, Variable (1)) as
+                                         ///< first element
+               const CFList& recResult, ///< [in] result of above equation
+                                         ///< modulo M
+               const CFList& M,          ///< [in] a list of powers of Variables
+                                         ///< of level higher than 1
+               const int d               ///< [in] precision
+                    );
+
+/// Hensel lifting from bivariate to trivariate.
+///
+/// @return @a henselLift23 returns a list of polynomials lifted to precision
+///          Variable (3)^l[1]
+/// @sa henselLift12(), henselLiftResume12(), henselLiftResume(), henselLift()
+CFList 
+henselLift23 (const CFList& eval,    ///< [in] contains compressed, bivariate
+                                     ///< as first element and trivariate one as
+                                     ///< second element 
+              const CFList& factors, ///< [in] monic bivariate factors, 
+                                     ///< including leading coefficient 
+                                     ///< as first element.
+              const int* l,          ///< [in] l[0]: precision of bivariate 
+                                     ///< lifting, l[1] as above
+              CFList& diophant,      ///< [in,out] returns the result of
+                                     ///< biDiophantine()
+              CFArray& Pi,           ///< [in,out] stores intermediate results
+              CFMatrix& M            ///< [in,out] stores intermediate results
+             );
+
+/// resume Hensel lifting.
+///
+/// @sa henselLift12(), henselLiftResume12(), henselLift23(), henselLift()
+void 
+henselLiftResume (
+              const CanonicalForm& F, ///< [in] compressed, multivariate poly 
+              CFList& factors,        ///< [in,out] monic multivariate factors
+                                      ///< lifted to precision F.mvar()^start, 
+                                      ///< including leading coefficient 
+                                      ///< as first element. Returns factors
+                                      ///< lifted to precision F.mvar()^end
+              int start,              ///< [in] starting precision 
+              int end,                ///< [in] end precision
+              CFArray& Pi,            ///< [in,out] stores intermediate results
+              const CFList& diophant, ///< [in] result of multiRecDiophantine()
+              CFMatrix& M,            ///< [in, out] stores intermediate results
+              const CFList& MOD       ///< [in] a list of powers of Variables
+                                      ///< of level higher than 1
+                 );
+
+/// Hensel lifting
+///
+/// @return @a henselLift returns a list of polynomials lifted to 
+///          precision F.getLast().mvar()^l_new
+/// @sa henselLift12(), henselLiftResume12(), henselLift23(), henselLiftResume()
+CFList
+henselLift (const CFList& F,       ///< [in] two compressed, multivariate
+                                   ///< polys F and G
+            const CFList& factors, ///< [in] monic multivariate factors                                      
+                                   ///< including leading coefficient 
+                                   ///< as first element. 
+            const CFList& MOD,     ///< [in] a list of powers of Variables
+                                   ///< of level higher than 1
+            CFList& diophant,      ///< [in,out] result of multiRecDiophantine()
+            CFArray& Pi,           ///< [in,out] stores intermediate results
+            CFMatrix& M,           ///< [in,out] stores intermediate results
+            const int lOld,       ///< [in] lifting precision of F.mvar() 
+            const int lNew        ///< [in] lifting precision of G.mvar()
+           );
+
+/// Hensel lifting from bivariate to multivariate 
+///
+/// @return @a henselLift returns a list of lifted factors 
+/// @sa henselLift12(), henselLiftResume12(), henselLift23(), henselLiftResume()         
+CFList
+henselLift (const CFList& eval,    ///< [in] a list of polynomials the last 
+                                   ///< element is a compressed multivariate 
+                                   ///< poly, last but one element equals the 
+                                   ///< last elements modulo its main variable
+                                   ///< and so on. The first element is a 
+                                   ///< compressed bivariate poly.
+            const CFList& factors, ///< [in] bivariate factors, including 
+                                   ///< leading coefficient
+            const int* l,          ///< [in] lifting bounds
+            const int lLength     ///< [in] length of l
+           );
+
+#endif
+/* FAC_HENSEL_H */
+
