Index: factory/cf_gcd_smallp.cc
===================================================================
--- factory/cf_gcd_smallp.cc	(revision 0e2e23228417d64209e60f3d18fa1bd11233839d)
+++ factory/cf_gcd_smallp.cc	(revision 81d96cef80a045a0f0adf2d6e851d67d9d72074e)
@@ -26,4 +26,5 @@
 
 #include "canonicalform.h"
+#include "algext.h"
 #include "cf_map.h"
 #include "cf_util.h"
@@ -4032,6 +4033,6 @@
   lcs [0]= shiftedLCsEval1.getFirst();
   lcs [1]= shiftedLCsEval2.getFirst();
-  henselLift122 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
-                 lcs, false);
+  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
+                        lcs, false);
 
   for (CFListIterator i= factors; i.hasItem(); i++)
@@ -4049,7 +4050,7 @@
     for (int i= 1; i < U.level() - 1; i++)
       liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
-    factors= henselLift2 (UEval, factors, liftBounds, U.level() - 1, false,
-                          shiftedLCsEval1, shiftedLCsEval2, Pi, diophant,
-                          noOneToOne);
+    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
+                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
+                                  diophant, noOneToOne);
     delete [] liftBounds;
     if (noOneToOne)
Index: factory/facFactorize.cc
===================================================================
--- factory/facFactorize.cc	(revision 0e2e23228417d64209e60f3d18fa1bd11233839d)
+++ factory/facFactorize.cc	(revision 81d96cef80a045a0f0adf2d6e851d67d9d72074e)
@@ -481,6 +481,6 @@
       CFArray Pi;
       CFList diophant;
-      henselLift122 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
-                     leadingCoeffs, false);
+      nonMonicHenselLift12 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
+                            leadingCoeffs, false);
 
       if (sqrfPartF.level() > 2)
Index: factory/facFqBivarUtil.cc
===================================================================
--- factory/facFqBivarUtil.cc	(revision 0e2e23228417d64209e60f3d18fa1bd11233839d)
+++ factory/facFqBivarUtil.cc	(revision 81d96cef80a045a0f0adf2d6e851d67d9d72074e)
@@ -16,4 +16,5 @@
 
 #include "cf_map.h"
+#include "algext.h"
 #include "cf_map_ext.h"
 #include "templates/ftmpl_functions.h"
Index: factory/facFqFactorize.cc
===================================================================
--- factory/facFqFactorize.cc	(revision 0e2e23228417d64209e60f3d18fa1bd11233839d)
+++ factory/facFqFactorize.cc	(revision 81d96cef80a045a0f0adf2d6e851d67d9d72074e)
@@ -1658,6 +1658,6 @@
       CFArray Pi;
       CFList diophant;
-      henselLift122 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
-                     leadingCoeffs, false);
+      nonMonicHenselLift12 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
+                            leadingCoeffs, false);
 
       if (sqrfPartF.level() > 2)
Index: factory/facHensel.cc
===================================================================
--- factory/facHensel.cc	(revision 0e2e23228417d64209e60f3d18fa1bd11233839d)
+++ factory/facHensel.cc	(revision 81d96cef80a045a0f0adf2d6e851d67d9d72074e)
@@ -25,7 +25,7 @@
 #include "timing.h"
 
+#include "algext.h"
 #include "facHensel.h"
 #include "facMul.h"
-#include "cf_util.h"
 #include "fac_util.h"
 #include "cf_algorithm.h"
@@ -125,5 +125,5 @@
 }
 
-static inline
+static
 CFList mapinto (const CFList& L)
 {
@@ -134,5 +134,5 @@
 }
 
-static inline
+static
 int mod (const CFList& L, const CanonicalForm& p)
 {
@@ -146,5 +146,5 @@
 
 
-static inline void
+static void
 chineseRemainder (const CFList & x1, const CanonicalForm & q1,
                   const CFList & x2, const CanonicalForm & q2,
@@ -162,5 +162,5 @@
 }
 
-static inline
+static
 CFList Farey (const CFList& L, const CanonicalForm& q)
 {
@@ -171,5 +171,5 @@
 }
 
-static inline
+static
 CFList replacevar (const CFList& L, const Variable& a, const Variable& b)
 {
@@ -359,5 +359,4 @@
 }
 
-static inline
 CFList diophantine (const CanonicalForm& F, const CFList& factors)
 {
@@ -669,7 +668,6 @@
 }
 
-static inline
 CFList
-biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
+biDiophantine (const CanonicalForm& F, const CFList& factors, int d)
 {
   Variable y= F.mvar();
@@ -770,8 +768,7 @@
 }
 
-static inline
 CFList
 multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
-                     const CFList& recResult, const CFList& M, const int d)
+                     const CFList& recResult, const CFList& M, int d)
 {
   Variable y= F.mvar();
@@ -859,5 +856,5 @@
 }
 
-static inline
+static
 void
 henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
@@ -1101,5 +1098,5 @@
 
 CFList
-henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
+henselLift23 (const CFList& eval, const CFList& factors, int* l, CFList&
               diophant, CFArray& Pi, CFMatrix& M)
 {
@@ -1168,6 +1165,5 @@
 CFList
 henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
-            diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
-            lNew)
+            diophant, CFArray& Pi, CFMatrix& M, int lOld, int lNew)
 {
   diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
@@ -1209,6 +1205,6 @@
 
 CFList
-henselLift (const CFList& eval, const CFList& factors, const int* l, const int
-            lLength, bool sort)
+henselLift (const CFList& eval, const CFList& factors, int* l, int lLength,
+            bool sort)
 {
   CFList diophant;
@@ -1243,8 +1239,10 @@
 }
 
+// nonmonic
+
 void
-henselStep122 (const CanonicalForm& F, const CFList& factors,
-              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
-              CFArray& Pi, int j, const CFArray& /*LCs*/)
+nonMonicHenselStep12 (const CanonicalForm& F, const CFList& factors,
+                      CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
+                      CFArray& Pi, int j, const CFArray& /*LCs*/)
 {
   Variable x= F.mvar();
@@ -1455,6 +1453,7 @@
 
 void
-henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
-              CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
+nonMonicHenselLift12 (const CanonicalForm& F, CFList& factors, int l,
+                      CFArray& Pi, CFList& diophant, CFMatrix& M,
+                      const CFArray& LCs, bool sort)
 {
   if (sort)
@@ -1522,5 +1521,5 @@
 
   for (i= 1; i < l; i++)
-    henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
+    nonMonicHenselStep12 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
 
   factors= CFList();
@@ -1533,5 +1532,4 @@
 /// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
 /// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
-static inline
 CFList
 diophantine (const CFList& recResult, const CFList& factors,
@@ -1612,9 +1610,9 @@
 }
 
-static inline
 void
-henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
-            const CFList& diophant, CFMatrix& M, CFArray& Pi,
-            const CFList& products, int j, const CFList& MOD, bool& noOneToOne)
+nonMonicHenselStep (const CanonicalForm& F, const CFList& factors,
+                    CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
+                    CFArray& Pi, const CFList& products, int j,
+                    const CFList& MOD, bool& noOneToOne)
 {
   CanonicalForm E;
@@ -1852,7 +1850,7 @@
 
 CFList
-henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
-              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
-              const CFList& LCs2, bool& bad)
+nonMonicHenselLift232(const CFList& eval, const CFList& factors, int* l, CFList&
+                      diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
+                      const CFList& LCs2, bool& bad)
 {
   CFList buf= factors;
@@ -1902,5 +1900,6 @@
   for (int d= 1; d < l[1]; d++)
   {
-    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
+    nonMonicHenselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, products,
+                        d, MOD, bad);
     if (bad)
       return CFList();
@@ -1914,7 +1913,8 @@
 
 CFList
-henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
-            diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
-            lNew, const CFList& LCs1, const CFList& LCs2, bool& bad)
+nonMonicHenselLift2 (const CFList& F, const CFList& factors, const CFList& MOD,
+                    CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
+                    int& lNew, const CFList& LCs1, const CFList& LCs2, bool& bad
+                   )
 {
   int k= 0;
@@ -1963,5 +1963,6 @@
   for (int d= 1; d < lNew; d++)
   {
-    henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
+    nonMonicHenselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, products,
+                        d, MOD, bad);
     if (bad)
       return CFList();
@@ -1975,7 +1976,7 @@
 
 CFList
-henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
-             lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
-             const CFArray& Pi, const CFList& diophant, bool& bad)
+nonMonicHenselLift2 (const CFList& eval, const CFList& factors, int* l, int
+                    lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
+                    const CFArray& Pi, const CFList& diophant, bool& bad)
 {
   CFList bufDiophant= diophant;
@@ -1985,6 +1986,6 @@
   CFArray bufPi= Pi;
   CFMatrix M= CFMatrix (l[1], factors.length());
-  CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2,
-                               bad);
+  CFList result= 
+    nonMonicHenselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2, bad);
   if (bad)
     return CFList();
@@ -2014,6 +2015,6 @@
     bufLCs2.append (jjj.getItem());
     M= CFMatrix (l[i], factors.length());
-    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
-                         l[i], bufLCs1, bufLCs2, bad);
+    result= nonMonicHenselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M,
+                                 l[i - 1], l[i], bufLCs1, bufLCs2, bad);
     if (bad)
       return CFList();
@@ -2112,5 +2113,6 @@
   for (int d= 1; d < liftBound; d++)
   {
-    henselStep2 (F, factors, bufFactors, diophant, M, Pi, products, d, MOD, bad);
+    nonMonicHenselStep (F, factors, bufFactors, diophant, M, Pi, products, d,
+                        MOD, bad);
     if (bad)
       return CFList();
@@ -2125,5 +2127,5 @@
 CFList
 nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
-                    CFList& diophant, CFArray& Pi, CFMatrix& M, const int lOld,
+                    CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
                     int& lNew, const CFList& MOD, bool& noOneToOne
                    )
@@ -2192,6 +2194,6 @@
   for (int d= 1; d < lNew; d++)
   {
-    henselStep2 (F.getLast(), factors, bufFactors, diophant, M, Pi, products, d,
-                 MOD, noOneToOne);
+    nonMonicHenselStep (F.getLast(), factors, bufFactors, diophant, M, Pi,
+                        products, d, MOD, noOneToOne);
     if (noOneToOne)
       return CFList();
Index: factory/facHensel.h
===================================================================
--- factory/facHensel.h	(revision 0e2e23228417d64209e60f3d18fa1bd11233839d)
+++ factory/facHensel.h	(revision 81d96cef80a045a0f0adf2d6e851d67d9d72074e)
@@ -25,7 +25,4 @@
 
 #include "canonicalform.h"
-#include "cf_iter.h"
-#include "templates/ftmpl_functions.h"
-#include "algext.h"
 
 #ifdef HAVE_NTL
@@ -87,5 +84,5 @@
                                      ///< including leading coefficient
                                      ///< as first element.
-              const int* l,          ///< [in] l[0]: precision of bivariate
+              int* l,          ///< [in] l[0]: precision of bivariate
                                      ///< lifting, l[1] as above
               CFList& diophant,      ///< [in,out] returns the result of
@@ -131,6 +128,6 @@
             CFArray& Pi,           ///< [in,out] stores intermediate results
             CFMatrix& M,           ///< [in,out] stores intermediate results
-            const int lOld,       ///< [in] lifting precision of F.mvar()
-            const int lNew        ///< [in] lifting precision of G.mvar()
+            int lOld,       ///< [in] lifting precision of F.mvar()
+            int lNew        ///< [in] lifting precision of G.mvar()
            );
 
@@ -148,14 +145,13 @@
             const CFList& factors, ///< [in] bivariate factors, including
                                    ///< leading coefficient
-            const int* l,          ///< [in] lifting bounds
-            const int lLength,     ///< [in] length of l
+            int* l,          ///< [in] lifting bounds
+            int lLength,     ///< [in] length of l
             bool sort= true        ///< [in] sort factors by degree in
                                    ///< Variable(1)
            );
 
