Index: factory/cf_gcd_smallp.cc
===================================================================
--- factory/cf_gcd_smallp.cc	(revision ce35eb186c60fc2922d2d106228a4d3a2cb87a5b)
+++ factory/cf_gcd_smallp.cc	(revision 139b671c6c11958a1f538c95dc0ade26520e76b0)
@@ -904,7 +904,5 @@
   bool inextensionextension= false;
   top_level= false;
-  CanonicalForm CF_buf;
   CFList source, dest;
-  CanonicalForm gcdcheck;
   do 
   {
@@ -942,5 +940,5 @@
       CFList list;
       CanonicalForm mipo;
-      int deg= 3;
+      int deg= 2;
       do {
         mipo= randomIrredpoly (deg, x); 
Index: factory/cf_map_ext.cc
===================================================================
--- factory/cf_map_ext.cc	(revision ce35eb186c60fc2922d2d106228a4d3a2cb87a5b)
+++ factory/cf_map_ext.cc	(revision 139b671c6c11958a1f538c95dc0ade26520e76b0)
@@ -340,7 +340,9 @@
       return 0;
   } while (1);
-  zz_pE::init (NTL_mipo);
-  zz_pEX NTL_alpha_mipo= convertFacCF2NTLzz_pEX (mipo, NTL_mipo);
-  zz_pE root= FindRoot (NTL_alpha_mipo);
+  
+  zz_pX alpha_mipo= convertFacCF2NTLzzpX (mipo);
+  zz_pE::init (alpha_mipo);
+  zz_pEX NTL_beta_mipo= to_zz_pEX (NTL_mipo);
+  zz_pE root= FindRoot (NTL_beta_mipo);
   return convertNTLzzpE2CF (root, alpha);
 }
Index: factory/facHensel.cc
===================================================================
--- factory/facHensel.cc	(revision ce35eb186c60fc2922d2d106228a4d3a2cb87a5b)
+++ factory/facHensel.cc	(revision 139b671c6c11958a1f538c95dc0ade26520e76b0)
@@ -26,8 +26,12 @@
 #include "facHensel.h"
 #include "cf_util.h"
+#include "fac_util.h"
+#include "cf_algorithm.h"
 
 #ifdef HAVE_NTL
 #include <NTL/lzz_pEX.h>
 #include "NTLconvert.h"
+
+CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c);
 
 static inline
@@ -920,7 +924,8 @@
 void 
 henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi, 
-              CFList& diophant, CFMatrix& M) 
-{
-  sortList (factors, Variable (1));
+              CFList& diophant, CFMatrix& M, bool sort) 
+{
+  if (sort)
+    sortList (factors, Variable (1));
   Pi= CFArray (factors.length() - 1);
   CFListIterator j= factors;
@@ -933,5 +938,5 @@
   if (j.hasItem())
     j++;
-  for (j; j.hasItem(); j++, i++) 
+  for (; j.hasItem(); j++, i++) 
   {
     Pi [i]= mulNTL (Pi [i - 1], j.getItem());
@@ -954,4 +959,254 @@
     k.getItem()= bufFactors[i];
   factors.removeFirst();
+  return;
+}
+
+void 
+henselStep122 (const CanonicalForm& F, const CFList& factors, 
+              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
+              CFArray& Pi, int j, const CFArray& LCs)
+{
+  Variable x= F.mvar();
+  CanonicalForm xToJ= power (x, j);
+  
+  CanonicalForm E;
+  // compute the error
+  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());
+
+  int k= 0;
+  CanonicalForm remainder;
+  // actual lifting
+  for (CFListIterator i= diophant; i.hasItem(); i++, k++) 
+  {
+    if (degree (bufFactors[k], x) > 0)
+      remainder= modNTL (E, bufFactors[k] [0]);
+    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= 0; 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]);
+    if (j + 2 <= M.rows())
+      M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
+  }
+  CanonicalForm uIZeroJ;
+  if (degBuf0 > 0 && degBuf1 > 0)
+    uIZeroJ= mulNTL (bufFactors[0] [0], buf [1]) + mulNTL (bufFactors[1] [0], buf[0]);
+  else if (degBuf0 > 0)
+    uIZeroJ= mulNTL (buf[0], 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)
+  {
+    while (one.exp() > j) one++;
+    while (two.exp() > j) two++;
+    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);
+    }
+
+    if (j + 2 <= M.rows())
+    {
+      if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
+        tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]), (bufFactors [1] [j + 1] + bufFactors [1] [0]))
+               - M(1,1) - M (j + 2,1);
+      else if (degBuf0 >= j + 1)
+        tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
+      else if (degBuf1 >= j + 1)
+        tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 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 
+henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi, 
+              CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort) 
+{
+  if (sort)
+    sortList (factors, Variable (1));
+  Pi= CFArray (factors.length() - 2);
+  CFList bufFactors2= factors;
+  bufFactors2.removeFirst();
+  diophant.removeFirst();
+  CFListIterator iter= diophant;
+  CanonicalForm test= 0, test2= 0;
+  CanonicalForm s,t;
+  extgcd (bufFactors2.getFirst(), bufFactors2.getLast(), s, t);
+  diophant= CFList();
+  diophant.append (t);
+  diophant.append (s);
+  DEBOUTLN (cerr, "diophant= " << diophant);
+
+  CFArray bufFactors= CFArray (bufFactors2.length());
+  int i= 0;
+  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
+    bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
+
+  Variable x= F.mvar();
+  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
+  {
+    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
+    Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) + 
+                        mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
+  }
+  else if (degree (bufFactors[0], x) > 0)
+  {
+    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
+    Pi [0]= M (1, 1) +
+            mulNTL (bufFactors [0] [1], bufFactors[1])*x;
+  }
+  else if (degree (bufFactors[1], x) > 0)
+  {
+    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
+    Pi [0]= M (1, 1) +
+            mulNTL (bufFactors [0], bufFactors[1] [1])*x;
+  }
+  else
+  {
+    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
+    Pi [0]= M (1, 1);
+  }          
+
+  for (i= 1; i < l; i++)
+    henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
+
+  factors= CFList();
+  for (i= 0; i < bufFactors.size(); i++)
+    factors.append (bufFactors[i]);
   return;
 }
