Index: factory/DegreePattern.cc
===================================================================
--- factory/DegreePattern.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/DegreePattern.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,152 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file DegreePattern.cc
+ * 
+ * This file provides functions for manipulating DegreePatterns      
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#include <config.h>
+
+#include "DegreePattern.h" 
+#include "cf_iter.h"
+#include "ftmpl_functions.h"
+#include "gfops.h"
+
+inline
+DegreePattern::DegreePattern (const CFList& l) 
+{ 
+  m_data = NULL;
+
+  if (l.length() == 0)
+    m_data = new Pattern();
+  else
+  {
+  
+  Variable x= Variable (1);
+  int p= getCharacteristic();
+  int d= getGFDegree();
+  char cGFName= gf_name;
+  setCharacteristic(0);
+  CanonicalForm buf= 1;
+  CFListIterator k= l;
+  for (int i= 0; i < l.length(); i++, k++)
+    buf *= (power (x, degree (k.getItem(), x)) + 1);
+
+  int j= 0;
+  for (CFIterator i= buf; i.hasTerms(); i++, j++)
+    ; 
+
+  ASSERT ( j > 1, "j > 1 expected" );
+  
+  m_data = new Pattern( j - 1 );
+
+  int i= 0;
+  for (CFIterator m = buf; i < getLength(); i++, m++)
+    (*this) [i]= m.exp();
+
+  if (d > 1)
+    setCharacteristic (p, d, cGFName);
+  else
+    setCharacteristic (p);
+  }
+}
+
+
+inline int
+DegreePattern::find (const int x) const
+{
+  if (getLength() == 0) return 0;
+  for (int i= 0; i < getLength(); i++) 
+    if ((*this)[i] == x) return i + 1;
+  return 0;
+}
+
+inline void
+DegreePattern::intersect (const DegreePattern& degPat)
+{
+  if (degPat.getLength() < getLength()) 
+  {
+    DegreePattern bufDeg= *this;
+    *this= degPat;
+    return (*this).intersect (bufDeg);
+  }
+
+  int count= 0;
+  int length= tmin (getLength(), degPat.getLength());
+  int* buf= new int [length];
+  for (int i= 0; i < length; i++) 
+  {
+    if (degPat.find ((*this)[i])) 
+    { 
+      buf[i]= (*this)[i];
+      count++;
+    }
+    else
+      buf[i]= -1;
+  }
+  ASSERT ( count > 0, "count > 0 expected" );
+
+  init (count);
+  count= 0;
+  for (int i= 0; i < length; i++) 
+  {
+    if (buf[i] != -1) 
+    {
+      (*this) [count]= buf[i];
+      count++;
+    }
+  }
+  delete[] buf;
+}
+   
+inline void
+DegreePattern::refine ()
+{
+  if (getLength() <= 1)
+    return;
+  int count= 0;
+  int* buf= new int [getLength()];
+  int d= (*this) [0];
+  int pos;
+  for (int i= 0; i < getLength(); i++)
+    buf[i]= -1;
+  for (int i= 1; i < getLength(); i++) 
+  {
+    pos= (*this).find (d - (*this)[i]);
+    if (pos) 
+    { 
+      buf[i]= (*this)[i];
+      count++;
+    }
+  }
+  buf[0]= d;
+  count++;
+  if (count == getLength())
+  {
+    return;
+  }
+  int length= getLength();
+
+  ASSERT ( count > 0, "count > 0 expected" );
+  init (count);
+  count= 0;
+  for (int i= 0; i < length; i++) 
+  {
+    if (buf[i] != -1) 
+    {
+      (*this)[count]= buf[i];
+      count++;
+    }
+  }
+
+  delete[] buf;
+  return;
+}
+
Index: factory/DegreePattern.h
===================================================================
--- factory/DegreePattern.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/DegreePattern.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,176 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file DegreePattern.h
+ * 
+ * This file provides a class to handle degree patterns.           
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+#ifndef DEGREE_PATTERN_H
+#define DEGREE_PATTERN_H
+
+#include <config.h>
+
+#include "assert.h"
+
+#include "canonicalform.h"
+#include "cf_iter.h"
+#include "ftmpl_functions.h"
+#include "gfops.h"
+
+/** @class DegreePattern DegreePattern.h "factory/DegreePattern.h"
+ * 
+ * DegreePattern provides a functionality to create, intersect and refine
+ * degree patterns.
+ *
+ * 
+ */
+class DegreePattern
+{
+private:
+  struct Pattern
+  {
+    int  m_refCounter; ///< reference counter
+    int  m_length;     ///< length of m_pattern
+    int* m_pattern;    ///< some array containing the degree pattern
+
+    /// construct a Pattern from an int
+    Pattern(int n): m_refCounter(1), m_length(n), m_pattern( new int[n]) {};
+    /// default contructor
+    Pattern(): m_refCounter(1), m_length(0), m_pattern(NULL) {};
+
+  }* m_data;
+
+  /// clear m_data
+  void release()
+  {
+    ASSERT ( m_data != NULL, "non-null pointer expected");
+    ASSERT ( m_data->m_refCounter == 0, "ref count of 0 expected");
+    if( m_data->m_pattern != NULL ) 
+      delete[] m_data->m_pattern;
+    m_data->m_pattern = NULL;
+    
+    delete m_data;
+    m_data = NULL;
+  }
+  /// initialise a DegreePattern
+  void init( int n )
+  {
+    ASSERT ( m_data != NULL, "non-null pointer expected" );
+    ASSERT( m_data->m_refCounter > 0, "ref count > 0 expected" );
+
+    if( (--m_data->m_refCounter) < 1 )
+      release();
+    
+    m_data = new Pattern(n);
+  }
+
+  /// getter
+  ///
+  /// @return @a getPattern returns a degree pattern
+  inline int* getPattern() const
+  {
+    ASSERT( m_data != NULL, "non-null pointer expected" );
+    ASSERT( m_data->m_pattern != NULL, "non-null pointer expected" );
+    return m_data->m_pattern;
+  }
+
+
+public:
+  /// getter
+  ///
+  /// @return @a getLength returns the length of the degree pattern
+  inline int getLength() const
+  {
+    ASSERT( m_data != NULL, "non-null pointer expected" );
+    return m_data->m_length;
+  }
+
+  /// operator []
+  ///
+  /// @return @a operator[] returns the element at @a index 
+  inline int operator[] (const int index ///< [in] some int >= 0, < getLength()
+                        ) const
+  {
+    ASSERT( m_data != NULL, "non-null pointer expected" );
+    ASSERT( index >= 0 && index < getLength(), "bad index" );
+    ASSERT( getPattern() != NULL, "non-null pointer expected" );
+    return getPattern()[index];
+  }
+
+  /// operator []
+  ///
+  /// @return @a operator[] sets the element at @a index  
+  inline int& operator[] (const int index ///< [in] some int >= 0, < getLength()
+                         )
+  {
+    ASSERT( m_data != NULL, "non-null pointer expected" );
+    ASSERT( index >= 0 && index < getLength(), "bad index" );
+    ASSERT( getPattern() != NULL, "non-null pointer expected" );
+    return getPattern()[index];
+  }
+
+  /// default constructor
+  DegreePattern(): m_data( new Pattern() ){}
+
+  /// copy constructor
+  DegreePattern (const DegreePattern& degPat ///< [in] some degree pattern
+                ): m_data( degPat.m_data ) 
+  {
+    ASSERT( degPat.m_data != NULL, "non-null pointer expected"  );
+    m_data->m_refCounter++;
+  };
+
+  /// construct a degree pattern from a list of (univariate) polys
+  DegreePattern (const CFList& l ///< [in] some list of (univariate) polys
+                );
+ 
+  /// assignment
+  DegreePattern& operator= (const DegreePattern& degPat ///< [in] some degree
+                                                        ///< pattern
+                           )
+  {
+    ASSERT( m_data != NULL, "non-null pointer expected"  );
+    ASSERT( degPat.m_data != NULL, "non-null pointer expected" );
+    if( m_data != degPat.m_data )
+    {
+      m_data = degPat.m_data;
+      m_data->m_refCounter++;
+    }
+
+    return *this;
+  }
+
+  /// destructor
+  ~DegreePattern () 
+  {
+    ASSERT( m_data !=  NULL, "non-null pointer expected"  );
+    if( (--m_data->m_refCounter) < 1 ) 
+      release();
+  }
+
+  /// find an element @a x 
+  ///
+  /// @return @a find returns the index + 1 of @a x, if @a x is an element of 
+  ///         the degree pattern, 0 otherwise
+  inline int find (const int x ///< [in] some int
+           ) const;
+  
+  /// intersect two degree patterns 
+  inline void intersect (const DegreePattern& degPat ///< [in] some degree pattern
+                 );
+  /// Refine a degree pattern. Assumes that (*this)[0]:= @a d is the degree
+  /// of the poly to be factored. Now for every other entry @a a there should be
+  /// some entry @a b such that @a a+b= d. Elements which do not satisfy this 
+  /// relation are removed.
+  inline void refine ();
+};
+
+#endif
+/* DEGREE_PATTERN_H */
+
Index: factory/ExtensionInfo.cc
===================================================================
--- factory/ExtensionInfo.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/ExtensionInfo.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,96 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file ExtensionInfo.cc
+ * 
+ * This file provides member functions for ExtensionInfo      
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#include <config.h>
+
+#include "ExtensionInfo.h"
+
+inline
+ExtensionInfo::ExtensionInfo (const bool extension)
+{
+  m_alpha= Variable (1);
+  m_beta= Variable (1);
+  m_gamma= CanonicalForm ();
+  m_delta= CanonicalForm ();
+  m_GFDegree= 1;
+  m_GFName= 'Z';
+  m_extension= extension;
+}
+
+inline
+ExtensionInfo::ExtensionInfo (const Variable& alpha, const Variable& beta, 
+                              const CanonicalForm& gamma, const CanonicalForm& 
+                              delta, const int nGFDegree, const char cGFName,
+                              const bool extension)
+{
+  m_alpha= alpha;
+  m_beta= beta;
+  m_gamma= gamma;
+  m_delta= delta;
+  m_GFDegree= nGFDegree;
+  m_GFName= cGFName;
+  m_extension= extension;
+}
+
+inline
+ExtensionInfo::ExtensionInfo (const Variable& alpha, const Variable& beta, 
+                              const CanonicalForm& gamma, const CanonicalForm& 
+                              delta)
+{
+  m_alpha= alpha;
+  m_beta= beta;
+  m_gamma= gamma;
+  m_delta= delta;
+  m_GFDegree= 0;
+  m_GFName= 'Z';
+  m_extension= true;
+}
+
+inline
+ExtensionInfo::ExtensionInfo (const Variable& alpha, const bool extension)
+{
+  m_alpha= alpha;
+  m_beta= Variable (1);
+  m_gamma= CanonicalForm ();
+  m_delta= CanonicalForm ();
+  m_GFDegree= 0;
+  m_GFName= 'Z';
+  m_extension= extension;
+}
+
+inline
+ExtensionInfo::ExtensionInfo (const Variable& alpha)
+{
+  m_alpha= alpha;
+  m_beta=  Variable (1);
+  m_gamma= CanonicalForm ();
+  m_delta= CanonicalForm ();
+  m_GFDegree= 1;
+  m_GFName= 'Z';
+  m_extension= true;
+}
+
+inline
+ExtensionInfo::ExtensionInfo (const int nGFDegree, const char cGFName, const
+                              bool extension)
+{
+  m_alpha= Variable (1);
+  m_beta= Variable (1);
+  m_gamma= CanonicalForm ();
+  m_delta= CanonicalForm ();
+  m_GFDegree= nGFDegree;
+  m_GFName= cGFName;
+  m_extension= extension;
+}
+
Index: factory/ExtensionInfo.h
===================================================================
--- factory/ExtensionInfo.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/ExtensionInfo.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,163 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file ExtensionInfo.h
+ * 
+ * This file provides a class to store information about finite fields and
+ * extensions thereof.   
+ *         
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#ifndef EXTENSION_INFO_H
+#define EXTENSION_INFO_H
+
+#include <config.h>
+
+#include "canonicalform.h"
+
+/** @class ExtensionInfo ExtensionInfo.h "factory/ExtensionInfo.h"
+ *  ExtensionInfo contains information about extension. 
+ *  If @a m_extension is true we are in an extension of some initial field. 
+ *  If the initial field is \f$ F_p \f$ and we pass to \f$ F_p (\alpha) \f$ 
+ *  then @a m_alpha is an algebraic variable, @a m_beta= Variable(1), 
+ *  @a m_gamma= @a m_delta= 1, @a m_GFDegree= 0, @a m_GFName= 'Z'. If we pass
+ *  to some GF (p^k) then @a m_alpha= Variable (1), @a m_beta= Variable(1), 
+ *  @a m_gamma= @a m_delta= 1, @a m_GFDegree= 1, @a m_GFName= 'Z'.  
+ *  @n If the initial field is \f$ F_p (\epsilon) \f$, then @a m_beta= 
+ *  \f$ \epsilon \f$, @a m_alpha an algebraic variable defining an extension of 
+ *  \f$ F_p (\epsilon) \f$, @a m_gamma is a primitive element of 
+ *  \f$ F_p (\alpha) \f$, @a m_delta is a primitive element of 
+ *  \f$ F_p (\beta) \f$, @a m_GFDegree= 0, @a m_GFName= 'Z'. 
+ *  @n If the initial field is GF(p^k), then @a m_alpha= Variable (1), 
+ *  @a m_beta= Variable (1), @a m_gamma= 1, @a m_delta= 1, @a m_GFDegree()= k, 
+ *  @a m_GFName= gf_name of the initial field. 
+ *  @n If @a m_extension is false and the current field is \f$ F_p \f$ then  
+ *  @a m_alpha= Variable (1), @a m_beta= Variable (1), @a m_gamma= 1, 
+ *  @a m_delta= 1, @a m_GFDegree= 1, @a m_GFName= 'Z'. 
+ *  @n If the current field is \f$ F_p (\alpha) \f$ then 
+ *  @a m_alpha is some algebraic variable, @a m_beta= Variable (1), 
+ *  @a m_gamma= 1, @a m_delta= 1, @a m_GFDegree= 0, @a m_GFName= 'Z'. 
+ *  @n If the current field is GF then @a m_alpha= Variable (1), 
+ *  @a m_beta= Variable (1),  @a m_gamma= 1, @a m_delta= 1, 
+ *  @a m_GFDegree= getGFDegree(), @a m_GFName= gf_name.
+ *
+ *  @sa facFqBivar.h, facFqFactorize.h 
+ */
+class ExtensionInfo
+{
+private:
+  /// an algebraic variable or Variable (1) 
+  Variable m_alpha;
+  /// an algebraic variable or Variable (1)
+  Variable m_beta;  
+  /// a primitive element of \f$ F_p (\alpha) \f$ or 1
+  CanonicalForm m_gamma;
+  /// a primitive element of \f$ F_p (\beta) \f$ or 1
+  CanonicalForm m_delta;
+  /// GF degree or 1
+  int m_GFDegree;
+  /// name of GF variable
+  char m_GFName;
+  /// indicates if we are in an extension of some initial field
+  bool m_extension;
+public:
+  /// \f$ F_p \f$ as initial field, if @a extension is true we are in some GF
+  ExtensionInfo (const bool extension  ///< [in] some bool 
+                );
+  /// Construct an @a ExtensionInfo 
+  ExtensionInfo (const Variable& alpha,      ///< [in] some algebraic variable 
+                 const Variable& beta,       ///< [in] some algebraic variable
+                 const CanonicalForm& gamma, ///< [in] some primitive element
+                                             ///< of \f$ F_p (\alpha) \f$
+                 const CanonicalForm& delta, ///< [in] some primitive element
+                                             ///< of \f$ F_p (\beta) \f$
+                 const int nGFDegree,        ///< [in] GFDegree of initial field
+                 const char cGFName,         ///< [in] name of GF variable of
+                                             ///< initial field
+                 const bool extension        ///< [in] some bool
+                );
+  /// \f$ F_p (\beta) \f$ as initial field and switch to an extension given by
+  /// @a alpha, needs primitive elements @a gamma and @a delta for maps
+  /// between \f$ F_p (\alpha) \subset F_p (\beta) \f$ 
+  ExtensionInfo (const Variable& alpha,      ///< [in] some algebraic variable
+                 const Variable& beta,       ///< [in] some algebraic variable
+                 const CanonicalForm& gamma, ///< [in] some primitive element
+                                             ///< of \f$ F_p (\alpha) \f$
+                 const CanonicalForm& delta  ///< [in] some primitive element
+                                             ///< of \f$ F_p (\beta) \f$
+                );
+  /// \f$ F_p (\alpha) \f$ as initial field, if @a extension is false.
+  /// Else initial field is \f$ F_p \f$ 
+  ExtensionInfo (const Variable& alpha, ///< [in] some algebraic variable
+                 const bool extension   ///< [in] some bool
+                ); 
+
+  ExtensionInfo (const Variable& alpha ///< [in] some algebraic variable
+                );
+ 
+  /// GF as initial field
+  ExtensionInfo (const int nGFDegree,   ///< [in] GF degree of initial field 
+                 const char cGFName,    ///< [in] name of GF variable 
+                 const bool extension   ///< [in] some bool
+                ); 
+
+  /// getter
+  ///
+  /// @return @a getAlpha() returns @a m_alpha
+  Variable getAlpha () const
+  {
+    return m_alpha;
+  }
+  /// getter
+  ///
+  /// @return @a getBeta() returns @a m_beta
+  Variable getBeta () const
+  {
+    return m_beta;
+  }
+  /// getter
+  ///
+  /// @return @a getGamma() returns @a m_gamma
+  CanonicalForm getGamma() const
+  {
+    return m_gamma;
+  }
+  /// getter
+  ///
+  /// @return @a getDelta() returns @a m_delta
+  CanonicalForm getDelta() const
+  {
+    return m_delta;
+  }
+  /// getter
+  ///
+  /// @return @a getGFDegree() returns @a m_GFDegree
+  int getGFDegree() const
+  {
+    return m_GFDegree;
+  }
+  /// getter
+  ///
+  /// @return @a getGFName() returns @a m_GFName
+  char getGFName() const
+  {
+    return m_GFName;
+  }
+  /// getter
+  ///
+  /// @return @a isInextension() returns @a m_extension
+  bool isInExtension() const
+  {
+    return m_extension;
+  }
+};
+
+#endif
+/* EXTENSION_INFO_H */
+
Index: factory/cf_factor.cc
===================================================================
--- factory/cf_factor.cc	(revision 9015fc56fb11f10a66eb0275771d88644708bac3)
+++ factory/cf_factor.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -27,5 +27,7 @@
 #include "fac_sqrfree.h"
 #include "cf_algorithm.h"
+#include "facFqFactorize.h"
 #include "cf_map.h"
+#include "algext.h"
 
 #include "int_int.h"