-/// two factor Hensel lifting from univariate to bivariate, factors need not to
-/// be monic
-void
-henselLift122 (const CanonicalForm& F,///< [in] a bivariate poly
+/// Hensel lifting from univariate to bivariate, factors need not to be monic
+void
+nonMonicHenselLift12 (const CanonicalForm& F,///< [in] a bivariate poly
                CFList& factors,       ///< [in, out] a list of univariate polys
                                       ///< lifted factors
@@ -174,5 +170,5 @@
 /// @return @a henselLift122 returns a list of lifted factors
 CFList
-henselLift2 (const CFList& eval,   ///< [in] a list of polynomials the last
+nonMonicHenselLift2 (const CFList& eval,///< [in] a list of polynomials the last
                                    ///< element is a compressed multivariate
                                    ///< poly, last but one element equals the
@@ -182,5 +178,5 @@
              const CFList& factors,///< [in] bivariate factors
              int* l,               ///< [in] lift bounds
-             const int lLength,    ///< [in] length of l
+             int lLength,          ///< [in] length of l
              bool sort,            ///< [in] if true factors are sorted by
                                    ///< their degree in Variable(1)
Index: factory/facMul.cc
===================================================================
--- factory/facMul.cc	(revision 0e2e23228417d64209e60f3d18fa1bd11233839d)
+++ factory/facMul.cc	(revision 81d96cef80a045a0f0adf2d6e851d67d9d72074e)
@@ -17,4 +17,5 @@
 #include "facMul.h"
 #include "algext.h"
+#include "cf_util.h"
 #include "templates/ftmpl_functions.h"
 
@@ -26,4 +27,6 @@
 #include "FLINTconvert.h"
 #endif
+
+// univariate polys
 
 #ifdef HAVE_FLINT
@@ -48,6 +51,6 @@
 
 CanonicalForm
-reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha,
-              const CanonicalForm& den)
+reverseSubstQa (const fmpz_poly_t F, int d, const Variable& alpha,
+                const CanonicalForm& den)
 {
   Variable x= Variable (1);
@@ -111,5 +114,5 @@
 
   denA *= denB;
-  A= reverseSubst (FLINTA, d, alpha, denA);
+  A= reverseSubstQa (FLINTA, d, alpha, denA);
 
   fmpz_poly_clear (FLINTA);
@@ -222,5 +225,5 @@
 
   denA *= denB;
-  A= reverseSubst (FLINTA, d, alpha, denA);
+  A= reverseSubstQa (FLINTA, d, alpha, denA);
   fmpz_poly_clear (FLINTA);
   fmpz_poly_clear (FLINTB);
@@ -570,4 +573,8 @@
 }
 
+// end univariate polys
+//*************************
+// bivariate polys
+
 #ifdef HAVE_FLINT
 void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
@@ -592,32 +599,4 @@
 }
 
-/*void kronSubQ (fmpz_poly_t result, const CanonicalForm& A, int d)
-{
-  int degAy= degree (A);
-  fmpz_poly_init2 (result, d*(degAy + 1));
-  _fmpz_poly_set_length (result, d*(degAy+1));
-
-  CFIterator j;
-  for (CFIterator i= A; i.hasTerms(); i++)
-  {
-    if (i.coeff().inBas
-    convertFacCF2Fmpz_poly_t (buf, i.coeff());
-
-    int k= i.exp()*d;
-    int bufRepLength= (int) fmpz_poly_length (buf);
-    for (int j= 0; j < bufRepLength; j++)
-    {
-      fmpz_poly_get_coeff_fmpz (coeff, buf, j);
-      fmpz_poly_set_coeff_fmpz (result, j + k, coeff);
-    }
-    fmpz_poly_clear (buf);
-  }
-  fmpz_clear (coeff);
-  _fmpz_poly_normalise (result);
-}*/
-
-// A is a bivariate poly over Qa!!!!
-// d2= 2*deg_alpha + 1
-// d1= 2*deg_x*d2+1
 void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
 {
@@ -664,593 +643,8 @@
   _fmpq_poly_normalise (result);
 }
-#endif
-
-zz_pX kronSubFp (const CanonicalForm& A, int d)
-{
-  int degAy= degree (A);
-  zz_pX result;
-  result.rep.SetLength (d*(degAy + 1));
-
-  zz_p *resultp;
-  resultp= result.rep.elts();
-  zz_pX buf;
-  zz_p *bufp;
-  int j, k, bufRepLength;
-
-  for (CFIterator i= A; i.hasTerms(); i++)
-  {
-    if (i.coeff().inCoeffDomain())
-      buf= convertFacCF2NTLzzpX (i.coeff());
-    else
-      buf= convertFacCF2NTLzzpX (i.coeff());
-
-    k= i.exp()*d;
-    bufp= buf.rep.elts();
-    bufRepLength= (int) buf.rep.length();
-    for (j= 0; j < bufRepLength; j++)
-      resultp [j + k]= bufp [j];
-  }
-  result.normalize();
-
-  return result;
-}
-
-zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
-{
-  int degAy= degree (A);
-  zz_pEX result;
-  result.rep.SetLength (d*(degAy + 1));
-
-  Variable v;
-  zz_pE *resultp;
-  resultp= result.rep.elts();
-  zz_pEX buf1;
-  zz_pE *buf1p;
-  zz_pX buf2;
-  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
-  int j, k, buf1RepLength;
-
-  for (CFIterator i= A; i.hasTerms(); i++)
-  {
-    if (i.coeff().inCoeffDomain())
-    {
-      buf2= convertFacCF2NTLzzpX (i.coeff());
-      buf1= to_zz_pEX (to_zz_pE (buf2));
-    }
-    else
-      buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
-
-    k= i.exp()*d;
-    buf1p= buf1.rep.elts();
-    buf1RepLength= (int) buf1.rep.length();
-    for (j= 0; j < buf1RepLength; j++)
-      resultp [j + k]= buf1p [j];
-  }
-  result.normalize();
-
-  return result;
-}
 
 void
-kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
-                const Variable& alpha)
-{
-  int degAy= degree (A);
-  subA1.rep.SetLength ((long) d*(degAy + 2));
-  subA2.rep.SetLength ((long) d*(degAy + 2));
-
-  Variable v;
-  zz_pE *subA1p;
-  zz_pE *subA2p;
-  subA1p= subA1.rep.elts();
-  subA2p= subA2.rep.elts();
-  zz_pEX buf;
-  zz_pE *bufp;
-  zz_pX buf2;
-  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
-  int j, k, kk, bufRepLength;
-
-  for (CFIterator i= A; i.hasTerms(); i++)
-  {
-    if (i.coeff().inCoeffDomain())
-    {
-      buf2= convertFacCF2NTLzzpX (i.coeff());
-      buf= to_zz_pEX (to_zz_pE (buf2));
-    }
-    else
-      buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
-
-    k= i.exp()*d;
-    kk= (degAy - i.exp())*d;
-    bufp= buf.rep.elts();
-    bufRepLength= (int) buf.rep.length();
-    for (j= 0; j < bufRepLength; j++)
-    {
-      subA1p [j + k] += bufp [j];
-      subA2p [j + kk] += bufp [j];
-    }
-  }
-  subA1.normalize();
-  subA2.normalize();
-}
-
-void
-kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
-{
-  int degAy= degree (A);
-  subA1.rep.SetLength ((long) d*(degAy + 2));
-  subA2.rep.SetLength ((long) d*(degAy + 2));
-
-  zz_p *subA1p;
-  zz_p *subA2p;
-  subA1p= subA1.rep.elts();
-  subA2p= subA2.rep.elts();
-  zz_pX buf;
-  zz_p *bufp;
-  int j, k, kk, bufRepLength;
-
-  for (CFIterator i= A; i.hasTerms(); i++)
-  {
-    buf= convertFacCF2NTLzzpX (i.coeff());
-
-    k= i.exp()*d;
-    kk= (degAy - i.exp())*d;
-    bufp= buf.rep.elts();
-    bufRepLength= (int) buf.rep.length();
-    for (j= 0; j < bufRepLength; j++)
-    {
-      subA1p [j + k] += bufp [j];
-      subA2p [j + kk] += bufp [j];
-    }
-  }
-  subA1.normalize();
-  subA2.normalize();
-}
-
-CanonicalForm
-reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
-              const Variable& alpha)
-{
-  Variable y= Variable (2);
-  Variable x= Variable (1);
-
-  zz_pEX f= F;
-  zz_pEX g= G;
-  int degf= deg(f);
-  int degg= deg(g);
-
-  zz_pEX buf1;
-  zz_pEX buf2;
-  zz_pEX buf3;
-
-  zz_pE *buf1p;
-  zz_pE *buf2p;
-  zz_pE *buf3p;
-  if (f.rep.length() < (long) d*(k+1)) //zero padding
-    f.rep.SetLength ((long)d*(k+1));
-
-  zz_pE *gp= g.rep.elts();
-  zz_pE *fp= f.rep.elts();
-  CanonicalForm result= 0;
-  int i= 0;
-  int lf= 0;
-  int lg= d*k;
-  int degfSubLf= degf;
-  int deggSubLg= degg-lg;
-  int repLengthBuf2, repLengthBuf1, ind, tmp;
-  zz_pE zzpEZero= zz_pE();
-
-  while (degf >= lf || lg >= 0)
-  {
-    if (degfSubLf >= d)
-      repLengthBuf1= d;
-    else if (degfSubLf < 0)
-      repLengthBuf1= 0;
-    else
-      repLengthBuf1= degfSubLf + 1;
-    buf1.rep.SetLength((long) repLengthBuf1);
-
-    buf1p= buf1.rep.elts();
-    for (ind= 0; ind < repLengthBuf1; ind++)
-      buf1p [ind]= fp [ind + lf];
-    buf1.normalize();
-
-    repLengthBuf1= buf1.rep.length();
-
-    if (deggSubLg >= d - 1)
-      repLengthBuf2= d - 1;
-    else if (deggSubLg < 0)
-      repLengthBuf2= 0;
-    else
-      repLengthBuf2= deggSubLg + 1;
-
-    buf2.rep.SetLength ((long) repLengthBuf2);
-    buf2p= buf2.rep.elts();
-    for (ind= 0; ind < repLengthBuf2; ind++)
-      buf2p [ind]= gp [ind + lg];
-    buf2.normalize();
-
-    repLengthBuf2= buf2.rep.length();
-
-    buf3.rep.SetLength((long) repLengthBuf2 + d);
-    buf3p= buf3.rep.elts();
-    buf2p= buf2.rep.elts();
-    buf1p= buf1.rep.elts();
-    for (ind= 0; ind < repLengthBuf1; ind++)
-      buf3p [ind]= buf1p [ind];
-    for (ind= repLengthBuf1; ind < d; ind++)
-      buf3p [ind]= zzpEZero;
-    for (ind= 0; ind < repLengthBuf2; ind++)
-      buf3p [ind + d]= buf2p [ind];
-    buf3.normalize();
-
-    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
-    i++;
-
-
-    lf= i*d;
-    degfSubLf= degf - lf;
-
-    lg= d*(k-i);
-    deggSubLg= degg - lg;
-
-    buf1p= buf1.rep.elts();
-
-    if (lg >= 0 && deggSubLg > 0)
-    {
-      if (repLengthBuf2 > degfSubLf + 1)
-        degfSubLf= repLengthBuf2 - 1;
-      tmp= tmin (repLengthBuf1, deggSubLg + 1);
-      for (ind= 0; ind < tmp; ind++)
-        gp [ind + lg] -= buf1p [ind];
-    }
-
-    if (lg < 0)
-      break;
-
-    buf2p= buf2.rep.elts();
-    if (degfSubLf >= 0)
-    {
-      for (ind= 0; ind < repLengthBuf2; ind++)
-        fp [ind + lf] -= buf2p [ind];
-    }
-  }
-
-  return result;
-}
-
-CanonicalForm
-reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
-{
-  Variable y= Variable (2);
-  Variable x= Variable (1);
-
-  zz_pX f= F;
-  zz_pX g= G;
-  int degf= deg(f);
-  int degg= deg(g);
-
-  zz_pX buf1;
-  zz_pX buf2;
-  zz_pX buf3;
-
-  zz_p *buf1p;
-  zz_p *buf2p;
-  zz_p *buf3p;
-
-  if (f.rep.length() < (long) d*(k+1)) //zero padding
-    f.rep.SetLength ((long)d*(k+1));
-
-  zz_p *gp= g.rep.elts();
-  zz_p *fp= f.rep.elts();
-  CanonicalForm result= 0;
-  int i= 0;
-  int lf= 0;
-  int lg= d*k;
-  int degfSubLf= degf;
-  int deggSubLg= degg-lg;
-  int repLengthBuf2, repLengthBuf1, ind, tmp;
-  zz_p zzpZero= zz_p();
-  while (degf >= lf || lg >= 0)
-  {
-    if (degfSubLf >= d)
-      repLengthBuf1= d;
-    else if (degfSubLf < 0)
-      repLengthBuf1= 0;
-    else
-      repLengthBuf1= degfSubLf + 1;
-    buf1.rep.SetLength((long) repLengthBuf1);
-
-    buf1p= buf1.rep.elts();
-    for (ind= 0; ind < repLengthBuf1; ind++)
-      buf1p [ind]= fp [ind + lf];
-    buf1.normalize();
-
-    repLengthBuf1= buf1.rep.length();
-
-    if (deggSubLg >= d - 1)
-      repLengthBuf2= d - 1;
-    else if (deggSubLg < 0)
-      repLengthBuf2= 0;
-    else
-      repLengthBuf2= deggSubLg + 1;
-
-    buf2.rep.SetLength ((long) repLengthBuf2);
-    buf2p= buf2.rep.elts();
-    for (ind= 0; ind < repLengthBuf2; ind++)
-      buf2p [ind]= gp [ind + lg];
-
-    buf2.normalize();
-
-    repLengthBuf2= buf2.rep.length();
-
-
-    buf3.rep.SetLength((long) repLengthBuf2 + d);
-    buf3p= buf3.rep.elts();
-    buf2p= buf2.rep.elts();
-    buf1p= buf1.rep.elts();
-    for (ind= 0; ind < repLengthBuf1; ind++)
-      buf3p [ind]= buf1p [ind];
-    for (ind= repLengthBuf1; ind < d; ind++)
-      buf3p [ind]= zzpZero;
-    for (ind= 0; ind < repLengthBuf2; ind++)
-      buf3p [ind + d]= buf2p [ind];
-    buf3.normalize();
-
-    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
-    i++;
-
-
-    lf= i*d;
-    degfSubLf= degf - lf;
-
-    lg= d*(k-i);
-    deggSubLg= degg - lg;
-
-    buf1p= buf1.rep.elts();
-
-    if (lg >= 0 && deggSubLg > 0)
-    {
-      if (repLengthBuf2 > degfSubLf + 1)
-        degfSubLf= repLengthBuf2 - 1;
-      tmp= tmin (repLengthBuf1, deggSubLg + 1);
-      for (ind= 0; ind < tmp; ind++)
-        gp [ind + lg] -= buf1p [ind];
-    }
-    if (lg < 0)
-      break;
-
-    buf2p= buf2.rep.elts();
-    if (degfSubLf >= 0)
-    {
-      for (ind= 0; ind < repLengthBuf2; ind++)
-        fp [ind + lf] -= buf2p [ind];
-    }
-  }
-
-  return result;
-}
-
-CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
-{
-  Variable y= Variable (2);
-  Variable x= Variable (1);
-
-  zz_pEX f= F;
-  zz_pE *fp= f.rep.elts();
-
-  zz_pEX buf;
-  zz_pE *bufp;
-  CanonicalForm result= 0;
-  int i= 0;
-  int degf= deg(f);
-  int k= 0;
-  int degfSubK, repLength, j;
-  while (degf >= k)
-  {
-    degfSubK= degf - k;
-    if (degfSubK >= d)
-      repLength= d;
-    else
-      repLength= degfSubK + 1;
-
-    buf.rep.SetLength ((long) repLength);
-    bufp= buf.rep.elts();
-    for (j= 0; j < repLength; j++)
-      bufp [j]= fp [j + k];
-    buf.normalize();
-
-    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
-    i++;
-    k= d*i;
-  }
-
-  return result;
-}
-
-CanonicalForm reverseSubstFp (const zz_pX& F, int d)
-{
-  Variable y= Variable (2);
-  Variable x= Variable (1);
-
-  zz_pX f= F;
-  zz_p *fp= f.rep.elts();
-
-  zz_pX buf;
-  zz_p *bufp;
-  CanonicalForm result= 0;
-  int i= 0;
-  int degf= deg(f);
-  int k= 0;
-  int degfSubK, repLength, j;
-  while (degf >= k)
-  {
-    degfSubK= degf - k;
-    if (degfSubK >= d)
-      repLength= d;
-    else
-      repLength= degfSubK + 1;
-
-    buf.rep.SetLength ((long) repLength);
-    bufp= buf.rep.elts();
-    for (j= 0; j < repLength; j++)
-      bufp [j]= fp [j + k];
-    buf.normalize();
-
-    result += convertNTLzzpX2CF (buf, x)*power (y, i);
-    i++;
-    k= d*i;
-  }
-
-  return result;
-}
-
-// assumes input to be reduced mod M and to be an element of Fq not Fp
-CanonicalForm
-mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
-                  CanonicalForm& M)
-{
-  int d1= degree (F, 1) + degree (G, 1) + 1;
-  d1 /= 2;
-  d1 += 1;
-
-  zz_pX F1, F2;
-  kronSubRecipro (F1, F2, F, d1);
-  zz_pX G1, G2;
-  kronSubRecipro (G1, G2, G, d1);
-
-  int k= d1*degree (M);
-  MulTrunc (F1, F1, G1, (long) k);
-
-  int degtailF= degree (tailcoeff (F), 1);
-  int degtailG= degree (tailcoeff (G), 1);
-  int taildegF= taildegree (F);
-  int taildegG= taildegree (G);
-  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
-
-  reverse (F2, F2);
-  reverse (G2, G2);
-  MulTrunc (F2, F2, G2, b + 1);
-  reverse (F2, F2, b);
-
-  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
-  return reverseSubst (F1, F2, d1, d2);
-}
-
-//Kronecker substitution
-CanonicalForm
-mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
-              CanonicalForm& M)
-{
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-
-  int degAx= degree (A, 1);
-  int degAy= degree (A, 2);
-  int degBx= degree (B, 1);
-  int degBy= degree (B, 2);
-  int d1= degAx + 1 + degBx;
-  int d2= tmax (degAy, degBy);
-
-  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
-    return mulMod2NTLFpReci (A, B, M);
-
-  zz_pX NTLA= kronSubFp (A, d1);
-  zz_pX NTLB= kronSubFp (B, d1);
-
-  int k= d1*degree (M);
-  MulTrunc (NTLA, NTLA, NTLB, (long) k);
-
-  A= reverseSubstFp (NTLA, d1);
-
-  return A;
-}
-
-// assumes input to be reduced mod M and to be an element of Fq not Fp
-CanonicalForm
-mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
-                  CanonicalForm& M, const Variable& alpha)
-{
-  int d1= degree (F, 1) + degree (G, 1) + 1;
-  d1 /= 2;
-  d1 += 1;
-
-  zz_pEX F1, F2;
-  kronSubRecipro (F1, F2, F, d1, alpha);
-  zz_pEX G1, G2;
-  kronSubRecipro (G1, G2, G, d1, alpha);
-
-  int k= d1*degree (M);
-  MulTrunc (F1, F1, G1, (long) k);
-
-  int degtailF= degree (tailcoeff (F), 1);
-  int degtailG= degree (tailcoeff (G), 1);
-  int taildegF= taildegree (F);
-  int taildegG= taildegree (G);
-  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
-
-  reverse (F2, F2);
-  reverse (G2, G2);
-  MulTrunc (F2, F2, G2, b + 1);
-  reverse (F2, F2, b);
-
-  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
-  return reverseSubst (F1, F2, d1, d2, alpha);
-}
-
-#ifdef HAVE_FLINT
-CanonicalForm
-mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
-                CanonicalForm& M);
-#endif
-
-CanonicalForm
-mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
-              CanonicalForm& M)
-{
-  Variable alpha;
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-
-  if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
-  {
-    int degAx= degree (A, 1);
-    int degAy= degree (A, 2);
-    int degBx= degree (B, 1);
-    int degBy= degree (B, 2);
-    int d1= degAx + degBx + 1;
-    int d2= tmax (degAy, degBy);
-    zz_p::init (getCharacteristic());
-    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
-    zz_pE::init (NTLMipo);
-
-    int degMipo= degree (getMipo (alpha));
-    if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
-        (2*degAy > degree (M)))
-      return mulMod2NTLFqReci (A, B, M, alpha);
-
-    zz_pEX NTLA= kronSub (A, d1, alpha);
-    zz_pEX NTLB= kronSub (B, d1, alpha);
-
-    int k= d1*degree (M);
-
-    MulTrunc (NTLA, NTLA, NTLB, (long) k);
-
-    A= reverseSubst (NTLA, d1, alpha);
-
-    return A;
-  }
-  else
-#ifdef HAVE_FLINT
-    return mulMod2FLINTFp (A, B, M);
-#else
-    return mulMod2NTLFp (A, B, M);
-#endif
-}
-
-#ifdef HAVE_FLINT
-void
-kronSubRecipro (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
-                int d)
+kronSubReciproFp (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
+                  int d)
 {
   int degAy= degree (A);
@@ -1291,6 +685,6 @@
 
 void
-kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
-                int d)
+kronSubReciproQ (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
+                 int d)
 {
   int degAy= degree (A);
@@ -1372,5 +766,5 @@
 
 CanonicalForm
-reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
+reverseSubstReciproFp (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
 {
   Variable y= Variable (2);
@@ -1490,5 +884,5 @@
 
 CanonicalForm
-reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
+reverseSubstReciproQ (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
 {
   Variable y= Variable (2);
@@ -1674,10 +1068,10 @@
   nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
   nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
-  kronSubRecipro (F1, F2, F, d1);
+  kronSubReciproFp (F1, F2, F, d1);
 
   nmod_poly_t G1, G2;
   nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
   nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
-  kronSubRecipro (G1, G2, G, d1);
+  kronSubReciproFp (G1, G2, G, d1);
 
   int k= d1*degree (M);
@@ -1696,5 +1090,5 @@
 
 
-  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
+  CanonicalForm result= reverseSubstReciproFp (F1, F2, d1, d2);
 
   nmod_poly_clear (F1);
@@ -1750,10 +1144,10 @@
   fmpz_poly_init (F1);
   fmpz_poly_init (F2);
-  kronSubRecipro (F1, F2, F, d1);
+  kronSubReciproQ (F1, F2, F, d1);
 
   fmpz_poly_t G1, G2;
   fmpz_poly_init (G1);
   fmpz_poly_init (G2);
-  kronSubRecipro (G1, G2, G, d1);
+  kronSubReciproQ (G1, G2, G, d1);
 
   int k= d1*degree (M);
@@ -1771,5 +1165,5 @@
   int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
 
-  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
+  CanonicalForm result= reverseSubstReciproQ (F1, F2, d1, d2);
 
   fmpz_poly_clear (F1);
@@ -1811,4 +1205,587 @@
 
 #endif
+
+zz_pX kronSubFp (const CanonicalForm& A, int d)
+{
+  int degAy= degree (A);
+  zz_pX result;
+  result.rep.SetLength (d*(degAy + 1));
+
+  zz_p *resultp;
+  resultp= result.rep.elts();
+  zz_pX buf;
+  zz_p *bufp;
+  int j, k, bufRepLength;
+
+  for (CFIterator i= A; i.hasTerms(); i++)
+  {
+    if (i.coeff().inCoeffDomain())
+      buf= convertFacCF2NTLzzpX (i.coeff());
+    else
+      buf= convertFacCF2NTLzzpX (i.coeff());
+
+    k= i.exp()*d;
+    bufp= buf.rep.elts();
+    bufRepLength= (int) buf.rep.length();
+    for (j= 0; j < bufRepLength; j++)
+      resultp [j + k]= bufp [j];
+  }
+  result.normalize();
+
+  return result;
+}
+
+zz_pEX kronSubFq (const CanonicalForm& A, int d, const Variable& alpha)
+{
+  int degAy= degree (A);
+  zz_pEX result;
+  result.rep.SetLength (d*(degAy + 1));
+
+  Variable v;
+  zz_pE *resultp;
+  resultp= result.rep.elts();
+  zz_pEX buf1;
+  zz_pE *buf1p;
+  zz_pX buf2;
+  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+  int j, k, buf1RepLength;
+
+  for (CFIterator i= A; i.hasTerms(); i++)
+  {
+    if (i.coeff().inCoeffDomain())
+    {
+      buf2= convertFacCF2NTLzzpX (i.coeff());
+      buf1= to_zz_pEX (to_zz_pE (buf2));
+    }
+    else
+      buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
+
+    k= i.exp()*d;
+    buf1p= buf1.rep.elts();
+    buf1RepLength= (int) buf1.rep.length();
+    for (j= 0; j < buf1RepLength; j++)
+      resultp [j + k]= buf1p [j];
+  }
+  result.normalize();
+
+  return result;
+}
+
+void
+kronSubReciproFq (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
+                  const Variable& alpha)
+{
+  int degAy= degree (A);
+  subA1.rep.SetLength ((long) d*(degAy + 2));
+  subA2.rep.SetLength ((long) d*(degAy + 2));
+
+  Variable v;
+  zz_pE *subA1p;
+  zz_pE *subA2p;
+  subA1p= subA1.rep.elts();
+  subA2p= subA2.rep.elts();
+  zz_pEX buf;
+  zz_pE *bufp;
+  zz_pX buf2;
+  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+  int j, k, kk, bufRepLength;
+
+  for (CFIterator i= A; i.hasTerms(); i++)
+  {
+    if (i.coeff().inCoeffDomain())
+    {
+      buf2= convertFacCF2NTLzzpX (i.coeff());
+      buf= to_zz_pEX (to_zz_pE (buf2));
+    }
+    else
+      buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
+
+    k= i.exp()*d;
+    kk= (degAy - i.exp())*d;
+    bufp= buf.rep.elts();
+    bufRepLength= (int) buf.rep.length();
+    for (j= 0; j < bufRepLength; j++)
+    {
+      subA1p [j + k] += bufp [j];
+      subA2p [j + kk] += bufp [j];
+    }
+  }
+  subA1.normalize();
+  subA2.normalize();
+}
+
+void
+kronSubReciproFp (zz_pX& subA1, zz_pX& subA2, const CanonicalForm& A, int d)
+{
+  int degAy= degree (A);
+  subA1.rep.SetLength ((long) d*(degAy + 2));
+  subA2.rep.SetLength ((long) d*(degAy + 2));
+
+  zz_p *subA1p;
+  zz_p *subA2p;
+  subA1p= subA1.rep.elts();
+  subA2p= subA2.rep.elts();
+  zz_pX buf;
+  zz_p *bufp;
+  int j, k, kk, bufRepLength;
+
+  for (CFIterator i= A; i.hasTerms(); i++)
+  {
+    buf= convertFacCF2NTLzzpX (i.coeff());
+
+    k= i.exp()*d;
+    kk= (degAy - i.exp())*d;
+    bufp= buf.rep.elts();
+    bufRepLength= (int) buf.rep.length();
+    for (j= 0; j < bufRepLength; j++)
+    {
+      subA1p [j + k] += bufp [j];
+      subA2p [j + kk] += bufp [j];
+    }
+  }
+  subA1.normalize();
+  subA2.normalize();
+}
+
+CanonicalForm
+reverseSubstReciproFq (const zz_pEX& F, const zz_pEX& G, int d, int k,
+                       const Variable& alpha)
+{
+  Variable y= Variable (2);
+  Variable x= Variable (1);
+
+  zz_pEX f= F;
+  zz_pEX g= G;
+  int degf= deg(f);
+  int degg= deg(g);
+
+  zz_pEX buf1;
+  zz_pEX buf2;
+  zz_pEX buf3;
+
+  zz_pE *buf1p;
+  zz_pE *buf2p;
+  zz_pE *buf3p;
+  if (f.rep.length() < (long) d*(k+1)) //zero padding
+    f.rep.SetLength ((long)d*(k+1));
+
+  zz_pE *gp= g.rep.elts();
+  zz_pE *fp= f.rep.elts();
+  CanonicalForm result= 0;
+  int i= 0;
+  int lf= 0;
+  int lg= d*k;
+  int degfSubLf= degf;
+  int deggSubLg= degg-lg;
+  int repLengthBuf2, repLengthBuf1, ind, tmp;
+  zz_pE zzpEZero= zz_pE();
+
+  while (degf >= lf || lg >= 0)
+  {
+    if (degfSubLf >= d)
+      repLengthBuf1= d;
+    else if (degfSubLf < 0)
+      repLengthBuf1= 0;
+    else
+      repLengthBuf1= degfSubLf + 1;
+    buf1.rep.SetLength((long) repLengthBuf1);
+
+    buf1p= buf1.rep.elts();
+    for (ind= 0; ind < repLengthBuf1; ind++)
+      buf1p [ind]= fp [ind + lf];
+    buf1.normalize();
+
+    repLengthBuf1= buf1.rep.length();
+
+    if (deggSubLg >= d - 1)
+      repLengthBuf2= d - 1;
+    else if (deggSubLg < 0)
+      repLengthBuf2= 0;
+    else
+      repLengthBuf2= deggSubLg + 1;
+
+    buf2.rep.SetLength ((long) repLengthBuf2);
+    buf2p= buf2.rep.elts();
+    for (ind= 0; ind < repLengthBuf2; ind++)
+      buf2p [ind]= gp [ind + lg];
+    buf2.normalize();
+
+    repLengthBuf2= buf2.rep.length();
+
+    buf3.rep.SetLength((long) repLengthBuf2 + d);
+    buf3p= buf3.rep.elts();
+    buf2p= buf2.rep.elts();
+    buf1p= buf1.rep.elts();
+    for (ind= 0; ind < repLengthBuf1; ind++)
+      buf3p [ind]= buf1p [ind];
+    for (ind= repLengthBuf1; ind < d; ind++)
+      buf3p [ind]= zzpEZero;
+    for (ind= 0; ind < repLengthBuf2; ind++)
+      buf3p [ind + d]= buf2p [ind];
+    buf3.normalize();
+
+    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
+    i++;
+
+
+    lf= i*d;
+    degfSubLf= degf - lf;
+
+    lg= d*(k-i);
+    deggSubLg= degg - lg;
+
+    buf1p= buf1.rep.elts();
+
+    if (lg >= 0 && deggSubLg > 0)
+    {
+      if (repLengthBuf2 > degfSubLf + 1)
+        degfSubLf= repLengthBuf2 - 1;
+      tmp= tmin (repLengthBuf1, deggSubLg + 1);
+      for (ind= 0; ind < tmp; ind++)
+        gp [ind + lg] -= buf1p [ind];
+    }
+
+    if (lg < 0)
+      break;
+
+    buf2p= buf2.rep.elts();
+    if (degfSubLf >= 0)
+    {
+      for (ind= 0; ind < repLengthBuf2; ind++)
+        fp [ind + lf] -= buf2p [ind];
+    }
+  }
+
+  return result;
+}
+
+CanonicalForm
+reverseSubstReciproFp (const zz_pX& F, const zz_pX& G, int d, int k)
+{
+  Variable y= Variable (2);
+  Variable x= Variable (1);
+
+  zz_pX f= F;
+  zz_pX g= G;
+  int degf= deg(f);
+  int degg= deg(g);
+
+  zz_pX buf1;
+  zz_pX buf2;
+  zz_pX buf3;
+
+  zz_p *buf1p;
+  zz_p *buf2p;
+  zz_p *buf3p;
+
+  if (f.rep.length() < (long) d*(k+1)) //zero padding
+    f.rep.SetLength ((long)d*(k+1));
+
+  zz_p *gp= g.rep.elts();
+  zz_p *fp= f.rep.elts();
+  CanonicalForm result= 0;
+  int i= 0;
+  int lf= 0;
+  int lg= d*k;
+  int degfSubLf= degf;
+  int deggSubLg= degg-lg;
+  int repLengthBuf2, repLengthBuf1, ind, tmp;
+  zz_p zzpZero= zz_p();
+  while (degf >= lf || lg >= 0)
+  {
+    if (degfSubLf >= d)
+      repLengthBuf1= d;
+    else if (degfSubLf < 0)
+      repLengthBuf1= 0;
+    else
+      repLengthBuf1= degfSubLf + 1;
+    buf1.rep.SetLength((long) repLengthBuf1);
+
+    buf1p= buf1.rep.elts();
+    for (ind= 0; ind < repLengthBuf1; ind++)
+      buf1p [ind]= fp [ind + lf];
+    buf1.normalize();
+
+    repLengthBuf1= buf1.rep.length();
+
+    if (deggSubLg >= d - 1)
+      repLengthBuf2= d - 1;
+    else if (deggSubLg < 0)
+      repLengthBuf2= 0;
+    else
+      repLengthBuf2= deggSubLg + 1;
+
+    buf2.rep.SetLength ((long) repLengthBuf2);
+    buf2p= buf2.rep.elts();
+    for (ind= 0; ind < repLengthBuf2; ind++)
+      buf2p [ind]= gp [ind + lg];
+
+    buf2.normalize();
+
+    repLengthBuf2= buf2.rep.length();
+
+
+    buf3.rep.SetLength((long) repLengthBuf2 + d);
+    buf3p= buf3.rep.elts();
+    buf2p= buf2.rep.elts();
+    buf1p= buf1.rep.elts();
+    for (ind= 0; ind < repLengthBuf1; ind++)
+      buf3p [ind]= buf1p [ind];
+    for (ind= repLengthBuf1; ind < d; ind++)
+      buf3p [ind]= zzpZero;
+    for (ind= 0; ind < repLengthBuf2; ind++)
+      buf3p [ind + d]= buf2p [ind];
+    buf3.normalize();
+
+    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
+    i++;
+
+
+    lf= i*d;
+    degfSubLf= degf - lf;
+
+    lg= d*(k-i);
+    deggSubLg= degg - lg;
+
+    buf1p= buf1.rep.elts();
+
+    if (lg >= 0 && deggSubLg > 0)
+    {
+      if (repLengthBuf2 > degfSubLf + 1)
+        degfSubLf= repLengthBuf2 - 1;
+      tmp= tmin (repLengthBuf1, deggSubLg + 1);
+      for (ind= 0; ind < tmp; ind++)
+        gp [ind + lg] -= buf1p [ind];
+    }
+    if (lg < 0)
+      break;
+
+    buf2p= buf2.rep.elts();
+    if (degfSubLf >= 0)
+    {
+      for (ind= 0; ind < repLengthBuf2; ind++)
+        fp [ind + lf] -= buf2p [ind];
+    }
+  }
+
+  return result;
+}
+
+CanonicalForm reverseSubstFq (const zz_pEX& F, int d, const Variable& alpha)
+{
+  Variable y= Variable (2);
+  Variable x= Variable (1);
+
+  zz_pEX f= F;
+  zz_pE *fp= f.rep.elts();
+
+  zz_pEX buf;
+  zz_pE *bufp;
+  CanonicalForm result= 0;
+  int i= 0;
+  int degf= deg(f);
+  int k= 0;
+  int degfSubK, repLength, j;
+  while (degf >= k)
+  {
+    degfSubK= degf - k;
+    if (degfSubK >= d)
+      repLength= d;
+    else
+      repLength= degfSubK + 1;
+
+    buf.rep.SetLength ((long) repLength);
+    bufp= buf.rep.elts();
+    for (j= 0; j < repLength; j++)
+      bufp [j]= fp [j + k];
+    buf.normalize();
+
+    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
+    i++;
+    k= d*i;
+  }
+
+  return result;
+}
+
+CanonicalForm reverseSubstFp (const zz_pX& F, int d)
+{
+  Variable y= Variable (2);
+  Variable x= Variable (1);
+
+  zz_pX f= F;
+  zz_p *fp= f.rep.elts();
+
+  zz_pX buf;
+  zz_p *bufp;
+  CanonicalForm result= 0;
+  int i= 0;
+  int degf= deg(f);
+  int k= 0;
+  int degfSubK, repLength, j;
+  while (degf >= k)
+  {
+    degfSubK= degf - k;
+    if (degfSubK >= d)
+      repLength= d;
+    else
+      repLength= degfSubK + 1;
+
+    buf.rep.SetLength ((long) repLength);
+    bufp= buf.rep.elts();
+    for (j= 0; j < repLength; j++)
+      bufp [j]= fp [j + k];
+    buf.normalize();
+
+    result += convertNTLzzpX2CF (buf, x)*power (y, i);
+    i++;
+    k= d*i;
+  }
+
+  return result;
+}
+
+// assumes input to be reduced mod M and to be an element of Fq not Fp
+CanonicalForm
+mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
+                  CanonicalForm& M)
+{
+  int d1= degree (F, 1) + degree (G, 1) + 1;
+  d1 /= 2;
+  d1 += 1;
+
+  zz_pX F1, F2;
+  kronSubReciproFp (F1, F2, F, d1);
+  zz_pX G1, G2;
+  kronSubReciproFp (G1, G2, G, d1);
+
+  int k= d1*degree (M);
+  MulTrunc (F1, F1, G1, (long) k);
+
+  int degtailF= degree (tailcoeff (F), 1);
+  int degtailG= degree (tailcoeff (G), 1);
+  int taildegF= taildegree (F);
+  int taildegG= taildegree (G);
+  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
+
+  reverse (F2, F2);
+  reverse (G2, G2);
+  MulTrunc (F2, F2, G2, b + 1);
+  reverse (F2, F2, b);
+
+  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
+  return reverseSubstReciproFp (F1, F2, d1, d2);
+}
+
+//Kronecker substitution
+CanonicalForm
+mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
+              CanonicalForm& M)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+
+  int degAx= degree (A, 1);
+  int degAy= degree (A, 2);
+  int degBx= degree (B, 1);
+  int degBy= degree (B, 2);
+  int d1= degAx + 1 + degBx;
+  int d2= tmax (degAy, degBy);
+
+  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
+    return mulMod2NTLFpReci (A, B, M);
+
+  zz_pX NTLA= kronSubFp (A, d1);
+  zz_pX NTLB= kronSubFp (B, d1);
+
+  int k= d1*degree (M);
+  MulTrunc (NTLA, NTLA, NTLB, (long) k);
+
+  A= reverseSubstFp (NTLA, d1);
+
+  return A;
+}
+
+// assumes input to be reduced mod M and to be an element of Fq not Fp
+CanonicalForm
+mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
+                  CanonicalForm& M, const Variable& alpha)
+{
+  int d1= degree (F, 1) + degree (G, 1) + 1;
+  d1 /= 2;
+  d1 += 1;
+
+  zz_pEX F1, F2;
+  kronSubReciproFq (F1, F2, F, d1, alpha);
+  zz_pEX G1, G2;
+  kronSubReciproFq (G1, G2, G, d1, alpha);
+
+  int k= d1*degree (M);
+  MulTrunc (F1, F1, G1, (long) k);
+
+  int degtailF= degree (tailcoeff (F), 1);
+  int degtailG= degree (tailcoeff (G), 1);
+  int taildegF= taildegree (F);
+  int taildegG= taildegree (G);
+  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
+
+  reverse (F2, F2);
+  reverse (G2, G2);
+  MulTrunc (F2, F2, G2, b + 1);
+  reverse (F2, F2, b);
+
+  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
+  return reverseSubstReciproFq (F1, F2, d1, d2, alpha);
+}
+
+#ifdef HAVE_FLINT
+CanonicalForm
+mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
+                CanonicalForm& M);
+#endif
+
+CanonicalForm
+mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
+              CanonicalForm& M)
+{
+  Variable alpha;
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+
+  if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
+  {
+    int degAx= degree (A, 1);
+    int degAy= degree (A, 2);
+    int degBx= degree (B, 1);
+    int degBy= degree (B, 2);
+    int d1= degAx + degBx + 1;
+    int d2= tmax (degAy, degBy);
+    zz_p::init (getCharacteristic());
+    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+    zz_pE::init (NTLMipo);
+
+    int degMipo= degree (getMipo (alpha));
+    if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
+        (2*degAy > degree (M)))
+      return mulMod2NTLFqReci (A, B, M, alpha);
+
+    zz_pEX NTLA= kronSubFq (A, d1, alpha);
+    zz_pEX NTLB= kronSubFq (B, d1, alpha);
+
+    int k= d1*degree (M);
+
+    MulTrunc (NTLA, NTLA, NTLB, (long) k);
+
+    A= reverseSubstFq (NTLA, d1, alpha);
+
+    return A;
+  }
+  else
+#ifdef HAVE_FLINT
+    return mulMod2FLINTFp (A, B, M);
+#else
+    return mulMod2NTLFp (A, B, M);
+#endif
+}
 
 CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
@@ -1886,4 +1863,8 @@
   DEBOUTLN (cerr, "fatal end in mulMod2");
 }
+
+// end bivariate polys
+//**********************
+// multivariate polys
 
 CanonicalForm mod (const CanonicalForm& F, const CFList& M)
@@ -2031,4 +2012,8 @@
 }
 
+// end multivariate polys
+//***************************
+// division
+
 CanonicalForm reverse (const CanonicalForm& F, int d)
 {
@@ -2510,3 +2495,5 @@
 }
 
+// end division
+
 #endif