@@ -1133,4 +1388,5 @@
   for (int i= 1; i < d; i++)
   {
+    CanonicalForm kk= 0;
     if (degree (e, y) > 0)
       coeffE= e[i];
@@ -1156,5 +1412,166 @@
       } 
     }
- 
+    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;
+}
+
+/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$ 
+/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
+static inline
+CFList 
+diophantine (const CFList& recResult, const CFList& factors, 
+             const CFList& products, const CFList& M, const CanonicalForm& E)
+{
+  if (M.isEmpty())
+  {
+    CFList result;
+    CFListIterator j= factors;
+    CanonicalForm buf;
+    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
+    {
+      ASSERT (E.isUnivariate() || E.inCoeffDomain(), 
+              "constant or univariate poly expected");
+      ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(), 
+              "constant or univariate poly expected");
+      ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(), 
+              "constant or univariate poly expected");
+      buf= mulNTL (E, i.getItem());
+      result.append (modNTL (buf, j.getItem()));
+    }
+    return result;
+  }
+  Variable y= M.getLast().mvar();
+  CFList bufFactors= factors;
+  for (CFListIterator i= bufFactors; i.hasItem(); i++)
+    i.getItem()= mod (i.getItem(), y);
+  CFList bufProducts= products;
+  for (CFListIterator i= bufProducts; i.hasItem(); i++)
+    i.getItem()= mod (i.getItem(), y);
+  CFList buf= M;
+  buf.removeLast();
+  CanonicalForm bufE= mod (E, y);
+  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf, 
+                                      bufE);
+
+  CanonicalForm e= E;
+  CFListIterator j= products;
+  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
+    e -= mulMod (i.getItem(), j.getItem(), M);
+
+  CFList result= recDiophantine;
+  int d= degree (M.getLast());
+  CanonicalForm coeffE;
+  for (int i= 1; i < d; i++)
+  {
+    if (degree (e, y) > 0)
+      coeffE= e[i];
+    else
+      coeffE= 0;
+    if (!coeffE.isZero())
+    {
+      CFListIterator k= result;
+      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf, 
+                                   coeffE);
+      CFListIterator l= products;
+      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
+      {
+        k.getItem() += j.getItem()*power (y, i);
+	e -= mulMod (j.getItem()*power (y, i), l.getItem(), M);
+      } 
+    }
+    if (e.isZero())
+      break;
+  }
+  return result;
+}
+
+// non monic case
+static inline 
+CFList
+multiRecDiophantine2 (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;
+  result= recResult;
+  CanonicalForm g;
+  buf.removeLast();
+  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;
+      CFList recDiophantine;
+      CFList buf2, buf3;
+      buf2= factors;
+      buf3= p;
+      for (CFListIterator iii= buf2; iii.hasItem(); iii++)
+        iii.getItem()= mod (iii.getItem(), y);
+      for (CFListIterator iii= buf3; iii.hasItem(); iii++)
+        iii.getItem()= mod (iii.getItem(), y);
+      recDiophantine= diophantine (recResult, buf2, buf3, buf, coeffE);
+      CFListIterator iter= recDiophantine;
+      for (; j.hasItem(); j++, k++, l++, ii++, iter++)
+      {
+	k.getItem() += iter.getItem()*power (y, i);
+	e -= mulMod (iter.getItem()*power (y, i), l.getItem(), M);
+      } 
+    }
     if (e.isZero())
       break;