@@ -427,110 +429,142 @@
   if ( getCharacteristic() > 0 )
   {
-    ASSERT( f.isUnivariate(), "multivariate factorization not implemented" );
-    #ifdef HAVE_NTL
-    if (isOn(SW_USE_NTL) && (isPurePoly(f)))
-    {
-      // USE NTL
-      if (getCharacteristic()!=2)
-      {
-        if (fac_NTL_char!=getCharacteristic())
-        {
-          fac_NTL_char=getCharacteristic();
-          #ifndef NTL_ZZ
-          if (fac_NTL_char >NTL_SP_BOUND)
-          {
-            ZZ r;
-            r=getCharacteristic();
-            ZZ_pContext ccc(r);
-            ccc.restore();
-            ZZ_p::init(r);
-          }
-          else
-          #endif
-          {
-            #ifdef NTL_ZZ
-            ZZ r;
-            r=getCharacteristic();
-            ZZ_pContext ccc(r);
-            #else
-            zz_pContext ccc(getCharacteristic());
-            #endif
-            ccc.restore();
-            #ifdef NTL_ZZ
-            ZZ_p::init(r);
-            #else
-            zz_p::init(getCharacteristic());
-            #endif
-          }
-        }
-       #ifndef NTL_ZZ
-       if (fac_NTL_char >NTL_SP_BOUND)
-       {
-          // convert to NTL
-          ZZ_pX f1=convertFacCF2NTLZZpX(f);
-          ZZ_p leadcoeff = LeadCoeff(f1);
-          //make monic
-          f1=f1 / LeadCoeff(f1);
-          // factorize
-          vec_pair_ZZ_pX_long factors;
-          CanZass(factors,f1);
-          // convert back to factory
-          F=convertNTLvec_pair_ZZpX_long2FacCFFList(factors,leadcoeff,f.mvar());
-        }
+    if (f.isUnivariate())
+    {
+      #ifdef HAVE_NTL
+      if (isOn(SW_USE_NTL) && (isPurePoly(f)))
+      {
+	// USE NTL
+	if (getCharacteristic()!=2)
+	{
+	  if (fac_NTL_char!=getCharacteristic())
+	  {
+	    fac_NTL_char=getCharacteristic();
+	    #ifndef NTL_ZZ
+	    if (fac_NTL_char >NTL_SP_BOUND)
+	    {
+	      ZZ r;
+	      r=getCharacteristic();
+	      ZZ_pContext ccc(r);
+	      ccc.restore();
+	      ZZ_p::init(r);
+	    }
+	    else
+	    #endif
+	    {
+	      #ifdef NTL_ZZ
+	      ZZ r;
+	      r=getCharacteristic();
+	      ZZ_pContext ccc(r);
+	      #else
+	      zz_pContext ccc(getCharacteristic());
+	      #endif
+	      ccc.restore();
+	      #ifdef NTL_ZZ
+	      ZZ_p::init(r);
+	      #else
+	      zz_p::init(getCharacteristic());
+	      #endif
+	    }
+	  }
+	#ifndef NTL_ZZ
+	if (fac_NTL_char >NTL_SP_BOUND)
+	{
+	    // convert to NTL
+	    ZZ_pX f1=convertFacCF2NTLZZpX(f);
+	    ZZ_p leadcoeff = LeadCoeff(f1);
+	    //make monic
+	    f1=f1 / LeadCoeff(f1);
+	    // factorize
+	    vec_pair_ZZ_pX_long factors;
+	    CanZass(factors,f1);
+	    // convert back to factory
+	    F=convertNTLvec_pair_ZZpX_long2FacCFFList(factors,leadcoeff,f.mvar());
+	  }
+	  else
+	  #endif
+	  {
+	    // convert to NTL
+	    #ifdef NTL_ZZ
+	    ZZ_pX f1=convertFacCF2NTLZZpX(f);
+	    ZZ_p leadcoeff = LeadCoeff(f1);
+	    #else
+	    zz_pX f1=convertFacCF2NTLzzpX(f);
+	    zz_p leadcoeff = LeadCoeff(f1);
+	    #endif
+	    //make monic
+	    f1=f1 / LeadCoeff(f1);
+	    // factorize
+	    #ifdef NTL_ZZ
+	    vec_pair_ZZ_pX_long factors;
+	    #else
+	    vec_pair_zz_pX_long factors;
+	    #endif
+	    CanZass(factors,f1);
+	    // convert back to factory
+	    #ifdef NTL_ZZ
+	    F=convertNTLvec_pair_ZZpX_long2FacCFFList(factors,leadcoeff,f.mvar());
+	    #else
+	    F=convertNTLvec_pair_zzpX_long2FacCFFList(factors,leadcoeff,f.mvar());
+	    #endif
+	  }
+	  //test_cff(F,f);
+	}
+	else
+	{
+	  // Specialcase characteristic==2
+	  if (fac_NTL_char!=2)
+	  {
+	    fac_NTL_char=2;
+	    zz_p::init(2);
+	  }
+	  // convert to NTL using the faster conversion routine for characteristic 2
+	  GF2X f1=convertFacCF2NTLGF2X(f);
+	  // no make monic necessary in GF2
+	  //factorize
+	  vec_pair_GF2X_long factors;
+	  CanZass(factors,f1);
+
+	  // convert back to factory again using the faster conversion routine for vectors over GF2X
+	  F=convertNTLvec_pair_GF2X_long2FacCFFList(factors,LeadCoeff(f1),f.mvar());
+	}
+      }
+      else
+      #endif
+      {  // Use Factory without NTL
+	if ( isOn( SW_BERLEKAMP ) )
+	  F=FpFactorizeUnivariateB( f, issqrfree );
+	else
+	  F=FpFactorizeUnivariateCZ( f, issqrfree, 0, Variable(), Variable() );
+      }
+    }
+    else
+    {
+      #ifdef HAVE_NTL
+      if (issqrfree)
+      {
+        CFList factors;
+        Variable alpha;
+        if (CFFactory::gettype() == GaloisFieldDomain)
+          factors= GFSqrfFactorize (f);
+        else if (hasFirstAlgVar (f, alpha))
+          factors= FqSqrfFactorize (f, alpha);
         else
-        #endif
-        {
-          // convert to NTL
-          #ifdef NTL_ZZ
-          ZZ_pX f1=convertFacCF2NTLZZpX(f);
-          ZZ_p leadcoeff = LeadCoeff(f1);
-          #else
-          zz_pX f1=convertFacCF2NTLzzpX(f);
-          zz_p leadcoeff = LeadCoeff(f1);
-          #endif
-          //make monic
-          f1=f1 / LeadCoeff(f1);
-          // factorize
-          #ifdef NTL_ZZ
-          vec_pair_ZZ_pX_long factors;
-          #else
-          vec_pair_zz_pX_long factors;
-          #endif
-          CanZass(factors,f1);
-          // convert back to factory
-          #ifdef NTL_ZZ
-          F=convertNTLvec_pair_ZZpX_long2FacCFFList(factors,leadcoeff,f.mvar());
-          #else
-          F=convertNTLvec_pair_zzpX_long2FacCFFList(factors,leadcoeff,f.mvar());
-          #endif
-        }
-        //test_cff(F,f);
+          factors= FpSqrfFactorize (f);
+        for (CFListIterator i= factors; i.hasItem(); i++)
+          F.append (CFFactor (i.getItem(), 1));
       }
       else
       {
-        // Specialcase characteristic==2
-        if (fac_NTL_char!=2)
-        {
-          fac_NTL_char=2;
-          zz_p::init(2);
-        }
-        // convert to NTL using the faster conversion routine for characteristic 2
-        GF2X f1=convertFacCF2NTLGF2X(f);
-        // no make monic necessary in GF2
-        //factorize
-        vec_pair_GF2X_long factors;
-        CanZass(factors,f1);
-
-        // convert back to factory again using the faster conversion routine for vectors over GF2X
-        F=convertNTLvec_pair_GF2X_long2FacCFFList(factors,LeadCoeff(f1),f.mvar());
-      }
-    }
-    else
-    #endif
-    {  // Use Factory without NTL
-      if ( isOn( SW_BERLEKAMP ) )
-         F=FpFactorizeUnivariateB( f, issqrfree );
-      else
-        F=FpFactorizeUnivariateCZ( f, issqrfree, 0, Variable(), Variable() );
+        Variable alpha;
+        if (CFFactory::gettype() == GaloisFieldDomain)
+          F= GFFactorize (f);
+        else if (hasFirstAlgVar (f, alpha))
+          F= FqFactorize (f, alpha);
+        else
+          F= FpFactorize (f); 
+      }
+      #else
+      ASSERT( f.isUnivariate(), "multivariate factorization not implemented" );
+      #endif
     }
   }
@@ -641,111 +675,122 @@
   CFFList F;
   ASSERT( alpha.level() < 0, "not an algebraic extension" );
-  ASSERT( f.isUnivariate(), "multivariate factorization not implemented" );
   ASSERT( getCharacteristic() > 0, "char 0 factorization not implemented" );
-  #ifdef HAVE_NTL
-  if  (isOn(SW_USE_NTL))
-  {
-    //USE NTL
-    if (getCharacteristic()!=2)
-    {
-      // First all cases with characteristic !=2
-      // set remainder
-      if (fac_NTL_char!=getCharacteristic())
-      {
-        fac_NTL_char=getCharacteristic();
-        #ifdef NTL_ZZ
-        ZZ r;
-        r=getCharacteristic();
-        ZZ_pContext ccc(r);
-        #else
-        zz_pContext ccc(getCharacteristic());
-        #endif
-        ccc.restore();
-        #ifdef NTL_ZZ
-        ZZ_p::init(r);
-        #else
-        zz_p::init(getCharacteristic());
-        #endif
-      }
-
-      // set minimal polynomial in NTL
-      #ifdef NTL_ZZ
-      ZZ_pX minPo=convertFacCF2NTLZZpX(getMipo(alpha));
-      ZZ_pEContext c(minPo);
-      #else
-      zz_pX minPo=convertFacCF2NTLzzpX(getMipo(alpha));
-      zz_pEContext c(minPo);
-      #endif
-
-      c.restore();
-
-      // convert to NTL
-      #ifdef NTL_ZZ
-      ZZ_pEX f1=convertFacCF2NTLZZ_pEX(f,minPo);
-      ZZ_pE leadcoeff= LeadCoeff(f1);
-      #else
-      zz_pEX f1=convertFacCF2NTLzz_pEX(f,minPo);
-      zz_pE leadcoeff= LeadCoeff(f1);
-      #endif
-
-      //make monic
-      f1=f1 / leadcoeff;
-
-      // factorize using NTL
-      #ifdef NTL_ZZ
-      vec_pair_ZZ_pEX_long factors;
-      #else
-      vec_pair_zz_pEX_long factors;
-      #endif
-      CanZass(factors,f1);
-
-      // return converted result
-      F=convertNTLvec_pair_zzpEX_long2FacCFFList(factors,leadcoeff,f.mvar(),alpha);
-    }
-    else
-    {
-      // special case : GF2
-
-      // remainder is two ==> nothing to do
-      // set remainder
-      ZZ r;
-      r=getCharacteristic();
-      ZZ_pContext ccc(r);
-      ccc.restore();
-
-      // set minimal polynomial in NTL using the optimized conversion routines for characteristic 2
-      GF2X minPo=convertFacCF2NTLGF2X(getMipo(alpha,f.mvar()));
-      GF2EContext c(minPo);
-      c.restore();
-
-      // convert to NTL again using the faster conversion routines
-      GF2EX f1;
-      if (isPurePoly(f))
-      {
-        GF2X f_tmp=convertFacCF2NTLGF2X(f);
-        f1=to_GF2EX(f_tmp);
+  if (f.isUnivariate())
+  {
+    #ifdef HAVE_NTL
+    if  (isOn(SW_USE_NTL))
+    {
+      //USE NTL
+      if (getCharacteristic()!=2)
+      {
+	// First all cases with characteristic !=2
+	// set remainder
+	if (fac_NTL_char!=getCharacteristic())
+	{
+	  fac_NTL_char=getCharacteristic();
+	  #ifdef NTL_ZZ
+	  ZZ r;
+	  r=getCharacteristic();
+	  ZZ_pContext ccc(r);
+	  #else
+	  zz_pContext ccc(getCharacteristic());
+	  #endif
+	  ccc.restore();
+	  #ifdef NTL_ZZ
+	  ZZ_p::init(r);
+	  #else
+	  zz_p::init(getCharacteristic());
+	  #endif
+	}
+
+	// set minimal polynomial in NTL
+	#ifdef NTL_ZZ
+	ZZ_pX minPo=convertFacCF2NTLZZpX(getMipo(alpha));
+	ZZ_pEContext c(minPo);
+	#else
+	zz_pX minPo=convertFacCF2NTLzzpX(getMipo(alpha));
+	zz_pEContext c(minPo);
+	#endif
+
+	c.restore();
+
+	// convert to NTL
+	#ifdef NTL_ZZ
+	ZZ_pEX f1=convertFacCF2NTLZZ_pEX(f,minPo);
+	ZZ_pE leadcoeff= LeadCoeff(f1);
+	#else
+	zz_pEX f1=convertFacCF2NTLzz_pEX(f,minPo);
+	zz_pE leadcoeff= LeadCoeff(f1);
+	#endif
+
+	//make monic
+	f1=f1 / leadcoeff;
+
+	// factorize using NTL
+	#ifdef NTL_ZZ
+	vec_pair_ZZ_pEX_long factors;
+	#else
+	vec_pair_zz_pEX_long factors;
+	#endif
+	CanZass(factors,f1);
+
+	// return converted result
+	F=convertNTLvec_pair_zzpEX_long2FacCFFList(factors,leadcoeff,f.mvar(),alpha);
       }
       else
       {
-        f1=convertFacCF2NTLGF2EX(f,minPo);
-      }
-
-      // make monic (in Z/2(a))
-      GF2E f1_coef=LeadCoeff(f1);
-      MakeMonic(f1);
-
-      // factorize using NTL
-      vec_pair_GF2EX_long factors;
-      CanZass(factors,f1);
-
-      // return converted result
-      F=convertNTLvec_pair_GF2EX_long2FacCFFList(factors,f1_coef,f.mvar(),alpha);
-    }
-
+	// special case : GF2
+
+	// remainder is two ==> nothing to do
+	// set remainder
+	ZZ r;
+	r=getCharacteristic();
+	ZZ_pContext ccc(r);
+	ccc.restore();
+
+	// set minimal polynomial in NTL using the optimized conversion routines for characteristic 2
+	GF2X minPo=convertFacCF2NTLGF2X(getMipo(alpha,f.mvar()));
+	GF2EContext c(minPo);
+	c.restore();
+
+	// convert to NTL again using the faster conversion routines
+	GF2EX f1;
+	if (isPurePoly(f))
+	{
+	  GF2X f_tmp=convertFacCF2NTLGF2X(f);
+	  f1=to_GF2EX(f_tmp);
+	}
+	else
+	{
+	  f1=convertFacCF2NTLGF2EX(f,minPo);
+	}
+
+	// make monic (in Z/2(a))
+	GF2E f1_coef=LeadCoeff(f1);
+	MakeMonic(f1);
+
+	// factorize using NTL
+	vec_pair_GF2EX_long factors;
+	CanZass(factors,f1);
+
+	// return converted result
+	F=convertNTLvec_pair_GF2EX_long2FacCFFList(factors,f1_coef,f.mvar(),alpha);
+      }
+
+    }
+    else
+    #endif
+    {
+      F=FpFactorizeUnivariateCZ( f, false, 1, alpha, Variable() );
+    }
   }
   else
-  #endif
-  {
-    F=FpFactorizeUnivariateCZ( f, false, 1, alpha, Variable() );
+  {
+    #ifdef HAVE_NTL
+    F= FqFactorize (f, alpha);
+    #else
+    ASSERT( f.isUnivariate(), "multivariate factorization not implemented" );
+    #endif
+    
   }
   return F;
Index: factory/facFqBivar.cc
===================================================================
--- factory/facFqBivar.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/facFqBivar.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,1301 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file facFqBivar.cc
+ * 
+ * This file provides functions for factorizing a bivariate polynomial over
+ * \f$ F_{p} \f$ , \f$ F_{p}(\alpha ) \f$ or GF, based on "Modern Computer 
+ * Algebra, Chapter 15" by J. von zur Gathen & J. Gerhard and "Factoring 
+ * multivariate polynomials over a finite field" by L. Bernardin.    
+ *
+ * ABSTRACT: In contrast to biFactorizer() in facFqFactorice.cc we evaluate and 
+ * factorize the polynomial in both variables. So far factor recombination is 
+ * done naive!         
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#include <config.h>
+
+#include "assert.h"
+#include "debug.h"
+#include "timing.h"
+
+#include "canonicalform.h"
+#include "cf_defs.h"
+#include "cf_map_ext.h"
+#include "cf_random.h"
+#include "facHensel.h"
+#include "cf_map.h"
+#include "cf_gcd_smallp.h"
+#include "DegreePattern.cc"
+#include "facFqBivarUtil.cc"
+#include "facFqBivar.h"
+
+
+#ifdef HAVE_NTL
+
+TIMING_DEFINE_PRINT(fac_uni_factorizer);
+TIMING_DEFINE_PRINT(fac_hensel_lift);
+TIMING_DEFINE_PRINT(fac_factor_recombination);
+
+CanonicalForm prodMod0 (const CFList& L, const CanonicalForm& M) 
+{
+  if (L.isEmpty())
+    return 1;
+  else if (L.length() == 1)
+    return mod (L.getFirst()(0, 1) , M);
+  else if (L.length() == 2)
+    return mod (L.getFirst()(0, 1)*L.getLast()(0, 1), 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= prodMod0 (tmp1, M);
+    buf2= prodMod0 (tmp2, M);
+    return mod (buf1*buf2, M);
+  }
+}
+
+CanonicalForm evalPoint (const CanonicalForm& F, CanonicalForm & eval, 
+                         const Variable& alpha, CFList& list, const bool& GF,
+                         bool& fail) 
+{
+  fail= false;
+  Variable x= Variable(2);
+  FFRandom genFF;
+  GFRandom genGF;
+  CanonicalForm random, mipo;
+  double bound;
+  int p= getCharacteristic ();
+  if (alpha != Variable(1)) 
+  {
+    mipo= getMipo (alpha);
+    int d= degree (mipo);
+    bound= ipower (p, d);
+  }
+  else if (GF) 
+  {
+    int d= getGFDegree();
+    bound= ipower (p, d);
+  }
+  else
+    bound= p;
+
+  do 
+  {
+    if (list.length() >= bound)
+    {
+      fail= true;
+      break;
+    }
+    if (GF)
+      random= genGF.generate();
+    else if (list.length() < p || alpha == x) 
+      random= genFF.generate();
+    else if (alpha != x && list.length() >= p) 
+    {
+      AlgExtRandomF genAlgExt (alpha);
+      random= genAlgExt.generate();
+    }
+    if (find (list, random)) continue;
+    eval= F (random, x);
+    if (degree (eval) != degree (F, Variable (1))) 
+    { //leading coeff vanishes
+      if (!find (list, random))
+        list.append (random);
+      continue;
+    }
+    if (degree (gcd (deriv (eval, eval.mvar()), eval), eval.mvar()) > 0) 
+    { //evaluated polynomial is not squarefree
+      if (!find (list, random))
+        list.append (random);
+      continue;
+    }
+  } while (find (list, random));
+
+  return random;
+}
+
+CFList 
+uniFactorizer (const CanonicalForm& A, const Variable& alpha, const bool& GF) 
+{
+  Variable x= A.mvar();
+  ASSERT (A.isUnivariate(), "univariate polynomial expected");
+  CFFList factorsA;
+  ZZ p= to_ZZ (getCharacteristic());
+  ZZ_p::init (p); 
+  if (GF) 
+  {
+    Variable beta= rootOf (gf_mipo); 
+    int k= getGFDegree();
+    char cGFName= gf_name;
+    setCharacteristic (getCharacteristic());     
+    CanonicalForm buf= GF2FalphaRep (A, beta); 
+    if (getCharacteristic() > 2)
+    {
+      ZZ_pX NTLMipo= convertFacCF2NTLZZpX (gf_mipo);
+      ZZ_pE::init (NTLMipo);
+      ZZ_pEX NTLA= convertFacCF2NTLZZ_pEX (buf, NTLMipo);
+      MakeMonic (NTLA);
+      vec_pair_ZZ_pEX_long NTLFactorsA= CanZass (NTLA);
+      ZZ_pE multi= to_ZZ_pE (1);
+      factorsA= convertNTLvec_pair_ZZpEX_long2FacCFFList (NTLFactorsA, multi,
+                                                         x, beta);
+    }
+    else
+    {
+      GF2X NTLMipo= convertFacCF2NTLGF2X (gf_mipo);
+      GF2E::init (NTLMipo);
+      GF2EX NTLA= convertFacCF2NTLGF2EX (buf, NTLMipo);
+      MakeMonic (NTLA);
+      vec_pair_GF2EX_long NTLFactorsA= CanZass (NTLA);
+      GF2E multi= to_GF2E (1);
+      factorsA= convertNTLvec_pair_GF2EX_long2FacCFFList (NTLFactorsA, multi,
+                                                           x, beta);
+    }
+    setCharacteristic (getCharacteristic(), k, cGFName);
+    for (CFFListIterator i= factorsA; i.hasItem(); i++) 
+    {
+      buf= i.getItem().factor();
+      buf= Falpha2GFRep (buf);
+      i.getItem()= CFFactor (buf, i.getItem().exp());
+    }
+  } 
+  else if (alpha != Variable(1)) 
+  {
+    if (getCharacteristic() > 2)
+    {
+      ZZ_pX NTLMipo= convertFacCF2NTLZZpX (getMipo (alpha));
+      ZZ_pE::init (NTLMipo);  
+      ZZ_pEX NTLA= convertFacCF2NTLZZ_pEX (A, NTLMipo);
+      MakeMonic (NTLA);
+      vec_pair_ZZ_pEX_long NTLFactorsA= CanZass (NTLA);
+      ZZ_pE multi= to_ZZ_pE (1);
+      factorsA= convertNTLvec_pair_ZZpEX_long2FacCFFList (NTLFactorsA, multi, 
+                                                           x, alpha);
+    }
+    else
+    {
+      GF2X NTLMipo= convertFacCF2NTLGF2X (getMipo (alpha));
+      GF2E::init (NTLMipo);  
+      GF2EX NTLA= convertFacCF2NTLGF2EX (A, NTLMipo);
+      MakeMonic (NTLA);
+      vec_pair_GF2EX_long NTLFactorsA= CanZass (NTLA);
+      GF2E multi= to_GF2E (1);
+      factorsA= convertNTLvec_pair_GF2EX_long2FacCFFList (NTLFactorsA, multi, 
+                                                           x, alpha);
+    } 
+  }
+  else 
+  {
+    if (getCharacteristic() > 2)
+    {
+      ZZ_pX NTLA= convertFacCF2NTLZZpX (A);
+      MakeMonic (NTLA);
+      vec_pair_ZZ_pX_long NTLFactorsA= CanZass (NTLA);
+      ZZ_p multi= to_ZZ_p (1);
+      factorsA= convertNTLvec_pair_ZZpX_long2FacCFFList (NTLFactorsA, multi, 
+                                                          x);
+    }
+    else
+    {
+      GF2X NTLA= convertFacCF2NTLGF2X (A);
+      vec_pair_GF2X_long NTLFactorsA= CanZass (NTLA);
+      GF2 multi= to_GF2 (1);
+      factorsA= convertNTLvec_pair_GF2X_long2FacCFFList (NTLFactorsA, multi, 
+                                                          x);
+    }
+  }
+  CFList uniFactors;
+  for (CFFListIterator i= factorsA; i.hasItem(); i++) 
+    uniFactors.append (i.getItem().factor());
+  return uniFactors;
+} 
+
+/// naive factor recombination as decribed in "Factoring 
+/// multivariate polynomials over a finite field" by L Bernardin. 
+CFList 
+extFactorRecombination (const CFList& factors, const CanonicalForm& F,  
+                        const CanonicalForm& M, const ExtensionInfo& info, 
+                        const DegreePattern& degs, const CanonicalForm& eval) 
+{ 
+  if (factors.length() == 0)
+    return CFList();
+
+  Variable alpha= info.getAlpha();
+  Variable beta= info.getBeta();
+  CanonicalForm gamma= info.getGamma();
+  CanonicalForm delta= info.getDelta();
+  int k= info.getGFDegree();
+
+  Variable y= F.mvar();
+  CFList source, dest;
+  if (degs.getLength() <= 1 || factors.length() == 1) 
+    return CFList(mapDown (F(y-eval, y), info, source, dest));  
+  
+  DEBOUTLN (cerr, "LC (F, 1)*prodMod (factors, M) == F " << 
+            (LC (F, 1)*prodMod (factors, M) == F));
+  int degMipoBeta;
+  if (!k && beta == Variable(1))
+    degMipoBeta= 1;
+  else if (!k && beta != Variable(1))
+    degMipoBeta= degree (getMipo (beta));
+
+  CFList T, S, Diff;
+  T= factors;
+
+  int s= 1;
+  CFList result;
+  CanonicalForm buf, buf2;
+  if (beta != Variable (1))
+    buf= mapDown (F, gamma, delta, alpha, source, dest); 
+  else
+    buf= F;
+
+  CanonicalForm g, LCBuf= LC (buf, Variable (1));
+  int v [T.length()];
+  for (int i= 0; i < T.length(); i++)
+    v[i]= 0;
+
+  CFArray TT;
+  DegreePattern bufDegs1, bufDegs2; 
+  bufDegs1= degs;
+  int subsetDeg;
+  TT= copy (factors); 
+  bool nosubset= false;
+  bool recombination= false;
+  while (T.length() >= 2*s) 
+  {
+    while (nosubset == false) 
+    {
+      if (T.length() == s) 
+      {
+        if (recombination) 
+        {
+          T.insert (LCBuf);
+          g= prodMod (T, M);
+          T.removeFirst();
+          g /= content(g);
+          g= g (y - eval, y);
+          g /= Lc (g);
+          appendTestMapDown (result, g, info, source, dest);
+          return result;
+        }
+        else
+        {
+          appendMapDown (result, F (y - eval, y), info, source, dest);
+          return result;
+        }
+      }
+      S= subset (v, s, TT, nosubset);
+      if (nosubset) break;
+      subsetDeg= subsetDegree (S);
+      // skip those combinations that are not possible
+      if (!degs.find (subsetDeg)) 
+        continue;
+      else 
+      { 
+        g= prodMod0 (S, M);
+        g= mod (g*LCBuf, M);
+        g /= content (g);
+        if (fdivides (LC (g), LCBuf))
+        {
+          S.insert (LCBuf);
+          g= prodMod (S, M);
+          S.removeFirst();
+          g /= content (g, Variable (1));
+          if (fdivides (g, buf)) 
+          {
+            buf2= g (y - eval, y);
+            appendTestMapDown (result, buf2, info, source, dest); 
+            if (!k && beta == Variable (1)) 
+            {
+              if (degree (buf2, alpha) < degMipoBeta)
+              {  
+                buf /= g;
+                LCBuf= LC (buf, Variable (1));
+                recombination= true;
+              }
+            }
+            else 
+            {
+              if (!isInExtension (buf2, delta, k))
+              {
+                buf /= g;
+                LCBuf= LC (buf, Variable (1));
+                recombination= true;
+              }
+            }
+            T= Difference (T, S);
+            // compute new possible degree pattern
+            bufDegs2= DegreePattern (T);
+            bufDegs1.intersect (bufDegs2);
+            bufDegs1.refine ();
+            if (T.length() < 2*s || T.length() == s || 
+                bufDegs1.getLength() == 1) 
+            {
+              if (recombination)
+              {
+                appendTestMapDown (result, buf (y - eval, y), info, source, 
+                                   dest);
+                return result; 
+              }
+              else
+              {
+                appendMapDown (result, F (y - eval, y), info, source, dest);
+                return result;                
+              }
+            }
+            TT= copy (T);
+            indexUpdate (v, s, T.length(), nosubset);
+            if (nosubset) break;        
+          }
+        }
+      }
+    }
+    s++;
+    if (T.length() < 2*s || T.length() == s) 
+    {
+      if (recombination)
+      {
+        appendTestMapDown (result, buf (y - eval, y), info, source, dest);
+        return result;
+      }
+      else
+      {
+        appendMapDown (result, F (y - eval, y), info, source, dest);
+        return result;
+      }
+    }
+    int v [T.length()];
+    for (int i= 0; i < T.length(); i++)
+      v[i]= 0;
+    nosubset= false;
+  }
+  if (T.length() < 2*s) 
+    appendMapDown (result, F (y - eval, y), info, source, dest);
+
+  return result; 
+}     
+
+/// naive factor recombination as decribed in "Factoring 
+/// multivariate polynomials over a finite field" by L Bernardin. 
+CFList 
+factorRecombination (const CFList& factors, const CanonicalForm& F, 
+                     const CanonicalForm& M, const DegreePattern& degs)
+{
+  if (factors.length() == 0)
+    return CFList ();
+  if (degs.getLength() <= 1 || factors.length() == 1)  
+    return CFList(F);
+  DEBOUTLN (cerr, "LC (F, 1)*prodMod (factors, M) == F " << 
+            (LC (F, 1)*prodMod (factors, M) == F));
+  CFList T, S;
+
+  T= factors;
+  int s= 1;
+  CFList result;
+  CanonicalForm LCBuf= LC (F, Variable (1));
+  CanonicalForm g, buf= F;
+  int v [T.length()];
+  for (int i= 0; i < T.length(); i++)
+    v[i]= 0;
+  bool nosubset= false;
+  CFArray TT;
+  DegreePattern bufDegs1, bufDegs2;
+  bufDegs1= degs;
+  int subsetDeg;
+  TT= copy (factors);
+  bool recombination= false;
+  while (T.length() >= 2*s) 
+  {
+    while (nosubset == false) 
+    {
+      if (T.length() == s) 
+      {
+        if (recombination)
+        {
+          T.insert (LCBuf);
+          g= prodMod (T, M);
+          T.removeFirst();
+          result.append (g/content (g, Variable (1)));
+          return result;
+        }
+        else
+          return CFList (F);
+      }
+      S= subset (v, s, TT, nosubset);
+      if (nosubset) break;
+      subsetDeg= subsetDegree (S);
+      // skip those combinations that are not possible
+      if (!degs.find (subsetDeg)) 
+        continue;
+      else 
+      { 
+        g= prodMod0 (S, M);
+        g= mod (g*LCBuf, M);
+        g /= content (g);
+        if (fdivides (LC(g), LCBuf)) 
+        {
+          S.insert (LCBuf);
+          g= prodMod (S, M);
+          S.removeFirst();
+          g /= content (g, Variable (1));
+          if (fdivides (g, buf))
+          {
+            recombination= true;
+            result.append (g);
+            buf /= g;
+            LCBuf= LC (buf, Variable(1));
+            T= Difference (T, S);
+
+            // compute new possible degree pattern
+            bufDegs2= DegreePattern (T);
+            bufDegs1.intersect (bufDegs2);
+            bufDegs1.refine ();
+            if (T.length() < 2*s || T.length() == s || 
+                bufDegs1.getLength() == 1) 
+            {
+              if (recombination)
+              {
+                result.append (buf);
+                return result;
+              }
+              else
+                return CFList (F); 
+            }
+            TT= copy (T);
+            indexUpdate (v, s, T.length(), nosubset);
+            if (nosubset) break; 
+          }       
+        }
+      }
+    }
+    s++;
+    if (T.length() < 2*s || T.length() == s) 
+    {
+      if (recombination)
+      {
+        result.append (buf);
+        return result;
+      }
+      else
+        return CFList (F);
+    }
+    int v [T.length()];
+    for (int i= 0; i < T.length(); i++)
+      v[i]= 0;
+    nosubset= false;
+  }
+  if (T.length() < 2*s)  
+    result.append (F);
+
+  return result; 
+} 
+
+Variable chooseExtension (const CanonicalForm & A, const Variable & alpha) 
+{
+  int p= getCharacteristic();
+  ZZ NTLp= to_ZZ (p);
+  ZZ_p::init (NTLp);
+  ZZ_pX NTLirredpoly; 
+  int d= degree (A); 
+  int m= 1;
+  if (alpha != Variable (1)) 
+    m= degree (getMipo (alpha));   
+  int i= (int) floor((double) d/(double) m) - 1;
+  if (i < 2) 
+    i= 2;
+  if (i > 8)
+    i= i/4;
+  BuildIrred (NTLirredpoly, i*m);
+  Variable x= A.mvar();
+  CanonicalForm newMipo= convertNTLZZpX2CF (NTLirredpoly, x);
+  return rootOf (newMipo);
+}
+
+CFList 
+earlyFactorDetection (CanonicalForm& F, CFList& factors,int& adaptedLiftBound, 
+                      DegreePattern& degs, bool& success, int deg)
+{
+  DegreePattern bufDegs1= degs;
+  DegreePattern bufDegs2;
+  CFList result;
+  CFList T= factors;
+  CanonicalForm buf= F;
+  CanonicalForm LCBuf= LC (buf, Variable (1));
+  CanonicalForm g;
+  CanonicalForm M= power (F.mvar(), deg);
+  adaptedLiftBound= 0;
+  int d;
+  if (degree (LCBuf) == degree (F))
+    d= degree (F);
+  else
+    d= degree (F) + degree (LCBuf);
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+  {
+    if (!bufDegs1.find (degree (i.getItem(), 1)))
+      continue;
+    else
+    {
+      g= i.getItem() (0, 1);
+      g *= LCBuf;
+      g= mod (g, M);
+      if (fdivides (LC (g), LCBuf))
+      {
+        g= mulMod2 (i.getItem(), LCBuf, M);
+        g /= content (g, Variable (1));
+        if (fdivides (g, buf))
+        {
+          result.append (g);
+          buf /= g;
+          d -= degree (g) + degree (LC (g, Variable (1)));
+          LCBuf= LC (buf, Variable (1)); 
+          T= Difference (T, CFList (i.getItem()));
+
+          // compute new possible degree pattern
+          bufDegs2= DegreePattern (T);
+          bufDegs1.intersect (bufDegs2);
+          bufDegs1.refine (); 
+          if (bufDegs1.getLength() <= 1)
+          {
+            result.append (buf);
+            break;
+          }
+        }
+      }
+    }
+  }
+  adaptedLiftBound= d + 1;
+  if (d < deg)
+  {
+    factors= T;
+    degs= bufDegs1;
+    F= buf;
+    success= true;
+  }
+  if (bufDegs1.getLength() <= 1)
+    degs= bufDegs1;
+  return result;
+}
+
+CFList 
+extEarlyFactorDetection (CanonicalForm& F, CFList& factors,  
+                         int& adaptedLiftBound, DegreePattern& degs, 
+                         bool& success, const ExtensionInfo& info,
+                         const CanonicalForm& eval, int deg)
+{
+  Variable alpha= info.getAlpha();
+  Variable beta= info.getBeta();
+  CanonicalForm gamma= info.getGamma();
+  CanonicalForm delta= info.getDelta();
+  int k= info.getGFDegree();
+  DegreePattern bufDegs1= degs, bufDegs2;
+  CFList result; 
+  CFList T= factors;
+  Variable y= F.mvar();
+  CanonicalForm buf= F, LCBuf= LC (buf, Variable (1)), g, buf2;
+  CanonicalForm M= power (y, deg);
+  adaptedLiftBound= 0;
+  int d;
+  if (degree (F) == degree (LCBuf))
+    d= degree (F);
+  else
+    d= degree (F) + degree (LCBuf);
+  CFList source, dest;
+  int degMipoBeta;
+  if (!k && beta == Variable(1))
+    degMipoBeta= 1;
+  else if (!k && beta != Variable(1))
+    degMipoBeta= degree (getMipo (beta));
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+  {
+    if (!bufDegs1.find (degree (i.getItem(), 1)))
+      continue;
+    else
+    {
+      g= i.getItem() (0, 1);
+      g *= LCBuf;
+      g= mod (g, M);
+      if (fdivides (LC (g), LCBuf))
+      {
+        g= mulMod2 (i.getItem(), LCBuf, M);
+        g /= content (g, Variable (1));
+        if (fdivides (g, buf))
+        {
+          buf2= g (y - eval, y);
+          buf2 /= Lc (buf2);
+           
+          if (!k && beta == Variable (1)) 
+          {
+            if (degree (buf2, alpha) < degMipoBeta)
+            {
+              appendTestMapDown (result, buf2, info, source, dest);
+              buf /= g;
+              d -= degree (g) + degree (LC (g, Variable (1)));
+              LCBuf= LC (buf, Variable (1));            
+            }
+          }
+          else 
+          {
+            if (!isInExtension (buf2, delta, k))
+            {
+              appendTestMapDown (result, buf2, info, source, dest);
+              buf /= g;
+              d -= degree (g) + degree (LC (g, Variable (1)));
+              LCBuf= LC (buf, Variable (1));
+            }
+          } 
+          T= Difference (T, CFList (i.getItem()));
+
+          // compute new possible degree pattern
+          bufDegs2= DegreePattern (T);
+          bufDegs1.intersect (bufDegs2);
+          bufDegs1.refine ();
+          if (bufDegs1.getLength() <= 1)
+          {
+            buf= buf (y - eval, y);
+            buf /= Lc (buf);
+            appendMapDown (result, buf, info, source, dest);
+            break;
+          }
+        }
+      }
+    }
+  }
+  adaptedLiftBound= d + 1;
+  if (adaptedLiftBound < deg)
+  {
+    success= true;
+    factors= T;
+    degs= bufDegs1;
+    F= buf;
+  }
+  if (bufDegs1.getLength() <= 1)
+    degs= bufDegs1;
+
+  return result;
+}
+
+CFList
+henselLiftAndEarly (CanonicalForm& A, bool& earlySuccess, CFList& 
+                    earlyFactors, DegreePattern& degs, int& liftBound,
+                    const CFList& uniFactors, const ExtensionInfo& info, 
+                    const CanonicalForm& eval)
+{
+  Variable alpha= info.getAlpha();
+  Variable beta= info.getBeta();
+  CanonicalForm gamma= info.getGamma();
+  CanonicalForm delta= info.getDelta();
+  int k= info.getGFDegree();
+  bool extension= info.isInExtension();
+
+  Variable x= Variable (1);
+  Variable y= Variable (2);
+  CFArray Pi;
+  CFList diophant;
+  CFList bufUniFactors= uniFactors;
+  bufUniFactors.insert (LC (A, x));
+  CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1);
+  earlySuccess= false;
+  int newLiftBound= 0;
+  int smallFactorDeg= 11; //this is a tunable parameter
+  if (smallFactorDeg >= liftBound)
+    henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M);
+  else if (smallFactorDeg >= degree (A, y) + 1)
+  {
+    henselLift12 (A, bufUniFactors, degree (A, y) + 1, Pi, diophant, M);
+    if (!extension)
+      earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
+                     degs, earlySuccess, degree (A, y) + 1);
+    else
+      earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
+                     newLiftBound, degs, earlySuccess, info, eval, 
+                     degree (A, y) + 1);
+    if (degs.getLength() > 1 && !earlySuccess)
+    {
+      if (newLiftBound > degree (A, y) + 1)
+      {
+        liftBound= newLiftBound;
+        bufUniFactors.insert (LC(A, x));
+        henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound, 
+                            Pi, diophant, M);
+      }
+    }
+    else if (earlySuccess)
+      liftBound= newLiftBound;
+  }
+  else if (smallFactorDeg < degree (A, y) + 1)
+  {  
+    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
+    if (!extension)
+      earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
+                                           degs, earlySuccess,
+                                           smallFactorDeg);
+    else
+      earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
+                     newLiftBound, degs, earlySuccess,
+                     info, eval, smallFactorDeg);
+    if (degs.getLength() > 1 && !earlySuccess)
+    {
+      bufUniFactors.insert (LC (A, x));
+      henselLiftResume12 (A, bufUniFactors, smallFactorDeg, degree (A, y) 
+                          + 1, Pi, diophant, M);
+      if (!extension)
+	earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
+                       degs, earlySuccess, degree (A, y) + 1);
+      else
+	earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
+                       newLiftBound, degs, earlySuccess,
+                       info, eval, degree(A,y) + 1);
+      if (degs.getLength() > 1 && !earlySuccess)
+      {
+	if (newLiftBound > degree (A, y) + 1)
+	{
+          if (newLiftBound < newLiftBound)
+	    liftBound= newLiftBound;
+	  bufUniFactors.insert (LC(A, x));
+	  henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
+                              Pi, diophant, M);
+	}
+      }
+      else if (earlySuccess)
+        liftBound= newLiftBound;
+    }
+    else if (earlySuccess)
+      liftBound= newLiftBound;
+  }
+  if (newLiftBound > 0)
+    liftBound= newLiftBound;
+  return bufUniFactors;
+}
+
+CFList 
+extBiFactorize (const CanonicalForm& F, const ExtensionInfo& info);
+
+/// bivariate factorization over finite fields as decribed in "Factoring 
+/// multivariate polynomials over a finite field" by L Bernardin.
+CFList 
+biFactorize (const CanonicalForm& F, const ExtensionInfo& info)
+{
+  if (F.inCoeffDomain())
+    return CFList(F);
+  
+  CanonicalForm A= F;
+  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+
+  Variable alpha= info.getAlpha();
+  Variable beta= info.getBeta();
+  CanonicalForm gamma= info.getGamma();
+  CanonicalForm delta= info.getDelta();
+  int k= info.getGFDegree();
+  bool extension= info.isInExtension();
+  if (A.isUnivariate()) 
+  {
+    if (extension == false) 
+      return uniFactorizer (F, alpha, GF);
+    else 
+    {
+      CFList source, dest;
+      A= mapDown (A, info, source, dest);
+      return uniFactorizer (A, beta, GF);
+    }
+  }
+
+  CFMap N;
+  A= compress (A, N);
+  Variable y= A.mvar();
+
+  if (y.level() > 2) return CFList (F);
+  Variable x= Variable (1);
+
+  //remove and factorize content
+  CanonicalForm contentAx= content (A, x);
+  CanonicalForm contentAy= content (A);
+
+  A= A/(contentAx*contentAy);
+  CFList contentAxFactors, contentAyFactors;
+
+  if (!extension) 
+  { 
+    contentAxFactors= uniFactorizer (contentAx, alpha, GF);
+    contentAyFactors= uniFactorizer (contentAy, alpha, GF);
+  }
+
+  //trivial case
+  CFList factors;
+  if (A.inCoeffDomain()) 
+  {
+    append (factors, contentAxFactors);
+    append (factors, contentAyFactors);
+    decompress (factors, N); 
+    return factors;
+  }
+  else if (A.isUnivariate()) 
+  {
+    factors= uniFactorizer (A, alpha, GF);
+    append (factors, contentAxFactors);
+    append (factors, contentAyFactors);
+    decompress (factors, N); 
+    return factors;
+  }
+
+  // check derivatives 
+  bool derivXZero= false;
+  CanonicalForm derivX= deriv (A, x);
+  CanonicalForm gcdDerivX;
+  if (derivX.isZero())
+    derivXZero= true;
+  else
+  {
+    if (GF)
+      gcdDerivX= GCD_GF (A, derivX);
+    else if (alpha == x)
+      gcdDerivX= GCD_small_p (A, derivX);
+    else 
+      gcdDerivX= GCD_Fp_extension (A, derivX, alpha);
+    if (degree (gcdDerivX) > 0)
+    {
+      CanonicalForm g= A/gcdDerivX;
+      CFList factorsG= 
+      Union (biFactorize (g, info), biFactorize (gcdDerivX, info));
+      append (factorsG, contentAxFactors);
+      append (factorsG, contentAyFactors);
+      decompress (factorsG, N);
+      normalize (factors);
+      return factorsG;    
+    }
+  }
+  bool derivYZero= false;
+  CanonicalForm derivY= deriv (A, y);
+  CanonicalForm gcdDerivY;
+  if (derivY.isZero())
+    derivYZero= true;
+  else
+  {
+    if (GF)
+      gcdDerivY= GCD_GF (A, derivY);
+    else if (alpha == x)
+      gcdDerivY= GCD_small_p (A, derivY);
+    else 
+      gcdDerivY= GCD_Fp_extension (A, derivY, alpha);
+    if (degree (gcdDerivY) > 0)
+    {
+      CanonicalForm g= A/gcdDerivY;
+      CFList factorsG= 
+      Union (biFactorize (g, info), biFactorize (gcdDerivY, info));
+      append (factorsG, contentAxFactors);
+      append (factorsG, contentAyFactors);
+      decompress (factorsG, N);
+      normalize (factors);
+      return factorsG;    
+    }
+  }
+  //main variable is chosen s.t. the degree in x is minimal
+  bool swap= false;
+  if ((degree (A) < degree (A, x)) || derivXZero) 
+  {
+    if (!derivYZero) 
+    { 
+      A= swapvar (A, y, x);
+      swap= derivXZero;  
+      derivXZero= derivYZero;
+      derivYZero= swap;
+      swap= true;
+    }
+  }
+
+  bool fail;
+  CanonicalForm Aeval, evaluation, bufAeval, bufEvaluation, buf;
+  CFList uniFactors, list, bufUniFactors;
+  DegreePattern degs;
+  DegreePattern bufDegs;
+
+  bool fail2= false;
+  CanonicalForm Aeval2, evaluation2, bufAeval2, bufEvaluation2;
+  CFList bufUniFactors2, list2, uniFactors2;
+  DegreePattern degs2;
+  DegreePattern bufDegs2;
+  bool swap2= false;
+
+  // several univariate factorizations to obtain more information about the 
+  // degree pattern therefore usually less combinations have to be tried during
+  // the recombination process
+  int factorNums= 5;
+  if (factorNums < (int) ceil (log (totaldegree (A)))) 
+    factorNums= (int) ceil (log (totaldegree (A)));
+  for (int i= 0; i < factorNums; i++) 
+  {
+    bufAeval= A; 
+    bufEvaluation= evalPoint (A, bufAeval, alpha, list, GF, fail);
+    if (!derivXZero && !fail2)
+    {
+      buf= swapvar (A, x, y);
+      bufAeval2= buf;
+      bufEvaluation2= evalPoint (buf, bufAeval2, alpha, list2, GF, fail2);
+    }
+    // first try to change main variable if there is no valid evaluation point
+    if (fail && (i == 0)) 
+    {
+      if (!derivXZero)
+        bufEvaluation= evalPoint (buf, bufAeval, alpha, list, GF, fail);
+      else
+        fail= true;
+
+      if (!fail) 
+      { 
+        A= buf;
+        swap2= true;
+      }
+    }
+
+    // if there is no valid evaluation point pass to a field extension
+    if (fail && (i == 0)) 
+    {
+      factors= extBiFactorize (A, info);
+      appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
+                            swap, swap2, N);
+      normalize (factors);
+      return factors;
+    }
+
+    // there is at least one valid evaluation point 
+    // but we do not compute more univariate factorization over an extension 
+    if (fail && (i != 0)) 
+      break;
+
+    // univariate factorization
+    TIMING_START (fac_uni_factorizer);
+    bufUniFactors= uniFactorizer (bufAeval, alpha, GF);
+    TIMING_END_AND_PRINT (fac_uni_factorizer, 
+                          "time for univariate factorization: ");
+    DEBOUTLN (cerr, "Lc (bufAeval)*prod (bufUniFactors)== bufAeval " << 
+              prod (bufUniFactors)*Lc (bufAeval) == bufAeval);
+
+    if (!derivXZero && !fail2)
+    {
+      TIMING_START (fac_uni_factorizer);
+      bufUniFactors2= uniFactorizer (bufAeval2, alpha, GF);
+      TIMING_END_AND_PRINT (fac_uni_factorizer, 
+                            "time for univariate factorization in y: ");
+      DEBOUTLN (cerr, "Lc (Aeval2)*prod (uniFactors2)== Aeval2 " << 
+                prod (bufUniFactors2)*Lc (bufAeval2) == bufAeval2);
+    }
+
+    if (bufUniFactors.length() == 1 || 
+        (!fail2 && !derivXZero && (bufUniFactors2.length() == 1))) 
+    {
+      if (extension)
+      {
+        CFList source, dest;
+        ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k, 
+                             info.getGFName(), info.isInExtension());
+        appendMapDown (factors, A, info2, source, dest);
+      }
+      else 
+        factors.append (A);
+
+      appendSwapDecompress (factors, contentAxFactors, contentAyFactors, 
+                            swap, swap2, N);
+
+      normalize (factors);
+      return factors;
+    }
+
+    // degree analysis    
+    bufDegs = DegreePattern (bufUniFactors); 
+    if (!derivXZero && !fail2)
+      bufDegs2= DegreePattern (bufUniFactors2);
+
+    if (i == 0) 
+    {
+      Aeval= bufAeval;
+      evaluation= bufEvaluation;
+      uniFactors= bufUniFactors;
+      degs= bufDegs;
+      if (!derivXZero && !fail2)
+      {
+        Aeval2= bufAeval2;
+        evaluation2= bufEvaluation2;
+        uniFactors2= bufUniFactors2;
+        degs2= bufDegs2;
+      }
+    }
+    else 
+    {
+      degs.intersect (bufDegs);
+      if (!derivXZero && !fail2) 
+      {
+        degs2.intersect (bufDegs2);
+        if (bufUniFactors2.length() < uniFactors2.length()) 
+        {  
+          uniFactors2= bufUniFactors2;
+          Aeval2= bufAeval2;
+          evaluation2= bufEvaluation2;
+        }
+      }   
+      if (bufUniFactors.length() < uniFactors.length()) 
+      { 
+        uniFactors= bufUniFactors;
+        Aeval= bufAeval;
+        evaluation= bufEvaluation;
+      }
+    }
+    list.append (bufEvaluation);
+    if (!derivXZero && !fail2)
+      list2.append (bufEvaluation2);
+  }
+  
+  if (!derivXZero && !fail2)
+  {
+    if (uniFactors.length() > uniFactors2.length() || 
+        (uniFactors.length() == uniFactors2.length() 
+         && degs.getLength() > degs2.getLength())) 
+    {
+      degs= degs2;
+      uniFactors= uniFactors2;
+      evaluation= evaluation2;
+      Aeval= Aeval2;
+      A= buf;
+      swap2= true;
+    }
+  }
+ 
+  if (degs.getLength() == 1) // A is irreducible 
+  {
+    if (extension) 
+    {
+      CFList source, dest;
+      ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k, 
+                           info.getGFName(), info.isInExtension());
+      appendMapDown (factors, A, info2, source, dest);
+    }
+    else 
+      factors.append (A);
+    appendSwapDecompress (factors, contentAxFactors, contentAyFactors, 
+                            swap, swap2, N);
+    normalize (factors);
+    return factors;
+  }
+
+  A= A (y + evaluation, y);
+
+  int liftBound;
+  if (degree (A, y) == degree (LC (A, x)))
+    liftBound= degree (LC (A, x)) + 1;
+  else
+    liftBound= degree (A, y) + 1 + degree (LC(A, x)); 
+
+  DEBOUTLN (cerr, "uniFactors= " << uniFactors);
+  bool earlySuccess= false;
+  CFList earlyFactors; 
+  TIMING_START (fac_hensel_lift);
+  uniFactors= henselLiftAndEarly 
+               (A, earlySuccess, earlyFactors, degs, liftBound,
+                uniFactors, info, evaluation); 
+  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
+  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
+
+  CanonicalForm MODl= power (y, liftBound);
+  TIMING_START (fac_factor_recombination);
+  if (!extension)
+    factors= factorRecombination (uniFactors, A, MODl, degs);
+  else 
+    factors= extFactorRecombination (uniFactors, A, MODl, info, degs, 
+                                     evaluation);
+  TIMING_END_AND_PRINT (fac_factor_recombination,
+                        "time for factor recombination: ");
+  if (earlySuccess)
+    factors= Union (earlyFactors, factors);
+  else if (!earlySuccess && degs.getLength() == 1)
+    factors= earlyFactors;
+  
+  if (!extension)
+  {
+    for (CFListIterator i= factors; i.hasItem(); i++)
+      i.getItem()= i.getItem() (y - evaluation, y);
+  }
+
+  appendSwapDecompress (factors, contentAxFactors, contentAyFactors, 
+                        swap, swap2, N);
+  normalize (factors); 
+
+  return factors;
+}
+
+CFList 
+extBiFactorize (const CanonicalForm& F, const ExtensionInfo& info) 
+{
+  
+  CanonicalForm A= F;
+  Variable alpha= info.getAlpha();
+  Variable beta= info.getBeta();
+  int k= info.getGFDegree();
+  char cGFName= info.getGFName();
+  CanonicalForm delta= info.getDelta();
+
+  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+  Variable x= Variable (1); 
+  CFList factors;
+  if (!GF && alpha == x)  // we are in F_p
+  {
+    bool extension= true;         
+    int p= getCharacteristic();
+    if (p*p < (1<<16)) // pass to GF if possible 
+    { 
+      setCharacteristic (getCharacteristic(), 2, 'Z');
+      A= A.mapinto();
+      ExtensionInfo info= ExtensionInfo (extension);
+      factors= biFactorize (A, info);
+
+      Variable vBuf= rootOf (gf_mipo);
+      setCharacteristic (getCharacteristic());
+      for (CFListIterator j= factors; j.hasItem(); j++) 
+	j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
+    }
+    else // not able to pass to GF, pass to F_p(\alpha) 
+    { 
+      CanonicalForm mipo= randomIrredpoly (3, Variable (1));
+      Variable v= rootOf (mipo);
+      ExtensionInfo info= ExtensionInfo (v);
+      factors= biFactorize (A, info);
+    } 
+    return factors;
+  }
+  else if (!GF && (alpha != x)) // we are in F_p(\alpha)
+  { 
+    if (k == 1) // need factorization over F_p
+    {
+      int extDeg= degree (getMipo (alpha));
+      extDeg++; 
+      CanonicalForm mipo= randomIrredpoly (extDeg + 1, Variable (1));
+      Variable v= rootOf (mipo);
+      ExtensionInfo info= ExtensionInfo (v);
+      factors= biFactorize (A, info);
+    }
+    else 
+    {
+      if (beta == Variable (1))
+      {
+        Variable v= chooseExtension (A, alpha);
+        CanonicalForm primElem, imPrimElem;
+        bool primFail= false;
+        Variable vBuf;
+        primElem= primitiveElement (alpha, vBuf, primFail);
+        ASSERT (!primFail, "failure in integer factorizer");
+        if (primFail)
+          ; //ERROR
+        else
+          imPrimElem= mapPrimElem (primElem, vBuf, v);
+
+        CFList source, dest;
+        CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem, 
+                                   source, dest);
+        ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
+        factors= biFactorize (bufA, info);
+      }
+      else
+      {
+        Variable v= chooseExtension (A, alpha);
+        CanonicalForm primElem, imPrimElem;
+        bool primFail= false;
+        Variable vBuf;
+        ASSERT (!primFail, "failure in integer factorizer");
+        if (primFail)
+          ; //ERROR
+        else
+          imPrimElem= mapPrimElem (delta, beta, v); //oder mapPrimElem (primElem, vBuf, v);
+
+        CFList source, dest;
+        CanonicalForm bufA= mapDown (A, info, source, dest);
+        source= CFList();
+        dest= CFList();
+        bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
+        ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
+        factors= biFactorize (bufA, info);
+      }
+    }
+    return factors;
+  }
+  else // we are in GF (p^k)
+  { 
+    int p= getCharacteristic();
+    int extensionDeg= getGFDegree();
+    bool extension= true;
+    if (k == 1) // need factorization over F_p
+    {
+      extensionDeg++;
+      if (ipower (p, extensionDeg) < (1<<16))  
+      // pass to GF(p^k+1) 
+      {
+        setCharacteristic (p);
+	Variable vBuf= rootOf (gf_mipo);
+	A= GF2FalphaRep (A, vBuf);
+	setCharacteristic (p, extensionDeg, 'Z');
+        ExtensionInfo info= ExtensionInfo (extension);
+	factors= biFactorize (A.mapinto(), info);
+      }
+      else // not able to pass to another GF, pass to F_p(\alpha)
+      { 
+	setCharacteristic (p);
+	Variable vBuf= rootOf (gf_mipo);
+	A= GF2FalphaRep (A, vBuf);
+	Variable v= chooseExtension (A, beta); 
+        ExtensionInfo info= ExtensionInfo (v, extension);
+	factors= biFactorize (A, info);
+      }
+    }
+    else // need factorization over GF (p^k)
+    {
+      if (ipower (p, 2*extensionDeg) < (1<<16)) 
+      // pass to GF (p^2k)
+      { 
+	setCharacteristic (p, 2*extensionDeg, 'Z');
+        ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
+	factors= biFactorize (GFMapUp (A, extensionDeg), info); 
+	setCharacteristic (p, extensionDeg, cGFName);
+      }
+      else // not able to pass to GF (p^2k), pass to F_p (\alpha)
+      { 
+	setCharacteristic (p);
+	Variable v1= rootOf (gf_mipo); 
+	A= GF2FalphaRep (A, v1);
+	Variable v2= chooseExtension (A, v1); 
+        CanonicalForm primElem, imPrimElem;
+        bool primFail= false;
+        Variable vBuf;
+        primElem= primitiveElement (v1, vBuf, primFail);
+        ASSERT (!primFail, "failure in integer factorizer");
+        if (primFail)
+          ; //ERROR
+        else
+          imPrimElem= mapPrimElem (primElem, vBuf, v2);
+
+        CFList source, dest;
+        CanonicalForm bufA= mapUp (A, v1, v2, primElem, imPrimElem, 
+                                     source, dest);
+        ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
+	factors= biFactorize (bufA, info); 
+	setCharacteristic (p, k, cGFName);
+	for (CFListIterator i= factors; i.hasItem(); i++) 
+	  i.getItem()= Falpha2GFRep (i.getItem());
+      }
+    }
+    return factors;
+  }
+}
+
+#endif
+/* HAVE_NTL */
+
+
Index: factory/facFqBivar.h
===================================================================
--- factory/facFqBivar.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/facFqBivar.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,391 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file facFqBivar.h
+ * 
+ * This file provides functions for factorizing a bivariate polynomial over
+ * \f$ F_{p} \f$ , \f$ F_{p}(\alpha ) \f$ or GF.   
+ *
+ * ABSTRACT: In contrast to biFactorizer() in facFqFactorice.cc we evaluate and 
+ * factorize the polynomial in both variables. So far factor recombination is 
+ * done naive!          
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#ifndef FAC_FQ_BIVAR_H
+#define FAC_FQ_BIVAR_H
+
+#include <config.h>
+
+#include "assert.h"
+
+#include "facFqBivarUtil.h"
+#include "DegreePattern.h"
+#include "ExtensionInfo.h"
+#include "ExtensionInfo.cc"
+#include "cf_util.h"
+#include "facFqSquarefree.h"
+
+/// Factorization of a squarefree bivariate polynomials over an arbitrary finite 
+/// field, information on the current field we work over is in @a info. @a info
+/// may also contain information about the initial field if initial and current
+/// field do not coincide. In this case the current field is an extension of the 
+/// initial field and the factors returned are factors of F over the initial
+/// field. 
+/// 
+/// @return @a biFactorize returns a list of factors of F. If F is not monic
+///         its leading coefficient is not outputted. 
+/// @sa extBifactorize()
+CFList 
+biFactorize (const CanonicalForm& F,       ///< [in] a bivariate poly
+             const ExtensionInfo& info     ///< [in] information about extension
+            );
+
+/// factorize a squarefree bivariate polynomial over \f$ F_{p} \f$.
+///
+/// @return @a FpBiSqrfFactorize returns a list of monic factors, the first 
+///         element is the leading coefficient.     
+/// @sa FqBiSqrfFactorize(), GFBiSqrfFactorize()
+inline
+CFList FpBiSqrfFactorize (const CanonicalForm & F ///< [in] a bivariate poly
+                         ) 
+{
+  ExtensionInfo info= ExtensionInfo (false);
+  CFList result= biFactorize (F, info);
+  result.insert (Lc(F));
+  return result;
+}
+
+/// factorize a squarefree bivariate polynomial over \f$ F_{p}(\alpha ) \f$.
+///
+/// @return @a FqBiSqrfFactorize returns a list of monic factors, the first 
+///         element is the leading coefficient.      
+/// @sa FpBiSqrfFactorize(), GFBiSqrfFactorize()
+inline
+CFList FqBiSqrfFactorize (const CanonicalForm & F, ///< [in] a bivariate poly
+                          const Variable& alpha    ///< [in] algebraic variable
+                         ) 
+{
+  ExtensionInfo info= ExtensionInfo (alpha, false);
+  CFList result= biFactorize (F, info);
+  result.insert (Lc(F));
+  return result;
+}
+
+/// factorize a squarefree bivariate polynomial over GF
+///
+/// @return @a GFBiSqrfFactorize returns a list of monic factors, the first 
+///         element is the leading coefficient.      
+/// @sa FpBiSqrfFactorize(), FqBiSqrfFactorize() 
+inline
+CFList GFBiSqrfFactorize (const CanonicalForm & F ///< [in] a bivariate poly
+                         ) 
+{
+  ASSERT (CFFactory::gettype() == GaloisFieldDomain, 
+          "GF as base field expected");
+  ExtensionInfo info= ExtensionInfo (getGFDegree(), gf_name, false);
+  CFList result= biFactorize (F, info);
+  result.insert (Lc(F));
+  return result;
+}
+
+/// factorize a bivariate polynomial over \f$ F_{p} \f$
+///
+/// @return @a FpBiFactorize returns a list of monic factors with 
+///         multiplicity, the first element is the leading coefficient.      
+/// @sa FqBiFactorize(), GFBiFactorize() 
+inline
+CFFList FpBiFactorize (const CanonicalForm & F ///< [in] a bivariate poly
+                      ) 
+{
+  ExtensionInfo info= ExtensionInfo (false);
+  bool GF= false;
+  CanonicalForm LcF= Lc (F); 
+  CanonicalForm pthRoot, A;
+  CanonicalForm sqrfP= sqrfPart (F/Lc(F), pthRoot, info.getAlpha());
+  CFList buf, bufRoot;
+  CFFList result, resultRoot;
+  int p= getCharacteristic();
+  int l;
+  if (degree (pthRoot) > 0)
+  {
+    pthRoot= maxpthRoot (pthRoot, p, l);
+    result= FpBiFactorize (pthRoot);
+    result.removeFirst();
+    for (CFFListIterator i= result; i.hasItem(); i++)
+      i.getItem()= CFFactor (i.getItem().factor(), i.getItem().exp()*l*p);
+    result.insert (CFFactor (LcF, 1));
+    return result;
+  }
+  else
+  {  
+    buf= biFactorize (sqrfP, info);
+    A= F/LcF;
+    result= multiplicity (A, buf);
+  }
+  if (degree (A) > 0)
+  {
+    resultRoot= FpBiFactorize (A);
+    resultRoot.removeFirst();
+    result= Union (result, resultRoot);
+  }
+  result.insert (CFFactor (LcF, 1));
+  return result;
+}
+
+/// factorize a bivariate polynomial over \f$ F_{p}(\alpha ) \f$
+///
+/// @return @a FqBiFactorize returns a list of monic factors with 
+///         multiplicity, the first element is the leading coefficient.      
+/// @sa FpBiFactorize(), FqBiFactorize() 
+inline
+CFFList FqBiFactorize (const CanonicalForm & F, ///< [in] a bivariate poly
+                       const Variable & alpha   ///< [in] algebraic variable
+                      )
+{
+  ExtensionInfo info= ExtensionInfo (alpha, false);
+  bool GF= false;
+  CanonicalForm LcF= Lc (F); 
+  CanonicalForm pthRoot, A;
+  CanonicalForm sqrfP= sqrfPart (F/Lc(F), pthRoot, alpha);
+  CFList buf, bufRoot;
+  CFFList result, resultRoot;
+  int p= getCharacteristic();
+  int q= ipower (p, degree (getMipo (alpha)));
+  int l;
+  if (degree (pthRoot) > 0)
+  {
+    pthRoot= maxpthRoot (pthRoot, q, l);
+    result= FqBiFactorize (pthRoot, alpha);
+    result.removeFirst();
+    for (CFFListIterator i= result; i.hasItem(); i++)
+      i.getItem()= CFFactor (i.getItem().factor(), i.getItem().exp()*l*p);
+    result.insert (CFFactor (LcF, 1));
+    return result;
+  }
+  else
+  {  
+    buf= biFactorize (sqrfP, info);
+    A= F/LcF;
+    result= multiplicity (A, buf);
+  }
+  if (degree (A) > 0)
+  {
+    resultRoot= FqBiFactorize (A, alpha);
+    resultRoot.removeFirst();
+    result= Union (result, resultRoot);
+  }
+  result.insert (CFFactor (LcF, 1));
+  return result;
+}
+
+/// factorize a bivariate polynomial over GF
+/// 
+/// @return @a GFBiFactorize returns a list of monic factors with 
+///         multiplicity, the first element is the leading coefficient.      
+/// @sa FpBiFactorize(), FqBiFactorize() 
+inline
+CFFList GFBiFactorize (const CanonicalForm & F ///< [in] a bivariate poly
+                      ) 
+{
+  ASSERT (CFFactory::gettype() == GaloisFieldDomain, 
+          "GF as base field expected");
+  ExtensionInfo info= ExtensionInfo (getGFDegree(), gf_name, false);
+  bool GF= true;
+  CanonicalForm LcF= Lc (F);
+  CanonicalForm pthRoot, A;
+  CanonicalForm sqrfP= sqrfPart (F/LcF, pthRoot, info.getAlpha());
+  CFList buf;
+  CFFList result, resultRoot;
+  int p= getCharacteristic();
+  int q= ipower (p, getGFDegree());
+  int l;
+  if (degree (pthRoot) > 0)
+  {
+    pthRoot= maxpthRoot (pthRoot, q, l);
+    result= GFBiFactorize (pthRoot);
+    result.removeFirst();
+    for (CFFListIterator i= result; i.hasItem(); i++)
+      i.getItem()= CFFactor (i.getItem().factor(), i.getItem().exp()*l*p);
+    result.insert (CFFactor (LcF, 1));
+    return result;
+  }
+  else
+  {
+    buf= biFactorize (sqrfP, info);
+    A= F/LcF;
+    result= multiplicity (A, buf);
+  }
+  if (degree (A) > 0)
+  {
+    resultRoot= GFBiFactorize (A);
+    resultRoot.removeFirst();
+    result= Union (result, resultRoot);
+  }
+  result.insert (CFFactor (LcF, 1));
+  return result;
+}
+
+/// \f$ \prod_{f\in L} {f (0, x)} \ mod\ M \f$ via divide-and-conquer
+///
+/// @return @a prodMod0 computes the above defined product
+/// @sa prodMod()
+CanonicalForm prodMod0 (const CFList& L,       ///< [in] a list of compressed, 
+                                               ///< bivariate polynomials
+                        const CanonicalForm& M ///< [in] a power of Variable (2)
+                       );
+
+/// find an evaluation point p, s.t. F(p,y) is squarefree and 
+/// \f$ deg_{y} (F(p,y))= deg_{y} (F(x,y)) \f$.  
+///
+/// @return @a evalPoint returns an evaluation point, which is valid if and only
+///         if fail == false.
+CanonicalForm 
+evalPoint (const CanonicalForm& F, ///< [in] compressed, bivariate poly 
+           CanonicalForm & eval,   ///< [in,out] F (p, y)
+           const Variable& alpha,  ///< [in] algebraic variable
+           CFList& list,           ///< [in] list of points already considered
+           const bool& GF,         ///< [in] GaloisFieldDomain?
+           bool& fail              ///< [in,out] equals true, if there is no
+                                   ///< valid evaluation point 
+          );
+
+/// Univariate factorization of squarefree monic polys over finite fields via
+/// NTL. If the characteristic is even special GF2 routines of NTL are used.  
+///
+/// @return @a uniFactorizer returns a list of monic factors
+inline CFList 
+uniFactorizer (const CanonicalForm& A, ///< [in] squarefree univariate poly
+               const Variable& alpha,  ///< [in] algebraic variable
+               const bool& GF          ///< [in] GaloisFieldDomain? 
+              );
+
+/// naive factor recombination over an extension of the initial field. 
+/// Uses precomputed data to exclude combinations that are not possible.
+///
+/// @return @a extFactorRecombination returns a list of factors over the initial
+///         field, whose shift to zero is reversed. 
+/// @sa factorRecombination(), extEarlyFactorDetection()   
+inline CFList 
+extFactorRecombination (
+         const CFList& factors,     ///< [in] list of lifted factors that are 
+                                    ///< monic wrt Variable (1)
+         const CanonicalForm& F,    ///< [in] poly to be factored 
+         const CanonicalForm& M,    ///< [in] Variable (2)^liftBound
+         const ExtensionInfo& info, ///< [in] contains information about 
+                                    ///< extension
+         const CanonicalForm& eval  ///< [in] evaluation point
+                       );
+
+/// naive factor recombination. 
+/// Uses precomputed data to exclude combinations that are not possible.
+///
+/// @return @a factorRecombination returns a list of factors of F.
+/// @sa extFactorRecombination(), earlyFactorDetectection()
+inline CFList 
+factorRecombination (
+                const CFList& factors,      ///< [in] list of lifted factors
+                                            ///< that are monic wrt Variable (1) 
+                const CanonicalForm& F,     ///< [in] poly to be factored
+                const CanonicalForm& M,     ///< [in] Variable (2)^liftBound
+                const DegreePattern& degs   ///< [in] degree pattern
+                    );
+
+/// chooses a field extension.
+///
+/// @return @a chooseExtension returns an extension specified by @a beta of
+///         appropiate size
+inline
+Variable chooseExtension (
+                      const CanonicalForm & A, ///< [in] some bivariate poly
+                      const Variable & beta    ///< [in] some algebraic
+                                               ///< variable 
+                         );
+
+/// detects factors of @a F at stage @a deg of Hensel lifting.
+/// No combinations of more than one factor are tested. Lift bound and possible
+/// degree pattern are updated. 
+/// 
+/// @return @a earlyFactorDetection returns a list of factors of F (possibly in-
+///         complete), in case of success. Otherwise an empty list.
+/// @sa factorRecombination(), extEarlyFactorDetection()
+inline CFList 
+earlyFactorDetection (
+           CanonicalForm& F,       ///< [in,out] poly to be factored, returns
+                                   ///< poly divided by detected factors in case
+                                   ///< of success
+           CFList& factors,        ///< [in,out] list of factors lifted up to 
+                                   ///< @a deg, returns a list of factors
+                                   ///< without detected factors
+           int& adaptedLiftBound,  ///< [in,out] adapted lift bound 
+           DegreePattern& degs,    ///< [in,out] degree pattern, is updated
+                                   ///< whenever we find a factor
+           bool& success,          ///< [in,out] indicating success
+           int deg                 ///< [in] stage of Hensel lifting
+                     );
+
+/// detects factors of @a F at stage @a deg of Hensel lifting.
+/// No combinations of more than one factor are tested. Lift bound and possible
+/// degree pattern are updated. 
+/// 
+/// @return @a extEarlyFactorDetection returns a list of factors of F (possibly 
+///         incomplete), whose shift to zero is reversed, in case of success.
+///         Otherwise an empty list.
+/// @sa factorRecombination(), earlyFactorDetection()
+inline CFList 
+extEarlyFactorDetection (
+        CanonicalForm& F,          ///< [in,out] poly to be factored, returns  
+                                   ///< poly divided by detected factors in case 
+                                   ///< of success
+        CFList& factors,           ///< [in,out] list of factors lifted up to 
+                                   ///< @a deg, returns a list of factors
+                                   ///< without detected factors
+        int& adaptedLiftBound,     ///< [in,out] adapted lift bound
+        DegreePattern& degs,       ///< [in,out] degree pattern, is updated
+                                   ///< whenever we find a factor
+        bool& success,             ///< [in,out] indicating success
+        const ExtensionInfo& info,  ///< [in] information about extension
+        const CanonicalForm& eval, ///< [in] evaluation point
+        int deg                    ///< [in] stage of Hensel lifting
+                      );
+
+/// hensel Lifting and early factor detection 
+///
+/// @return @a henselLiftAndEarly returns monic (wrt Variable (1)) lifted 
+///         factors without factors which have been detected at an early stage
+///         of Hensel lifting
+/// @sa earlyFactorDetection(), extEarlyFactorDetection()  
+
+inline CFList
+henselLiftAndEarly (
+        CanonicalForm& A,          ///< [in,out] poly to be factored,                                     
+                                   ///< returns poly divided by detected factors  
+                                   ///< in case of success
+        bool& earlySuccess,        ///< [in,out] indicating success
+        CFList& earlyFactors,      ///< [in,out] list of factors detected
+                                   ///< at early stage of Hensel lifting 
+        DegreePattern& degs,       ///< [in,out] degree pattern 
+        int& liftBound,            ///< [in,out] (adapted) lift bound
+        const CFList& uniFactors,  ///< [in] univariate factors
+        const ExtensionInfo& info, ///< [in] information about extension
+        const CanonicalForm& eval  ///< [in] evaluation point
+                  );
+
+/// Factorization over an extension of initial field
+///
+/// @return @a extBiFactorize returns factorization of F over initial field
+/// @sa biFactorize()
+inline CFList 
+extBiFactorize (const CanonicalForm& F,    ///< [in] poly to be factored 
+                const ExtensionInfo& info  ///< [in] info about extension 
+               );
+
+
+#endif
+/* FAC_FQ_BIVAR_H */
+
Index: factory/facFqBivarUtil.cc
===================================================================
--- factory/facFqBivarUtil.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/facFqBivarUtil.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,384 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file facFqBivarUtil.cc
+ * 
+ * This file provides utility functions for bivariate factorization       
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#include <config.h>
+
+#include "cf_map.h"
+#include "cf_map_ext.h"
+#include "ftmpl_functions.h"
+#include "ExtensionInfo.h"
+#include "facFqBivarUtil.h"
+
+
+void append (CFList& factors1, const CFList& factors2) 
+{
+  for (CFListIterator i= factors2; i.hasItem(); i++) 
+  {
+    if (!i.getItem().inCoeffDomain())
+      factors1.append (i.getItem());
+  }
+  return;
+}
+
+void decompress (CFList& factors, const CFMap& N) 
+{
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+    i.getItem()= N (i.getItem());
+  return;
+}
+
+void appendSwapDecompress (CFList& factors1, const CFList& factors2, 
+                           const CFList& factors3, const bool swap1, 
+                           const bool swap2, const CFMap& N) 
+{
+  Variable x= Variable (1);
+  Variable y= Variable (2);
+  for (CFListIterator i= factors1; i.hasItem(); i++) 
+  {
+    if (swap1) 
+    {
+      if (!swap2) 
+        i.getItem()= swapvar (i.getItem(), x, y);
+    }
+    else 
+    {
+      if (swap2) 
+        i.getItem()= swapvar (i.getItem(), y, x);
+    }
+    i.getItem()= N (i.getItem());
+  }
+  for (CFListIterator i= factors2; i.hasItem(); i++) 
+    factors1.append (N (i.getItem()));
+  for (CFListIterator i= factors3; i.hasItem(); i++) 
+    factors1.append (N (i.getItem()));
+  return;
+}
+
+void swapDecompress (CFList& factors, const bool swap, const CFMap& N) 
+{
+  Variable x= Variable (1);
+  Variable y= Variable (2);
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+  {
+    if (swap) 
+      i.getItem()= swapvar (i.getItem(), x, y);
+    i.getItem()= N (i.getItem());
+  }
+  return;
+}
+
+static inline 
+bool GFInExtensionHelper (const CanonicalForm& F, const int number) 
+{
+  if (F.isOne()) return false;
+  InternalCF* buf;
+  int exp;
+  bool result= false;
+  if (F.inBaseDomain()) 
+  {
+    buf= F.getval();
+    exp= imm2int (buf);
+    if (exp%number != 0) 
+      return true;
+    else 
+      return result;
+  }
+  else 
+  {
+    for (CFIterator i= F; i.hasTerms(); i++) 
+    {
+      result= GFInExtensionHelper (i.coeff(), number);
+      if (result == true)
+        return result;
+    }
+  }
+  return result;
+}
+
+static inline 
+bool FqInExtensionHelper (const CanonicalForm& F, const CanonicalForm& gamma) 
+{
+  bool result= false;
+  if (F.inBaseDomain()) 
+    return result;
+  else if (F.inCoeffDomain()) 
+  {
+    if (!fdivides (gamma, F)) 
+      return true;
+    else 
+      return result;
+  }
+  else 
+  {
+    for (CFIterator i= F; i.hasTerms(); i++) 
+    {
+      result= FqInExtensionHelper (i.coeff(), gamma);
+      if (result == true)
+        return result; 
+    }
+  }
+  return result;
+}
+
+bool isInExtension (const CanonicalForm& F, const CanonicalForm& gamma,
+                    const int k) 
+{
+  bool result;
+  if (CFFactory::gettype() == GaloisFieldDomain) 
+  {
+    int p= getCharacteristic();
+    int orderFieldExtension= ipower (p, getGFDegree()) - 1;
+    int order= ipower (p, k) - 1;
+    int number= orderFieldExtension/order;
+    result= GFInExtensionHelper (F, number);
+    return result;
+  }
+  else 
+  {
+    result= FqInExtensionHelper (F, gamma);
+    return result;
+  }
+}
+
+CanonicalForm 
+mapDown (const CanonicalForm& F, const ExtensionInfo& info, CFList& source, 
+         CFList& dest) 
+{
+  int k= info.getGFDegree();
+  Variable beta= info.getAlpha();
+  CanonicalForm primElem= info.getGamma();
+  CanonicalForm imPrimElem= info.getDelta();
+  if (k > 1)
+    return GFMapDown (F, k);
+  else if (k == 1)
+    return F;
+  else if (!k && beta == Variable (1))
+    return F;
+  else if (!k && beta != Variable (1)) 
+    return mapDown (F, primElem, imPrimElem, beta, source, dest);  
+}
+
+void appendTestMapDown (CFList& factors, const CanonicalForm& f, 
+                        const ExtensionInfo& info, CFList& source, CFList& dest) 
+{
+  int k= info.getGFDegree();
+  Variable beta= info.getBeta();
+  Variable alpha= info.getAlpha();
+  CanonicalForm delta= info.getDelta();
+  CanonicalForm gamma= info.getGamma();
+  CanonicalForm g= f;
+  int degMipoBeta;
+  if (!k && beta == Variable(1))
+    degMipoBeta= 1;
+  else if (!k && beta != Variable(1))
+    degMipoBeta= degree (getMipo (beta));
+  if (k > 1) 
+  {
+    if (!isInExtension (g, delta, k)) 
+    {
+      g= GFMapDown (g, k);
+      factors.append (g);
+    }
+  }
+  else if (k == 1) 
+  {           
+    if (!isInExtension (g, delta, k));
+      factors.append (g);
+  }
+  else if (!k && beta == Variable (1)) 
+  {
+    if (degree (g, alpha) < degMipoBeta)
+      factors.append (g);
+  }
+  else if (!k && beta != Variable (1)) 
+  {
+    if (!isInExtension (g, delta, k)) 
+    {
+      g= mapDown (g, gamma, delta, alpha, source, dest); 
+      factors.append (g);
+    }             
+  }
+  return;
+}
+
+void 
+appendMapDown (CFList& factors, const CanonicalForm& g, 
+               const ExtensionInfo& info, CFList& source, CFList& dest) 
+{
+  int k= info.getGFDegree();
+  Variable beta= info.getBeta();
+  Variable alpha= info.getAlpha();
+  CanonicalForm gamma= info.getGamma();
+  CanonicalForm delta= info.getDelta();
+  if (k > 1) 
+    factors.append (GFMapDown (g, k));
+  else if (k == 1)  
+    factors.append (g);
+  else if (!k && beta == Variable (1)) 
+    factors.append (g);
+  else if (!k && beta != Variable (1)) 
+    factors.append (mapDown (g, gamma, delta, alpha, source, dest)); 
+  return;
+}
+
+void normalize (CFList& factors) 
+{
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+    i.getItem() /= Lc(i.getItem());
+  return;
+}
+
+CFList subset (int index [], const int& s, const CFArray& elements, 
+               bool& noSubset) 
+{
+  int r= elements.size();
+  int i= 0;
+  CFList result;
+  noSubset= false;
+  if (index[s - 1] == 0) 
+  {
+    while (i < s) 
+    {
+      index[i]= i + 1;
+      result.append (elements[i]);
+      i++;
+    } 
+    return result;
+  }
+  int buf;
+  int k;
+  bool found= false;
+  if (index[s - 1] == r) 
+  {
+    if (index[0] == r - s + 1) 
+    {
+      noSubset= true;
+      return result;
+    }
+    else {
+      while (found == false) 
+      {
+        if (index[s - 2 - i] < r - i - 1)
+          found= true;
+        i++; 
+      }
+      buf= index[s - i - 1];
+      k= 0;
+      while (s - i - 1 + k < s) 
+      {
+        index[s - i - 1 + k]= buf + k + 1;
+        k++;
+      } 
+    }
+    for (int j= 0; j < s; j++) 
+      result.append (elements[index[j] - 1]);
+    return result;
+  }
+  else 
+  {
+    index[s - 1] += 1;
+    for (int j= 0; j < s; j++)
+      result.append (elements[index[j] - 1]);
+    return result;
+  }
+} 
+
+CFArray copy (const CFList& list) 
+{
+  CFArray array= CFArray (list.length());
+  int j= 0;
+  for (CFListIterator i= list; i.hasItem(); i++, j++) 
+    array[j]= i.getItem();
+  return array;
+}
+
+void indexUpdate (int index [], const int& subsetSize, const int& setSize, 
+                   bool& noSubset) 
+{
+  noSubset= false;
+  if (subsetSize > setSize) 
+  {
+    noSubset= true;
+    return;
+  }
+  int v [setSize];
+  for (int i= 0; i < setSize; i++) 
+    v[i]= index[i];
+  if (subsetSize == 1) 
+  {
+    v[0]= v[0] - 1;
+    if (v[0] >= setSize) 
+    {
+      noSubset= true;
+      return;
+    }
+  }
+  else 
+  {
+    if (v[subsetSize - 1] - v[0] + 1 == subsetSize && v[0] > 1) 
+    {
+      if (v[0] + subsetSize - 1 > setSize) 
+      { 
+        noSubset= true;
+        return;
+      }
+      v[0]= v[0] - 1;
+      for (int i= 1; i < subsetSize - 1; i++) 
+        v[i]= v[i - 1] + 1;
+      v[subsetSize - 1]= v[subsetSize - 2];
+    }
+    else 
+    {
+      if (v[0] + subsetSize - 1 > setSize) 
+      { 
+        noSubset= true;
+        return;
+      }
+      for (int i= 1; i < subsetSize - 1; i++)
+        v[i]= v[i - 1] + 1;
+      v[subsetSize - 1]= v[subsetSize - 2];
+    }
+  }
+  for (int i= 0; i < setSize; i++)
+    index[i]= v[i];
+  return;
+}
+
+int subsetDegree (const CFList& S) 
+{
+  int result= 0;
+  for (CFListIterator i= S; i.hasItem(); i++)
+    result += degree (i.getItem(), Variable (1));
+  return result;
+}
+
+CFFList multiplicity (CanonicalForm& F, const CFList& factors)
+{
+  if (F.inCoeffDomain())
+    return CFFList (CFFactor (F, 1));
+  CFFList result;
+  int multi= 0;
+  for (CFListIterator i= factors; i.hasItem(); i++)
+  {
+    while (fdivides (i.getItem(), F))
+    {
+      multi++;
+      F /= i.getItem();
+    }
+    result.append (CFFactor (i.getItem(), multi));
+    multi= 0;
+  }
+  return result;
+}
+
Index: factory/facFqBivarUtil.h
===================================================================
--- factory/facFqBivarUtil.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/facFqBivarUtil.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,162 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file facFqBivarUtil.h
+ * 
+ * This file provides utility functions for bivariate factorization       
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#ifndef FAC_FQ_BIVAR_UTIL_H
+#define FAC_FQ_BIVAR_UTIL_H
+
+#include <config.h>
+
+#include "cf_map.h"
+#include "ExtensionInfo.h"
+
+/// append @a factors2 on @a factors1
+void append (CFList& factors1,       ///< [in,out] a list of polys
+             const CFList& factors2  ///< [in] a list of polys
+            );
+
+/// decompress a list of polys @a factors using the map @a N
+void decompress (CFList& factors, ///< [in,out] a list of polys
+                 const CFMap& N   ///< [in] a map
+                );
+
+/// first swap Variables in @a factors1 if necessary, then append @a factors2 
+/// and @a factors3 on @a factors1 and finally decompress @a factors1 
+void appendSwapDecompress (CFList& factors1,       ///< [in,out] a list of polys 
+                           const CFList& factors2, ///< [in] a list of polys
+                           const CFList& factors3, ///< [in] a list of polys 
+                           const bool swap1,       ///< [in] indicates the need
+                                                   ///< to swap
+                           const bool swap2,       ///< [in] indicates the need
+                                                   ///< to swap
+                           const CFMap& N          ///< [in] a map
+                          );
+
+/// swap Variables in @a factors, then decompress @a factors
+void swapDecompress (CFList& factors, ///< [in,out] a list of polys
+                     const bool swap, ///< [in] indicates the need to swap
+                     const CFMap& N   ///< [in] a map
+                    );
+
+/// tests if F is not contained in a subfield defined by @a gamma (Fq case) or 
+/// @a k (GF case)
+///
+/// @return @a isInExtension returns true if @a F is not contained in a subfield 
+///         defined by @a gamma (Fq case) or @a k (GF case), false else
+/// @sa appendTestMapDown()
+bool isInExtension (const CanonicalForm& F,      ///< [in] a poly over 
+                                                 ///< F_p (alpha)=Fq or GF(p^l) 
+                    const CanonicalForm& gamma,  ///< [in] a primitive element
+                                                 ///< defining a subfield of
+                                                 ///< Fq if we are not over some
+                                                 ///< GF
+                    const int k                  ///< some int k such that k 
+                                                 ///< divides l if we are not
+                                                 ///< over some Fq
+                   );
+
+/// map a poly into a subfield of the current field, no testing is performed!
+///
+/// @return @a mapDown returns @a F mapped into a subfield of the current field
+/// @sa appendTestMapDown(), appendMapDown()
+CanonicalForm 
+mapDown (const CanonicalForm& F,    ///< [in] a poly
+         const ExtensionInfo& info, ///< [in] information about the sub- and
+                                    ///< current field
+         CFList& source,            ///< [in,out] in case we are over some 
+                                    ///< F_p (alpha) and want to map down into
+                                    ///< F_p (beta) source contains powers of
+                                    ///< the primitive element of F_p (alpha) 
+         CFList& dest               ///< [in,out] as source but contains
+                                    ///< the corresponding powers of the  
+                                    ///< primitive element of F_p (beta)
+        );
+
+/// test if @a g is in a subfield of the current field, if so map it down and 
+/// append it to @a factors
+///
+/// @sa mapDown(), isInExtension()
+void appendTestMapDown (CFList& factors,          ///< [in,out] a list of polys
+                        const CanonicalForm& f,   ///< [in] a poly
+                        const ExtensionInfo& info,///< [in] information about 
+                                                  ///< extension
+                        CFList& source,           ///< [in,out] @sa mapDown()
+                        CFList& dest              ///< [in,out] @sa mapDown()
+                       );
+
+/// map @a g down into a subfield of the current field and append it to @a 
+/// factors
+///
+/// @sa mapDown(), appendTestMapDown()  
+void 
+appendMapDown (CFList& factors,          ///< [in,out] a list of polys
+               const CanonicalForm& g,   ///< [in] a poly
+               const ExtensionInfo& info,///< [in] information about extension  
+               CFList& source,           ///< [in,out] @sa mapDown()
+               CFList& dest              ///< [in,out] @sa mapDown()
+              );
+
+/// normalize factors, i.e. make factors monic
+void normalize (CFList& factors ///< [in,out] a list of polys 
+               );
+
+/// extract a subset given by @a index of size @a s from @a elements, if there 
+/// is no subset we have not yet considered noSubset is set to true. @a index 
+/// encodes the next subset, e.g. if s= 3, elements= {a,b,c,d},
+/// index= {1, 2, 4, 0}, then subset= {a,c,d}. @a index is of size 
+/// @a elements.size(). 
+///
+/// @return @a subset returns a list of polys of length @a s if there is a 
+///         subset we have not yet considered.  
+CFList subset (int index [],            ///< [in,out] an array encoding the next
+                                        ///< subset
+               const int& s,            ///< [in] size of the subset
+               const CFArray& elements, ///< [in] an array of polys
+               bool& noSubset           ///< [in,out] if there is no subset we 
+                                        ///< have not yet considered @a noSubset
+                                        ///< is true
+              );
+
+/// write elements of @a list into an array  
+///
+/// @return an array of polys 
+CFArray copy (const CFList& list ///< [in] a list of polys
+             );
+
+/// update @a index 
+void indexUpdate (int index [],          ///< [in,out] an array encoding a 
+                                         ///< subset of size subsetSize  
+                  const int& subsetSize, ///< [in] size of the subset
+                  const int& setSize,    ///< [in] size of the given set
+                  bool& noSubset         ///< [in,out] if there is no subset we 
+                                         ///< have not yet considered @a 
+                                         ///< noSubset is true 
+                 );
+
+/// compute the sum of degrees in Variable(1) of elements in S
+///
+/// @return @a subsetDegree returns the sum of degrees in Variable(1) of 
+///         elements in S 
+int subsetDegree (const CFList& S ///< [in] a list of polys
+                 );
+
+/// determine multiplicity of the factors
+///
+/// @return a list of factors of F with their multiplicity 
+CFFList multiplicity (CanonicalForm& F,     ///< [in] a poly
+                      const CFList& factors ///< [in] a list of factors of F
+                     );
+
+#endif
+/* FAC_FQ_BIVAR_UTIL_H */
+
Index: factory/facFqFactorize.cc
===================================================================
--- factory/facFqFactorize.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/facFqFactorize.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,2087 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file facFqFactorize.cc
+ * 
+ * This file implements functions for factoring a multivariate polynomial over
+ * a finite field. 
+ * 
+ * ABSTRACT: "Efficient Multivariate Factorization over Finite Fields" by 
+ * L. Bernardin & M. Monagon.  
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#include <config.h>
+
+#include "assert.h"
+#include "debug.h"
+#include "timing.h"
+
+#include "facFqFactorizeUtil.h"
+#include "facFqFactorize.h"
+#include "cf_random.h"
+#include "facHensel.h"
+#include "cf_gcd_smallp.h"
+#include "DegreePattern.cc"
+#include "algext.h"
+
+#ifdef HAVE_NTL
+#include <NTL/ZZ_pEX.h>
+#include "NTLconvert.h"
+
+TIMING_DEFINE_PRINT(fac_bi_factorizer);
+TIMING_DEFINE_PRINT(fac_hensel_lift);
+TIMING_DEFINE_PRINT(fac_factor_recombination);
+
+static inline
+CanonicalForm 
+listGCD (const CFList& L);
+
+static inline
+CanonicalForm 
+myContent (const CanonicalForm& F)
+{
+  CanonicalForm G= swapvar (F, F.mvar(), Variable (1));
+  CFList L;
+  Variable alpha;
+  bool algExt= hasFirstAlgVar (G, alpha);
+  for (CFIterator i= G; i.hasTerms(); i++)
+    L.append (i.coeff());
+  if (L.length() == 2)
+  {
+    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+    if (GF)
+    {
+      return swapvar (GCD_GF (L.getFirst(), L.getLast()), F.mvar(), 
+                      Variable (1));
+    }
+    else if (!GF && algExt)
+    {
+      return swapvar (GCD_Fp_extension (L.getFirst(), L.getLast(), alpha), 
+                                       F.mvar(), Variable (1));
+    }
+    else
+      return swapvar (GCD_small_p (L.getFirst(), L.getLast()), F.mvar(), 
+                     Variable (1));
+  }
+  if (L.length() == 1)
+  {
+    return LC (F, Variable (1));
+  }
+  return swapvar (listGCD (L), F.mvar(), Variable (1));
+}
+
+static inline
+CanonicalForm 
+listGCD (const CFList& L)
+{
+  if (L.length() == 1)
+    return L.getFirst();
+  Variable alpha;
+  if (L.length() == 2)
+  {
+    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+    if (GF)
+    {
+      return GCD_GF (L.getFirst(), L.getLast());
+    }
+    else if (!GF && (hasFirstAlgVar (L.getFirst(), alpha) || 
+                     hasFirstAlgVar (L.getLast(), alpha)))
+    {
+      return GCD_Fp_extension (L.getFirst(), L.getLast(), alpha);
+    }
+    else
+      return GCD_small_p (L.getFirst(), L.getLast());
+  }
+  else
+  {
+    CFList lHi, lLo;
+    CanonicalForm resultHi, resultLo;
+    int length= L.length()/2;
+    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+    int j= 0;
+    for (CFListIterator i= L; j < length; i++, j++)
+      lHi.append (i.getItem());    
+    lLo= Difference (L, lHi);
+    resultHi= listGCD (lHi);
+    resultLo= listGCD (lLo);
+    if (resultHi.isOne() || resultLo.isOne())
+      return 1;
+    if (GF)
+    {
+      return GCD_GF (resultHi, resultLo);
+    }
+    else if (!GF && (hasFirstAlgVar (resultHi, alpha) || 
+                     hasFirstAlgVar (resultLo, alpha)))
+    {
+      return GCD_Fp_extension (resultHi, resultLo, alpha);
+    }
+    else
+      return GCD_small_p (resultHi, resultLo);
+  }
+}
+
+static inline
+CanonicalForm
+myContent (const CanonicalForm& F, const Variable& x)
+{
+  if (degree (F, x) <= 0)
+    return 1;
+  CanonicalForm G= F;
+  bool swap= false;
+  if (x != F.mvar())
+  {
+    swap= true;
+    G= swapvar (F, x, F.mvar());
+  }
+  CFList L;
+  Variable alpha;
+  bool algExt= hasFirstAlgVar (G, alpha);
+  for (CFIterator i= G; i.hasTerms(); i++)
+    L.append (i.coeff());
+  if (L.length() == 2)
+  {
+    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+    if (GF)
+    {
+      if (swap)
+        return swapvar(GCD_GF (L.getFirst(), L.getLast()), F.mvar(), x);
+      else
+        return GCD_GF (L.getFirst(), L.getLast());
+    }
+    else if (!GF && algExt)
+    {
+      if (swap)
+        return swapvar(GCD_Fp_extension (L.getFirst(), L.getLast(), alpha), 
+                                         F.mvar(), x);
+      else
+        return GCD_Fp_extension (L.getFirst(), L.getLast(), alpha);
+    }
+    else
+    {
+      if (swap)
+        return swapvar (GCD_small_p (L.getFirst(), L.getLast()), F.mvar(), 
+                        x);
+      else
+        return GCD_small_p (L.getFirst(), L.getLast());
+    }
+  }
+  if (L.length() == 1)
+  {
+    return LC (F, x);
+  }
+  if (swap)
+    return swapvar (listGCD (L), F.mvar(), x);
+  else
+    return listGCD (L);
+}
+
+static inline
+CanonicalForm
+myLcm (const CanonicalForm& F, const CanonicalForm& G)
+{
+  if ( F.isZero() || G.isZero() )
+    return 0;
+  else
+  {
+    Variable alpha;
+    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+    if (GF)
+      return (F/GCD_GF (F, G))*G;
+    else if (!GF && (hasFirstAlgVar (F, alpha) || 
+                     hasFirstAlgVar (G, alpha)))
+      return (F/GCD_Fp_extension (F, G, alpha))*G;
+    else
+      return (F/GCD_small_p (F, G))*G;
+  }
+}
+
+
+static inline 
+CanonicalForm myCompress (const CanonicalForm& F, CFMap& N)
+{
+  int n= F.level();
+  int * degsf= new int [n + 1];
+  int ** swap;
+  swap= new int* [n + 1];
+  for (int i= 0; i <= n; i++) 
+  {
+    degsf[i]= 0;
+    swap [i]= new int [2];
+    swap [i] [1]= 0;
+    swap [i] [2]= 0;
+  }
+  int i= 1;
+  n= 1;  
+  degsf= degrees (F, degsf);
+
+  CanonicalForm result= F;
+  while ( i <= F.level() ) 
+  {
+    while( degsf[i] == 0 ) i++;
+    swap[n][1]= i;
+    swap[n][2]= degsf[i]; 
+    if (i != n)
+      result= swapvar (result, Variable (n), Variable(i));
+    n++; i++;
+  }
+
+  int buf1, buf2;
+  n--;
+
+  for (i= 1; i < n; i++) 
+  {
+    for (int j= 1; j < n - i + 1; j++) 
+    {
+      if (swap[j][2] < swap[j + 1][2]) 
+      {
+        buf1= swap [j + 1] [1];
+        buf2= swap [j + 1] [2];
+        swap[j + 1] [1]= swap[j] [1];
+        swap[j + 1] [2]= swap[j] [2];
+        swap[j][1]= buf1;
+        swap[j][2]= buf2;
+        result= swapvar (result, Variable (j + 1), Variable (j));
+      }     
+    }
+  }
+
+  for (i= n; i > 0; i--) 
+  {
+    if (i != swap[i] [1])
+      N.newpair (Variable (i), Variable (swap[i] [1]));
+  }
+
+  for (i= 0; i <= n; i++)
+    delete [] swap[i];
+  delete [] swap;
+
+  delete [] degsf;
+
+  return result;
+}
+ 
+CFList 
+monicFactorRecombi (const CFList& factors,const CanonicalForm& F, const
+                    CanonicalForm& M, const DegreePattern& degs) 
+{
+  if (degs.getLength() == 1) 
+    return CFList(F); 
+
+  CFList T, S;
+  
+  T= factors;
+
+  int s= 1;
+  CFList result;
+  CanonicalForm buf= F;
+  CanonicalForm LCBuf= LC (buf, Variable (1));
+  CanonicalForm g, gg;
+  int v [T.length()];
+  for (int i= 0; i < T.length(); i++)
+    v[i]= 0;
+  bool noSubset= false;
+  CFArray TT;
+  int subsetDeg;
+  DegreePattern bufDegs1= degs, bufDegs2;
+  TT= copy (factors); 
+  while (T.length() >= 2*s) 
+  {
+    while (noSubset == false) 
+    {
+      if (T.length() == s) 
+      {
+        result.append (prodMod (T, M));
+        return result;
+      }
+      S= subset (v, s, TT, noSubset);
+      if (noSubset) break;
+      subsetDeg= subsetDegree (S);
+      if (!degs.find (subsetDeg)) 
+        continue;
+      else 
+      { 
+        g= prodMod0 (S, M);
+        gg= mod (g*LCBuf, M); //univariate polys
+        gg /= content (gg);
+        if (fdivides (LC (gg), LCBuf))
+        {
+          g= prodMod (S, M);
+          gg= mulMod2 (LCBuf, g, M);
+          gg /= content (gg, Variable (1));
+          if (fdivides (gg, buf))
+          {
+            result.append (g);
+            buf /= gg;
+            LCBuf= LC (buf, Variable(1));
+            T= Difference (T, S);
+            // compute new possible degree pattern
+            bufDegs2= DegreePattern (T);
+            bufDegs1.intersect (bufDegs2);
+            bufDegs1.refine (); 
+            if (T.length() < 2*s || T.length() == s || 
+                bufDegs1.getLength() == 1) 
+            {
+              result.append (prodMod (T, M)); 
+              return result; 
+            }
+            TT= copy (T);
+            indexUpdate (v, s, T.length(), noSubset);
+            if (noSubset) break; 
+          }
+        }
+      }
+    }
+    s++;
+    if (T.length() < 2*s || T.length() == s) 
+    {
+      result.append (prodMod (T, M));
+      return result;
+    }
+    int v [T.length()];
+    for (int i= 0; i < T.length(); i++)
+      v[i]= 0;
+    noSubset= false;
+  }
+  if (T.length() < 2*s) 
+    result.append (prodMod (T, M)); 
+
+  return result; 
+} 
+
+CFList 
+earlyMonicFactorDetect (CanonicalForm& F, CFList& factors, 
+                        int& adaptedLiftBound, DegreePattern& degs,
+                        bool& success, int deg, const int bound)
+{
+  DegreePattern bufDegs1= degs;
+  DegreePattern bufDegs2;
+  CFList result;
+  CFList T= factors;
+  CanonicalForm buf= F;
+  CanonicalForm LCBuf= LC (buf, Variable (1));
+  CanonicalForm g, gg;
+  CanonicalForm M= power (F.mvar(), deg);
+  int d= bound;
+  int e= 0;
+  adaptedLiftBound= 0;
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+  {
+    if (!bufDegs1.find (degree (i.getItem(), 1)))
+      continue;
+    else
+    {
+      gg= i.getItem() (0, 1);
+      gg *= LCBuf;
+      gg= mod (gg, M);
+      if (fdivides (LC (gg), LCBuf))
+      {
+        g= i.getItem();
+        gg= mulMod2 (g, LCBuf, M);
+        gg /= content (gg, Variable (1));
+        if (fdivides (gg, buf))
+        {
+          result.append (g);
+          CanonicalForm bufbuf= buf;
+          buf /= gg;
+          d -= degree (gg) + degree (LC (gg, 1));
+          e= tmax (e, degree (gg) + degree (LC (gg, 1)));
+          LCBuf= LC (buf, Variable (1)); 
+          T= Difference (T, CFList (i.getItem()));
+
+          // compute new possible degree pattern
+          bufDegs2= DegreePattern (T);
+          bufDegs1.intersect (bufDegs2);
+          bufDegs1.refine (); 
+          if (bufDegs1.getLength() <= 1)
+          {
+            result.append (prodMod (T, M)); 
+            break;
+          }
+        }
+      }
+    }
+  }
+  adaptedLiftBound= d;
+
+  if (adaptedLiftBound < deg) 
+  {
+    factors= T;
+    degs= bufDegs1;
+    if (adaptedLiftBound < degree (F) + 1)
+    {
+      if (d == 1)
+      {
+        if (e + 1 > deg)
+        {
+          success= false;
+          adaptedLiftBound= deg;
+        }
+        else
+        {
+          F= buf;
+          success= true;
+          if (e + 1 < degree (F) + 1)
+            adaptedLiftBound= deg;
+          else
+            adaptedLiftBound= e + 1;
+        }
+      }
+      else
+      {
+        success= true;
+        F= buf;
+        adaptedLiftBound= deg;
+      }
+    }
+    else
+    {
+      F= buf;
+      success= true;
+    }
+  }
+  if (bufDegs1.getLength() <= 1)
+  {
+    degs= bufDegs1;
+    if (!success)
+      adaptedLiftBound= deg;
+  }
+
+  return result;
+}
+
+CFList biFactorizer (const CanonicalForm& F, const Variable& alpha,
+                     CanonicalForm& bivarEval, int& liftBound) 
+{
+  CanonicalForm A= F;
+  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+
+  if (A.isUnivariate())  
+    return uniFactorizer (F, alpha, GF);
+  
+
+  CFMap N;
+  A= compress (A, N);
+  Variable y= A.mvar();
+  A /= Lc(A);
+
+  if (y.level() > 2) return CFList (F);
+  Variable x= Variable (1);
+
+  //remove content and factorize content
+  CanonicalForm contentAy= content (A, y);
+  CanonicalForm contentAx= content (A, x);
+  A= A/(contentAy*contentAx);
+  CFList contentAyFactors= uniFactorizer (contentAy, alpha, GF);
+
+  //trivial case
+  CFList factors;
+  if (A.inCoeffDomain()) 
+  {
+    append (factors, contentAyFactors);
+    decompress (factors, N);
+    bivarEval= N (y - bivarEval);
+    return factors;
+  }
+  else if (A.isUnivariate()) 
+  {
+    if (A.mvar() == x) 
+      factors= uniFactorizer (A, alpha, GF);
+    append (factors, contentAyFactors);
+    decompress (factors, N);
+    bivarEval= N (y - bivarEval);
+    return factors;
+  }
+  bool fail;
+  CanonicalForm Aeval, evaluation, bufAeval, bufEvaluation, buf;
+  CFList uniFactors, list, bufUniFactors;
+  DegreePattern degs;
+  DegreePattern bufDegs;
+ 
+  int factorNums= 5;
+  if (factorNums < (int) ceil (log (totaldegree (A)))) 
+    factorNums= (int) ceil (log (totaldegree (A)));
+  for (int i= 0; i < factorNums; i++) 
+  {
+    if (i == 0) 
+      Aeval= A (bivarEval, y);
+    else if (i > 0 && contentAx.inCoeffDomain())
+    {
+      Aeval= A; 
+      evaluation= evalPoint (A, Aeval, alpha, list, GF, fail);
+    }
+ 
+    if (fail && (i != 0)) 
+      break;
+
+    // univariate factorization
+    uniFactors= uniFactorizer (Aeval, alpha, GF);
+
+    if (uniFactors.length() == 1) 
+    { 
+      append (factors, contentAyFactors);
+      if (contentAyFactors.isEmpty())
+        factors.append (F/lc(F));
+      else 
+      {
+        A= A (y - bivarEval, y);
+        A= A/lc (A);
+        if (!LC(A, 1).isOne())
+        {
+          CanonicalForm s,t;
+          (void) extgcd (LC (A, 1), power (y, liftBound), s, t);
+          A *= s;
+          A= mod (A, power (y, liftBound)); 
+        }
+        factors.append (A);
+      }
+      decompress (factors, N);
+      bivarEval= N (bivarEval); 
+      return factors;
+    }
+
+    // degree analysis
+    degs= DegreePattern (uniFactors);
+
+    if (i == 0) 
+    {
+      bufAeval= Aeval;
+      bufEvaluation= bivarEval;
+      bufUniFactors= uniFactors;
+      bufDegs= degs;
+      if (!contentAx.inCoeffDomain())
+        break;
+    }
+    else 
+    {
+      bufDegs.intersect (degs); 
+      if (uniFactors.length() < bufUniFactors.length()) 
+      { 
+        bufUniFactors= uniFactors;
+        bufAeval= Aeval;
+        bufEvaluation= evaluation;
+      }
+    }
+
+    if (bufDegs.getLength() == 1)
+    {
+      append (factors, contentAyFactors);
+      if (contentAyFactors.isEmpty())
+        factors.append (F/lc(F));
+      else 
+      {
+        A= A (y - bivarEval, y);
+        A= A/lc (A);
+        if (!LC(A, 1).isOne())
+        {
+          CanonicalForm s,t;
+          (void) extgcd (LC (A, 1), power (y, liftBound), s, t);
+          A *= s;
+          A= mod (A, power (y, liftBound)); 
+        }
+        factors.append (A);
+      }
+      decompress (factors, N);
+      bivarEval= N (bivarEval); 
+      return factors;
+    }
+    list.append (evaluation);
+  }
+
+  bivarEval= y - bufEvaluation;
+  A= A (y + bufEvaluation, y);
+  bufUniFactors.insert (LC (A, x));
+
+  // Hensel lifting
+  CFList diophant;
+  CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1); 
+  CFArray Pi;
+  bool earlySuccess= false;
+  int newLiftBound= 0;
+  CFList earlyFactors;
+  int smallFactorDeg= 11; //tunable parameter
+  if (smallFactorDeg >= liftBound)
+    henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M);
+  else if (smallFactorDeg >= degree (A, y) + 1)
+  {
+    henselLift12 (A, bufUniFactors, degree (A, y) + 1, Pi, diophant, M);
+    earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
+                                           bufDegs, earlySuccess,
+                                           degree (A, y) + 1, liftBound);
+    if (bufDegs.getLength() > 1 && !earlySuccess)
+    {
+      if (newLiftBound > degree (A, y) + 1)
+      {
+        liftBound= newLiftBound;
+        bufUniFactors.insert (LC(A, x));
+        henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
+                            Pi, diophant, M);
+      }
+    }
+    else if (earlySuccess)
+      liftBound= newLiftBound;
+  }
+  else if (smallFactorDeg < degree (A, y) + 1)
+  {
+    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
+    earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
+                                           bufDegs, earlySuccess,
+                                           smallFactorDeg, liftBound);
+    if (bufDegs.getLength() > 1 && !earlySuccess)
+    {
+      bufUniFactors.insert (LC (A, x));
+      henselLiftResume12 (A, bufUniFactors, smallFactorDeg, degree (A, y) +
+                          1, Pi, diophant, M);
+      earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
+                                             bufDegs, earlySuccess, 
+                                             degree (A, y) + 1, liftBound);
+      if (bufDegs.getLength() > 1 && !earlySuccess)
+      {
+	if (newLiftBound > degree (A, y) + 1)
+	{
+          if (newLiftBound < newLiftBound)
+	    liftBound= newLiftBound;
+	  bufUniFactors.insert (LC(A, x));
+	  henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
+                              Pi, diophant, M);
+	}
+      }
+      else if (earlySuccess)
+        liftBound= newLiftBound;
+    }
+    else if (earlySuccess)
+      liftBound= newLiftBound;
+  }
+
+  if (newLiftBound > 0)
+    liftBound= newLiftBound;
+
+  CanonicalForm MODl= power (y, liftBound);
+
+  if (bufDegs.getLength() > 1)
+    factors= monicFactorRecombi (bufUniFactors, A, MODl, bufDegs);
+
+  if (earlySuccess) 
+    factors= Union (earlyFactors, factors);
+  else if (!earlySuccess && bufDegs.getLength() == 1)
+    factors= earlyFactors;
+
+  decompressAppend (factors, contentAyFactors, N);
+
+  bivarEval= N (bivarEval);
+
+  return factors;
+}
+
+CFList 
+extFactorRecombination (const CFList& factors, const CanonicalForm& F, 
+                        const CFList& M, const ExtensionInfo& info,  
+                        const CFList& evaluation)
+{ 
+  Variable alpha= info.getAlpha();
+  Variable beta= info.getBeta();
+  CanonicalForm gamma= info.getGamma();
+  CanonicalForm delta= info.getDelta();
+  int k= info.getGFDegree();
+  CFList source, dest;
+  if (factors.length() == 1) 
+  {
+    CanonicalForm buf= reverseShift (F, evaluation);
+    return CFList (mapDown (buf, info, source, dest)); 
+  } 
+  if (factors.length() < 1)
+    return CFList();
+
+  int degMipoBeta;
+  if (!k && beta == Variable(1))
+    degMipoBeta= 1;
+  else if (!k && beta != Variable(1))
+    degMipoBeta= degree (getMipo (beta));
+
+  CFList T, S;
+  T= factors;
+
+  int s= 1;
+  CFList result;
+  CanonicalForm buf;
+  if (beta != Variable (1))
+    buf= mapDown (F, gamma, delta, alpha, source, dest); 
+  else
+    buf= F;
+
+  CanonicalForm g, LCBuf= LC (buf, Variable (1)); 
+  CanonicalForm buf2;
+  int v [T.length()];
+  for (int i= 0; i < T.length(); i++)
+    v[i]= 0;
+  bool noSubset= false;
+  CFArray TT;
+  int subsetDeg;
+  TT= copy (factors); 
+  bool recombination= false;
+  while (T.length() >= 2*s) 
+  {
+    while (noSubset == false) 
+    {
+      if (T.length() == s) 
+      {
+        if (recombination)
+        {
+          T.insert (LCBuf);
+          g= prodMod (T, M);
+          T.removeFirst();
+          result.append (g/myContent (g));
+          g= reverseShift (g, evaluation);
+          g /= Lc (g);
+          appendTestMapDown (result, g, info, source, dest);
+          return result;
+        }
+        else
+        {
+          buf= reverseShift (buf, evaluation);
+          return CFList (buf);
+        }
+        return result;
+      }
+
+      S= subset (v, s, TT, noSubset);
+      if (noSubset) break;
+
+      S.insert (LCBuf);
+      g= prodMod (S, M);
+      S.removeFirst();
+      g /= myContent (g);
+      if (fdivides (g, buf)) 
+      {
+	buf2= reverseShift (g, evaluation); 
+        buf2 /= Lc (buf2);
+	if (!k && beta == Variable (1)) 
+	{
+	  if (degree (buf2, alpha) < degMipoBeta)
+	  {  
+	    appendTestMapDown (result, buf2, info, source, dest);
+	    buf /= g;
+	    LCBuf= LC (buf, Variable (1));
+	    recombination= true;
+	  }
+	}
+	else 
+	{
+	  if (!isInExtension (buf2, delta, k))
+	  {
+	    appendTestMapDown (result, buf2, info, source, dest);
+	    buf /= g;
+	    LCBuf= LC (buf, Variable (1));
+	    recombination= true;
+	  }
+	}
+	T= Difference (T, S);
+
+	if (T.length() < 2*s || T.length() == s) 
+	{
+	  buf= reverseShift (buf, evaluation);
+          buf /= Lc (buf);
+	  appendTestMapDown (result, buf, info, source, dest);
+	  return result; 
+	}
+	TT= copy (T);
+	indexUpdate (v, s, T.length(), noSubset);
+	if (noSubset) break;        
+      }
+    }
+    s++;
+    if (T.length() < 2*s || T.length() == s) 
+    {
+      buf= reverseShift (buf, evaluation);
+      appendTestMapDown (result, buf, info, source, dest);
+      return result;
+    }
+    int v [T.length()];
+    for (int i= 0; i < T.length(); i++)
+      v[i]= 0;
+    noSubset= false;
+  }
+  if (T.length() < 2*s)
+  {
+    buf= reverseShift (F, evaluation); 
+    appendMapDown (result, buf, info, source, dest); 
+  }
+
+  return result; 
+} 
+
+CFList 
+factorRecombination (const CanonicalForm& F, const CFList& factors, 
+                     const CFList& M)
+{
+  if (factors.length() == 1) 
+    return CFList(F);
+  if (factors.length() < 1)
+    return CFList();
+
+  CFList T, S;
+
+  T= factors;
+
+  int s= 1;
+  CFList result;
+  CanonicalForm LCBuf= LC (F, Variable (1));
+  CanonicalForm g, buf= F;
+  int v [T.length()];
+  for (int i= 0; i < T.length(); i++)
+    v[i]= 0;
+  bool noSubset= false;
+  CFArray TT;
+  int subsetDeg;
+  TT= copy (factors); 
+  Variable y= F.level() - 1;
+  bool recombination= false;
+  CanonicalForm h;
+  while (T.length() >= 2*s) 
+  {
+    while (noSubset == false) 
+    {
+      if (T.length() == s) 
+      {
+        if (recombination)
+        {
+          T.insert (LC (buf));
+          g= prodMod (T, M);
+          result.append (g/myContent (g));
+          return result;
+        }
+        else
+          return CFList (F);
+      }
+      S= subset (v, s, TT, noSubset);
+      if (noSubset) break;
+      S.insert (LCBuf);
+      g= prodMod (S, M);
+      S.removeFirst();
+      g /= myContent (g);
+      if (fdivides (g, buf))
+      {
+        recombination= true;
+        result.append (g);
+        buf /= g;
+        LCBuf= LC (buf, Variable(1));
+        T= Difference (T, S);
+	if (T.length() < 2*s || T.length() == s) 
+        {
+	  result.append (buf);
+	  return result; 
+	}
+	TT= copy (T);
+	indexUpdate (v, s, T.length(), noSubset);
+	if (noSubset) break;        
+      }
+    }
+    s++;
+    if (T.length() < 2*s || T.length() == s) 
+    {
+      result.append (buf);
+      return result;
+    }
+    int v [T.length()];
+    for (int i= 0; i < T.length(); i++)
+      v[i]= 0;
+    noSubset= false;
+  }
+  if (T.length() < 2*s) 
+    result.append (F);
+
+  return result; 
+}
+
+int
+liftBoundAdaption (const CanonicalForm& F, const CFList& factors, bool&
+                   success, const int deg, const CFList& MOD, const int bound)
+{
+  int adaptedLiftBound= 0;
+  CanonicalForm buf= F;
+  Variable y= F.mvar();
+  CanonicalForm LCBuf= LC (buf, Variable (1));
+  CanonicalForm g;
+  CFList M= MOD;
+  M.append (power (y, deg));
+  int d= bound;
+  int e= 0;
+  int nBuf;
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+  {
+    g= mulMod (i.getItem(), LCBuf, M); 
+    g /= myContent (g);
+    if (fdivides (g, buf))
+    {
+      nBuf= degree (g, y) + degree (LC (g, 1), y);
+      d -= nBuf;
+      e= tmax (e, nBuf);
+      buf /= g;
+      LCBuf= LC (buf, Variable (1)); 
+    }
+  }
+  adaptedLiftBound= d;
+
+  if (adaptedLiftBound < deg) 
+  {  
+    if (adaptedLiftBound < degree (F) + 1)
+    {
+      if (d == 1)
+      {
+        if (e + 1 > deg)
+        {
+          adaptedLiftBound= deg;
+          success= false;
+        }
+        else 
+        {
+          success= true;
+          if (e + 1 < degree (F) + 1)
+            adaptedLiftBound= deg;
+          else
+            adaptedLiftBound= e + 1;
+        }
+      }
+      else
+      {
+        success= true;
+        adaptedLiftBound= deg; 
+      }
+    }
+    else
+    {
+      success= true;
+    }
+  }
+  return adaptedLiftBound;
+}
+
+int 
+extLiftBoundAdaption (const CanonicalForm& F, const CFList& factors, bool& 
+                      success, const ExtensionInfo& info, const CFList& eval, 
+                      const int deg, const CFList& MOD, const int bound)
+{
+  Variable alpha= info.getAlpha();
+  Variable beta= info.getBeta();
+  CanonicalForm gamma= info.getGamma();
+  CanonicalForm delta= info.getDelta();
+  int k= info.getGFDegree();
+  int adaptedLiftBound= 0;
+  CanonicalForm buf= F;
+  Variable y= F.mvar();
+  CanonicalForm LCBuf= LC (buf, Variable (1));
+  CanonicalForm g, gg;
+  CFList M= MOD;
+  M.append (power (y, deg));
+  adaptedLiftBound= 0;
+  int d= bound;
+  int e= 0; 
+  int nBuf;
+  int degMipoBeta;
+  if (!k && beta == Variable(1))
+    degMipoBeta= 1;
+  else if (!k && beta != Variable(1))
+    degMipoBeta= degree (getMipo (beta));
+
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+  {
+    g= mulMod (i.getItem(), LCBuf, M);
+    g /= myContent (g);
+    if (fdivides (g, buf))
+    {
+      gg= reverseShift (g, eval);
+      gg /= Lc (gg);
+      if (!k && beta == Variable (1)) 
+      {
+        if (degree (gg, alpha) < degMipoBeta)
+        {  
+          buf /= g;
+          nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
+          d -= nBuf;
+          e= tmax (e, nBuf);
+          LCBuf= LC (buf, Variable (1));            
+        }
+      }
+      else 
+      {
+        if (!isInExtension (gg, delta, k))
+        {
+          buf /= g;
+          nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
+          d -= nBuf;
+          e= tmax (e, nBuf);
+          LCBuf= LC (buf, Variable (1));
+        }
+      }
+    } 
+  }
+  adaptedLiftBound= d;
+
+  if (adaptedLiftBound < deg)
+  {
+    if (adaptedLiftBound < degree (F) + 1)
+    {
+      if (d == 1)
+      {
+        if (e + 1 > deg)
+        {
+          adaptedLiftBound= deg;
+          success= false;
+        }
+        else 
+        {
+          success= true;
+          if (e + 1 < degree (F) + 1)
+            adaptedLiftBound= deg;
+          else
+            adaptedLiftBound= e + 1;
+        }
+      }
+      else 
+      {
+        success= true;
+        adaptedLiftBound= deg; 
+      }
+    }
+    else
+    {
+      success= true;
+    }
+  }
+
+  return adaptedLiftBound;
+}
+
+CFList 
+earlyFactorDetect (CanonicalForm& F, CFList& factors, int& adaptedLiftBound,
+                   bool& success, const int deg, const CFList& MOD, 
+                   const int bound)
+{
+  CFList result;
+  CFList T= factors;
+  CanonicalForm buf= F;
+  Variable y= F.mvar();
+  CanonicalForm LCBuf= LC (buf, Variable (1));
+  CanonicalForm g;
+  CFList M= MOD;
+  M.append (power (y, deg));
+  adaptedLiftBound= 0;
+  int d= bound;
+  int e= 0;
+  int nBuf;
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+  {
+    g= mulMod (i.getItem(), LCBuf, M);
+    g /= myContent (g);
+    if (fdivides (g, buf))
+    {
+      result.append (g);
+      nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
+      d -= nBuf;
+      e= tmax (e, nBuf);
+      buf /= g;
+      LCBuf= LC (buf, Variable (1)); 
+      T= Difference (T, CFList (i.getItem()));
+    }
+  }
+  adaptedLiftBound= d;
+
+  if (adaptedLiftBound < deg) 
+  {
+    if (adaptedLiftBound < degree (F) + 1)
+    {
+      if (d == 1)
+        adaptedLiftBound= tmin (e + 1, deg);
+      else
+        adaptedLiftBound= deg; 
+    }
+    factors= T;
+    F= buf;
+    success= true;
+  }
+  return result;
+}
+
+CFList 
+extEarlyFactorDetect (CanonicalForm& F, CFList& factors, int& adaptedLiftBound,
+                      bool& success, const ExtensionInfo& info, const CFList& 
+                      eval, const int deg, const CFList& MOD, const int bound)
+{
+  Variable alpha= info.getAlpha();
+  Variable beta= info.getBeta();
+  CanonicalForm gamma= info.getGamma();
+  CanonicalForm delta= info.getDelta();
+  int k= info.getGFDegree();
+  CFList result;
+  CFList T= factors;
+  CanonicalForm buf= F;
+  Variable y= F.mvar();
+  CanonicalForm LCBuf= LC (buf, Variable (1));
+  CanonicalForm g, gg;
+  CFList M= MOD;
+  M.append (power (y, deg));
+  adaptedLiftBound= 0;
+  int d= bound;
+  int e= 0;
+  int nBuf;
+  CFList source, dest;
+
+  int degMipoBeta;
+  if (!k && beta == Variable(1))
+    degMipoBeta= 1;
+  else if (!k && beta != Variable(1))
+    degMipoBeta= degree (getMipo (beta));
+
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+  {
+    g= mulMod (i.getItem(), LCBuf, M);
+    g /= myContent (g);
+    if (fdivides (g, buf))
+    {
+      gg= reverseShift (g, eval);
+      gg /= Lc (gg);
+          if (!k && beta == Variable (1)) 
+          {
+            if (degree (gg, alpha) < degMipoBeta)
+            {  
+              appendTestMapDown (result, gg, info, source, dest);
+              buf /= g;
+              nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
+              d -= nBuf;
+              e= tmax (e, nBuf);
+              LCBuf= LC (buf, Variable (1));            
+            }
+          }
+          else 
+          {
+            if (!isInExtension (gg, delta, k))
+            {
+              appendTestMapDown (result, gg, info, source, dest);
+              buf /= g;
+              nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
+              d -= nBuf;
+              e= tmax (e, nBuf);
+              LCBuf= LC (buf, Variable (1));
+            }
+          } 
+      T= Difference (T, CFList (i.getItem()));
+    }
+  }
+  adaptedLiftBound= d;
+
+  if (adaptedLiftBound < deg)
+  {
+    if (adaptedLiftBound < degree (F) + 1)
+    {
+      if (d == 1)
+        adaptedLiftBound= tmin (e + 1, deg);
+      else
+        adaptedLiftBound= deg; 
+    }
+    success= true;
+    factors= T;
+    F= buf;
+  }
+  return result;
+}
+
+CFList 
+evalPoints (const CanonicalForm& F, CFList & eval, const Variable& alpha,
+            CFList& list, const bool& GF, bool& fail) 
+{  
+  int k= F.level() - 1;
+  Variable x= Variable (1);
+  CFList result;
+  FFRandom genFF;
+  GFRandom genGF;
+  int p= getCharacteristic ();
+  double bound; 
+  if (alpha != Variable (1)) 
+  {
+    bound= pow ((double) p, (double) degree (getMipo(alpha)));
+    bound= pow ((double) bound, (double) k);
+  }
+  else if (GF) 
+  {
+    bound= pow ((double) p, (double) getGFDegree());
+    bound= pow ((double) bound, (double) k);
+  }
+  else
+    bound= pow ((double) p, (double) k);
+
+  CanonicalForm random;
+  CanonicalForm deriv_x, gcd_deriv;
+  do 
+  {
+    random= 0;
+    // possible overflow if list.length() does not fit into a int
+    if (list.length() >= bound) 
+    { 
+      fail= true;
+      break;
+    }
+    for (int i= 0; i < k; i++) 
+    { 
+      if (GF) 
+      {
+        result.append (genGF.generate());
+        random += result.getLast()*power (x, i);
+      }
+      else if (alpha != Variable(1)) 
+      {    
+        AlgExtRandomF genAlgExt (alpha);
+        result.append (genAlgExt.generate());
+        random += result.getLast()*power (x, i); 
+      }
+      else 
+      { 
+        result.append (genFF.generate());
+        random += result.getLast()*power (x, i);
+      } 
+    }
+    if (find (list, random)) 
+    {
+      result= CFList();
+      continue;
+    }
+    int l= F.level();
+    eval.insert (F);
+    for (CFListIterator i= result; i.hasItem(); i++, l--) 
+      eval.insert (eval.getFirst()(i.getItem(), l));
+
+    if (degree (eval.getFirst()) != degree (F, 1)) 
+    {
+      if (!find (list, random))
+        list.append (random);
+      result= CFList();
+      eval= CFList();
+      continue;
+    }     
+   
+    deriv_x= deriv (eval.getFirst(), x);
+    gcd_deriv= gcd (eval.getFirst(), deriv_x);
+    if (degree (gcd_deriv) > 0) 
+    { 
+      if (!find (list, random))
+        list.append (random);
+      result= CFList();
+      eval= CFList();
+      continue;
+    }
+    CFListIterator i= eval;
+    i++;
+    CanonicalForm contentx= content (i.getItem(), x);
+    if (degree (contentx) > 0)
+    {
+      if (!find (list, random))
+        list.append (random);
+      result= CFList();
+      eval= CFList();
+      continue;
+    }
+
+    if (list.length() >= bound) 
+    {
+      fail= true;
+      break;
+    }
+  } while (find (list, random));
+
+  if (!eval.isEmpty()) 
+    eval.removeFirst();
+
+  return result;
+}
+
+static inline  
+int newMainVariableSearch (CanonicalForm& A, CFList& Aeval, CFList&
+                           evaluation, const Variable& alpha, const int lev) 
+{ 
+  Variable x= Variable (1);
+  CanonicalForm derivI, buf;
+  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+  int swapLevel= 0;
+  CFList list;
+  bool fail= false;
+  buf= A;
+  Aeval= CFList();
+  evaluation= CFList();
+  for (int i= lev; i <= A.level(); i++) 
+  {
+    derivI= deriv (buf, Variable (i));
+    if (!derivI.isZero())
+    {
+      buf= swapvar (buf, x, Variable (i)); 
+      Aeval= CFList();
+      evaluation= CFList();
+      fail= false;
+      evaluation= evalPoints (buf, Aeval, alpha, list, GF, fail);
+      if (!fail) 
+      {
+	A= buf;
+	swapLevel= i;
+	break;
+      }
+      else 
+        buf= A;
+    }
+  }
+  return swapLevel;
+}
+
+static inline 
+CanonicalForm lcmContent (const CanonicalForm& A, CFList& contentAi) 
+{
+  int i= A.level();
+  contentAi.append (myContent (A, i));
+  contentAi.append (myContent (A, i - 1));
+  CanonicalForm result= myLcm (contentAi.getFirst(), contentAi.getLast()); 
+  for (i= i - 2; i > 0; i--) 
+  {
+    contentAi.append (content (A, i));
+    result= myLcm (result, contentAi.getLast());
+  }
+  return result;
+} 
+ 
+CFList
+henselLiftAndEarly (CanonicalForm& A, CFList& MOD, int*& liftBounds, bool&
+                    earlySuccess, CFList& earlyFactors, const CFList& Aeval,
+                    const CFList& biFactors, const CFList& evaluation, 
+                    const ExtensionInfo& info)
+{
+  bool extension= info.isInExtension();
+  CFList bufFactors= biFactors; 
+  bufFactors.insert (LC (Aeval.getFirst(), 1));
+
+  sortList (bufFactors, Variable (1));
+
+  CFList diophant;
+  CFArray Pi;
+  int smallFactorDeg= 11; //tunable parameter
+  CFList result; 
+  int newLiftBound= 0;
+  int adaptedLiftBound= 0;
+  int liftBound= liftBounds[1];
+
+  earlySuccess= false;
+  bool earlyReconst= false;
+  CFList earlyReconstFactors;
+  CFListIterator j= Aeval;
+  j++;
+  CanonicalForm buf= j.getItem();
+  CFMatrix Mat= CFMatrix (liftBound, bufFactors.length() - 1);
+  MOD= CFList (power (Variable (2), liftBounds[0]));
+  if (smallFactorDeg >= liftBound)
+  {
+    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
+  }
+  else if (smallFactorDeg >= degree (buf) + 1)
+  {
+    liftBounds[1]= degree (buf) + 1;
+    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
+    if (Aeval.length() == 2)
+    {
+      if (!extension)
+        earlyFactors= earlyFactorDetect 
+                       (buf, result, adaptedLiftBound, earlySuccess,
+                        degree (buf) + 1, MOD, liftBound);
+      else
+        earlyFactors= extEarlyFactorDetect 
+                       (buf, result, adaptedLiftBound, earlySuccess,
+                        info, evaluation, degree       
+                        (buf) + 1, MOD, liftBound);
+    }
+    else
+    {
+      if (!extension)
+        adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
+                                             degree (buf) + 1, MOD, liftBound);
+      else
+        adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, info,
+                                                evaluation, degree (buf) + 1,
+                                                MOD, liftBound);       
+    }
+    if (!earlySuccess)
+    {
+      result.insert (LC (buf, 1));
+      liftBounds[1]= adaptedLiftBound;
+      liftBound= adaptedLiftBound;
+      henselLiftResume (buf, result, degree (buf) + 1, liftBound, 
+                        Pi, diophant, Mat, MOD);
+    }  
+    else
+      liftBounds[1]= adaptedLiftBound;
+  }
+  else if (smallFactorDeg < degree (buf) + 1)
+  {
+    liftBounds[1]= smallFactorDeg;
+    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
+    if (Aeval.length() == 2)
+    {
+      if (!extension)
+        earlyFactors= earlyFactorDetect (buf, result, adaptedLiftBound,  
+                                         earlySuccess, smallFactorDeg, MOD,
+                                         liftBound);
+      else 
+        earlyFactors= extEarlyFactorDetect (buf, result, adaptedLiftBound,
+                                            earlySuccess, info, evaluation,
+                                            smallFactorDeg, MOD, liftBound);
+    }
+    else
+    {
+      if (!extension)
+        adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
+                                             smallFactorDeg, MOD, liftBound);
+      else
+        adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, info,
+                                                evaluation, smallFactorDeg, MOD, 
+                                                liftBound);
+    }
+
+    if (!earlySuccess)
+    {
+      result.insert (LC (buf, 1));
+      henselLiftResume (buf, result, smallFactorDeg, degree (buf) + 1, 
+                        Pi, diophant, Mat, MOD);
+      if (Aeval.length() == 2)
+      {
+         if (!extension)
+           earlyFactors= earlyFactorDetect (buf, result, adaptedLiftBound,
+                                            earlySuccess, degree (buf) + 1,
+                                            MOD, liftBound);
+         else
+           earlyFactors= extEarlyFactorDetect (buf, result, adaptedLiftBound,
+                                               earlySuccess, info, evaluation,
+                                               degree (buf) + 1, MOD, 
+                                               liftBound);
+      }
+      else 
+      {
+        if (!extension)
+          adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
+                                               degree (buf) + 1, MOD,liftBound); 
+        else
+          adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, 
+                                                  info, evaluation,
+                                                  degree (buf) + 1, MOD,
+                                                  liftBound);
+      }
+      if (!earlySuccess)
+      {
+        result.insert (LC (buf, 1));
+        liftBounds[1]= adaptedLiftBound;
+        liftBound= adaptedLiftBound;
+        henselLiftResume (buf, result, degree (buf) + 1, liftBound, 
+                          Pi, diophant, Mat, MOD);
+      }
+      else
+        liftBounds[1]= adaptedLiftBound;
+    }
+    else
+      liftBounds[1]= adaptedLiftBound;
+  }
+
+  MOD.append (power (Variable (3), liftBounds[1]));
+
+  if (Aeval.length() > 2)
+  {
+    CFListIterator j= Aeval;
+    j++;
+    CFList bufEval;
+    bufEval.append (j.getItem());
+    j++;
+    int liftBoundsLength= Aeval.getLast().level() - 1;
+    for (int i= 2; i <= liftBoundsLength && j.hasItem(); i++, j++)
+    { 
+      earlySuccess= false;
+      result.insert (LC (bufEval.getFirst(), 1));
+      bufEval.append (j.getItem());
+      liftBound= liftBounds[i];
+      Mat= CFMatrix (liftBounds[i], result.length() - 1);
+
+      buf= j.getItem();
+      if (smallFactorDeg >= liftBound)
+        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
+			    liftBounds[i -  1], liftBounds[i]);
+      else if (smallFactorDeg >= degree (buf) + 1)
+      {
+	result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
+			    liftBounds[i -  1], degree (buf) + 1);
+
+	if (Aeval.length() == i + 1)
+	{
+	  if (!extension)
+	    earlyFactors= earlyFactorDetect 
+                           (buf, result, adaptedLiftBound, earlySuccess,
+                            degree (buf) + 1, MOD, liftBound);
+	  else
+	    earlyFactors= extEarlyFactorDetect 
+                           (buf, result, adaptedLiftBound, earlySuccess,
+                            info, evaluation, degree (buf) + 1, MOD, liftBound);
+	}
+	else
+	{
+	  if (!extension)
+	    adaptedLiftBound= liftBoundAdaption 
+                                (buf, result, earlySuccess, degree (buf)
+                                 + 1,  MOD, liftBound);
+	  else
+	    adaptedLiftBound= extLiftBoundAdaption 
+                                (buf, result, earlySuccess, info, evaluation, 
+                                 degree (buf) + 1, MOD, liftBound);	  
+	}
+
+	if (!earlySuccess)
+	{
+	  result.insert (LC (buf, 1));
+          liftBounds[i]= adaptedLiftBound;
+          liftBound= adaptedLiftBound;
+	  henselLiftResume (buf, result, degree (buf) + 1, liftBound, 
+			    Pi, diophant, Mat, MOD);
+	}  
+	else
+	{
+	  liftBounds[i]= adaptedLiftBound;
+	}
+      }
+      else if (smallFactorDeg < degree (buf) + 1)
+      {
+	result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
+			    liftBounds[i -  1], smallFactorDeg);
+
+	if (Aeval.length() == i + 1)
+	{
+	  if (!extension)
+	    earlyFactors= earlyFactorDetect 
+                           (buf, result, adaptedLiftBound, earlySuccess,
+                            smallFactorDeg, MOD, liftBound);
+	  else
+	    earlyFactors= extEarlyFactorDetect 
+                           (buf, result, adaptedLiftBound, earlySuccess,
+                            info, evaluation, smallFactorDeg, MOD, liftBound);
+	}
+	else
+	{
+	  if (!extension)
+	    adaptedLiftBound= liftBoundAdaption 
+                                (buf, result, earlySuccess,
+                                 smallFactorDeg, MOD, liftBound);
+	  else
+	    adaptedLiftBound= extLiftBoundAdaption 
+                                (buf, result, earlySuccess, info, evaluation,
+                                 smallFactorDeg, MOD, liftBound);	  
+	}
+
+	if (!earlySuccess)
+	{
+	  result.insert (LC (buf, 1));
+	  henselLiftResume (buf, result, smallFactorDeg, 
+                            degree (buf) + 1, Pi, diophant, Mat, MOD);
+	  if (Aeval.length() == i + 1)
+	  {
+	    if (!extension)
+	      earlyFactors= earlyFactorDetect 
+                             (buf, result, adaptedLiftBound, earlySuccess,
+                              degree (buf) +  1,  MOD, liftBound);
+	    else
+	      earlyFactors= extEarlyFactorDetect 
+                             (buf, result, adaptedLiftBound, earlySuccess,
+                              info, evaluation, degree (buf) + 1, MOD, 
+                              liftBound);
+	  }
+	  else
+	  {
+	    if (!extension)
+	      adaptedLiftBound= liftBoundAdaption 
+                                  (buf, result, earlySuccess, degree
+                                   (buf) +  1,  MOD, liftBound);
+	    else
+	      adaptedLiftBound= extLiftBoundAdaption 
+                                  (buf, result, earlySuccess, info, evaluation, 
+                                   degree (buf) + 1,  MOD, liftBound);    
+	  }
+
+	  if (!earlySuccess)
+	  {
+	    result.insert (LC (buf, 1));
+            liftBounds[i]= adaptedLiftBound;
+            liftBound= adaptedLiftBound;
+	    henselLiftResume (buf, result, degree (buf) + 1, liftBound, 
+			      Pi, diophant, Mat, MOD);
+	  }
+	  else
+	    liftBounds[i]= adaptedLiftBound;
+	}
+	else
+	  liftBounds[i]= adaptedLiftBound;
+      }
+      MOD.append (power (Variable (i + 2), liftBounds[i]));
+      bufEval.removeFirst();
+    }
+    bufFactors= result;
+  }
+  else
+    bufFactors= result;
+
+  if (earlySuccess)
+    A= buf;
+  return result;
+}
+
+CFList 
+extFactorize (const CanonicalForm& F, const ExtensionInfo& info); 
+
+CFList 
+multiFactorize (const CanonicalForm& F, const ExtensionInfo& info) 
+{
+
+  if (F.inCoeffDomain())
+    return CFList (F);
+
+  // compress and find main Variable
+  CFMap N;
+  CanonicalForm A= myCompress (F, N);
+
+  A /= Lc (A); // make monic
+
+  Variable alpha= info.getAlpha();
+  Variable beta= info.getBeta();
+  CanonicalForm gamma= info.getGamma();
+  CanonicalForm delta= info.getDelta();
+  int k= info.getGFDegree();
+  bool extension= info.isInExtension();
+  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+  //univariate case
+  if (F.isUnivariate()) 
+  {
+    if (extension == false)
+      return uniFactorizer (F, alpha, GF);
+    else 
+    {  
+      CFList source, dest;
+      A= mapDown (F, info, source, dest);
+      return uniFactorizer (A, beta, GF);
+    }
+  }
+
+  //bivariate case
+  if (A.level() == 2) 
+  {
+    CFList buf= biFactorize (F, info); 
+    return buf;
+  }
+ 
+  Variable x= Variable (1);
+  Variable y= Variable (2);
+ 
+  // remove content
+  CFList contentAi;
+  CanonicalForm lcmCont= lcmContent (A, contentAi);
+  A /= lcmCont; 
+
+  // trivial after content removal
+  CFList contentAFactors;
+  if (A.inCoeffDomain()) 
+  { 
+    for (CFListIterator i= contentAi; i.hasItem(); i++) 
+    {
+      if (i.getItem().inCoeffDomain())
+        continue;
+      else 
+      {
+        lcmCont /= i.getItem();
+        contentAFactors= 
+        Union (multiFactorize (lcmCont, info), 
+               multiFactorize (i.getItem(), info));
+        break;
+      }
+    }
+    decompress (contentAFactors, N);
+    normalize (contentAFactors);
+    return contentAFactors;
+  }
+
+  // factorize content
+  contentAFactors= multiFactorize (lcmCont, info); 
+
+  // univariate after content removal
+  CFList factors;
+  if (A.isUnivariate ()) 
+  {
+    factors= uniFactorizer (A, alpha, GF);
+    append (factors, contentAFactors);
+    decompress (factors, N);
+    return factors;
+  }
+
+  // check main variable
+  int swapLevel= 0;
+  CanonicalForm derivZ;
+  CanonicalForm gcdDerivZ;
+  CanonicalForm bufA= A;
+  Variable z;
+  for (int i= 1; i <= A.level(); i++)
+  {  
+    z= Variable (i);
+    derivZ= deriv (bufA, z);
+    if (derivZ.isZero())
+    {
+      if (i == 1)
+        swapLevel= 1;
+      else
+        continue;
+    }
+    else
+    {
+      if (swapLevel == 1)
+      {
+        swapLevel= i;
+        A= swapvar (A, x, z);
+      }
+      if (GF)
+        gcdDerivZ= GCD_GF (bufA, derivZ); 
+      else if (alpha == Variable (1))
+        gcdDerivZ= GCD_small_p (bufA, derivZ);
+      else 
+        gcdDerivZ= GCD_Fp_extension (bufA, derivZ, alpha);
+      if (degree (gcdDerivZ) > 0 && !derivZ.isZero()) 
+      {
+        CanonicalForm g= bufA/gcdDerivZ;
+        CFList factorsG= 
+        Union (multiFactorize (g, info),
+               multiFactorize (gcdDerivZ, info));
+        appendSwapDecompress (factorsG, contentAFactors, N, swapLevel, x);
+        normalize (factorsG);
+        return factorsG;
+      }
+    }
+  }
+
+
+  CFList Aeval, list, evaluation, bufEvaluation, bufAeval;
+  bool fail= false;
+  int swapLevel2= 0;
+  int level;
+  int factorNums= 3;
+  CanonicalForm bivarEval;
+  CFList biFactors, bufBiFactors;
+  CanonicalForm evalPoly;
+  int lift, bufLift;
+  if (factorNums < (int) ceil (log (totaldegree (A)))) 
+    factorNums= (int) ceil (log (totaldegree (A)));
+  // several bivariate factorizations 
+  for (int i= 0; i < factorNums; i++) 
+  {
+    bufA= A;
+    bufAeval= CFList();
+    bufEvaluation= evalPoints (bufA, bufAeval, alpha, list, GF, fail);
+    evalPoly= 0;  
+
+    if (fail && (i == 0)) 
+    {
+      if (!swapLevel)
+	level= 2;
+      else 
+	level= swapLevel + 1;
+
+      swapLevel2= newMainVariableSearch (A, Aeval, evaluation, alpha, level);
+
+      if (!swapLevel2) // need to pass to an extension
+      { 
+	factors= extFactorize (A, info);
+	appendSwapDecompress (factors, contentAFactors, N, swapLevel, x);
+	normalize (factors);
+	return factors;
+      }
+      else
+      {
+        fail= false;
+        bufAeval= Aeval;
+        bufA= A;
+        bufEvaluation= evaluation; 
+      }
+    }
+    else if (fail && (i > 0))
+      break;
+
+    bivarEval= bufEvaluation.getLast();
+    bufEvaluation.removeLast(); 
+
+    bufLift= degree (A, y) + 1 + degree (LC(A, x), y);
+
+    TIMING_START (fac_bi_factorizer);
+    bufBiFactors= biFactorizer (bufAeval.getFirst(), alpha, bivarEval, bufLift); 
+    TIMING_END_AND_PRINT (fac_bi_factorizer, 
+			  "time for bivariate factorization: ");
+
+    if (bufBiFactors.length() == 1) 
+    {
+      if (extension) 
+      {
+	CFList source, dest;
+	A= mapDown (A, info, source, dest);
+      }
+      factors.append (A);
+      appendSwapDecompress (factors, contentAFactors, N, swapLevel,
+			    swapLevel2, x);
+      normalize (factors);
+      return factors;
+    }
+
+    bufEvaluation.append (-bivarEval[0]);
+    if (i == 0)
+    {
+      Aeval= bufAeval;
+      evaluation= bufEvaluation;
+      biFactors= bufBiFactors;
+      lift= bufLift;
+    }
+    else
+    {
+      if (bufBiFactors.length() < biFactors.length() || 
+          ((bufLift < lift) && (bufBiFactors.length() == biFactors.length())))
+      {
+        Aeval= bufAeval;
+        evaluation= bufEvaluation;
+        biFactors= bufBiFactors;
+        lift= bufLift;
+      }
+    }
+    int k= 0;
+    for (CFListIterator j= bufEvaluation; j.hasItem(); j++, k++)
+      evalPoly += j.getItem()*power (x, k);
+    list.append (evalPoly);  
+  }
+
+  //shifting to zero
+  A= shift2Zero (A, Aeval, evaluation);
+
+  int* liftBounds;
+  int liftBoundsLength= F.level() - 1;
+  liftBounds= liftingBounds (A, lift);
+  
+  CFList MOD;
+  bool earlySuccess;
+  CFList earlyFactors;
+  TIMING_START (fac_hensel_lift);
+  CFList liftedFactors= henselLiftAndEarly 
+                        (A, MOD, liftBounds, earlySuccess, earlyFactors,
+                         Aeval, biFactors, evaluation, info);
+  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
+
+  if (!extension) 
+  {
+    TIMING_START (fac_factor_recombination);
+    factors= factorRecombination (A, liftedFactors, MOD);
+    TIMING_END_AND_PRINT (fac_factor_recombination, 
+                          "time for factor recombination: ");
+  }
+  else 
+  {
+    TIMING_START (fac_factor_recombination);
+    factors= extFactorRecombination (liftedFactors, A, MOD, info, evaluation);
+    TIMING_END_AND_PRINT (fac_factor_recombination, 
+                          "time for factor recombination: "); 
+  }
+
+  if (earlySuccess)
+    factors= Union (factors, earlyFactors);
+
+  if (!extension) 
+  {
+    for (CFListIterator i= factors; i.hasItem(); i++)
+    {
+      int kk= Aeval.getLast().level();
+      for (CFListIterator j= evaluation; j.hasItem(); j++, kk--)
+      {
+        if (i.getItem().level() < kk)
+          continue;
+        i.getItem()= i.getItem() (Variable (kk) - j.getItem(), kk);
+      }
+    }
+  }
+
+  swap (factors, swapLevel, swapLevel2, x); 
+  append (factors, contentAFactors);
+  decompress (factors, N);
+  normalize (factors);
+
+  delete[] liftBounds;
+
+  return factors;
+}
+
+/// multivariate factorization over an extension of the initial field
+CFList 
+extFactorize (const CanonicalForm& F, const ExtensionInfo& info) 
+{
+  CanonicalForm A= F;
+
+  Variable alpha= info.getAlpha();
+  Variable beta= info.getBeta();
+  int k= info.getGFDegree();
+  char cGFName= info.getGFName();
+  CanonicalForm delta= info.getDelta();
+  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
+  Variable w= Variable (1);
+
+  CFList factors;
+  if (!GF && alpha == w)  // we are in F_p
+  {
+    CFList factors;
+    bool extension= true;         
+    int p= getCharacteristic();
+    if (p*p < (1<<16)) // pass to GF if possible
+    { 
+      setCharacteristic (getCharacteristic(), 2, 'Z');
+      ExtensionInfo info= ExtensionInfo (extension);
+      A= A.mapinto();
+      factors= multiFactorize (A, info);
+
+      Variable vBuf= rootOf (gf_mipo);
+      setCharacteristic (getCharacteristic());
+      for (CFListIterator j= factors; j.hasItem(); j++) 
+	j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
+    }
+    else  // not able to pass to GF, pass to F_p(\alpha)
+    {
+      Variable v= chooseExtension (A, beta);
+      ExtensionInfo info= ExtensionInfo (v, extension);
+      factors= multiFactorize (A, info);
+    }
+    return factors;
+  }
+  else if (!GF && (alpha != w)) // we are in F_p(\alpha)
+  { 
+    if (k == 1) // need factorization over F_p
+    {
+      int extDeg= degree (getMipo (alpha));
+      extDeg++; 
+      CanonicalForm mipo= randomIrredpoly (extDeg + 1, Variable (1));
+      Variable v= rootOf (mipo);
+      ExtensionInfo info= ExtensionInfo (v);
+      factors= biFactorize (A, info);
+    }
+    else 
+    {
+      if (beta == Variable (1))
+      {
+        Variable v= chooseExtension (A, alpha);
+        CanonicalForm primElem, imPrimElem;
+        bool primFail= false;
+        Variable vBuf;
+        primElem= primitiveElement (alpha, vBuf, primFail);
+        ASSERT (!primFail, "failure in integer factorizer");
+        if (primFail)
+          ; //ERROR
+        else
+          imPrimElem= mapPrimElem (primElem, vBuf, v);
+
+        CFList source, dest;
+        CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem, 
+                                   source, dest);
+        ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
+        factors= biFactorize (bufA, info);
+      }
+      else
+      {
+        Variable v= chooseExtension (A, alpha);
+        CanonicalForm primElem, imPrimElem;
+        bool primFail= false;
+        Variable vBuf;
+        ASSERT (!primFail, "failure in integer factorizer");
+        if (primFail)
+          ; //ERROR
+        else
+          imPrimElem= mapPrimElem (delta, beta, v); //oder mapPrimElem (primElem, vBuf, v);
+
+        CFList source, dest;
+        CanonicalForm bufA= mapDown (A, info, source, dest);
+        source= CFList();
+        dest= CFList();
+        bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
+        ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
+        factors= biFactorize (bufA, info);
+      }
+    }
+    return factors;
+  }
+  else // we are in GF (p^k)
+  { 
+    int p= getCharacteristic();
+    int extensionDeg= getGFDegree();
+    bool extension= true;
+    if (k == 1) // need factorization over F_p
+    { 
+      extensionDeg++;
+      if (pow ((double) p, (double) extensionDeg) < (1<<16)) 
+      // pass to GF(p^k+1)
+      { 
+        setCharacteristic (p);
+	Variable vBuf= rootOf (gf_mipo);
+	A= GF2FalphaRep (A, vBuf);
+	setCharacteristic (p, extensionDeg, 'Z');
+        ExtensionInfo info= ExtensionInfo (extension);
+	factors= multiFactorize (A.mapinto(), info);
+      }
+      else // not able to pass to another GF, pass to F_p(\alpha)
+      { 
+	setCharacteristic (p);
+	Variable vBuf= rootOf (gf_mipo);
+	A= GF2FalphaRep (A, vBuf);
+	Variable v= chooseExtension (A, beta);
+        ExtensionInfo info= ExtensionInfo (v, extension);
+	factors= multiFactorize (A, info);
+      }
+    }
+    else // need factorization over GF (p^k)
+    {
+      if (pow ((double) p, (double) 2*extensionDeg) < (1<<16)) 
+      // pass to GF(p^2k)
+      { 
+	setCharacteristic (p, 2*extensionDeg, 'Z');
+        ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
+	factors= multiFactorize (GFMapUp (A, extensionDeg), info); 
+	setCharacteristic (p, extensionDeg, cGFName);
+      }
+      else // not able to pass to GF (p^2k), pass to F_p (\alpha)
+      { 
+	setCharacteristic (p);
+	Variable v1= rootOf (gf_mipo);
+	A= GF2FalphaRep (A, v1);
+	Variable v2= chooseExtension (A, v1); 
+        CanonicalForm primElem, imPrimElem;
+        bool primFail= false;
+        Variable vBuf;
+        primElem= primitiveElement (v1, vBuf, primFail);
+        if (primFail)
+        {
+          ; //ERROR
+        }
+        else
+        {
+          imPrimElem= mapPrimElem (primElem, vBuf, v2);
+        }
+        CFList source, dest;
+        CanonicalForm bufA= mapUp (A, alpha, beta, primElem, imPrimElem, 
+                                     source, dest);
+        ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
+	factors= multiFactorize (bufA, info); 
+	setCharacteristic (p, k, cGFName);
+	for (CFListIterator i= factors; i.hasItem(); i++) 
+	  i.getItem()= Falpha2GFRep (i.getItem());
+      }
+    }
+    return factors;
+  }
+}
+
+#endif
+/* HAVE_NTL */
+
Index: factory/facFqFactorize.h
===================================================================
--- factory/facFqFactorize.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/facFqFactorize.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,532 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file facFqFactorize.h
+ * 
+ * This file provides functions for factorizing a multivariate polynomial over
+ * \f$ F_{p} \f$ , \f$ F_{p}(\alpha ) \f$ or GF.   
+ *
+ * ABSTRACT: So far factor recombination is done naive!        
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#ifndef FAC_FQ_FACTORIZE_H
+#define FAC_FQ_FACTORIZE_H
+
+#include <config.h>
+
+#include "facFqBivar.h"
+#include "DegreePattern.h"
+#include "ExtensionInfo.h"
+#include "cf_util.h"
+#include "facFqSquarefree.h"
+#include "facFqBivarUtil.h"
+
+/// Factorization over a finite field
+///  
+/// @return @a multiFactorize returns a factorization of F
+/// @sa biFactorize(), extFactorize()
+CFList 
+multiFactorize (const CanonicalForm& F,    ///< [in] poly to be factored
+                const ExtensionInfo& info  ///< [in] info about extension
+               );
+
+/// factorize a squarefree multivariate polynomial over \f$ F_{p} \f$
+///
+/// @return @a FpSqrfFactorize returns a list of monic factors, the first 
+///         element is the leading coefficient.     
+/// @sa FqSqrfFactorize(), GFSqrfFactorize()
+inline
+CFList FpSqrfFactorize (const CanonicalForm & F ///< [in] a multivariate poly
+                       ) 
+{
+  ExtensionInfo info= ExtensionInfo (false);
+  CFList result= multiFactorize (F, info);
+  result.insert (Lc(F));
+  return result;
+}
+
+/// factorize a squarefree multivariate polynomial over \f$ F_{p} (\alpha ) \f$
+///
+/// @return @a FqSqrfFactorize returns a list of monic factors, the first 
+///         element is the leading coefficient.     
+/// @sa FpSqrfFactorize(), GFSqrfFactorize() 
+inline
+CFList FqSqrfFactorize (const CanonicalForm & F, ///< [in] a multivariate poly
+                        const Variable& alpha    ///< [in] algebraic variable
+                       ) 
+{
+  ExtensionInfo info= ExtensionInfo (alpha, false);
+  CFList result= multiFactorize (F, info);
+  result.insert (Lc(F));
+  return result;
+}
+
+/// factorize a squarefree multivariate polynomial over GF 
+///
+/// @return @a GFSqrfFactorize returns a list of monic factors, the first 
+///         element is the leading coefficient.     
+/// @sa FpSqrfFactorize(), FqSqrfFactorize()
+inline
+CFList GFSqrfFactorize (const CanonicalForm & F ///< [in] a multivariate poly
+                       ) 
+{
+  ASSERT (CFFactory::gettype() == GaloisFieldDomain, 
+          "GF as base field expected");
+  ExtensionInfo info= ExtensionInfo (getGFDegree(), gf_name, false);
+  CFList result= multiFactorize (F, info);
+  result.insert (Lc(F));
+  return result;
+}
+
+/// factorize a multivariate polynomial over \f$ F_{p} \f$ 
+///
+/// @return @a FpFactorize returns a list of monic factors with 
+///         multiplicity, the first element is the leading coefficient.      
+/// @sa FqFactorize(), GFFactorize()
+inline
+CFFList FpFactorize (const CanonicalForm& F ///< [in] a multivariate poly
+                    ) 
+{
+  ExtensionInfo info= ExtensionInfo (false);
+  Variable a= Variable (1);
+  CanonicalForm LcF= Lc (F); 
+  CanonicalForm pthRoot, A;
+  CanonicalForm sqrfP= sqrfPart (F/Lc(F), pthRoot, a);
+  CFList buf, bufRoot;
+  CFFList result, resultRoot;
+  int p= getCharacteristic();
+  int l;
+  if (degree (pthRoot) > 0)
+  {
+    pthRoot= maxpthRoot (pthRoot, p, l);
+    result= FpFactorize (pthRoot);
+    result.removeFirst();
+    for (CFFListIterator i= result; i.hasItem(); i++)
+      i.getItem()= CFFactor (i.getItem().factor(), i.getItem().exp()*l*p);
+    result.insert (CFFactor (LcF, 1));
+    return result;
+  }
+  else
+  {  
+    buf= multiFactorize (sqrfP, info);
+    A= F/LcF;
+    result= multiplicity (A, buf);
+  }
+  if (degree (A) > 0)
+  {
+    resultRoot= FpFactorize (A);
+    resultRoot.removeFirst();
+    result= Union (result, resultRoot);
+  }
+  result.insert (CFFactor (LcF, 1));
+  return result;
+}
+
+/// factorize a multivariate polynomial over \f$ F_{p} (\alpha ) \f$
+///
+/// @return @a FqFactorize returns a list of monic factors with 
+///         multiplicity, the first element is the leading coefficient.      
+/// @sa FpFactorize(), GFFactorize()
+inline
+CFFList FqFactorize (const CanonicalForm& F, ///< [in] a multivariate poly
+                     const Variable& alpha   ///< [in] algebraic variable
+                    ) 
+{
+  ExtensionInfo info= ExtensionInfo (alpha, false); 
+  CanonicalForm LcF= Lc (F); 
+  CanonicalForm pthRoot, A;
+  CanonicalForm sqrfP= sqrfPart (F/Lc(F), pthRoot, alpha);
+  CFList buf, bufRoot;
+  CFFList result, resultRoot;
+  int p= getCharacteristic();
+  int q= ipower (p, degree (getMipo (alpha)));
+  int l;
+  if (degree (pthRoot) > 0)
+  {
+    pthRoot= maxpthRoot (pthRoot, q, l);
+    result= FqFactorize (pthRoot, alpha);
+    result.removeFirst();
+    for (CFFListIterator i= result; i.hasItem(); i++)
+      i.getItem()= CFFactor (i.getItem().factor(), i.getItem().exp()*l*p);
+    result.insert (CFFactor (LcF, 1));
+    return result;
+  }
+  else
+  {  
+    buf= multiFactorize (sqrfP, info);
+    A= F/LcF;
+    result= multiplicity (A, buf);
+  }
+  if (degree (A) > 0)
+  {
+    resultRoot= FqFactorize (A, alpha);
+    resultRoot.removeFirst();
+    result= Union (result, resultRoot);
+  }
+  result.insert (CFFactor (LcF, 1));
+  return result;
+}
+
+/// factorize a multivariate polynomial over GF
+///
+/// @return @a GFFactorize returns a list of monic factors with 
+///         multiplicity, the first element is the leading coefficient.      
+/// @sa FpFactorize(), FqFactorize()
+inline
+CFFList GFFactorize (const CanonicalForm& F ///< [in] a multivariate poly
+                    ) 
+{
+  ASSERT (CFFactory::gettype() == GaloisFieldDomain, 
+          "GF as base field expected");
+  Variable a= Variable (1);
+  ExtensionInfo info= ExtensionInfo (getGFDegree(), gf_name, false); 
+  CanonicalForm LcF= Lc (F);
+  CanonicalForm pthRoot, A;
+  CanonicalForm sqrfP= sqrfPart (F/LcF, pthRoot, a);
+  CFList buf;
+  CFFList result, resultRoot;
+  int p= getCharacteristic();
+  int q= ipower (p, getGFDegree());
+  int l;
+  if (degree (pthRoot) > 0)
+  {
+    pthRoot= maxpthRoot (pthRoot, q, l);
+    result= GFFactorize (pthRoot);
+    result.removeFirst();
+    for (CFFListIterator i= result; i.hasItem(); i++)
+      i.getItem()= CFFactor (i.getItem().factor(), i.getItem().exp()*(l-1)*p);
+    result.insert (CFFactor (LcF, 1));
+    return result;
+  }
+  else
+  {
+    buf= multiFactorize (sqrfP, info);
+    A= F/LcF;
+    result= multiplicity (A, buf);
+  }
+  if (degree (A) > 0)
+  {
+    resultRoot= GFFactorize (A);
+    resultRoot.removeFirst();
+    result= Union (result, resultRoot);
+  }
+  result.insert (CFFactor (LcF, 1));
+  return result;
+}
+
+/// compute the content of @a F wrt Variable (1) using routines from 
+/// @a cf_gcd_smallp.h
+///
+/// @return @a myContent returns the content of F wrt Variable (1)
+static inline
+CanonicalForm 
+myContent (const CanonicalForm& F ///< [in] a poly
+          );
+
+/// compute the content of @a F wrt @a x using routines from 
+/// @a cf_gcd_smallp.h
+///
+/// @return @a myContent returns the content of F wrt x
+static inline
+CanonicalForm
+myContent (const CanonicalForm& F, ///< [in] a poly 
+           const Variable& x       ///< [in] a variable
+          );
+
+/// compute the GCD of all element in @a L using routines from 
+/// @a cf_gcd_smallp.h
+///
+/// @return @a listGCD returns the GCD of all elements in @a L
+static inline
+CanonicalForm 
+listGCD (const CFList& L ///< [in] a list of polys
+        );
+
+/// compute the LCM of @a F and @a G using routines from 
+/// @a cf_gcd_smallp.h
+///
+/// @return @a myLcm returns the LCM of @a F and @a G
+static inline
+CanonicalForm
+myLcm (const CanonicalForm& F,   ///< [in] a poly
+       const CanonicalForm& G    ///< [in] a poly
+      );
+
+/// compute the LCM of the contents of @a A wrt to each variable occuring in @a
+/// A.
+///
+/// @return @a lcmContent returns the LCM of the contents of @a A wrt to each 
+///         variable occuring in @a A.
+static inline 
+CanonicalForm 
+lcmContent (const CanonicalForm& A, ///< [in] a compressed multivariate poly 
+            CFList& contentAi       ///< [in,out] an empty list, returns a list
+                                    ///< of the contents of @a A wrt to each 
+                                    ///< variable occuring in @a A starting from
+                                    ///< @a A.mvar().
+           );
+
+/// compress a polynomial s.t. \f$ deg_{x_{i}} (F) >= deg_{x_{i+1}} (F) \f$ and 
+/// no gaps between the variables occur
+///
+/// @return a compressed poly with the above properties
+static inline
+CanonicalForm myCompress (const CanonicalForm& F, ///< [in] a poly
+                          CFMap& N                ///< [in,out] a map to 
+                                                  ///< decompress
+                         );
+
+/// naive factor recombination for bivariate factorization.
+/// Uses precomputed data to exclude combinations that are not possible.
+///
+/// @return @a monicFactorRecombi returns a list of factors of F, that are 
+///         monic wrt Variable (1).
+/// @sa extFactorRecombination(), factorRecombination(), biFactorizer()
+CFList 
+monicFactorRecombi (
+                  const CFList& factors,    ///< [in] list of lifted factors
+                                            ///< that are monic wrt Variable (1)
+                  const CanonicalForm& F,   ///< [in] bivariate poly 
+                  const CanonicalForm& M,   ///< [in] Variable (2)^liftBound
+                  const DegreePattern& degs ///< [in] degree pattern
+                   );
+
+/// detects factors of @a F at stage @a deg of Hensel lifting.
+/// No combinations of more than one factor are tested. Lift bound and possible
+/// degree pattern are updated. 
+/// 
+/// @return @a earlyMonicFactorDetect returns a list of factors of F (possibly 
+///         incomplete) that are monic wrt Variable (1)
+/// @sa monicFactorRecombi(), earlyFactorDetection(), monicFactorDetect(),
+///     biFactorizer()
+CFList 
+earlyMonicFactorDetect (
+           CanonicalForm& F,       ///< [in,out] poly to be factored, returns
+                                   ///< poly divided by detected factors in case
+                                   ///< of success
+           CFList& factors,        ///< [in,out] list of factors lifted up to 
+                                   ///< @a deg, returns a list of factors
+                                   ///< without detected factors
+           int& adaptedLiftBound,  ///< [in,out] adapted lift bound
+           DegreePattern& degs,    ///< [in,out] degree pattern, is updated
+                                   ///< whenever we find a factor
+           bool& success,          ///< [in,out] indicating success
+           int deg,                ///< [in] stage of Hensel lifting
+           const int bound         ///< [in] degree (A, 2) + 1 + 
+                                   ///< degree (LC (A, 1), 2), where A is the
+                                   ///< multivariate polynomial corresponding to
+                                   ///< F. 
+                       );
+
+/// Bivariate factorization. In contrast to biFactorize() the factors returned 
+/// are monic wrt Variable (1), if @a F is not irreducible. And no factorization
+/// wrt Variable (2) are performed. However, 
+///
+/// @return @a biFactorizer returns a list of factors that are monic wrt 
+///         Variable (1), if @a F is irreducible @a F is returned
+/// @sa multiFactorize(), biFactorize() 
+CFList biFactorizer (const CanonicalForm& F,   ///< [in] a bivariate poly 
+                     const Variable& alpha,    ///< [in] algebraic variable
+                     CanonicalForm& bivarEval, ///< [in,out] a valid evaluation 
+                                               ///< point  
+                     int& liftBound            ///< [in,out] lift bound, may be
+                                               ///< adapted during Hensel 
+                                               ///< lifting 
+                    );
+
+/// Naive factor recombination for multivariate factorization over an extension 
+/// of the initial field. No precomputed is used to exclude combinations. 
+///
+/// @return @a extFactorRecombination returns a list of factors of @a F, whose
+///         shift to zero is reversed.
+/// @sa factorRecombination()   
+CFList 
+extFactorRecombination (
+                 const CFList& factors,     ///< [in] list of lifted factors
+                                            ///< that are monic wrt Variable (1) 
+                 const CanonicalForm& F,    ///< [in] poly to be factored
+                 const CFList& M,           ///< [in] a list of powers of
+                                            ///< Variables 
+                 const ExtensionInfo& info, ///< [in] info about extension 
+                 const CFList& evaluation   ///< [in] evaluation point
+                       );
+
+/// Naive factor recombination for multivariate factorization.
+/// No precomputed is used to exclude combinations.
+///
+/// @return @a factorRecombination returns a list of factors of @a F 
+/// @sa extFactorRecombination() 
+CFList 
+factorRecombination (const CanonicalForm& F,///< [in] poly to be factored
+                     const CFList& factors, ///< [in] list of lifted factors
+                                            ///< that are monic wrt Variable (1)
+                     const CFList& M        ///< [in] a list of powers of
+                                            ///< Variables
+                    );
+
+/// Lift bound adaption. Essentially an early factor detection but only the lift
+/// bound is adapted.
+///
+/// @return @a liftBoundAdaption returns an adapted lift bound.
+/// @sa earlyFactorDetect(), earlyFactorDetection()    
+int
+liftBoundAdaption (const CanonicalForm& F, ///< [in] a poly
+                   const CFList& factors,  ///< [in] list of list of lifted 
+                                           ///< factors that are monic wrt 
+                                           ///< Variable (1)
+                   bool& success,          ///< [in,out] indicates that no
+                                           ///< further lifting is necessary
+                   const int deg,          ///< [in] stage of Hensel lifting
+                   const CFList& MOD,      ///< [in] a list of powers of
+                                           ///< Variables
+                   const int bound         ///< [in] initial lift bound
+                  );
+
+/// Lift bound adaption over an extension of the initial field. Essentially an
+///early factor detection but only the lift bound is adapted.
+///
+/// @return @a liftBoundAdaption returns an adapted lift bound.
+/// @sa earlyFactorDetect(), earlyFactorDetection() 
+int 
+extLiftBoundAdaption (
+	    const CanonicalForm& F,    ///< [in] a poly
+	    const CFList& factors,     ///< [in] list of list of lifted 
+                                       ///< factors that are monic wrt 
+	    bool& success,             ///< [in,out] indicates that no further 
+                                       ///< lifting is necessary
+	    const ExtensionInfo& info, ///< [in] info about extension
+	    const CFList& eval,        ///< [in] evaluation point
+	    const int deg,             ///< [in] stage of Hensel lifting
+	    const CFList& MOD,         ///< [in] a list of powers of
+                                       ///< Variables 
+	    const int bound            ///< [in] initial lift bound
+                     );
+
+/// detects factors of @a F at stage @a deg of Hensel lifting.
+/// No combinations of more than one factor are tested. Lift bound is adapted. 
+/// 
+/// @return @a earlyFactorDetect returns a list of factors of F (possibly 
+///         incomplete), in case of success. Otherwise an empty list.
+/// @sa factorRecombination(), extEarlyFactorDetect()
+CFList 
+earlyFactorDetect (
+		CanonicalForm& F,      ///< [in,out] poly to be factored, 
+                                       ///< returns poly divided by detected
+                                       ///< factors in case of success
+		CFList& factors,       ///< [in,out] list of factors lifted up
+                                       ///< to @a deg, returns a list of factors
+                                       ///< without detected factors
+		int& adaptedLiftBound, ///< [in,out] adapted lift bound
+		bool& success,         ///< [in,out] indicating success
+		const int deg,         ///< [in] stage of Hensel lifting
+		const CFList& MOD,     ///< [in] a list of powers of
+                                       ///< Variables
+		const int bound        ///< [in] initial lift bound
+                  );
+
+/// detects factors of @a F at stage @a deg of Hensel lifting.
+/// No combinations of more than one factor are tested. Lift bound is adapted. 
+/// 
+/// @return @a extEarlyFactorDetect returns a list of factors of F (possibly 
+///         incomplete), whose shift to zero is reversed, in case of success.
+///         Otherwise an empty list. 
+/// @sa factorRecombination(), earlyFactorDetection()
+CFList 
+extEarlyFactorDetect (
+	    CanonicalForm& F,          ///< [in,out] poly to be factored, 
+                                       ///< returns poly divided by detected
+                                       ///< factors in case of success
+	    CFList& factors,           ///< [in,out] list of factors lifted up
+                                       ///< to @a deg, returns a list of factors
+                                       ///< without detected factors
+	    int& adaptedLiftBound,     ///< [in,out] adapted lift bound
+	    bool& success,             ///< [in,out] indicating succes
+	    const ExtensionInfo& info, ///< [in] info about extension
+	    const CFList& eval,        ///< [in] evaluation point
+	    const int deg,             ///< [in] stage of Hensel lifting
+	    const CFList& MOD,         ///< [in] a list of powers of Variables
+	    const int bound            ///< [in] initial lift bound
+                     );
+
+/// evaluation point search for multivariate factorization, 
+/// looks for a (F.level() - 1)-tuple such that the resulting univariate
+/// polynomial has main variable F.mvar(), is squarefree and its degree
+/// coincides with degree(F) and the bivariate one is primitive wrt. 
+/// Variable(1), fails if there are no valid evaluation points, eval contains 
+/// the intermediate evaluated polynomials.  
+///
+/// @return @a evalPoints returns an evaluation point, which is valid if and
+///         only if fail == false. 
+CFList 
+evalPoints (const CanonicalForm& F,  ///< [in] a compressed poly 
+            CFList & eval,           ///< [in,out] an empty list, returns @a F
+                                     ///< successive evaluated
+            const Variable& alpha,   ///< [in] algebraic variable
+            CFList& list,            ///< [in,out] a list of points already
+                                     ///< considered, a point is encoded as a
+                                     ///< poly of degree F.level()-1 in
+                                     ///< Variable(1) 
+            const bool& GF,          ///< [in] GF?
+            bool& fail               ///< [in,out] indicates failure
+           );
+
+/// looks for a new main variable of level higher than lev which omitts a valid
+/// evaluation, and returns its level. If there is no such variable 0 is
+/// returned
+///
+/// @return @a newMainVariableSearch returns the level of the new main variable,
+///         if there no variable which omitts a valid evaluation 0 is returned 
+static inline
+int newMainVariableSearch (
+              CanonicalForm& A,     ///< [in,out] a compressed poly, returns
+                                    ///< the swapped initial poly or the 
+                                    ///< initial poly if 0 is returned 
+	      CFList& Aeval,        ///< [in,out] @a A successively evaluated,
+                                    ///< empty list if 0 is returned
+	      CFList& evaluation,   ///< [in,out] evaluation point, empty list
+                                    ///< if 0 is returned 
+	      const Variable& alpha,///< [in] algebraic variable
+	      const int lev         ///< [in] some int <= A.level()
+                          );
+/// hensel Lifting and early factor detection 
+///
+/// @return @a henselLiftAndEarly returns monic (wrt Variable (1)) lifted 
+///         factors without factors which have been detected at an early stage
+///         of Hensel lifting
+/// @sa earlyFactorDetectn(), extEarlyFactorDetect()
+CFList
+henselLiftAndEarly (
+	    CanonicalForm& A,         ///< [in,out] poly to be factored,  
+                                      ///< returns poly divided by detected
+                                      ///< factors, in case of success 
+	    CFList& MOD,              ///< [in,out] a list of powers of
+                                      ///< Variables
+	    int*& liftBounds,         ///< [in,out] initial lift bounds, returns
+                                      ///< adapted lift bounds
+	    bool& earlySuccess,       ///< [in,out] indicating success
+	    CFList& earlyFactors,     ///< [in,out] early factors
+	    const CFList& Aeval,      ///< [in] @a A successively evaluated at
+                                      ///< elements of @a evaluation
+	    const CFList& biFactors,  ///< [in] bivariate factors
+	    const CFList& evaluation, ///< [in] evaluation point
+	    const ExtensionInfo& info ///< [in] info about extension
+                   );
+
+/// Factorization over an extension of initial field
+///
+/// @return @a extFactorize returns factorization of F over initial field
+/// @sa extBiFactorize(), multiFactorize()
+CFList 
+extFactorize (const CanonicalForm& F,   ///< [in] poly to be factored 
+              const ExtensionInfo& info ///< [in] info about extension 
+             ); 
+
+#endif
+/* FAC_FQ_FACTORIZE_H */
+
Index: factory/facFqFactorizeUtil.cc
===================================================================
--- factory/facFqFactorizeUtil.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/facFqFactorizeUtil.cc	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,177 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file facFqFactorizeUtil.cc
+ * 
+ * This file provides utility functions for multivariate factorization
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#include <config.h>
+
+#include "canonicalform.h"
+#include "cf_map.h"
+
+static inline 
+void decompressAppend (CFList& factors1, const CFList& factors2, const CFMap& N) 
+{
+  for (CFListIterator i= factors1; i.hasItem(); i++) 
+    i.getItem()= N (i.getItem());
+  for (CFListIterator i= factors2; i.hasItem(); i++)
+    factors1.append (N (i.getItem()));
+  return; 
+}
+
+static inline 
+void appendSwap (CFList& factors1, const CFList& factors2, const int
+                  swapLevel1, const int swapLevel2, const Variable& x) 
+{
+  for (CFListIterator i= factors2; i.hasItem(); i++) 
+  { 
+    if (swapLevel1) 
+    {
+      if (swapLevel2)
+	factors1.append (swapvar (swapvar (i.getItem(), x, 
+                         Variable (swapLevel2)), Variable (swapLevel1), x));
+      else
+	factors1.append (swapvar (i.getItem(), Variable (swapLevel1), x));
+    }
+    else 
+    {
+      if (swapLevel2)
+	factors1.append (swapvar (i.getItem(), x, Variable (swapLevel2)));
+      else 
+	factors1.append (i.getItem());
+    }
+  }
+  return;
+}
+
+static inline 
+void swap (CFList& factors, const int swapLevel1, const int swapLevel2, const
+           Variable& x) 
+{
+  for (CFListIterator i= factors; i.hasItem(); i++) 
+  { 
+    if (swapLevel1) 
+    {
+      if (swapLevel2)
+	i.getItem()= swapvar (swapvar (i.getItem(), x, Variable (swapLevel2)),
+                              Variable (swapLevel1), x);
+      else
+	i.getItem()= swapvar (i.getItem(), Variable (swapLevel1), x);
+    }
+    else 
+    {
+      if (swapLevel2)
+	i.getItem()= swapvar (i.getItem(), x, Variable (swapLevel2));
+    }
+  }
+  return;
+}
+
+static inline 
+void appendSwapDecompress (CFList& factors1, const CFList& factors2, 
+                             const CFMap& N, const int swapLevel, const
+                             Variable& x) 
+{
+  for (CFListIterator i= factors1; i.hasItem(); i++) 
+  {
+    if (swapLevel)   
+      i.getItem()= swapvar (i.getItem(), Variable (swapLevel), x);
+    i.getItem()= N(i.getItem());
+  }
+  for (CFListIterator i= factors2; i.hasItem(); i++) 
+  {
+    if (!i.getItem().inCoeffDomain()) 
+      factors1.append (N (i.getItem()));
+  }
+  return;
+}
+
+static inline 
+void appendSwapDecompress (CFList& factors1, const CFList& factors2, 
+                             const CFMap& N, const int swapLevel1, 
+                             const int swapLevel2, const Variable& x) 
+{
+  for (CFListIterator i= factors1; i.hasItem(); i++) 
+  {
+    if (swapLevel1) 
+    {
+      if (swapLevel2) 
+        i.getItem()=
+        N (swapvar (swapvar (i.getItem(), Variable (swapLevel2), x), x,
+                    Variable (swapLevel1)));
+      else 
+        i.getItem()= N (swapvar (i.getItem(), x, Variable (swapLevel1)));
+    }
+    else 
+    {
+      if (swapLevel2)
+        i.getItem()= N (swapvar (i.getItem(), Variable (swapLevel2), x));
+      else
+        i.getItem()= N (i.getItem());
+    }
+  }
+  for (CFListIterator i= factors2; i.hasItem(); i++) 
+  { 
+    if (!i.getItem().inCoeffDomain())  
+      factors1.append (N (i.getItem()));
+  }
+  return;
+}
+
+static inline
+int* liftingBounds (const CanonicalForm& A, const int& bivarLiftBound) 
+{
+  int j= A.level() - 1;
+  int* liftBounds= new int [j];
+  liftBounds[0]= bivarLiftBound;
+  for (int i= 1; i < j; i++) 
+  {
+    liftBounds[i]= degree (A, Variable (i + 2)) + 1 + 
+                            degree (LC (A, 1), Variable (i + 2));
+  }
+  return liftBounds;
+}
+
+static inline
+CanonicalForm
+shift2Zero (const CanonicalForm& F, CFList& Feval, const CFList& evaluation)
+{
+  CanonicalForm A= F;
+  int k= A.level();
+  for (CFListIterator i= evaluation; i.hasItem(); i++, k--)
+    A= A (Variable (k) + i.getItem(), k);
+  CanonicalForm buf= A;
+  Feval= CFList();
+  Feval.append (buf);
+  for (k= A.level(); k > 2; k--)
+  {
+    buf= mod (buf, Variable (k));
+    Feval.insert (buf);
+  }
+  return A;
+}
+
+static inline
+CanonicalForm reverseShift (const CanonicalForm& F, const CFList& evaluation)
+{
+  int l= evaluation.length() + 1;
+  CanonicalForm result= F;
+  CFListIterator j= evaluation;
+  for (int i= l; i > 1; i--, j++)
+  {
+    if (F.level() < i)
+      continue;
+    result= result (Variable (i) - j.getItem(), i);
+  }
+  return result;
+}
+
+
Index: factory/facFqFactorizeUtil.h
===================================================================
--- factory/facFqFactorizeUtil.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
+++ factory/facFqFactorizeUtil.h	(revision 7bf14550d613edf557228ca502041b4a91b07ce3)
@@ -0,0 +1,116 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR    
+\*****************************************************************************/
+/** @file facFqFactorizeUtil.h
+ * 
+ * This file provides utility functions for multivariate factorization       
+ *
+ * @author Martin Lee
+ *
+ * @internal @version \$Id$
+ *
+ **/
+/*****************************************************************************/
+
+#ifndef FAC_FQ_FACTORIZE_UTIL_H
+#define FAC_FQ_FACTORIZE_UTIL_H
+
+#include <config.h>
+
+#include "canonicalform.h"
+#include "facFqFactorizeUtil.cc"
+
+/// append @a factors2 to @a factors1 and decompress
+///
+void decompressAppend (CFList& factors1,       ///< [in,out] a list of polys,
+                                               ///< returns @a factors2 appended
+                                               ///< to it and everything is 
+                                               ///< decompressed
+                       const CFList& factors2, ///< [in] a list of polys
+                       const CFMap& N          ///< [in] a map
+                      );
+
+/// swap elements in @a factors2 and append them to @a factors1
+/// 
+void appendSwap (CFList& factors1,       ///< [in,out] a list of polys,
+                                         ///< returns swapped elements of @a 
+                                         ///< factors2 appended to it
+                 const CFList& factors2, ///< [in] a list of polys
+                 const int swapLevel1,   ///< [in] level of variable to be 
+                                         ///< swapped with x, 0 if no swapping
+                 const int swapLevel2,   ///< [in] level of variable to be 
+                                         ///< swapped with x, 0 if no swapping
+                 const Variable& x       ///< [in] a variable
+                );
+
+/// swap elements in @a factors
+void swap (CFList& factors,       ///< [in] a list of polys, returns swapped 
+                                  ///< elements of factors
+           const int swapLevel1, ///< [in] level of variable to be 
+                                  ///< swapped with x, 0 if no swapping
+           const int swapLevel2, ///< [in] level of variable to be 
+                                  ///< swapped with x, 0 if no swapping
+           const Variable& x      ///< [in] a variable
+          );
+
+/// swap elements of @a factors2, append them to @a factors1 and decompress
+void appendSwapDecompress ( 
+                  CFList& factors1,       ///< [in,out] a list of polys,
+                                          ///< returns swapped elements of @a 
+                                          ///< factors2 appended to it and 
+                                          ///< everything is decompressed
+                  const CFList& factors2, ///< [in] a list of polys
+                  const CFMap& N,         ///< [in] a map
+                  const int swapLevel,    ///< [in] level of variable to be 
+                                          ///< swapped with x, 0 if no swapping
+                  const Variable& x       ///< [in] a variable
+                          );
+
+/// swap elements of @a factors2, append them to @a factors1 and decompress
+void appendSwapDecompress (
+                  CFList& factors1,       ///< [in,out] a list of polys,
+                                          ///< returns swapped elements of @a 
+                                          ///< factors2 appended to it and 
+                                          ///< everything is decompressed
+                  const CFList& factors2, ///< [in] a list of polys
+                  const CFMap& N,         ///< [in] a map
+                  const int swapLevel1,   ///< [in] level of variable to be 
+                                          ///< swapped with x, 0 if no swapping
+                  const int swapLevel2,   ///< [in] level of variable to be 
+                                          ///< swapped with x, 0 if no swapping
+                  const Variable& x       ///< [in] a variable
+                          ); 
+
+
+/// compute lifting bounds
+///
+/// @return @a liftingBounds returns an array containing the lift bounds for A 
+int* liftingBounds (const CanonicalForm& A,    ///< [in] a compressed poly    
+                    const int& bivarLiftBound  ///< [in] lift bound for 
+                                               ///< @a biFactorizer()
+                   );
+
+/// shift evaluation point to zero
+/// @return @a shift2Zero returns @a F shifted by @a evaluation s.t. 0 is a 
+///         valid evaluation point
+/// @sa evalPoints(), reverseShift() 
+CanonicalForm
+shift2Zero (const CanonicalForm& F,  ///< [in] a compressed poly
+            CFList& Feval,           ///< [in,out] an empty list, returns
+                                     ///< @a F successively evaluated
+                                     ///< at 0
+            const CFList& evaluation ///< [in] a valid evaluation point
+           );
+
+/// reverse shifting the evaluation point to zero
+///
+/// @return @a reverseShift returns a poly whose shift to zero is reversed
+/// @sa shift2Zero(), evalPoints()
+CanonicalForm reverseShift (const CanonicalForm& F,  ///< [in] a compressed poly
+                            const CFList& evaluation ///< [in] a valid 
+                                                     ///< evaluation point
+                           );
+
+#endif
+/* FAC_FQ_FACTORIZE_UTIL_H */
+