@@ -1410,4 +1827,212 @@
 }
 
+// so far only tested for two! factor Hensel lifting
+static inline
+void 
+henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors, 
+            const CFList& diophant, CFMatrix& M, CFArray& Pi, 
+            const CFList& products, int j, const CFList& MOD) 
+{
+  CanonicalForm E;
+  CanonicalForm xToJ= power (F.mvar(), j);
+  Variable x= F.mvar();
+
+  // compute the error
+#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());
+
+  // actual lifting
+  CFList diophantine2= diophantine (diophant, factors, products, MOD, E);
+
+  int k= 0;
+  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
+  {
+    buf[k]= i.getItem();
+    bufFactors[k] += xToJ*i.getItem();
+  }
+
+  // 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);
+    if (j + 2 <= M.rows())
+      M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
+  }
+  CanonicalForm uIZeroJ;
+
+    if (degBuf0 > 0 && degBuf1 > 0)
+      uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) + 
+               mulMod (bufFactors[1] [0], buf[0], MOD);
+    else if (degBuf0 > 0)
+      uIZeroJ= mulMod (buf[0], 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)
+  {
+    while (one.exp() > j) one++;
+    while (two.exp() > j) two++;
+    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);
+      }
+    }
+    
+    if (j + 2 <= M.rows())
+    {
+      if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
+        tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]), 
+                           (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
+                   - M(1,1) - M (j + 2,1);
+      else if (degBuf0 >= j + 1)
+        tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
+      else if (degBuf1 >= j + 1)
+        tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
+    }
+  }
+  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&
@@ -1416,5 +2041,4 @@
   CFList buf= factors;
   int k= 0;
-  int liftBound;
   int liftBoundBivar= l[k];
   diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
@@ -1438,5 +2062,5 @@
   if (i.hasItem())
     i++;
-  for (i; i.hasItem(); i++, k++) 
+  for (; i.hasItem(); i++, k++) 
   {
     Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
@@ -1505,5 +2129,5 @@
   if (i.hasItem())
     i++;
-  for (i; i.hasItem(); i++, k++) 
+  for (; i.hasItem(); i++, k++) 
   {
     Pi [k]= mod (Pi [k], xToLOld);
@@ -1521,10 +2145,11 @@
 CFList
 henselLift (const CFList& eval, const CFList& factors, const int* l, const int
-            lLength)
+            lLength, bool sort)
 {
   CFList diophant; 
   CFList buf= factors; 
   buf.insert (LC (eval.getFirst(), 1));
-  sortList (buf, Variable (1));
+  if (sort)
+    sortList (buf, Variable (1));
   CFArray Pi;
   CFMatrix M= CFMatrix (l[1], factors.length());
@@ -1553,4 +2178,171 @@
 }
 
+// wrt. Variable (1)
+CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
+{
+  if (degree (F, 1) <= 0)
+   return c;
+  else
+  {
+    CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
+    result += (swapvar (c, Variable (F.level() + 1), Variable (1)) 
+              - LC (result))*power (result.mvar(), degree (result));
+    return swapvar (result, Variable (F.level() + 1), Variable (1));
+  } 
+}
+
+// so far only for two factor Hensel lifting
+CFList 
+henselLift232 (const CFList& eval, const CFList& factors, const int* l, CFList&
+              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1, 
+              const CFList& LCs2)
+{
+  CFList buf= factors;
+  int k= 0;
+  int liftBoundBivar= l[k];
+  CFList bufbuf= factors;
+  Variable v= Variable (2);
+
+  CFList MOD;
+  MOD.append (power (Variable (2), liftBoundBivar));
+  CFArray bufFactors= CFArray (factors.length());
+  k= 0;
+  CFListIterator j= eval;
+  j++;
+  CFListIterator iter1= LCs1;
+  CFListIterator iter2= LCs2;
+  iter1++;
+  iter2++;
+  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
+  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
+
+  CFListIterator i= buf;
+  i++;
+  Variable y= j.getItem().mvar();
+  if (y.level() != 3)
+    y= Variable (3);
+
+  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
+  M (1, 1)= Pi[0];
+  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
+    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) + 
+                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
+  else if (degree (bufFactors[0], y) > 0)
+    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
+  else if (degree (bufFactors[1], y) > 0)
+    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
+  
+  CFList products;
+  for (int i= 0; i < bufFactors.size(); i++)
+  {
+    if (degree (bufFactors[i], y) > 0)
+      products.append (M (1, 1)/bufFactors[i] [0]);
+    else
+      products.append (M (1, 1)/bufFactors[i]);
+  }
+  
+  for (int d= 1; d < l[1]; d++)
+    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD);
+  CFList result;
+  for (k= 0; k < factors.length(); k++)
+    result.append (bufFactors[k]);
+  return result;
+} 
+
+
+CFList
+henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
+            diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
+            lNew, const CFList& LCs1, const CFList& LCs2)
+{
+  int k= 0;
+  CFArray bufFactors= CFArray (factors.length());
+  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
+  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
+  CFList buf= factors;
+  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];
+
+  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
+    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) + 
+                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
+  else if (degree (bufFactors[0], y) > 0)
+    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
+  else if (degree (bufFactors[1], y) > 0)
+    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
+  
+  CFList products;
+  for (int i= 0; i < bufFactors.size(); i++)
+  {
+    if (degree (bufFactors[i], y) > 0)
+    {
+      ASSERT (fdivides (bufFactors[i] [0], M (1, 1)), "expected exact division");
+      products.append (M (1, 1)/bufFactors[i] [0]);
+    }
+    else
+    {
+      ASSERT (fdivides (bufFactors[i], M (1, 1)), "expected exact division");
+      products.append (M (1, 1)/bufFactors[i]);
+    }
+  }
+  
+  for (int d= 1; d < lNew; d++)
+    henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD);
+  CFList result;
+  for (k= 0; k < factors.length(); k++)
+    result.append (bufFactors[k]);
+  return result;
+}
+
+// so far only for two! factor Hensel lifting
+CFList
+henselLift2 (const CFList& eval, const CFList& factors, const int* l, const int
+             lLength, bool sort, const CFList& LCs1, const CFList& LCs2, 
+             const CFArray& Pi, const CFList& diophant)
+{
+  CFList bufDiophant= diophant;
+  CFList buf= factors;
+  if (sort)
+    sortList (buf, Variable (1));
+  CFArray bufPi= Pi;
+  CFMatrix M= CFMatrix (l[1], factors.length());
+  CFList result= henselLift232 (eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2);
+  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++;
+  CFListIterator jj= LCs1;
+  CFListIterator jjj= LCs2;
+  CFList bufLCs1, bufLCs2;
+  jj++, jjj++;
+  bufLCs1.append (jj.getItem());
+  bufLCs2.append (jjj.getItem());
+  jj++, jjj++;
+ 
+  for (int i= 2; i <= lLength && j.hasItem(); i++, j++, jj++, jjj++)
+  {
+    bufEval.append (j.getItem());
+    bufLCs1.append (jj.getItem());
+    bufLCs2.append (jjj.getItem());
+    M= CFMatrix (l[i], factors.length());
+    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1], 
+                         l[i], bufLCs1, bufLCs2);
+    MOD.append (power (Variable (i + 2), l[i]));
+    bufEval.removeFirst();
+    bufLCs1.removeFirst();
+    bufLCs2.removeFirst();
+  }
+  return result;
+}
+
 #endif
 /* HAVE_NTL */
Index: factory/facHensel.h
===================================================================
--- factory/facHensel.h	(revision ce35eb186c60fc2922d2d106228a4d3a2cb87a5b)
+++ factory/facHensel.h	(revision 139b671c6c11958a1f538c95dc0ade26520e76b0)
@@ -364,4 +364,43 @@
            );
 
+/// two factor Hensel lifting from univariate to bivariate, factors need not to
+/// be monic
+void 
+henselLift122 (const CanonicalForm& F,///< [in] a bivariate poly
+               CFList& factors,       ///< [in, out] a list of univariate polys
+                                      ///< lifted factors
+               int l,                 ///< [in] lift bound
+               CFArray& Pi,           ///< [in, out] stores intermediate results
+               CFList& diophant,      ///< [in, out] result of diophantine
+               CFMatrix& M,           ///< [in, out] stores intermediate results
+               const CFArray& LCs,    ///< [in] leading coefficients
+               bool sort              ///< [in] if true factors are sorted by
+                                      ///< their degree
+              );
+     
+
+/// two factor Hensel lifting from bivariate to multivariate, factors need not
+/// to be monic
+///
+/// @return @a henselLift122 returns a list of lifted factors   
+CFList
+henselLift2 (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
+             const int* l,         ///< [in] lift bounds
+             const int lLength,    ///< [in] length of l
+             bool sort,            ///< [in] if true factors are sorted by
+                                   ///< their degree in Variable(1)
+             const CFList& LCs1,   ///< [in] a list of evaluated LC of first 
+                                   ///< factor
+             const CFList& LCs2,   ///< [in] a list of evaluated LC of second 
+                                   ///< factor
+             const CFArray& Pi,    ///< [in] intermediate result
+             const CFList& diophant///< [in] result of diophantine
+            );
 #endif
 /* FAC_HENSEL_H */
