[SCM] an open source computer algebra system branch, cleanedupstream, updated. 6125e540ca6d66c307958938a9d53b245507c323

Bernhard R. Link brlink at debian.org
Tue Apr 24 15:54:47 UTC 2012


The following commit has been merged in the cleanedupstream branch:
commit 5b64ae6261d338d3ec706cdf793dd3def49ea6ce
Author: Martin Lee <martinlee84 at web.de>
Date:   Thu Feb 9 16:47:39 2012 +0100

    chg: separated multiplication and Hensel lifting functions

diff --git a/factory/GNUmakefile.in b/factory/GNUmakefile.in
index 75ba67e..da2c7c4 100644
--- a/factory/GNUmakefile.in
+++ b/factory/GNUmakefile.in
@@ -180,6 +180,7 @@ basefactorysrc := \
 		facFqSquarefree.cc \
 		facHensel.cc \
 		facIrredTest.cc \
+		facMul.cc \
 		facSparseHensel.cc \
 		fieldGCD.cc \
 		ffops.cc \
@@ -268,6 +269,7 @@ basefactoryincl := \
 		facFqSquarefree.h \
 		facHensel.h \
 		facIrredTest.h \
+		facMul.h \
 		facSparseHensel.h \
 		fieldGCD.h \
 		ffops.h \
diff --git a/factory/facBivar.cc b/factory/facBivar.cc
index eedf9c7..f049c91 100644
--- a/factory/facBivar.cc
+++ b/factory/facBivar.cc
@@ -19,6 +19,7 @@
 #include "facFqBivar.h"
 #include "facBivar.h"
 #include "facHensel.h"
+#include "facMul.h"
 
 #ifdef HAVE_NTL
 TIMING_DEFINE_PRINT(fac_uni_factorizer)
diff --git a/factory/facFqBivar.cc b/factory/facFqBivar.cc
index ce3a055..543c475 100644
--- a/factory/facFqBivar.cc
+++ b/factory/facFqBivar.cc
@@ -27,6 +27,7 @@
 #include "cf_map_ext.h"
 #include "cf_random.h"
 #include "facHensel.h"
+#include "facMul.h"
 #include "cf_map.h"
 #include "cf_gcd_smallp.h"
 #include "facFqBivarUtil.h"
diff --git a/factory/facFqBivarUtil.cc b/factory/facFqBivarUtil.cc
index 11ee282..9f48bd4 100644
--- a/factory/facFqBivarUtil.cc
+++ b/factory/facFqBivarUtil.cc
@@ -24,6 +24,7 @@
 #include "facFqBivarUtil.h"
 #include "cfNewtonPolygon.h"
 #include "facHensel.h"
+#include "facMul.h"
 
 
 void append (CFList& factors1, const CFList& factors2)
diff --git a/factory/facFqFactorize.cc b/factory/facFqFactorize.cc
index 807fa6b..e6351d0 100644
--- a/factory/facFqFactorize.cc
+++ b/factory/facFqFactorize.cc
@@ -29,6 +29,7 @@
 #include "cf_map_ext.h"
 #include "algext.h"
 #include "facSparseHensel.h"
+#include "facMul.h"
 
 #ifdef HAVE_NTL
 #include <NTL/ZZ_pEX.h>
diff --git a/factory/facHensel.cc b/factory/facHensel.cc
index 12509b9..bde2c1c 100644
--- a/factory/facHensel.cc
+++ b/factory/facHensel.cc
@@ -3,8 +3,7 @@
 \*****************************************************************************/
 /** @file facHensel.cc
  *
- * This file implements functions to lift factors via Hensel lifting and
- * functions for modular multiplication and division with remainder.
+ * This file implements functions to lift factors via Hensel lifting.
  *
  * ABSTRACT: Hensel lifting is described in "Efficient Multivariate
  * Factorization over Finite Fields" by L. Bernardin & M. Monagon. Division with
@@ -22,6 +21,7 @@
 #include "timing.h"
 
 #include "facHensel.h"
+#include "facMul.h"
 #include "cf_util.h"
 #include "fac_util.h"
 #include "cf_algorithm.h"
@@ -31,373 +31,6 @@
 #include <NTL/lzz_pEX.h>
 #include "NTLconvert.h"
 
-#ifdef HAVE_FLINT
-#include "FLINTconvert.h"
-#endif
-
-#ifdef HAVE_FLINT
-void kronSub (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().inBaseDomain())
-      convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d), i.coeff());
-    else
-      for (j= i.coeff(); j.hasTerms(); j++)
-        convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d+j.exp()),
-                        j.coeff());
-  }
-  _fmpz_poly_normalise(result);
-}
-
-
-CanonicalForm
-reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha,
-              const CanonicalForm& den)
-{
-  Variable x= Variable (1);
-
-  CanonicalForm result= 0;
-  int i= 0;
-  int degf= fmpz_poly_degree (F);
-  int k= 0;
-  int degfSubK;
-  int repLength, j;
-  CanonicalForm coeff;
-  fmpz* tmp;
-  while (degf >= k)
-  {
-    coeff= 0;
-    degfSubK= degf - k;
-    if (degfSubK >= d)
-      repLength= d;
-    else
-      repLength= degfSubK + 1;
-
-    for (j= 0; j < repLength; j++)
-    {
-      tmp= fmpz_poly_get_coeff_ptr (F, j+k);
-      if (!fmpz_is_zero (tmp))
-      {
-        CanonicalForm ff= convertFmpz2CF (tmp)/den;
-        coeff += ff*power (alpha, j);
-      }
-    }
-    result += coeff*power (x, i);
-    i++;
-    k= d*i;
-  }
-  return result;
-}
-
-CanonicalForm
-mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G,
-            const Variable& alpha)
-{
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-
-  CanonicalForm denA= bCommonDen (A);
-  CanonicalForm denB= bCommonDen (B);
-
-  A *= denA;
-  B *= denB;
-  int degAa= degree (A, alpha);
-  int degBa= degree (B, alpha);
-  int d= degAa + 1 + degBa;
-
-  fmpz_poly_t FLINTA,FLINTB;
-  fmpz_poly_init (FLINTA);
-  fmpz_poly_init (FLINTB);
-  kronSub (FLINTA, A, d);
-  kronSub (FLINTB, B, d);
-
-  fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
-
-  denA *= denB;
-  A= reverseSubst (FLINTA, d, alpha, denA);
-
-  fmpz_poly_clear (FLINTA);
-  fmpz_poly_clear (FLINTB);
-  return A;
-}
-
-CanonicalForm
-mulFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
-{
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-
-  CanonicalForm denA= bCommonDen (A);
-  CanonicalForm denB= bCommonDen (B);
-
-  A *= denA;
-  B *= denB;
-  fmpz_poly_t FLINTA,FLINTB;
-  convertFacCF2Fmpz_poly_t (FLINTA, A);
-  convertFacCF2Fmpz_poly_t (FLINTB, B);
-  fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
-  denA *= denB;
-  A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
-  A /= denA;
-  fmpz_poly_clear (FLINTA);
-  fmpz_poly_clear (FLINTB);
-
-  return A;
-}
-
-/*CanonicalForm
-mulFLINTQ2 (const CanonicalForm& F, const CanonicalForm& G)
-{
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-
-  fmpq_poly_t FLINTA,FLINTB;
-  convertFacCF2Fmpq_poly_t (FLINTA, A);
-  convertFacCF2Fmpq_poly_t (FLINTB, B);
-
-  fmpq_poly_mul (FLINTA, FLINTA, FLINTB);
-  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
-  fmpq_poly_clear (FLINTA);
-  fmpq_poly_clear (FLINTB);
-  return A;
-}*/
-
-CanonicalForm
-divFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
-{
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-
-  fmpq_poly_t FLINTA,FLINTB;
-  convertFacCF2Fmpq_poly_t (FLINTA, A);
-  convertFacCF2Fmpq_poly_t (FLINTB, B);
-
-  fmpq_poly_div (FLINTA, FLINTA, FLINTB);
-  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
-
-  fmpq_poly_clear (FLINTA);
-  fmpq_poly_clear (FLINTB);
-  return A;
-}
-
-CanonicalForm
-modFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
-{
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-
-  fmpq_poly_t FLINTA,FLINTB;
-  convertFacCF2Fmpq_poly_t (FLINTA, A);
-  convertFacCF2Fmpq_poly_t (FLINTB, B);
-
-  fmpq_poly_rem (FLINTA, FLINTA, FLINTB);
-  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
-
-  fmpq_poly_clear (FLINTA);
-  fmpq_poly_clear (FLINTB);
-  return A;
-}
-
-CanonicalForm
-mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G,
-                 const Variable& alpha, int m)
-{
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-
-  CanonicalForm denA= bCommonDen (A);
-  CanonicalForm denB= bCommonDen (B);
-
-  A *= denA;
-  B *= denB;
-
-  int degAa= degree (A, alpha);
-  int degBa= degree (B, alpha);
-  int d= degAa + 1 + degBa;
-
-  fmpz_poly_t FLINTA,FLINTB;
-  fmpz_poly_init (FLINTA);
-  fmpz_poly_init (FLINTB);
-  kronSub (FLINTA, A, d);
-  kronSub (FLINTB, B, d);
-
-  int k= d*m;
-  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
-
-  denA *= denB;
-  A= reverseSubst (FLINTA, d, alpha, denA);
-  fmpz_poly_clear (FLINTA);
-  fmpz_poly_clear (FLINTB);
-  return A;
-}
-
-CanonicalForm
-mulFLINTQTrunc (const CanonicalForm& F, const CanonicalForm& G, int m)
-{
-  if (F.inCoeffDomain() || G.inCoeffDomain())
-    return mod (F*G, power (Variable (1), m));
-  Variable alpha;
-  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
-    return mulFLINTQaTrunc (F, G, alpha, m);
-
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-
-  CanonicalForm denA= bCommonDen (A);
-  CanonicalForm denB= bCommonDen (B);
-
-  A *= denA;
-  B *= denB;
-  fmpz_poly_t FLINTA,FLINTB;
-  convertFacCF2Fmpz_poly_t (FLINTA, A);
-  convertFacCF2Fmpz_poly_t (FLINTB, B);
-  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, m);
-  denA *= denB;
-  A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
-  A /= denA;
-  fmpz_poly_clear (FLINTA);
-  fmpz_poly_clear (FLINTB);
-
-  return A;
-}
-
-CanonicalForm uniReverse (const CanonicalForm& F, int d)
-{
-  if (d == 0)
-    return F;
-  if (F.inCoeffDomain())
-    return F*power (Variable (1),d);
-  Variable x= Variable (1);
-  CanonicalForm result= 0;
-  CFIterator i= F;
-  while (d - i.exp() < 0)
-    i++;
-
-  for (; i.hasTerms() && (d - i.exp() >= 0); i++)
-    result += i.coeff()*power (x, d - i.exp());
-  return result;
-}
-
-CanonicalForm
-newtonInverse (const CanonicalForm& F, const int n)
-{
-  int l= ilog2(n);
-
-  CanonicalForm g= F [0];
-
-  ASSERT (!g.isZero(), "expected a unit");
-
-  if (!g.isOne())
-    g = 1/g;
-  Variable x= Variable (1);
-  CanonicalForm result;
-  int exp= 0;
-  if (n & 1)
-  {
-    result= g;
-    exp= 1;
-  }
-  CanonicalForm h;
-
-  for (int i= 1; i <= l; i++)
-  {
-    h= mulNTL (g, mod (F, power (x, (1 << i))));
-    h= mod (h, power (x, (1 << i)) - 1);
-    h= div (h, power (x, (1 << (i - 1))));
-    g -= power (x, (1 << (i - 1)))*
-         mulFLINTQTrunc (g, h, 1 << (i-1));
-
-    if (n & (1 << i))
-    {
-      if (exp)
-      {
-        h= mulNTL (result, mod (F, power (x, exp + (1 << i))));
-        h= mod (h, power (x, exp + (1 << i)) - 1);
-        h= div (h, power (x, exp));
-        result -= power(x, exp)*mulFLINTQTrunc (g, h, 1 << i);
-        exp += (1 << i);
-      }
-      else
-      {
-        exp= (1 << i);
-        result= g;
-      }
-    }
-  }
-
-  return result;
-}
-
-void
-newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
-              CanonicalForm& R)
-{
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-  Variable x= Variable (1);
-  int degA= degree (A, x);
-  int degB= degree (B, x);
-  int m= degA - degB;
-
-  if (m < 0)
-  {
-    R= A;
-    Q= 0;
-    return;
-  }
-
-  if (degB <= 1)
-    divrem (A, B, Q, R);
-  else
-  {
-    R= uniReverse (A, degA);
-
-    CanonicalForm revB= uniReverse (B, degB);
-    CanonicalForm buf= revB;
-    revB= newtonInverse (revB, m + 1);
-    Q= mulFLINTQTrunc (R, revB, m + 1);
-    Q= uniReverse (Q, m);
-
-    R= A - mulNTL (Q, B);
-  }
-}
-
-void
-newtonDiv (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q)
-{
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-  Variable x= Variable (1);
-  int degA= degree (A, x);
-  int degB= degree (B, x);
-  int m= degA - degB;
-
-  if (m < 0)
-  {
-    Q= 0;
-    return;
-  }
-
-  if (degB <= 1)
-    Q= div (A, B);
-  else
-  {
-    CanonicalForm R= uniReverse (A, degA);
-
-    CanonicalForm revB= uniReverse (B, degB);
-    revB= newtonInverse (revB, m + 1);
-    Q= mulFLINTQTrunc (R, revB, m + 1);
-    Q= uniReverse (Q, m);
-  }
-}
-
-#endif
-
 static
 CFList productsNTL (const CFList& factors, const CanonicalForm& M)
 {
@@ -694,2184 +327,6 @@ modularDiophant (const CanonicalForm& f, const CFList& factors,
   return result;
 }
 
-CanonicalForm
-mulNTL (const CanonicalForm& F, const CanonicalForm& G)
-{
-  if (F.inCoeffDomain() || G.inCoeffDomain() || getCharacteristic() == 0)
-  {
-    Variable alpha;
-#ifdef HAVE_FLINT
-    if ((!F.inCoeffDomain() && !G.inCoeffDomain()) &&
-        (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha)))
-    {
-      CanonicalForm result= mulFLINTQa (F, G, alpha);
-      return result;
-    }
-    else if (!F.inCoeffDomain() && !G.inCoeffDomain())
-      return mulFLINTQ (F, G);
-#endif
-    return F*G;
-  }
-  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
-  ASSERT (F.level() == G.level(), "expected polys of same level");
-  if (CFFactory::gettype() == GaloisFieldDomain)
-    return F*G;
-  zz_p::init (getCharacteristic());
-  Variable alpha;
-  CanonicalForm result;
-  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
-  {
-    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
-    zz_pE::init (NTLMipo);
-    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
-    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
-    mul (NTLF, NTLF, NTLG);
-    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
-  }
-  else
-  {
-#ifdef HAVE_FLINT
-    nmod_poly_t FLINTF, FLINTG;
-    convertFacCF2nmod_poly_t (FLINTF, F);
-    convertFacCF2nmod_poly_t (FLINTG, G);
-    nmod_poly_mul (FLINTF, FLINTF, FLINTG);
-    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
-    nmod_poly_clear (FLINTF);
-    nmod_poly_clear (FLINTG);
-#else
-    zz_pX NTLF= convertFacCF2NTLzzpX (F);
-    zz_pX NTLG= convertFacCF2NTLzzpX (G);
-    mul (NTLF, NTLF, NTLG);
-    result= convertNTLzzpX2CF(NTLF, F.mvar());
-#endif
-  }
-  return result;
-}
-
-CanonicalForm
-modNTL (const CanonicalForm& F, const CanonicalForm& G)
-{
-  if (F.inCoeffDomain() && G.isUnivariate())
-    return F;
-  else if (F.inCoeffDomain() && G.inCoeffDomain())
-    return mod (F, G);
-  else if (F.isUnivariate() && G.inCoeffDomain())
-    return mod (F,G);
-
-  if (getCharacteristic() == 0)
-  {
-#ifdef HAVE_FLINT
-    Variable alpha;
-    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
-      return modFLINTQ (F, G);
-    else
-    {
-      CanonicalForm Q, R;
-      newtonDivrem (F, G, Q, R);
-      return R;
-    }
-#else
-    return mod (F, G);
-#endif
-  }
-
-  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
-  ASSERT (F.level() == G.level(), "expected polys of same level");
-  if (CFFactory::gettype() == GaloisFieldDomain)
-    return mod (F, G);
-  zz_p::init (getCharacteristic());
-  Variable alpha;
-  CanonicalForm result;
-  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
-  {
-    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
-    zz_pE::init (NTLMipo);
-    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
-    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
-    rem (NTLF, NTLF, NTLG);
-    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
-  }
-  else
-  {
-#ifdef HAVE_FLINT
-    nmod_poly_t FLINTF, FLINTG;
-    convertFacCF2nmod_poly_t (FLINTF, F);
-    convertFacCF2nmod_poly_t (FLINTG, G);
-    nmod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
-    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
-    nmod_poly_clear (FLINTF);
-    nmod_poly_clear (FLINTG);
-#else
-    zz_pX NTLF= convertFacCF2NTLzzpX (F);
-    zz_pX NTLG= convertFacCF2NTLzzpX (G);
-    rem (NTLF, NTLF, NTLG);
-    result= convertNTLzzpX2CF(NTLF, F.mvar());
-#endif
-  }
-  return result;
-}
-
-CanonicalForm
-divNTL (const CanonicalForm& F, const CanonicalForm& G)
-{
-  if (F.inCoeffDomain() && G.isUnivariate())
-    return F;
-  else if (F.inCoeffDomain() && G.inCoeffDomain())
-    return div (F, G);
-  else if (F.isUnivariate() && G.inCoeffDomain())
-    return div (F,G);
-
-  if (getCharacteristic() == 0)
-  {
-#ifdef HAVE_FLINT
-    Variable alpha;
-    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
-      return divFLINTQ (F,G);
-    else
-    {
-      CanonicalForm Q;
-      newtonDiv (F, G, Q);
-      return Q;
-    }
-#else
-    return div (F, G);
-#endif
-  }
-
-  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
-  ASSERT (F.level() == G.level(), "expected polys of same level");
-  if (CFFactory::gettype() == GaloisFieldDomain)
-    return div (F, G);
-  zz_p::init (getCharacteristic());
-  Variable alpha;
-  CanonicalForm result;
-  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
-  {
-    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
-    zz_pE::init (NTLMipo);
-    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
-    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
-    div (NTLF, NTLF, NTLG);
-    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
-  }
-  else
-  {
-#ifdef HAVE_FLINT
-    nmod_poly_t FLINTF, FLINTG;
-    convertFacCF2nmod_poly_t (FLINTF, F);
-    convertFacCF2nmod_poly_t (FLINTG, G);
-    nmod_poly_div (FLINTF, FLINTF, FLINTG);
-    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
-    nmod_poly_clear (FLINTF);
-    nmod_poly_clear (FLINTG);
-#else
-    zz_pX NTLF= convertFacCF2NTLzzpX (F);
-    zz_pX NTLG= convertFacCF2NTLzzpX (G);
-    div (NTLF, NTLF, NTLG);
-    result= convertNTLzzpX2CF(NTLF, F.mvar());
-#endif
-  }
-  return result;
-}
-
-/*
-void
-divremNTL (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
-           CanonicalForm& R)
-{
-  if (F.inCoeffDomain() && G.isUnivariate())
-  {
-    R= F;
-    Q= 0;
-  }
-  else if (F.inCoeffDomain() && G.inCoeffDomain())
-  {
-    divrem (F, G, Q, R);
-    return;
-  }
-  else if (F.isUnivariate() && G.inCoeffDomain())
-  {
-    divrem (F, G, Q, R);
-    return;
-  }
-
-  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
-  ASSERT (F.level() == G.level(), "expected polys of same level");
-  if (CFFactory::gettype() == GaloisFieldDomain)
-  {
-    divrem (F, G, Q, R);
-    return;
-  }
-  zz_p::init (getCharacteristic());
-  Variable alpha;
-  CanonicalForm result;
-  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
-  {
-    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
-    zz_pE::init (NTLMipo);
-    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
-    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
-    zz_pEX NTLQ;
-    zz_pEX NTLR;
-    DivRem (NTLQ, NTLR, NTLF, NTLG);
-    Q= convertNTLzz_pEX2CF(NTLQ, F.mvar(), alpha);
-    R= convertNTLzz_pEX2CF(NTLR, F.mvar(), alpha);
-    return;
-  }
-  else
-  {
-    zz_pX NTLF= convertFacCF2NTLzzpX (F);
-    zz_pX NTLG= convertFacCF2NTLzzpX (G);
-    zz_pX NTLQ;
-    zz_pX NTLR;
-    DivRem (NTLQ, NTLR, NTLF, NTLG);
-    Q= convertNTLzzpX2CF(NTLQ, F.mvar());
-    R= convertNTLzzpX2CF(NTLR, F.mvar());
-    return;
-  }
-}*/
-
-CanonicalForm mod (const CanonicalForm& F, const CFList& M)
-{
-  CanonicalForm A= F;
-  for (CFListIterator i= M; i.hasItem(); i++)
-    A= mod (A, i.getItem());
-  return A;
-}
-
-#ifdef HAVE_FLINT
-void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
-{
-  int degAy= degree (A);
-  nmod_poly_init2 (result, getCharacteristic(), d*(degAy + 1));
-
-  nmod_poly_t buf;
-
-  int j, k, bufRepLength;
-  for (CFIterator i= A; i.hasTerms(); i++)
-  {
-    convertFacCF2nmod_poly_t (buf, i.coeff());
-
-    k= i.exp()*d;
-    bufRepLength= (int) nmod_poly_length (buf);
-    for (j= 0; j < bufRepLength; j++)
-      nmod_poly_set_coeff_ui (result, j + k, nmod_poly_get_coeff_ui (buf, j));
-    nmod_poly_clear (buf);
-  }
-  _nmod_poly_normalise (result);
-}
-
-/*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)
-{
-  int degAy= degree (A);
-  fmpq_poly_init2 (result, d1*(degAy + 1));
-
-  fmpq_poly_t buf;
-  fmpq_t coeff;
-
-  int k, l, bufRepLength;
-  CFIterator j;
-  for (CFIterator i= A; i.hasTerms(); i++)
-  {
-    if (i.coeff().inCoeffDomain())
-    {
-      k= d1*i.exp();
-      convertFacCF2Fmpq_poly_t (buf, i.coeff());
-      bufRepLength= (int) fmpq_poly_length(buf);
-      for (l= 0; l < bufRepLength; l++)
-      {
-        fmpq_poly_get_coeff_fmpq (coeff, buf, l);
-        fmpq_poly_set_coeff_fmpq (result, l + k, coeff);
-      }
-      fmpq_poly_clear (buf);
-    }
-    else
-    {
-      for (j= i.coeff(); j.hasTerms(); j++)
-      {
-        k= d1*i.exp();
-        k += d2*j.exp();
-        convertFacCF2Fmpq_poly_t (buf, j.coeff());
-        bufRepLength= (int) fmpq_poly_length(buf);
-        for (l= 0; l < bufRepLength; l++)
-        {
-          fmpq_poly_get_coeff_fmpq (coeff, buf, l);
-          fmpq_poly_set_coeff_fmpq (result, k + l, coeff);
-        }
-        fmpq_poly_clear (buf);
-      }
-    }
-  }
-  fmpq_clear (coeff);
-  _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)
-{
-  int degAy= degree (A);
-  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
-  nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
-  nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
-
-  nmod_poly_t buf;
-
-  int k, kk, j, bufRepLength;
-  for (CFIterator i= A; i.hasTerms(); i++)
-  {
-    convertFacCF2nmod_poly_t (buf, i.coeff());
-
-    k= i.exp()*d;
-    kk= (degAy - i.exp())*d;
-    bufRepLength= (int) nmod_poly_length (buf);
-    for (j= 0; j < bufRepLength; j++)
-    {
-      nmod_poly_set_coeff_ui (subA1, j + k,
-                              n_addmod (nmod_poly_get_coeff_ui (subA1, j+k),
-                                        nmod_poly_get_coeff_ui (buf, j),
-                                        getCharacteristic()
-                                       )
-                             );
-      nmod_poly_set_coeff_ui (subA2, j + kk,
-                              n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),
-                                        nmod_poly_get_coeff_ui (buf, j),
-                                        getCharacteristic()
-                                       )
-                             );
-    }
-    nmod_poly_clear (buf);
-  }
-  _nmod_poly_normalise (subA1);
-  _nmod_poly_normalise (subA2);
-}
-
-void
-kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
-                int d)
-{
-  int degAy= degree (A);
-  fmpz_poly_init2 (subA1, d*(degAy + 2));
-  fmpz_poly_init2 (subA2, d*(degAy + 2));
-
-  fmpz_poly_t buf;
-  fmpz_t coeff1, coeff2;
-
-  int k, kk, j, bufRepLength;
-  for (CFIterator i= A; i.hasTerms(); i++)
-  {
-    convertFacCF2Fmpz_poly_t (buf, i.coeff());
-
-    k= i.exp()*d;
-    kk= (degAy - i.exp())*d;
-    bufRepLength= (int) fmpz_poly_length (buf);
-    for (j= 0; j < bufRepLength; j++)
-    {
-      fmpz_poly_get_coeff_fmpz (coeff1, subA1, j+k);
-      fmpz_poly_get_coeff_fmpz (coeff2, buf, j);
-      fmpz_add (coeff1, coeff1, coeff2);
-      fmpz_poly_set_coeff_fmpz (subA1, j + k, coeff1);
-      fmpz_poly_get_coeff_fmpz (coeff1, subA2, j + kk);
-      fmpz_add (coeff1, coeff1, coeff2);
-      fmpz_poly_set_coeff_fmpz (subA2, j + kk, coeff1);
-    }
-    fmpz_poly_clear (buf);
-  }
-  fmpz_clear (coeff1);
-  fmpz_clear (coeff2);
-  _fmpz_poly_normalise (subA1);
-  _fmpz_poly_normalise (subA2);
-}
-
-CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
-{
-  Variable y= Variable (2);
-  Variable x= Variable (1);
-
-  fmpz_poly_t f;
-  fmpz_poly_init (f);
-  fmpz_poly_set (f, F);
-
-  fmpz_poly_t buf;
-  CanonicalForm result= 0;
-  int i= 0;
-  int degf= fmpz_poly_degree(f);
-  int k= 0;
-  int degfSubK, repLength, j;
-  fmpz_t coeff;
-  while (degf >= k)
-  {
-    degfSubK= degf - k;
-    if (degfSubK >= d)
-      repLength= d;
-    else
-      repLength= degfSubK + 1;
-
-    fmpz_poly_init2 (buf, repLength);
-    fmpz_init (coeff);
-    for (j= 0; j < repLength; j++)
-    {
-      fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
-      fmpz_poly_set_coeff_fmpz (buf, j, coeff);
-    }
-    _fmpz_poly_normalise (buf);
-
-    result += convertFmpz_poly_t2FacCF (buf, x)*power (y, i);
-    i++;
-    k= d*i;
-    fmpz_poly_clear (buf);
-    fmpz_clear (coeff);
-  }
-  fmpz_poly_clear (f);
-
-  return result;
-}
-
-CanonicalForm
-reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
-{
-  Variable y= Variable (2);
-  Variable x= Variable (1);
-
-  nmod_poly_t f, g;
-  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
-  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
-  nmod_poly_init_preinv (g, getCharacteristic(), ninv);
-  nmod_poly_set (f, F);
-  nmod_poly_set (g, G);
-  int degf= nmod_poly_degree(f);
-  int degg= nmod_poly_degree(g);
-
-
-  nmod_poly_t buf1,buf2, buf3;
-
-  if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
-    nmod_poly_fit_length (f,(long)d*(k+1));
-
-  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;
-  while (degf >= lf || lg >= 0)
-  {
-    if (degfSubLf >= d)
-      repLengthBuf1= d;
-    else if (degfSubLf < 0)
-      repLengthBuf1= 0;
-    else
-      repLengthBuf1= degfSubLf + 1;
-    nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
-
-    for (ind= 0; ind < repLengthBuf1; ind++)
-      nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
-    _nmod_poly_normalise (buf1);
-
-    repLengthBuf1= nmod_poly_length (buf1);
-
-    if (deggSubLg >= d - 1)
-      repLengthBuf2= d - 1;
-    else if (deggSubLg < 0)
-      repLengthBuf2= 0;
-    else
-      repLengthBuf2= deggSubLg + 1;
-
-    nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
-    for (ind= 0; ind < repLengthBuf2; ind++)
-      nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
-
-    _nmod_poly_normalise (buf2);
-    repLengthBuf2= nmod_poly_length (buf2);
-
-    nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
-    for (ind= 0; ind < repLengthBuf1; ind++)
-      nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
-    for (ind= repLengthBuf1; ind < d; ind++)
-      nmod_poly_set_coeff_ui (buf3, ind, 0);
-    for (ind= 0; ind < repLengthBuf2; ind++)
-      nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
-    _nmod_poly_normalise (buf3);
-
-    result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
-    i++;
-
-
-    lf= i*d;
-    degfSubLf= degf - lf;
-
-    lg= d*(k-i);
-    deggSubLg= degg - lg;
-
-    if (lg >= 0 && deggSubLg > 0)
-    {
-      if (repLengthBuf2 > degfSubLf + 1)
-        degfSubLf= repLengthBuf2 - 1;
-      tmp= tmin (repLengthBuf1, deggSubLg + 1);
-      for (ind= 0; ind < tmp; ind++)
-        nmod_poly_set_coeff_ui (g, ind + lg,
-                                n_submod (nmod_poly_get_coeff_ui (g, ind + lg),
-                                          nmod_poly_get_coeff_ui (buf1, ind),
-                                          getCharacteristic()
-                                         )
-                               );
-    }
-    if (lg < 0)
-    {
-      nmod_poly_clear (buf1);
-      nmod_poly_clear (buf2);
-      nmod_poly_clear (buf3);
-      break;
-    }
-    if (degfSubLf >= 0)
-    {
-      for (ind= 0; ind < repLengthBuf2; ind++)
-        nmod_poly_set_coeff_ui (f, ind + lf,
-                                n_submod (nmod_poly_get_coeff_ui (f, ind + lf),
-                                          nmod_poly_get_coeff_ui (buf2, ind),
-                                          getCharacteristic()
-                                         )
-                               );
-    }
-    nmod_poly_clear (buf1);
-    nmod_poly_clear (buf2);
-    nmod_poly_clear (buf3);
-  }
-
-  nmod_poly_clear (f);
-  nmod_poly_clear (g);
-
-  return result;
-}
-
-CanonicalForm
-reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
-{
-  Variable y= Variable (2);
-  Variable x= Variable (1);
-
-  fmpz_poly_t f, g;
-  fmpz_poly_init (f);
-  fmpz_poly_init (g);
-  fmpz_poly_set (f, F);
-  fmpz_poly_set (g, G);
-  int degf= fmpz_poly_degree(f);
-  int degg= fmpz_poly_degree(g);
-
-
-  fmpz_poly_t buf1,buf2, buf3;
-
-  if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
-    fmpz_poly_fit_length (f,(long)d*(k+1));
-
-  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;
-  fmpz_t tmp1, tmp2;
-  while (degf >= lf || lg >= 0)
-  {
-    if (degfSubLf >= d)
-      repLengthBuf1= d;
-    else if (degfSubLf < 0)
-      repLengthBuf1= 0;
-    else
-      repLengthBuf1= degfSubLf + 1;
-
-    fmpz_poly_init2 (buf1, repLengthBuf1);
-
-    for (ind= 0; ind < repLengthBuf1; ind++)
-    {
-      fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
-      fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
-    }
-    _fmpz_poly_normalise (buf1);
-
-    repLengthBuf1= fmpz_poly_length (buf1);
-
-    if (deggSubLg >= d - 1)
-      repLengthBuf2= d - 1;
-    else if (deggSubLg < 0)
-      repLengthBuf2= 0;
-    else
-      repLengthBuf2= deggSubLg + 1;
-
-    fmpz_poly_init2 (buf2, repLengthBuf2);
-
-    for (ind= 0; ind < repLengthBuf2; ind++)
-    {
-      fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
-      fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
-    }
-
-    _fmpz_poly_normalise (buf2);
-    repLengthBuf2= fmpz_poly_length (buf2);
-
-    fmpz_poly_init2 (buf3, repLengthBuf2 + d);
-    for (ind= 0; ind < repLengthBuf1; ind++)
-    {
-      fmpz_poly_get_coeff_fmpz (tmp1, buf1, ind);  //oder fmpz_set (fmpz_poly_get_coeff_ptr (buf3, ind),fmpz_poly_get_coeff_ptr (buf1, ind))
-      fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
-    }
-    for (ind= repLengthBuf1; ind < d; ind++)
-      fmpz_poly_set_coeff_ui (buf3, ind, 0);
-    for (ind= 0; ind < repLengthBuf2; ind++)
-    {
-      fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
-      fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
-    }
-    _fmpz_poly_normalise (buf3);
-
-    result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
-    i++;
-
-
-    lf= i*d;
-    degfSubLf= degf - lf;
-
-    lg= d*(k-i);
-    deggSubLg= degg - lg;
-
-    if (lg >= 0 && deggSubLg > 0)
-    {
-      if (repLengthBuf2 > degfSubLf + 1)
-        degfSubLf= repLengthBuf2 - 1;
-      tmp= tmin (repLengthBuf1, deggSubLg + 1);
-      for (ind= 0; ind < tmp; ind++)
-      {
-        fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
-        fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
-        fmpz_sub (tmp1, tmp1, tmp2);
-        fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
-      }
-    }
-    if (lg < 0)
-    {
-      fmpz_poly_clear (buf1);
-      fmpz_poly_clear (buf2);
-      fmpz_poly_clear (buf3);
-      break;
-    }
-    if (degfSubLf >= 0)
-    {
-      for (ind= 0; ind < repLengthBuf2; ind++)
-      {
-        fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
-        fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
-        fmpz_sub (tmp1, tmp1, tmp2);
-        fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
-      }
-    }
-    fmpz_poly_clear (buf1);
-    fmpz_poly_clear (buf2);
-    fmpz_poly_clear (buf3);
-  }
-
-  fmpz_poly_clear (f);
-  fmpz_poly_clear (g);
-  fmpz_clear (tmp1);
-  fmpz_clear (tmp2);
-
-  return result;
-}
-
-CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
-{
-  Variable y= Variable (2);
-  Variable x= Variable (1);
-
-  nmod_poly_t f;
-  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
-  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
-  nmod_poly_set (f, F);
-
-  nmod_poly_t buf;
-  CanonicalForm result= 0;
-  int i= 0;
-  int degf= nmod_poly_degree(f);
-  int k= 0;
-  int degfSubK, repLength, j;
-  while (degf >= k)
-  {
-    degfSubK= degf - k;
-    if (degfSubK >= d)
-      repLength= d;
-    else
-      repLength= degfSubK + 1;
-
-    nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
-    for (j= 0; j < repLength; j++)
-      nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
-    _nmod_poly_normalise (buf);
-
-    result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
-    i++;
-    k= d*i;
-    nmod_poly_clear (buf);
-  }
-  nmod_poly_clear (f);
-
-  return result;
-}
-
-CanonicalForm
-mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
-                    CanonicalForm& M)
-{
-  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
-  d1 /= 2;
-  d1 += 1;
-
-  nmod_poly_t F1, F2;
-  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
-  nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
-  nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
-  kronSubRecipro (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);
-
-  int k= d1*degree (M);
-  nmod_poly_mullow (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= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG
-         + d1*(2+taildegF + taildegG);
-  nmod_poly_mulhigh (F2, F2, G2, b);
-  nmod_poly_shift_right (F2, F2, b);
-  int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
-
-
-  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
-
-  nmod_poly_clear (F1);
-  nmod_poly_clear (F2);
-  nmod_poly_clear (G1);
-  nmod_poly_clear (G2);
-  return result;
-}
-
-CanonicalForm
-mulMod2FLINTFp (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 mulMod2FLINTFpReci (A, B, M);
-
-  nmod_poly_t FLINTA, FLINTB;
-  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
-  nmod_poly_init_preinv (FLINTA, getCharacteristic(), ninv);
-  nmod_poly_init_preinv (FLINTB, getCharacteristic(), ninv);
-  kronSubFp (FLINTA, A, d1);
-  kronSubFp (FLINTB, B, d1);
-
-  int k= d1*degree (M);
-  nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
-
-  A= reverseSubstFp (FLINTA, d1);
-
-  nmod_poly_clear (FLINTA);
-  nmod_poly_clear (FLINTB);
-  return A;
-}
-
-CanonicalForm
-mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
-                    CanonicalForm& M)
-{
-  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
-  d1 /= 2;
-  d1 += 1;
-
-  fmpz_poly_t F1, F2;
-  fmpz_poly_init (F1);
-  fmpz_poly_init (F2);
-  kronSubRecipro (F1, F2, F, d1);
-
-  fmpz_poly_t G1, G2;
-  fmpz_poly_init (G1);
-  fmpz_poly_init (G2);
-  kronSubRecipro (G1, G2, G, d1);
-
-  int k= d1*degree (M);
-  fmpz_poly_mullow (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= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG
-         + d1*(2+taildegF + taildegG);
-  fmpz_poly_mulhigh_n (F2, F2, G2, b);
-  fmpz_poly_shift_right (F2, F2, b);
-  int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
-
-  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
-
-  fmpz_poly_clear (F1);
-  fmpz_poly_clear (F2);
-  fmpz_poly_clear (G1);
-  fmpz_poly_clear (G2);
-  return result;
-}
-
-CanonicalForm
-mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
-                CanonicalForm& M)
-{
-  CanonicalForm A= F;
-  CanonicalForm B= G;
-
-  int degAx= degree (A, 1);
-  int degBx= degree (B, 1);
-  int d1= degAx + 1 + degBx;
-
-  CanonicalForm f= bCommonDen (F);
-  CanonicalForm g= bCommonDen (G);
-  A *= f;
-  B *= g;
-
-  fmpz_poly_t FLINTA, FLINTB;
-  fmpz_poly_init (FLINTA);
-  fmpz_poly_init (FLINTB);
-  kronSub (FLINTA, A, d1);
-  kronSub (FLINTB, B, d1);
-  int k= d1*degree (M);
-
-  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
-  A= reverseSubstQ (FLINTA, d1);
-  fmpz_poly_clear (FLINTA);
-  fmpz_poly_clear (FLINTB);
-  return A/(f*g);
-}
-
-#endif
-
-CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
-                       const CanonicalForm& M)
-{
-  if (A.isZero() || B.isZero())
-    return 0;
-
-  ASSERT (M.isUnivariate(), "M must be univariate");
-
-  CanonicalForm F= mod (A, M);
-  CanonicalForm G= mod (B, M);
-  if (F.inCoeffDomain() || G.inCoeffDomain())
-    return F*G;
-  Variable y= M.mvar();
-  int degF= degree (F, y);
-  int degG= degree (G, y);
-
-  if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
-      (F.level() == G.level()))
-  {
-    CanonicalForm result= mulNTL (F, G);
-    return mod (result, M);
-  }
-  else if (degF <= 1 && degG <= 1)
-  {
-    CanonicalForm result= F*G;
-    return mod (result, M);
-  }
-
-  int sizeF= size (F);
-  int sizeG= size (G);
-
-  int fallBackToNaive= 50;
-  if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
-    return mod (F*G, M);
-
-#ifdef HAVE_FLINT
-  Variable alpha;
-  if (getCharacteristic() == 0 && !hasFirstAlgVar (F, alpha) && ! hasFirstAlgVar (G, alpha))
-    return mulMod2FLINTQ (F, G, M);
-#endif
-
-  if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
-      (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
-    return mulMod2NTLFq (F, G, M);
-
-  int m= (int) ceil (degree (M)/2.0);
-  if (degF >= m || degG >= m)
-  {
-    CanonicalForm MLo= power (y, m);
-    CanonicalForm MHi= power (y, degree (M) - m);
-    CanonicalForm F0= mod (F, MLo);
-    CanonicalForm F1= div (F, MLo);
-    CanonicalForm G0= mod (G, MLo);
-    CanonicalForm G1= div (G, MLo);
-    CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
-    CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
-    CanonicalForm F0G0= mulMod2 (F0, G0, M);
-    return F0G0 + MLo*(F0G1 + F1G0);
-  }
-  else
-  {
-    m= (int) ceil (tmax (degF, degG)/2.0);
-    CanonicalForm yToM= power (y, m);
-    CanonicalForm F0= mod (F, yToM);
-    CanonicalForm F1= div (F, yToM);
-    CanonicalForm G0= mod (G, yToM);
-    CanonicalForm G1= div (G, yToM);
-    CanonicalForm H00= mulMod2 (F0, G0, M);
-    CanonicalForm H11= mulMod2 (F1, G1, M);
-    CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
-    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
-  }
-  DEBOUTLN (cerr, "fatal end in mulMod2");
-}
-
-CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
-                      const CFList& MOD)
-{
-  if (A.isZero() || B.isZero())
-    return 0;
-
-  if (MOD.length() == 1)
-    return mulMod2 (A, B, MOD.getLast());
-
-  CanonicalForm M= MOD.getLast();
-  CanonicalForm F= mod (A, M);
-  CanonicalForm G= mod (B, M);
-  if (F.inCoeffDomain() || G.inCoeffDomain())
-    return F*G;
-  Variable y= M.mvar();
-  int degF= degree (F, y);
-  int degG= degree (G, y);
-
-  if ((degF <= 1 && F.level() <= M.level()) &&
-      (degG <= 1 && G.level() <= M.level()))
-  {
-    CFList buf= MOD;
-    buf.removeLast();
-    if (degF == 1 && degG == 1)
-    {
-      CanonicalForm F0= mod (F, y);
-      CanonicalForm F1= div (F, y);
-      CanonicalForm G0= mod (G, y);
-      CanonicalForm G1= div (G, y);
-      if (degree (M) > 2)
-      {
-        CanonicalForm H00= mulMod (F0, G0, buf);
-        CanonicalForm H11= mulMod (F1, G1, buf);
-        CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
-        return H11*y*y + (H01 - H00 - H11)*y + H00;
-      }
-      else //here degree (M) == 2
-      {
-        buf.append (y);
-        CanonicalForm F0G1= mulMod (F0, G1, buf);
-        CanonicalForm F1G0= mulMod (F1, G0, buf);
-        CanonicalForm F0G0= mulMod (F0, G0, MOD);
-        CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
-        return result;
-      }
-    }
-    else if (degF == 1 && degG == 0)
-      return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
-    else if (degF == 0 && degG == 1)
-      return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
-    else
-      return mulMod (F, G, buf);
-  }
-  int m= (int) ceil (degree (M)/2.0);
-  if (degF >= m || degG >= m)
-  {
-    CanonicalForm MLo= power (y, m);
-    CanonicalForm MHi= power (y, degree (M) - m);
-    CanonicalForm F0= mod (F, MLo);
-    CanonicalForm F1= div (F, MLo);
-    CanonicalForm G0= mod (G, MLo);
-    CanonicalForm G1= div (G, MLo);
-    CFList buf= MOD;
-    buf.removeLast();
-    buf.append (MHi);
-    CanonicalForm F0G1= mulMod (F0, G1, buf);
-    CanonicalForm F1G0= mulMod (F1, G0, buf);
-    CanonicalForm F0G0= mulMod (F0, G0, MOD);
-    return F0G0 + MLo*(F0G1 + F1G0);
-  }
-  else
-  {
-    m= (int) ceil (tmax (degF, degG)/2.0);
-    CanonicalForm yToM= power (y, m);
-    CanonicalForm F0= mod (F, yToM);
-    CanonicalForm F1= div (F, yToM);
-    CanonicalForm G0= mod (G, yToM);
-    CanonicalForm G1= div (G, yToM);
-    CanonicalForm H00= mulMod (F0, G0, MOD);
-    CanonicalForm H11= mulMod (F1, G1, MOD);
-    CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
-    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
-  }
-  DEBOUTLN (cerr, "fatal end in mulMod");
-}
-
-CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
-{
-  if (L.isEmpty())
-    return 1;
-  int l= L.length();
-  if (l == 1)
-    return mod (L.getFirst(), M);
-  else if (l == 2) {
-    CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
-    return result;
-  }
-  else
-  {
-    l /= 2;
-    CFList tmp1, tmp2;
-    CFListIterator i= L;
-    CanonicalForm buf1, buf2;
-    for (int j= 1; j <= l; j++, i++)
-      tmp1.append (i.getItem());
-    tmp2= Difference (L, tmp1);
-    buf1= prodMod (tmp1, M);
-    buf2= prodMod (tmp2, M);
-    CanonicalForm result= mulMod2 (buf1, buf2, M);
-    return result;
-  }
-}
-
-CanonicalForm prodMod (const CFList& L, const CFList& M)
-{
-  if (L.isEmpty())
-    return 1;
-  else if (L.length() == 1)
-    return L.getFirst();
-  else if (L.length() == 2)
-    return mulMod (L.getFirst(), L.getLast(), M);
-  else
-  {
-    int l= L.length()/2;
-    CFListIterator i= L;
-    CFList tmp1, tmp2;
-    CanonicalForm buf1, buf2;
-    for (int j= 1; j <= l; j++, i++)
-      tmp1.append (i.getItem());
-    tmp2= Difference (L, tmp1);
-    buf1= prodMod (tmp1, M);
-    buf2= prodMod (tmp2, M);
-    return mulMod (buf1, buf2, M);
-  }
-}
-
-
-CanonicalForm reverse (const CanonicalForm& F, int d)
-{
-  if (d == 0)
-    return F;
-  CanonicalForm A= F;
-  Variable y= Variable (2);
-  Variable x= Variable (1);
-  if (degree (A, x) > 0)
-  {
-    A= swapvar (A, x, y);
-    CanonicalForm result= 0;
-    CFIterator i= A;
-    while (d - i.exp() < 0)
-      i++;
-
-    for (; i.hasTerms() && (d - i.exp() >= 0); i++)
-      result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
-    return result;
-  }
-  else
-    return A*power (x, d);
-}
-
-CanonicalForm
-newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
-{
-  int l= ilog2(n);
-
-  CanonicalForm g= mod (F, M)[0] [0];
-
-  ASSERT (!g.isZero(), "expected a unit");
-
-  Variable alpha;
-
-  if (!g.isOne())
-    g = 1/g;
-  Variable x= Variable (1);
-  CanonicalForm result;
-  int exp= 0;
-  if (n & 1)
-  {
-    result= g;
-    exp= 1;
-  }
-  CanonicalForm h;
-
-  for (int i= 1; i <= l; i++)
-  {
-    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
-    h= mod (h, power (x, (1 << i)) - 1);
-    h= div (h, power (x, (1 << (i - 1))));
-    h= mod (h, M);
-    g -= power (x, (1 << (i - 1)))*
-         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
-
-    if (n & (1 << i))
-    {
-      if (exp)
-      {
-        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
-        h= mod (h, power (x, exp + (1 << i)) - 1);
-        h= div (h, power (x, exp));
-        h= mod (h, M);
-        result -= power(x, exp)*mod (mulMod2 (g, h, M),
-                                       power (x, (1 << i)));
-        exp += (1 << i);
-      }
-      else
-      {
-        exp= (1 << i);
-        result= g;
-      }
-    }
-  }
-
-  return result;
-}
-
-CanonicalForm
-newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
-           M)
-{
-  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
-  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
-
-  CanonicalForm A= mod (F, M);
-  CanonicalForm B= mod (G, M);
-
-  Variable x= Variable (1);
-  int degA= degree (A, x);
-  int degB= degree (B, x);
-  int m= degA - degB;
-  if (m < 0)
-    return 0;
-
-  Variable v;
-  CanonicalForm Q;
-  if (degB < 1 || CFFactory::gettype() == GaloisFieldDomain)
-  {
-    CanonicalForm R;
-    divrem2 (A, B, Q, R, M);
-  }
-  else
-  {
-    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
-    {
-      CanonicalForm R= reverse (A, degA);
-      CanonicalForm revB= reverse (B, degB);
-      revB= newtonInverse (revB, m + 1, M);
-      Q= mulMod2 (R, revB, M);
-      Q= mod (Q, power (x, m + 1));
-      Q= reverse (Q, m);
-    }
-    else
-    {
-      zz_pX mipo= convertFacCF2NTLzzpX (M);
-      Variable y= Variable (2);
-      zz_pEX NTLA, NTLB;
-      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
-      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
-      div (NTLA, NTLA, NTLB);
-      Q= convertNTLzz_pEX2CF (NTLA, x, y);
-    }
-  }
-
-  return Q;
-}
-
-void
-newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
-              CanonicalForm& R, const CanonicalForm& M)
-{
-  CanonicalForm A= mod (F, M);
-  CanonicalForm B= mod (G, M);
-  Variable x= Variable (1);
-  int degA= degree (A, x);
-  int degB= degree (B, x);
-  int m= degA - degB;
-
-  if (m < 0)
-  {
-    R= A;
-    Q= 0;
-    return;
-  }
-
-  Variable v;
-  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
-  {
-     divrem2 (A, B, Q, R, M);
-  }
-  else
-  {
-    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
-    {
-      R= reverse (A, degA);
-
-      CanonicalForm revB= reverse (B, degB);
-      revB= newtonInverse (revB, m + 1, M);
-      Q= mulMod2 (R, revB, M);
-
-      Q= mod (Q, power (x, m + 1));
-      Q= reverse (Q, m);
-
-      R= A - mulMod2 (Q, B, M);
-    }
-    else
-    {
-      zz_pX mipo= convertFacCF2NTLzzpX (M);
-      Variable y= Variable (2);
-      zz_pEX NTLA, NTLB;
-      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
-      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
-      zz_pEX NTLQ, NTLR;
-      DivRem (NTLQ, NTLR, NTLA, NTLB);
-      Q= convertNTLzz_pEX2CF (NTLQ, x, y);
-      R= convertNTLzz_pEX2CF (NTLR, x, y);
-    }
-  }
-}
-
-static inline
-CFList split (const CanonicalForm& F, const int m, const Variable& x)
-{
-  CanonicalForm A= F;
-  CanonicalForm buf= 0;
-  bool swap= false;
-  if (degree (A, x) <= 0)
-    return CFList(A);
-  else if (x.level() != A.level())
-  {
-    swap= true;
-    A= swapvar (A, x, A.mvar());
-  }
-
-  int j= (int) floor ((double) degree (A)/ m);
-  CFList result;
-  CFIterator i= A;
-  for (; j >= 0; j--)
-  {
-    while (i.hasTerms() && i.exp() - j*m >= 0)
-    {
-      if (swap)
-        buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
-      else
-        buf += i.coeff()*power (x, i.exp() - j*m);
-      i++;
-    }
-    if (swap)
-      result.append (swapvar (buf, x, F.mvar()));
-    else
-      result.append (buf);
-    buf= 0;
-  }
-  return result;
-}
-
-static inline
-void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
-               CanonicalForm& R, const CFList& M);
-
-static inline
-void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
-               CanonicalForm& R, const CFList& M)
-{
-  CanonicalForm A= mod (F, M);
-  CanonicalForm B= mod (G, M);
-  Variable x= Variable (1);
-  int degB= degree (B, x);
-  int degA= degree (A, x);
-  if (degA < degB)
-  {
-    Q= 0;
-    R= A;
-    return;
-  }
-  ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
-  if (degB < 1)
-  {
-    divrem (A, B, Q, R);
-    Q= mod (Q, M);
-    R= mod (R, M);
-    return;
-  }
-
-  int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
-  CFList splitA= split (A, m, x);
-  if (splitA.length() == 3)
-    splitA.insert (0);
-  if (splitA.length() == 2)
-  {
-    splitA.insert (0);
-    splitA.insert (0);
-  }
-  if (splitA.length() == 1)
-  {
-    splitA.insert (0);
-    splitA.insert (0);
-    splitA.insert (0);
-  }
-
-  CanonicalForm xToM= power (x, m);
-
-  CFListIterator i= splitA;
-  CanonicalForm H= i.getItem();
-  i++;
-  H *= xToM;
-  H += i.getItem();
-  i++;
-  H *= xToM;
-  H += i.getItem();
-  i++;
-
-  divrem32 (H, B, Q, R, M);
-
-  CFList splitR= split (R, m, x);
-  if (splitR.length() == 1)
-    splitR.insert (0);
-
-  H= splitR.getFirst();
-  H *= xToM;
-  H += splitR.getLast();
-  H *= xToM;
-  H += i.getItem();
-
-  CanonicalForm bufQ;
-  divrem32 (H, B, bufQ, R, M);
-
-  Q *= xToM;
-  Q += bufQ;
-  return;
-}
-
-static inline
-void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
-               CanonicalForm& R, const CFList& M)
-{
-  CanonicalForm A= mod (F, M);
-  CanonicalForm B= mod (G, M);
-  Variable x= Variable (1);
-  int degB= degree (B, x);
-  int degA= degree (A, x);
-  if (degA < degB)
-  {
-    Q= 0;
-    R= A;
-    return;
-  }
-  ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
-  if (degB < 1)
-  {
-    divrem (A, B, Q, R);
-    Q= mod (Q, M);
-    R= mod (R, M);
-    return;
-  }
-  int m= (int) ceil ((double) (degB + 1)/ 2.0);
-
-  CFList splitA= split (A, m, x);
-  CFList splitB= split (B, m, x);
-
-  if (splitA.length() == 2)
-  {
-    splitA.insert (0);
-  }
-  if (splitA.length() == 1)
-  {
-    splitA.insert (0);
-    splitA.insert (0);
-  }
-  CanonicalForm xToM= power (x, m);
-
-  CanonicalForm H;
-  CFListIterator i= splitA;
-  i++;
-
-  if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
-  {
-    H= splitA.getFirst()*xToM + i.getItem();
-    divrem21 (H, splitB.getFirst(), Q, R, M);
-  }
-  else
-  {
-    R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
-       splitB.getFirst()*xToM;
-    Q= xToM - 1;
-  }
-
-  H= mulMod (Q, splitB.getLast(), M);
-
-  R= R*xToM + splitA.getLast() - H;
-
-  while (degree (R, x) >= degB)
-  {
-    xToM= power (x, degree (R, x) - degB);
-    Q += LC (R, x)*xToM;
-    R -= mulMod (LC (R, x), B, M)*xToM;
-    Q= mod (Q, M);
-    R= mod (R, M);
-  }
-
-  return;
-}
-
-void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
-              CanonicalForm& R, const CanonicalForm& M)
-{
-  CanonicalForm A= mod (F, M);
-  CanonicalForm B= mod (G, M);
-
-  if (B.inCoeffDomain())
-  {
-    divrem (A, B, Q, R);
-    return;
-  }
-  if (A.inCoeffDomain() && !B.inCoeffDomain())
-  {
-    Q= 0;
-    R= A;
-    return;
-  }
-
-  if (B.level() < A.level())
-  {
-    divrem (A, B, Q, R);
-    return;
-  }
-  if (A.level() > B.level())
-  {
-    R= A;
-    Q= 0;
-    return;
-  }
-  if (B.level() == 1 && B.isUnivariate())
-  {
-    divrem (A, B, Q, R);
-    return;
-  }
-  if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
-  {
-    Q= 0;
-    R= A;
-    return;
-  }
-
-  Variable x= Variable (1);
-  int degB= degree (B, x);
-  if (degB > degree (A, x))
-  {
-    Q= 0;
-    R= A;
-    return;
-  }
-
-  CFList splitA= split (A, degB, x);
-
-  CanonicalForm xToDegB= power (x, degB);
-  CanonicalForm H, bufQ;
-  Q= 0;
-  CFListIterator i= splitA;
-  H= i.getItem()*xToDegB;
-  i++;
-  H += i.getItem();
-  CFList buf;
-  while (i.hasItem())
-  {
-    buf= CFList (M);
-    divrem21 (H, B, bufQ, R, buf);
-    i++;
-    if (i.hasItem())
-      H= R*xToDegB + i.getItem();
-    Q *= xToDegB;
-    Q += bufQ;
-  }
-  return;
-}
-
-void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
-             CanonicalForm& R, const CFList& MOD)
-{
-  CanonicalForm A= mod (F, MOD);
-  CanonicalForm B= mod (G, MOD);
-  Variable x= Variable (1);
-  int degB= degree (B, x);
-  if (degB > degree (A, x))
-  {
-    Q= 0;
-    R= A;
-    return;
-  }
-
-  if (degB <= 0)
-  {
-    divrem (A, B, Q, R);
-    Q= mod (Q, MOD);
-    R= mod (R, MOD);
-    return;
-  }
-  CFList splitA= split (A, degB, x);
-
-  CanonicalForm xToDegB= power (x, degB);
-  CanonicalForm H, bufQ;
-  Q= 0;
-  CFListIterator i= splitA;
-  H= i.getItem()*xToDegB;
-  i++;
-  H += i.getItem();
-  while (i.hasItem())
-  {
-    divrem21 (H, B, bufQ, R, MOD);
-    i++;
-    if (i.hasItem())
-      H= R*xToDegB + i.getItem();
-    Q *= xToDegB;
-    Q += bufQ;
-  }
-  return;
-}
-
 void sortList (CFList& list, const Variable& x)
 {
   int l= 1;
diff --git a/factory/facHensel.h b/factory/facHensel.h
index bf495d4..e3bf585 100644
--- a/factory/facHensel.h
+++ b/factory/facHensel.h
@@ -3,7 +3,7 @@
 \*****************************************************************************/
 /** @file facHensel.h
  *
- * This file defines functions for Hensel lifting and modular multiplication.
+ * This file defines functions for Hensel lifting.
  *
  * ABSTRACT: function are used for bi- and multivariate factorization over
  * finite fields
@@ -26,142 +26,6 @@
 #include "templates/ftmpl_functions.h"
 #include "algext.h"
 
-/// multiplication of univariate polys over a finite field using NTL, if we are
-/// in GF factory's default multiplication is used.
-///
-/// @return @a mulNTL returns F*G
-CanonicalForm
-mulNTL (const CanonicalForm& F, ///< [in] a univariate poly
-        const CanonicalForm& G  ///< [in] a univariate poly
-       );
-
-/// mod of univariate polys over a finite field using NTL, if we are
-/// in GF factory's default mod is used.
-///
-/// @return @a modNTL returns F mod G
-CanonicalForm
-modNTL (const CanonicalForm& F, ///< [in] a univariate poly
-        const CanonicalForm& G  ///< [in] a univariate poly
-       );
-
-/// division of univariate polys over a finite field using NTL, if we are
-/// in GF factory's default division is used.
-///
-/// @return @a divNTL returns F/G
-CanonicalForm
-divNTL (const CanonicalForm& F, ///< [in] a univariate poly
-        const CanonicalForm& G  ///< [in] a univariate poly
-       );
-
-/*/// division with remainder of univariate polys over a finite field using NTL,
-/// if we are in GF factory's default division with remainder is used.
-void
-divremNTL (const CanonicalForm& F, ///< [in] a univariate poly
-           const CanonicalForm& G, ///< [in] a univariate poly
-           CanonicalForm& Q,       ///< [in,out] dividend
-           CanonicalForm& R        ///< [in,out] remainder
-          );*/
-
-/// division with remainder of @a F by
-/// @a G wrt Variable (1) modulo @a M.
-///
-/// @return @a Q returns the dividend, @a R returns the remainder.
-/// @sa divrem()
-void divrem2 (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
-              const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
-              CanonicalForm& Q,       ///< [in,out] dividend
-              CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
-                                      ///< degree (G, 1)
-              const CanonicalForm& M  ///< [in] power of Variable (2)
-             );
-
-/// division with remainder of @a F by
-/// @a G wrt Variable (1) modulo @a MOD.
-///
-/// @sa divrem2()
-void divrem (
-           const CanonicalForm& F, ///< [in] multivariate, compressed polynomial
-           const CanonicalForm& G, ///< [in] multivariate, compressed polynomial
-           CanonicalForm& Q,       ///< [in,out] dividend
-           CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
-                                   ///< degree (G, 1)
-           const CFList& MOD       ///< [in] only contains powers of
-                                   ///< Variables of level higher than 1
-            );
-
-
-/// division with remainder of @a F by
-/// @a G wrt Variable (1) modulo @a M using Newton inversion
-///
-/// @return @a Q returns the dividend, @a R returns the remainder.
-/// @sa divrem2(), newtonDiv()
-void
-newtonDivrem (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
-              const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
-                                      ///< which is monic in Variable (1)
-              CanonicalForm& Q,       ///< [in,out] dividend
-              CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
-                                      ///< degree (G, 1)
-              const CanonicalForm& M  ///< [in] power of Variable (2)
-             );
-
-/// division of @a F by
-/// @a G wrt Variable (1) modulo @a M using Newton inversion
-///
-/// @return @a newtonDiv returns the dividend
-/// @sa divrem2(), newtonDivrem()
-CanonicalForm
-newtonDiv (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
-           const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
-                                   ///< which is monic in Variable (1)
-           const CanonicalForm& M  ///< [in] power of Variable (2)
-          );
-
-/// reduce @a F modulo elements in @a M.
-///
-/// @return @a mod returns @a F modulo @a M
-CanonicalForm mod (const CanonicalForm& F, ///< [in] compressed polynomial
-                   const CFList& M         ///< [in] list containing only
-                                           ///< univariate polynomials
-                  );
-
-/// Karatsuba style modular multiplication for bivariate polynomials.
-///
-/// @return @a mulMod2 returns @a A * @a B mod @a M.
-CanonicalForm
-mulMod2 (const CanonicalForm& A, ///< [in] bivariate, compressed polynomial
-         const CanonicalForm& B, ///< [in] bivariate, compressed polynomial
-         const CanonicalForm& M  ///< [in] power of Variable (2)
-        );
-
-/// Karatsuba style modular multiplication for multivariate polynomials.
-///
-/// @return @a mulMod2 returns @a A * @a B mod @a MOD.
-CanonicalForm
-mulMod (const CanonicalForm& A, ///< [in] multivariate, compressed polynomial
-        const CanonicalForm& B, ///< [in] multivariate, compressed polynomial
-        const CFList& MOD       ///< [in] only contains powers of
-                                ///< Variables of level higher than 1
-       );
-
-/// product of all elements in @a L modulo @a M via divide-and-conquer.
-///
-/// @return @a prodMod returns product of all elements in @a L modulo @a M.
-CanonicalForm
-prodMod (const CFList& L,       ///< [in] contains only bivariate, compressed
-                                ///< polynomials
-         const CanonicalForm& M ///< [in] power of Variable (2)
-        );
-
-/// product of all elements in @a L modulo @a M via divide-and-conquer.
-///
-/// @return @a prodMod returns product of all elements in @a L modulo @a M.
-CanonicalForm
-prodMod (const CFList& L, ///< [in] contains multivariate, compressed
-                          ///< polynomials
-         const CFList& M  ///< [in] contains only powers of Variables
-        );
-
 /// sort a list of polynomials by their degree in @a x.
 ///
 void sortList (CFList& list,     ///< [in, out] list of polys, sorted list
diff --git a/factory/facMul.cc b/factory/facMul.cc
new file mode 100644
index 0000000..f022ac5
--- /dev/null
+++ b/factory/facMul.cc
@@ -0,0 +1,2512 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR
+\*****************************************************************************/
+/** @file facHensel.cc
+ *
+ * This file implements functions for fast multiplication and division with
+ * remainder
+ *
+ * @author Martin Lee
+ *
+ **/
+/*****************************************************************************/
+
+#include "debug.h"
+
+#include "canonicalform.h"
+#include "facMul.h"
+#include "algext.h"
+#include "templates/ftmpl_functions.h"
+
+#ifdef HAVE_NTL
+#include <NTL/lzz_pEX.h>
+#include "NTLconvert.h"
+
+#ifdef HAVE_FLINT
+#include "FLINTconvert.h"
+#endif
+
+#ifdef HAVE_FLINT
+void kronSub (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().inBaseDomain())
+      convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d), i.coeff());
+    else
+      for (j= i.coeff(); j.hasTerms(); j++)
+        convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d+j.exp()),
+                        j.coeff());
+  }
+  _fmpz_poly_normalise(result);
+}
+
+
+CanonicalForm
+reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha,
+              const CanonicalForm& den)
+{
+  Variable x= Variable (1);
+
+  CanonicalForm result= 0;
+  int i= 0;
+  int degf= fmpz_poly_degree (F);
+  int k= 0;
+  int degfSubK;
+  int repLength, j;
+  CanonicalForm coeff;
+  fmpz* tmp;
+  while (degf >= k)
+  {
+    coeff= 0;
+    degfSubK= degf - k;
+    if (degfSubK >= d)
+      repLength= d;
+    else
+      repLength= degfSubK + 1;
+
+    for (j= 0; j < repLength; j++)
+    {
+      tmp= fmpz_poly_get_coeff_ptr (F, j+k);
+      if (!fmpz_is_zero (tmp))
+      {
+        CanonicalForm ff= convertFmpz2CF (tmp)/den;
+        coeff += ff*power (alpha, j);
+      }
+    }
+    result += coeff*power (x, i);
+    i++;
+    k= d*i;
+  }
+  return result;
+}
+
+CanonicalForm
+mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G,
+            const Variable& alpha)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+
+  CanonicalForm denA= bCommonDen (A);
+  CanonicalForm denB= bCommonDen (B);
+
+  A *= denA;
+  B *= denB;
+  int degAa= degree (A, alpha);
+  int degBa= degree (B, alpha);
+  int d= degAa + 1 + degBa;
+
+  fmpz_poly_t FLINTA,FLINTB;
+  fmpz_poly_init (FLINTA);
+  fmpz_poly_init (FLINTB);
+  kronSub (FLINTA, A, d);
+  kronSub (FLINTB, B, d);
+
+  fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
+
+  denA *= denB;
+  A= reverseSubst (FLINTA, d, alpha, denA);
+
+  fmpz_poly_clear (FLINTA);
+  fmpz_poly_clear (FLINTB);
+  return A;
+}
+
+CanonicalForm
+mulFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+
+  CanonicalForm denA= bCommonDen (A);
+  CanonicalForm denB= bCommonDen (B);
+
+  A *= denA;
+  B *= denB;
+  fmpz_poly_t FLINTA,FLINTB;
+  convertFacCF2Fmpz_poly_t (FLINTA, A);
+  convertFacCF2Fmpz_poly_t (FLINTB, B);
+  fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
+  denA *= denB;
+  A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
+  A /= denA;
+  fmpz_poly_clear (FLINTA);
+  fmpz_poly_clear (FLINTB);
+
+  return A;
+}
+
+/*CanonicalForm
+mulFLINTQ2 (const CanonicalForm& F, const CanonicalForm& G)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+
+  fmpq_poly_t FLINTA,FLINTB;
+  convertFacCF2Fmpq_poly_t (FLINTA, A);
+  convertFacCF2Fmpq_poly_t (FLINTB, B);
+
+  fmpq_poly_mul (FLINTA, FLINTA, FLINTB);
+  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
+  fmpq_poly_clear (FLINTA);
+  fmpq_poly_clear (FLINTB);
+  return A;
+}*/
+
+CanonicalForm
+divFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+
+  fmpq_poly_t FLINTA,FLINTB;
+  convertFacCF2Fmpq_poly_t (FLINTA, A);
+  convertFacCF2Fmpq_poly_t (FLINTB, B);
+
+  fmpq_poly_div (FLINTA, FLINTA, FLINTB);
+  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
+
+  fmpq_poly_clear (FLINTA);
+  fmpq_poly_clear (FLINTB);
+  return A;
+}
+
+CanonicalForm
+modFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+
+  fmpq_poly_t FLINTA,FLINTB;
+  convertFacCF2Fmpq_poly_t (FLINTA, A);
+  convertFacCF2Fmpq_poly_t (FLINTB, B);
+
+  fmpq_poly_rem (FLINTA, FLINTA, FLINTB);
+  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
+
+  fmpq_poly_clear (FLINTA);
+  fmpq_poly_clear (FLINTB);
+  return A;
+}
+
+CanonicalForm
+mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G,
+                 const Variable& alpha, int m)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+
+  CanonicalForm denA= bCommonDen (A);
+  CanonicalForm denB= bCommonDen (B);
+
+  A *= denA;
+  B *= denB;
+
+  int degAa= degree (A, alpha);
+  int degBa= degree (B, alpha);
+  int d= degAa + 1 + degBa;
+
+  fmpz_poly_t FLINTA,FLINTB;
+  fmpz_poly_init (FLINTA);
+  fmpz_poly_init (FLINTB);
+  kronSub (FLINTA, A, d);
+  kronSub (FLINTB, B, d);
+
+  int k= d*m;
+  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
+
+  denA *= denB;
+  A= reverseSubst (FLINTA, d, alpha, denA);
+  fmpz_poly_clear (FLINTA);
+  fmpz_poly_clear (FLINTB);
+  return A;
+}
+
+CanonicalForm
+mulFLINTQTrunc (const CanonicalForm& F, const CanonicalForm& G, int m)
+{
+  if (F.inCoeffDomain() || G.inCoeffDomain())
+    return mod (F*G, power (Variable (1), m));
+  Variable alpha;
+  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+    return mulFLINTQaTrunc (F, G, alpha, m);
+
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+
+  CanonicalForm denA= bCommonDen (A);
+  CanonicalForm denB= bCommonDen (B);
+
+  A *= denA;
+  B *= denB;
+  fmpz_poly_t FLINTA,FLINTB;
+  convertFacCF2Fmpz_poly_t (FLINTA, A);
+  convertFacCF2Fmpz_poly_t (FLINTB, B);
+  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, m);
+  denA *= denB;
+  A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
+  A /= denA;
+  fmpz_poly_clear (FLINTA);
+  fmpz_poly_clear (FLINTB);
+
+  return A;
+}
+
+CanonicalForm uniReverse (const CanonicalForm& F, int d)
+{
+  if (d == 0)
+    return F;
+  if (F.inCoeffDomain())
+    return F*power (Variable (1),d);
+  Variable x= Variable (1);
+  CanonicalForm result= 0;
+  CFIterator i= F;
+  while (d - i.exp() < 0)
+    i++;
+
+  for (; i.hasTerms() && (d - i.exp() >= 0); i++)
+    result += i.coeff()*power (x, d - i.exp());
+  return result;
+}
+
+CanonicalForm
+newtonInverse (const CanonicalForm& F, const int n)
+{
+  int l= ilog2(n);
+
+  CanonicalForm g= F [0];
+
+  ASSERT (!g.isZero(), "expected a unit");
+
+  if (!g.isOne())
+    g = 1/g;
+  Variable x= Variable (1);
+  CanonicalForm result;
+  int exp= 0;
+  if (n & 1)
+  {
+    result= g;
+    exp= 1;
+  }
+  CanonicalForm h;
+
+  for (int i= 1; i <= l; i++)
+  {
+    h= mulNTL (g, mod (F, power (x, (1 << i))));
+    h= mod (h, power (x, (1 << i)) - 1);
+    h= div (h, power (x, (1 << (i - 1))));
+    g -= power (x, (1 << (i - 1)))*
+         mulFLINTQTrunc (g, h, 1 << (i-1));
+
+    if (n & (1 << i))
+    {
+      if (exp)
+      {
+        h= mulNTL (result, mod (F, power (x, exp + (1 << i))));
+        h= mod (h, power (x, exp + (1 << i)) - 1);
+        h= div (h, power (x, exp));
+        result -= power(x, exp)*mulFLINTQTrunc (g, h, 1 << i);
+        exp += (1 << i);
+      }
+      else
+      {
+        exp= (1 << i);
+        result= g;
+      }
+    }
+  }
+
+  return result;
+}
+
+void
+newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+              CanonicalForm& R)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+  Variable x= Variable (1);
+  int degA= degree (A, x);
+  int degB= degree (B, x);
+  int m= degA - degB;
+
+  if (m < 0)
+  {
+    R= A;
+    Q= 0;
+    return;
+  }
+
+  if (degB <= 1)
+    divrem (A, B, Q, R);
+  else
+  {
+    R= uniReverse (A, degA);
+
+    CanonicalForm revB= uniReverse (B, degB);
+    CanonicalForm buf= revB;
+    revB= newtonInverse (revB, m + 1);
+    Q= mulFLINTQTrunc (R, revB, m + 1);
+    Q= uniReverse (Q, m);
+
+    R= A - mulNTL (Q, B);
+  }
+}
+
+void
+newtonDiv (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+  Variable x= Variable (1);
+  int degA= degree (A, x);
+  int degB= degree (B, x);
+  int m= degA - degB;
+
+  if (m < 0)
+  {
+    Q= 0;
+    return;
+  }
+
+  if (degB <= 1)
+    Q= div (A, B);
+  else
+  {
+    CanonicalForm R= uniReverse (A, degA);
+
+    CanonicalForm revB= uniReverse (B, degB);
+    revB= newtonInverse (revB, m + 1);
+    Q= mulFLINTQTrunc (R, revB, m + 1);
+    Q= uniReverse (Q, m);
+  }
+}
+
+#endif
+
+CanonicalForm
+mulNTL (const CanonicalForm& F, const CanonicalForm& G)
+{
+  if (F.inCoeffDomain() || G.inCoeffDomain() || getCharacteristic() == 0)
+  {
+    Variable alpha;
+#ifdef HAVE_FLINT
+    if ((!F.inCoeffDomain() && !G.inCoeffDomain()) &&
+        (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha)))
+    {
+      CanonicalForm result= mulFLINTQa (F, G, alpha);
+      return result;
+    }
+    else if (!F.inCoeffDomain() && !G.inCoeffDomain())
+      return mulFLINTQ (F, G);
+#endif
+    return F*G;
+  }
+  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
+  ASSERT (F.level() == G.level(), "expected polys of same level");
+  if (CFFactory::gettype() == GaloisFieldDomain)
+    return F*G;
+  zz_p::init (getCharacteristic());
+  Variable alpha;
+  CanonicalForm result;
+  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+  {
+    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+    zz_pE::init (NTLMipo);
+    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
+    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
+    mul (NTLF, NTLF, NTLG);
+    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
+  }
+  else
+  {
+#ifdef HAVE_FLINT
+    nmod_poly_t FLINTF, FLINTG;
+    convertFacCF2nmod_poly_t (FLINTF, F);
+    convertFacCF2nmod_poly_t (FLINTG, G);
+    nmod_poly_mul (FLINTF, FLINTF, FLINTG);
+    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
+    nmod_poly_clear (FLINTF);
+    nmod_poly_clear (FLINTG);
+#else
+    zz_pX NTLF= convertFacCF2NTLzzpX (F);
+    zz_pX NTLG= convertFacCF2NTLzzpX (G);
+    mul (NTLF, NTLF, NTLG);
+    result= convertNTLzzpX2CF(NTLF, F.mvar());
+#endif
+  }
+  return result;
+}
+
+CanonicalForm
+modNTL (const CanonicalForm& F, const CanonicalForm& G)
+{
+  if (F.inCoeffDomain() && G.isUnivariate())
+    return F;
+  else if (F.inCoeffDomain() && G.inCoeffDomain())
+    return mod (F, G);
+  else if (F.isUnivariate() && G.inCoeffDomain())
+    return mod (F,G);
+
+  if (getCharacteristic() == 0)
+  {
+#ifdef HAVE_FLINT
+    Variable alpha;
+    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
+      return modFLINTQ (F, G);
+    else
+    {
+      CanonicalForm Q, R;
+      newtonDivrem (F, G, Q, R);
+      return R;
+    }
+#else
+    return mod (F, G);
+#endif
+  }
+
+  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
+  ASSERT (F.level() == G.level(), "expected polys of same level");
+  if (CFFactory::gettype() == GaloisFieldDomain)
+    return mod (F, G);
+  zz_p::init (getCharacteristic());
+  Variable alpha;
+  CanonicalForm result;
+  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+  {
+    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
+    zz_pE::init (NTLMipo);
+    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
+    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
+    rem (NTLF, NTLF, NTLG);
+    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
+  }
+  else
+  {
+#ifdef HAVE_FLINT
+    nmod_poly_t FLINTF, FLINTG;
+    convertFacCF2nmod_poly_t (FLINTF, F);
+    convertFacCF2nmod_poly_t (FLINTG, G);
+    nmod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
+    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
+    nmod_poly_clear (FLINTF);
+    nmod_poly_clear (FLINTG);
+#else
+    zz_pX NTLF= convertFacCF2NTLzzpX (F);
+    zz_pX NTLG= convertFacCF2NTLzzpX (G);
+    rem (NTLF, NTLF, NTLG);
+    result= convertNTLzzpX2CF(NTLF, F.mvar());
+#endif
+  }
+  return result;
+}
+
+CanonicalForm
+divNTL (const CanonicalForm& F, const CanonicalForm& G)
+{
+  if (F.inCoeffDomain() && G.isUnivariate())
+    return F;
+  else if (F.inCoeffDomain() && G.inCoeffDomain())
+    return div (F, G);
+  else if (F.isUnivariate() && G.inCoeffDomain())
+    return div (F,G);
+
+  if (getCharacteristic() == 0)
+  {
+#ifdef HAVE_FLINT
+    Variable alpha;
+    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
+      return divFLINTQ (F,G);
+    else
+    {
+      CanonicalForm Q;
+      newtonDiv (F, G, Q);
+      return Q;
+    }
+#else
+    return div (F, G);
+#endif
+  }
+
+  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
+  ASSERT (F.level() == G.level(), "expected polys of same level");
+  if (CFFactory::gettype() == GaloisFieldDomain)
+    return div (F, G);
+  zz_p::init (getCharacteristic());
+  Variable alpha;
+  CanonicalForm result;
+  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+  {
+    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
+    zz_pE::init (NTLMipo);
+    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
+    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
+    div (NTLF, NTLF, NTLG);
+    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
+  }
+  else
+  {
+#ifdef HAVE_FLINT
+    nmod_poly_t FLINTF, FLINTG;
+    convertFacCF2nmod_poly_t (FLINTF, F);
+    convertFacCF2nmod_poly_t (FLINTG, G);
+    nmod_poly_div (FLINTF, FLINTF, FLINTG);
+    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
+    nmod_poly_clear (FLINTF);
+    nmod_poly_clear (FLINTG);
+#else
+    zz_pX NTLF= convertFacCF2NTLzzpX (F);
+    zz_pX NTLG= convertFacCF2NTLzzpX (G);
+    div (NTLF, NTLF, NTLG);
+    result= convertNTLzzpX2CF(NTLF, F.mvar());
+#endif
+  }
+  return result;
+}
+
+#ifdef HAVE_FLINT
+void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
+{
+  int degAy= degree (A);
+  nmod_poly_init2 (result, getCharacteristic(), d*(degAy + 1));
+
+  nmod_poly_t buf;
+
+  int j, k, bufRepLength;
+  for (CFIterator i= A; i.hasTerms(); i++)
+  {
+    convertFacCF2nmod_poly_t (buf, i.coeff());
+
+    k= i.exp()*d;
+    bufRepLength= (int) nmod_poly_length (buf);
+    for (j= 0; j < bufRepLength; j++)
+      nmod_poly_set_coeff_ui (result, j + k, nmod_poly_get_coeff_ui (buf, j));
+    nmod_poly_clear (buf);
+  }
+  _nmod_poly_normalise (result);
+}
+
+/*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)
+{
+  int degAy= degree (A);
+  fmpq_poly_init2 (result, d1*(degAy + 1));
+
+  fmpq_poly_t buf;
+  fmpq_t coeff;
+
+  int k, l, bufRepLength;
+  CFIterator j;
+  for (CFIterator i= A; i.hasTerms(); i++)
+  {
+    if (i.coeff().inCoeffDomain())
+    {
+      k= d1*i.exp();
+      convertFacCF2Fmpq_poly_t (buf, i.coeff());
+      bufRepLength= (int) fmpq_poly_length(buf);
+      for (l= 0; l < bufRepLength; l++)
+      {
+        fmpq_poly_get_coeff_fmpq (coeff, buf, l);
+        fmpq_poly_set_coeff_fmpq (result, l + k, coeff);
+      }
+      fmpq_poly_clear (buf);
+    }
+    else
+    {
+      for (j= i.coeff(); j.hasTerms(); j++)
+      {
+        k= d1*i.exp();
+        k += d2*j.exp();
+        convertFacCF2Fmpq_poly_t (buf, j.coeff());
+        bufRepLength= (int) fmpq_poly_length(buf);
+        for (l= 0; l < bufRepLength; l++)
+        {
+          fmpq_poly_get_coeff_fmpq (coeff, buf, l);
+          fmpq_poly_set_coeff_fmpq (result, k + l, coeff);
+        }
+        fmpq_poly_clear (buf);
+      }
+    }
+  }
+  fmpq_clear (coeff);
+  _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)
+{
+  int degAy= degree (A);
+  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+  nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
+  nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
+
+  nmod_poly_t buf;
+
+  int k, kk, j, bufRepLength;
+  for (CFIterator i= A; i.hasTerms(); i++)
+  {
+    convertFacCF2nmod_poly_t (buf, i.coeff());
+
+    k= i.exp()*d;
+    kk= (degAy - i.exp())*d;
+    bufRepLength= (int) nmod_poly_length (buf);
+    for (j= 0; j < bufRepLength; j++)
+    {
+      nmod_poly_set_coeff_ui (subA1, j + k,
+                              n_addmod (nmod_poly_get_coeff_ui (subA1, j+k),
+                                        nmod_poly_get_coeff_ui (buf, j),
+                                        getCharacteristic()
+                                       )
+                             );
+      nmod_poly_set_coeff_ui (subA2, j + kk,
+                              n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),
+                                        nmod_poly_get_coeff_ui (buf, j),
+                                        getCharacteristic()
+                                       )
+                             );
+    }
+    nmod_poly_clear (buf);
+  }
+  _nmod_poly_normalise (subA1);
+  _nmod_poly_normalise (subA2);
+}
+
+void
+kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
+                int d)
+{
+  int degAy= degree (A);
+  fmpz_poly_init2 (subA1, d*(degAy + 2));
+  fmpz_poly_init2 (subA2, d*(degAy + 2));
+
+  fmpz_poly_t buf;
+  fmpz_t coeff1, coeff2;
+
+  int k, kk, j, bufRepLength;
+  for (CFIterator i= A; i.hasTerms(); i++)
+  {
+    convertFacCF2Fmpz_poly_t (buf, i.coeff());
+
+    k= i.exp()*d;
+    kk= (degAy - i.exp())*d;
+    bufRepLength= (int) fmpz_poly_length (buf);
+    for (j= 0; j < bufRepLength; j++)
+    {
+      fmpz_poly_get_coeff_fmpz (coeff1, subA1, j+k);
+      fmpz_poly_get_coeff_fmpz (coeff2, buf, j);
+      fmpz_add (coeff1, coeff1, coeff2);
+      fmpz_poly_set_coeff_fmpz (subA1, j + k, coeff1);
+      fmpz_poly_get_coeff_fmpz (coeff1, subA2, j + kk);
+      fmpz_add (coeff1, coeff1, coeff2);
+      fmpz_poly_set_coeff_fmpz (subA2, j + kk, coeff1);
+    }
+    fmpz_poly_clear (buf);
+  }
+  fmpz_clear (coeff1);
+  fmpz_clear (coeff2);
+  _fmpz_poly_normalise (subA1);
+  _fmpz_poly_normalise (subA2);
+}
+
+CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
+{
+  Variable y= Variable (2);
+  Variable x= Variable (1);
+
+  fmpz_poly_t f;
+  fmpz_poly_init (f);
+  fmpz_poly_set (f, F);
+
+  fmpz_poly_t buf;
+  CanonicalForm result= 0;
+  int i= 0;
+  int degf= fmpz_poly_degree(f);
+  int k= 0;
+  int degfSubK, repLength, j;
+  fmpz_t coeff;
+  while (degf >= k)
+  {
+    degfSubK= degf - k;
+    if (degfSubK >= d)
+      repLength= d;
+    else
+      repLength= degfSubK + 1;
+
+    fmpz_poly_init2 (buf, repLength);
+    fmpz_init (coeff);
+    for (j= 0; j < repLength; j++)
+    {
+      fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
+      fmpz_poly_set_coeff_fmpz (buf, j, coeff);
+    }
+    _fmpz_poly_normalise (buf);
+
+    result += convertFmpz_poly_t2FacCF (buf, x)*power (y, i);
+    i++;
+    k= d*i;
+    fmpz_poly_clear (buf);
+    fmpz_clear (coeff);
+  }
+  fmpz_poly_clear (f);
+
+  return result;
+}
+
+CanonicalForm
+reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
+{
+  Variable y= Variable (2);
+  Variable x= Variable (1);
+
+  nmod_poly_t f, g;
+  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
+  nmod_poly_init_preinv (g, getCharacteristic(), ninv);
+  nmod_poly_set (f, F);
+  nmod_poly_set (g, G);
+  int degf= nmod_poly_degree(f);
+  int degg= nmod_poly_degree(g);
+
+
+  nmod_poly_t buf1,buf2, buf3;
+
+  if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
+    nmod_poly_fit_length (f,(long)d*(k+1));
+
+  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;
+  while (degf >= lf || lg >= 0)
+  {
+    if (degfSubLf >= d)
+      repLengthBuf1= d;
+    else if (degfSubLf < 0)
+      repLengthBuf1= 0;
+    else
+      repLengthBuf1= degfSubLf + 1;
+    nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
+
+    for (ind= 0; ind < repLengthBuf1; ind++)
+      nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
+    _nmod_poly_normalise (buf1);
+
+    repLengthBuf1= nmod_poly_length (buf1);
+
+    if (deggSubLg >= d - 1)
+      repLengthBuf2= d - 1;
+    else if (deggSubLg < 0)
+      repLengthBuf2= 0;
+    else
+      repLengthBuf2= deggSubLg + 1;
+
+    nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
+    for (ind= 0; ind < repLengthBuf2; ind++)
+      nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
+
+    _nmod_poly_normalise (buf2);
+    repLengthBuf2= nmod_poly_length (buf2);
+
+    nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
+    for (ind= 0; ind < repLengthBuf1; ind++)
+      nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
+    for (ind= repLengthBuf1; ind < d; ind++)
+      nmod_poly_set_coeff_ui (buf3, ind, 0);
+    for (ind= 0; ind < repLengthBuf2; ind++)
+      nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
+    _nmod_poly_normalise (buf3);
+
+    result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
+    i++;
+
+
+    lf= i*d;
+    degfSubLf= degf - lf;
+
+    lg= d*(k-i);
+    deggSubLg= degg - lg;
+
+    if (lg >= 0 && deggSubLg > 0)
+    {
+      if (repLengthBuf2 > degfSubLf + 1)
+        degfSubLf= repLengthBuf2 - 1;
+      tmp= tmin (repLengthBuf1, deggSubLg + 1);
+      for (ind= 0; ind < tmp; ind++)
+        nmod_poly_set_coeff_ui (g, ind + lg,
+                                n_submod (nmod_poly_get_coeff_ui (g, ind + lg),
+                                          nmod_poly_get_coeff_ui (buf1, ind),
+                                          getCharacteristic()
+                                         )
+                               );
+    }
+    if (lg < 0)
+    {
+      nmod_poly_clear (buf1);
+      nmod_poly_clear (buf2);
+      nmod_poly_clear (buf3);
+      break;
+    }
+    if (degfSubLf >= 0)
+    {
+      for (ind= 0; ind < repLengthBuf2; ind++)
+        nmod_poly_set_coeff_ui (f, ind + lf,
+                                n_submod (nmod_poly_get_coeff_ui (f, ind + lf),
+                                          nmod_poly_get_coeff_ui (buf2, ind),
+                                          getCharacteristic()
+                                         )
+                               );
+    }
+    nmod_poly_clear (buf1);
+    nmod_poly_clear (buf2);
+    nmod_poly_clear (buf3);
+  }
+
+  nmod_poly_clear (f);
+  nmod_poly_clear (g);
+
+  return result;
+}
+
+CanonicalForm
+reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
+{
+  Variable y= Variable (2);
+  Variable x= Variable (1);
+
+  fmpz_poly_t f, g;
+  fmpz_poly_init (f);
+  fmpz_poly_init (g);
+  fmpz_poly_set (f, F);
+  fmpz_poly_set (g, G);
+  int degf= fmpz_poly_degree(f);
+  int degg= fmpz_poly_degree(g);
+
+
+  fmpz_poly_t buf1,buf2, buf3;
+
+  if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
+    fmpz_poly_fit_length (f,(long)d*(k+1));
+
+  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;
+  fmpz_t tmp1, tmp2;
+  while (degf >= lf || lg >= 0)
+  {
+    if (degfSubLf >= d)
+      repLengthBuf1= d;
+    else if (degfSubLf < 0)
+      repLengthBuf1= 0;
+    else
+      repLengthBuf1= degfSubLf + 1;
+
+    fmpz_poly_init2 (buf1, repLengthBuf1);
+
+    for (ind= 0; ind < repLengthBuf1; ind++)
+    {
+      fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
+      fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
+    }
+    _fmpz_poly_normalise (buf1);
+
+    repLengthBuf1= fmpz_poly_length (buf1);
+
+    if (deggSubLg >= d - 1)
+      repLengthBuf2= d - 1;
+    else if (deggSubLg < 0)
+      repLengthBuf2= 0;
+    else
+      repLengthBuf2= deggSubLg + 1;
+
+    fmpz_poly_init2 (buf2, repLengthBuf2);
+
+    for (ind= 0; ind < repLengthBuf2; ind++)
+    {
+      fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
+      fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
+    }
+
+    _fmpz_poly_normalise (buf2);
+    repLengthBuf2= fmpz_poly_length (buf2);
+
+    fmpz_poly_init2 (buf3, repLengthBuf2 + d);
+    for (ind= 0; ind < repLengthBuf1; ind++)
+    {
+      fmpz_poly_get_coeff_fmpz (tmp1, buf1, ind);  //oder fmpz_set (fmpz_poly_get_coeff_ptr (buf3, ind),fmpz_poly_get_coeff_ptr (buf1, ind))
+      fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
+    }
+    for (ind= repLengthBuf1; ind < d; ind++)
+      fmpz_poly_set_coeff_ui (buf3, ind, 0);
+    for (ind= 0; ind < repLengthBuf2; ind++)
+    {
+      fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
+      fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
+    }
+    _fmpz_poly_normalise (buf3);
+
+    result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
+    i++;
+
+
+    lf= i*d;
+    degfSubLf= degf - lf;
+
+    lg= d*(k-i);
+    deggSubLg= degg - lg;
+
+    if (lg >= 0 && deggSubLg > 0)
+    {
+      if (repLengthBuf2 > degfSubLf + 1)
+        degfSubLf= repLengthBuf2 - 1;
+      tmp= tmin (repLengthBuf1, deggSubLg + 1);
+      for (ind= 0; ind < tmp; ind++)
+      {
+        fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
+        fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
+        fmpz_sub (tmp1, tmp1, tmp2);
+        fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
+      }
+    }
+    if (lg < 0)
+    {
+      fmpz_poly_clear (buf1);
+      fmpz_poly_clear (buf2);
+      fmpz_poly_clear (buf3);
+      break;
+    }
+    if (degfSubLf >= 0)
+    {
+      for (ind= 0; ind < repLengthBuf2; ind++)
+      {
+        fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
+        fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
+        fmpz_sub (tmp1, tmp1, tmp2);
+        fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
+      }
+    }
+    fmpz_poly_clear (buf1);
+    fmpz_poly_clear (buf2);
+    fmpz_poly_clear (buf3);
+  }
+
+  fmpz_poly_clear (f);
+  fmpz_poly_clear (g);
+  fmpz_clear (tmp1);
+  fmpz_clear (tmp2);
+
+  return result;
+}
+
+CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
+{
+  Variable y= Variable (2);
+  Variable x= Variable (1);
+
+  nmod_poly_t f;
+  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
+  nmod_poly_set (f, F);
+
+  nmod_poly_t buf;
+  CanonicalForm result= 0;
+  int i= 0;
+  int degf= nmod_poly_degree(f);
+  int k= 0;
+  int degfSubK, repLength, j;
+  while (degf >= k)
+  {
+    degfSubK= degf - k;
+    if (degfSubK >= d)
+      repLength= d;
+    else
+      repLength= degfSubK + 1;
+
+    nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
+    for (j= 0; j < repLength; j++)
+      nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
+    _nmod_poly_normalise (buf);
+
+    result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
+    i++;
+    k= d*i;
+    nmod_poly_clear (buf);
+  }
+  nmod_poly_clear (f);
+
+  return result;
+}
+
+CanonicalForm
+mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
+                    CanonicalForm& M)
+{
+  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
+  d1 /= 2;
+  d1 += 1;
+
+  nmod_poly_t F1, F2;
+  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+  nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
+  nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
+  kronSubRecipro (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);
+
+  int k= d1*degree (M);
+  nmod_poly_mullow (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= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG
+         + d1*(2+taildegF + taildegG);
+  nmod_poly_mulhigh (F2, F2, G2, b);
+  nmod_poly_shift_right (F2, F2, b);
+  int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
+
+
+  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
+
+  nmod_poly_clear (F1);
+  nmod_poly_clear (F2);
+  nmod_poly_clear (G1);
+  nmod_poly_clear (G2);
+  return result;
+}
+
+CanonicalForm
+mulMod2FLINTFp (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 mulMod2FLINTFpReci (A, B, M);
+
+  nmod_poly_t FLINTA, FLINTB;
+  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+  nmod_poly_init_preinv (FLINTA, getCharacteristic(), ninv);
+  nmod_poly_init_preinv (FLINTB, getCharacteristic(), ninv);
+  kronSubFp (FLINTA, A, d1);
+  kronSubFp (FLINTB, B, d1);
+
+  int k= d1*degree (M);
+  nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
+
+  A= reverseSubstFp (FLINTA, d1);
+
+  nmod_poly_clear (FLINTA);
+  nmod_poly_clear (FLINTB);
+  return A;
+}
+
+CanonicalForm
+mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
+                    CanonicalForm& M)
+{
+  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
+  d1 /= 2;
+  d1 += 1;
+
+  fmpz_poly_t F1, F2;
+  fmpz_poly_init (F1);
+  fmpz_poly_init (F2);
+  kronSubRecipro (F1, F2, F, d1);
+
+  fmpz_poly_t G1, G2;
+  fmpz_poly_init (G1);
+  fmpz_poly_init (G2);
+  kronSubRecipro (G1, G2, G, d1);
+
+  int k= d1*degree (M);
+  fmpz_poly_mullow (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= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG
+         + d1*(2+taildegF + taildegG);
+  fmpz_poly_mulhigh_n (F2, F2, G2, b);
+  fmpz_poly_shift_right (F2, F2, b);
+  int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
+
+  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
+
+  fmpz_poly_clear (F1);
+  fmpz_poly_clear (F2);
+  fmpz_poly_clear (G1);
+  fmpz_poly_clear (G2);
+  return result;
+}
+
+CanonicalForm
+mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
+                CanonicalForm& M)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+
+  int degAx= degree (A, 1);
+  int degBx= degree (B, 1);
+  int d1= degAx + 1 + degBx;
+
+  CanonicalForm f= bCommonDen (F);
+  CanonicalForm g= bCommonDen (G);
+  A *= f;
+  B *= g;
+
+  fmpz_poly_t FLINTA, FLINTB;
+  fmpz_poly_init (FLINTA);
+  fmpz_poly_init (FLINTB);
+  kronSub (FLINTA, A, d1);
+  kronSub (FLINTB, B, d1);
+  int k= d1*degree (M);
+
+  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
+  A= reverseSubstQ (FLINTA, d1);
+  fmpz_poly_clear (FLINTA);
+  fmpz_poly_clear (FLINTB);
+  return A/(f*g);
+}
+
+#endif
+
+CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
+                       const CanonicalForm& M)
+{
+  if (A.isZero() || B.isZero())
+    return 0;
+
+  ASSERT (M.isUnivariate(), "M must be univariate");
+
+  CanonicalForm F= mod (A, M);
+  CanonicalForm G= mod (B, M);
+  if (F.inCoeffDomain() || G.inCoeffDomain())
+    return F*G;
+  Variable y= M.mvar();
+  int degF= degree (F, y);
+  int degG= degree (G, y);
+
+  if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
+      (F.level() == G.level()))
+  {
+    CanonicalForm result= mulNTL (F, G);
+    return mod (result, M);
+  }
+  else if (degF <= 1 && degG <= 1)
+  {
+    CanonicalForm result= F*G;
+    return mod (result, M);
+  }
+
+  int sizeF= size (F);
+  int sizeG= size (G);
+
+  int fallBackToNaive= 50;
+  if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
+    return mod (F*G, M);
+
+#ifdef HAVE_FLINT
+  Variable alpha;
+  if (getCharacteristic() == 0 && !hasFirstAlgVar (F, alpha) && ! hasFirstAlgVar (G, alpha))
+    return mulMod2FLINTQ (F, G, M);
+#endif
+
+  if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
+      (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
+    return mulMod2NTLFq (F, G, M);
+
+  int m= (int) ceil (degree (M)/2.0);
+  if (degF >= m || degG >= m)
+  {
+    CanonicalForm MLo= power (y, m);
+    CanonicalForm MHi= power (y, degree (M) - m);
+    CanonicalForm F0= mod (F, MLo);
+    CanonicalForm F1= div (F, MLo);
+    CanonicalForm G0= mod (G, MLo);
+    CanonicalForm G1= div (G, MLo);
+    CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
+    CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
+    CanonicalForm F0G0= mulMod2 (F0, G0, M);
+    return F0G0 + MLo*(F0G1 + F1G0);
+  }
+  else
+  {
+    m= (int) ceil (tmax (degF, degG)/2.0);
+    CanonicalForm yToM= power (y, m);
+    CanonicalForm F0= mod (F, yToM);
+    CanonicalForm F1= div (F, yToM);
+    CanonicalForm G0= mod (G, yToM);
+    CanonicalForm G1= div (G, yToM);
+    CanonicalForm H00= mulMod2 (F0, G0, M);
+    CanonicalForm H11= mulMod2 (F1, G1, M);
+    CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
+    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
+  }
+  DEBOUTLN (cerr, "fatal end in mulMod2");
+}
+
+CanonicalForm mod (const CanonicalForm& F, const CFList& M)
+{
+  CanonicalForm A= F;
+  for (CFListIterator i= M; i.hasItem(); i++)
+    A= mod (A, i.getItem());
+  return A;
+}
+
+CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
+                      const CFList& MOD)
+{
+  if (A.isZero() || B.isZero())
+    return 0;
+
+  if (MOD.length() == 1)
+    return mulMod2 (A, B, MOD.getLast());
+
+  CanonicalForm M= MOD.getLast();
+  CanonicalForm F= mod (A, M);
+  CanonicalForm G= mod (B, M);
+  if (F.inCoeffDomain() || G.inCoeffDomain())
+    return F*G;
+  Variable y= M.mvar();
+  int degF= degree (F, y);
+  int degG= degree (G, y);
+
+  if ((degF <= 1 && F.level() <= M.level()) &&
+      (degG <= 1 && G.level() <= M.level()))
+  {
+    CFList buf= MOD;
+    buf.removeLast();
+    if (degF == 1 && degG == 1)
+    {
+      CanonicalForm F0= mod (F, y);
+      CanonicalForm F1= div (F, y);
+      CanonicalForm G0= mod (G, y);
+      CanonicalForm G1= div (G, y);
+      if (degree (M) > 2)
+      {
+        CanonicalForm H00= mulMod (F0, G0, buf);
+        CanonicalForm H11= mulMod (F1, G1, buf);
+        CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
+        return H11*y*y + (H01 - H00 - H11)*y + H00;
+      }
+      else //here degree (M) == 2
+      {
+        buf.append (y);
+        CanonicalForm F0G1= mulMod (F0, G1, buf);
+        CanonicalForm F1G0= mulMod (F1, G0, buf);
+        CanonicalForm F0G0= mulMod (F0, G0, MOD);
+        CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
+        return result;
+      }
+    }
+    else if (degF == 1 && degG == 0)
+      return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
+    else if (degF == 0 && degG == 1)
+      return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
+    else
+      return mulMod (F, G, buf);
+  }
+  int m= (int) ceil (degree (M)/2.0);
+  if (degF >= m || degG >= m)
+  {
+    CanonicalForm MLo= power (y, m);
+    CanonicalForm MHi= power (y, degree (M) - m);
+    CanonicalForm F0= mod (F, MLo);
+    CanonicalForm F1= div (F, MLo);
+    CanonicalForm G0= mod (G, MLo);
+    CanonicalForm G1= div (G, MLo);
+    CFList buf= MOD;
+    buf.removeLast();
+    buf.append (MHi);
+    CanonicalForm F0G1= mulMod (F0, G1, buf);
+    CanonicalForm F1G0= mulMod (F1, G0, buf);
+    CanonicalForm F0G0= mulMod (F0, G0, MOD);
+    return F0G0 + MLo*(F0G1 + F1G0);
+  }
+  else
+  {
+    m= (int) ceil (tmax (degF, degG)/2.0);
+    CanonicalForm yToM= power (y, m);
+    CanonicalForm F0= mod (F, yToM);
+    CanonicalForm F1= div (F, yToM);
+    CanonicalForm G0= mod (G, yToM);
+    CanonicalForm G1= div (G, yToM);
+    CanonicalForm H00= mulMod (F0, G0, MOD);
+    CanonicalForm H11= mulMod (F1, G1, MOD);
+    CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
+    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
+  }
+  DEBOUTLN (cerr, "fatal end in mulMod");
+}
+
+CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
+{
+  if (L.isEmpty())
+    return 1;
+  int l= L.length();
+  if (l == 1)
+    return mod (L.getFirst(), M);
+  else if (l == 2) {
+    CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
+    return result;
+  }
+  else
+  {
+    l /= 2;
+    CFList tmp1, tmp2;
+    CFListIterator i= L;
+    CanonicalForm buf1, buf2;
+    for (int j= 1; j <= l; j++, i++)
+      tmp1.append (i.getItem());
+    tmp2= Difference (L, tmp1);
+    buf1= prodMod (tmp1, M);
+    buf2= prodMod (tmp2, M);
+    CanonicalForm result= mulMod2 (buf1, buf2, M);
+    return result;
+  }
+}
+
+CanonicalForm prodMod (const CFList& L, const CFList& M)
+{
+  if (L.isEmpty())
+    return 1;
+  else if (L.length() == 1)
+    return L.getFirst();
+  else if (L.length() == 2)
+    return mulMod (L.getFirst(), L.getLast(), M);
+  else
+  {
+    int l= L.length()/2;
+    CFListIterator i= L;
+    CFList tmp1, tmp2;
+    CanonicalForm buf1, buf2;
+    for (int j= 1; j <= l; j++, i++)
+      tmp1.append (i.getItem());
+    tmp2= Difference (L, tmp1);
+    buf1= prodMod (tmp1, M);
+    buf2= prodMod (tmp2, M);
+    return mulMod (buf1, buf2, M);
+  }
+}
+
+CanonicalForm reverse (const CanonicalForm& F, int d)
+{
+  if (d == 0)
+    return F;
+  CanonicalForm A= F;
+  Variable y= Variable (2);
+  Variable x= Variable (1);
+  if (degree (A, x) > 0)
+  {
+    A= swapvar (A, x, y);
+    CanonicalForm result= 0;
+    CFIterator i= A;
+    while (d - i.exp() < 0)
+      i++;
+
+    for (; i.hasTerms() && (d - i.exp() >= 0); i++)
+      result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
+    return result;
+  }
+  else
+    return A*power (x, d);
+}
+
+CanonicalForm
+newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
+{
+  int l= ilog2(n);
+
+  CanonicalForm g= mod (F, M)[0] [0];
+
+  ASSERT (!g.isZero(), "expected a unit");
+
+  Variable alpha;
+
+  if (!g.isOne())
+    g = 1/g;
+  Variable x= Variable (1);
+  CanonicalForm result;
+  int exp= 0;
+  if (n & 1)
+  {
+    result= g;
+    exp= 1;
+  }
+  CanonicalForm h;
+
+  for (int i= 1; i <= l; i++)
+  {
+    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
+    h= mod (h, power (x, (1 << i)) - 1);
+    h= div (h, power (x, (1 << (i - 1))));
+    h= mod (h, M);
+    g -= power (x, (1 << (i - 1)))*
+         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
+
+    if (n & (1 << i))
+    {
+      if (exp)
+      {
+        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
+        h= mod (h, power (x, exp + (1 << i)) - 1);
+        h= div (h, power (x, exp));
+        h= mod (h, M);
+        result -= power(x, exp)*mod (mulMod2 (g, h, M),
+                                       power (x, (1 << i)));
+        exp += (1 << i);
+      }
+      else
+      {
+        exp= (1 << i);
+        result= g;
+      }
+    }
+  }
+
+  return result;
+}
+
+CanonicalForm
+newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
+           M)
+{
+  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
+  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
+
+  CanonicalForm A= mod (F, M);
+  CanonicalForm B= mod (G, M);
+
+  Variable x= Variable (1);
+  int degA= degree (A, x);
+  int degB= degree (B, x);
+  int m= degA - degB;
+  if (m < 0)
+    return 0;
+
+  Variable v;
+  CanonicalForm Q;
+  if (degB < 1 || CFFactory::gettype() == GaloisFieldDomain)
+  {
+    CanonicalForm R;
+    divrem2 (A, B, Q, R, M);
+  }
+  else
+  {
+    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
+    {
+      CanonicalForm R= reverse (A, degA);
+      CanonicalForm revB= reverse (B, degB);
+      revB= newtonInverse (revB, m + 1, M);
+      Q= mulMod2 (R, revB, M);
+      Q= mod (Q, power (x, m + 1));
+      Q= reverse (Q, m);
+    }
+    else
+    {
+      zz_pX mipo= convertFacCF2NTLzzpX (M);
+      Variable y= Variable (2);
+      zz_pEX NTLA, NTLB;
+      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
+      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
+      div (NTLA, NTLA, NTLB);
+      Q= convertNTLzz_pEX2CF (NTLA, x, y);
+    }
+  }
+
+  return Q;
+}
+
+void
+newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+              CanonicalForm& R, const CanonicalForm& M)
+{
+  CanonicalForm A= mod (F, M);
+  CanonicalForm B= mod (G, M);
+  Variable x= Variable (1);
+  int degA= degree (A, x);
+  int degB= degree (B, x);
+  int m= degA - degB;
+
+  if (m < 0)
+  {
+    R= A;
+    Q= 0;
+    return;
+  }
+
+  Variable v;
+  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
+  {
+     divrem2 (A, B, Q, R, M);
+  }
+  else
+  {
+    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
+    {
+      R= reverse (A, degA);
+
+      CanonicalForm revB= reverse (B, degB);
+      revB= newtonInverse (revB, m + 1, M);
+      Q= mulMod2 (R, revB, M);
+
+      Q= mod (Q, power (x, m + 1));
+      Q= reverse (Q, m);
+
+      R= A - mulMod2 (Q, B, M);
+    }
+    else
+    {
+      zz_pX mipo= convertFacCF2NTLzzpX (M);
+      Variable y= Variable (2);
+      zz_pEX NTLA, NTLB;
+      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
+      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
+      zz_pEX NTLQ, NTLR;
+      DivRem (NTLQ, NTLR, NTLA, NTLB);
+      Q= convertNTLzz_pEX2CF (NTLQ, x, y);
+      R= convertNTLzz_pEX2CF (NTLR, x, y);
+    }
+  }
+}
+
+static inline
+CFList split (const CanonicalForm& F, const int m, const Variable& x)
+{
+  CanonicalForm A= F;
+  CanonicalForm buf= 0;
+  bool swap= false;
+  if (degree (A, x) <= 0)
+    return CFList(A);
+  else if (x.level() != A.level())
+  {
+    swap= true;
+    A= swapvar (A, x, A.mvar());
+  }
+
+  int j= (int) floor ((double) degree (A)/ m);
+  CFList result;
+  CFIterator i= A;
+  for (; j >= 0; j--)
+  {
+    while (i.hasTerms() && i.exp() - j*m >= 0)
+    {
+      if (swap)
+        buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
+      else
+        buf += i.coeff()*power (x, i.exp() - j*m);
+      i++;
+    }
+    if (swap)
+      result.append (swapvar (buf, x, F.mvar()));
+    else
+      result.append (buf);
+    buf= 0;
+  }
+  return result;
+}
+
+static inline
+void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+               CanonicalForm& R, const CFList& M);
+
+static inline
+void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+               CanonicalForm& R, const CFList& M)
+{
+  CanonicalForm A= mod (F, M);
+  CanonicalForm B= mod (G, M);
+  Variable x= Variable (1);
+  int degB= degree (B, x);
+  int degA= degree (A, x);
+  if (degA < degB)
+  {
+    Q= 0;
+    R= A;
+    return;
+  }
+  ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
+  if (degB < 1)
+  {
+    divrem (A, B, Q, R);
+    Q= mod (Q, M);
+    R= mod (R, M);
+    return;
+  }
+
+  int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
+  CFList splitA= split (A, m, x);
+  if (splitA.length() == 3)
+    splitA.insert (0);
+  if (splitA.length() == 2)
+  {
+    splitA.insert (0);
+    splitA.insert (0);
+  }
+  if (splitA.length() == 1)
+  {
+    splitA.insert (0);
+    splitA.insert (0);
+    splitA.insert (0);
+  }
+
+  CanonicalForm xToM= power (x, m);
+
+  CFListIterator i= splitA;
+  CanonicalForm H= i.getItem();
+  i++;
+  H *= xToM;
+  H += i.getItem();
+  i++;
+  H *= xToM;
+  H += i.getItem();
+  i++;
+
+  divrem32 (H, B, Q, R, M);
+
+  CFList splitR= split (R, m, x);
+  if (splitR.length() == 1)
+    splitR.insert (0);
+
+  H= splitR.getFirst();
+  H *= xToM;
+  H += splitR.getLast();
+  H *= xToM;
+  H += i.getItem();
+
+  CanonicalForm bufQ;
+  divrem32 (H, B, bufQ, R, M);
+
+  Q *= xToM;
+  Q += bufQ;
+  return;
+}
+
+static inline
+void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+               CanonicalForm& R, const CFList& M)
+{
+  CanonicalForm A= mod (F, M);
+  CanonicalForm B= mod (G, M);
+  Variable x= Variable (1);
+  int degB= degree (B, x);
+  int degA= degree (A, x);
+  if (degA < degB)
+  {
+    Q= 0;
+    R= A;
+    return;
+  }
+  ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
+  if (degB < 1)
+  {
+    divrem (A, B, Q, R);
+    Q= mod (Q, M);
+    R= mod (R, M);
+    return;
+  }
+  int m= (int) ceil ((double) (degB + 1)/ 2.0);
+
+  CFList splitA= split (A, m, x);
+  CFList splitB= split (B, m, x);
+
+  if (splitA.length() == 2)
+  {
+    splitA.insert (0);
+  }
+  if (splitA.length() == 1)
+  {
+    splitA.insert (0);
+    splitA.insert (0);
+  }
+  CanonicalForm xToM= power (x, m);
+
+  CanonicalForm H;
+  CFListIterator i= splitA;
+  i++;
+
+  if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
+  {
+    H= splitA.getFirst()*xToM + i.getItem();
+    divrem21 (H, splitB.getFirst(), Q, R, M);
+  }
+  else
+  {
+    R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
+       splitB.getFirst()*xToM;
+    Q= xToM - 1;
+  }
+
+  H= mulMod (Q, splitB.getLast(), M);
+
+  R= R*xToM + splitA.getLast() - H;
+
+  while (degree (R, x) >= degB)
+  {
+    xToM= power (x, degree (R, x) - degB);
+    Q += LC (R, x)*xToM;
+    R -= mulMod (LC (R, x), B, M)*xToM;
+    Q= mod (Q, M);
+    R= mod (R, M);
+  }
+
+  return;
+}
+
+void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+              CanonicalForm& R, const CanonicalForm& M)
+{
+  CanonicalForm A= mod (F, M);
+  CanonicalForm B= mod (G, M);
+
+  if (B.inCoeffDomain())
+  {
+    divrem (A, B, Q, R);
+    return;
+  }
+  if (A.inCoeffDomain() && !B.inCoeffDomain())
+  {
+    Q= 0;
+    R= A;
+    return;
+  }
+
+  if (B.level() < A.level())
+  {
+    divrem (A, B, Q, R);
+    return;
+  }
+  if (A.level() > B.level())
+  {
+    R= A;
+    Q= 0;
+    return;
+  }
+  if (B.level() == 1 && B.isUnivariate())
+  {
+    divrem (A, B, Q, R);
+    return;
+  }
+  if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
+  {
+    Q= 0;
+    R= A;
+    return;
+  }
+
+  Variable x= Variable (1);
+  int degB= degree (B, x);
+  if (degB > degree (A, x))
+  {
+    Q= 0;
+    R= A;
+    return;
+  }
+
+  CFList splitA= split (A, degB, x);
+
+  CanonicalForm xToDegB= power (x, degB);
+  CanonicalForm H, bufQ;
+  Q= 0;
+  CFListIterator i= splitA;
+  H= i.getItem()*xToDegB;
+  i++;
+  H += i.getItem();
+  CFList buf;
+  while (i.hasItem())
+  {
+    buf= CFList (M);
+    divrem21 (H, B, bufQ, R, buf);
+    i++;
+    if (i.hasItem())
+      H= R*xToDegB + i.getItem();
+    Q *= xToDegB;
+    Q += bufQ;
+  }
+  return;
+}
+
+void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+             CanonicalForm& R, const CFList& MOD)
+{
+  CanonicalForm A= mod (F, MOD);
+  CanonicalForm B= mod (G, MOD);
+  Variable x= Variable (1);
+  int degB= degree (B, x);
+  if (degB > degree (A, x))
+  {
+    Q= 0;
+    R= A;
+    return;
+  }
+
+  if (degB <= 0)
+  {
+    divrem (A, B, Q, R);
+    Q= mod (Q, MOD);
+    R= mod (R, MOD);
+    return;
+  }
+  CFList splitA= split (A, degB, x);
+
+  CanonicalForm xToDegB= power (x, degB);
+  CanonicalForm H, bufQ;
+  Q= 0;
+  CFListIterator i= splitA;
+  H= i.getItem()*xToDegB;
+  i++;
+  H += i.getItem();
+  while (i.hasItem())
+  {
+    divrem21 (H, B, bufQ, R, MOD);
+    i++;
+    if (i.hasItem())
+      H= R*xToDegB + i.getItem();
+    Q *= xToDegB;
+    Q += bufQ;
+  }
+  return;
+}
+
+#endif
diff --git a/factory/facMul.h b/factory/facMul.h
new file mode 100644
index 0000000..76337e0
--- /dev/null
+++ b/factory/facMul.h
@@ -0,0 +1,147 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR
+\*****************************************************************************/
+/** @file facHensel.h
+ *
+ * This file defines functions for fast multiplication and division with
+ * remainder
+ *
+ * @author Martin Lee
+ *
+ **/
+/*****************************************************************************/
+
+#ifndef FAC_MUL_H
+#define FAC_MUL_H
+
+#include "canonicalform.h"
+
+/// multiplication of univariate polys over a finite field using NTL, if we are
+/// in GF factory's default multiplication is used.
+///
+/// @return @a mulNTL returns F*G
+CanonicalForm
+mulNTL (const CanonicalForm& F, ///< [in] a univariate poly
+        const CanonicalForm& G  ///< [in] a univariate poly
+       );
+
+/// mod of univariate polys over a finite field using NTL, if we are
+/// in GF factory's default mod is used.
+///
+/// @return @a modNTL returns F mod G
+CanonicalForm
+modNTL (const CanonicalForm& F, ///< [in] a univariate poly
+        const CanonicalForm& G  ///< [in] a univariate poly
+       );
+
+/// division of univariate polys over a finite field using NTL, if we are
+/// in GF factory's default division is used.
+///
+/// @return @a divNTL returns F/G
+CanonicalForm
+divNTL (const CanonicalForm& F, ///< [in] a univariate poly
+        const CanonicalForm& G  ///< [in] a univariate poly
+       );
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a M.
+///
+/// @return @a Q returns the dividend, @a R returns the remainder.
+/// @sa divrem()
+void divrem2 (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
+              const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
+              CanonicalForm& Q,       ///< [in,out] dividend
+              CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
+                                      ///< degree (G, 1)
+              const CanonicalForm& M  ///< [in] power of Variable (2)
+             );
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a MOD.
+///
+/// @sa divrem2()
+void divrem (
+           const CanonicalForm& F, ///< [in] multivariate, compressed polynomial
+           const CanonicalForm& G, ///< [in] multivariate, compressed polynomial
+           CanonicalForm& Q,       ///< [in,out] dividend
+           CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
+                                   ///< degree (G, 1)
+           const CFList& MOD       ///< [in] only contains powers of
+                                   ///< Variables of level higher than 1
+            );
+
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a M using Newton inversion
+///
+/// @return @a Q returns the dividend, @a R returns the remainder.
+/// @sa divrem2(), newtonDiv()
+void
+newtonDivrem (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
+              const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
+                                      ///< which is monic in Variable (1)
+              CanonicalForm& Q,       ///< [in,out] dividend
+              CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
+                                      ///< degree (G, 1)
+              const CanonicalForm& M  ///< [in] power of Variable (2)
+             );
+
+/// division of @a F by
+/// @a G wrt Variable (1) modulo @a M using Newton inversion
+///
+/// @return @a newtonDiv returns the dividend
+/// @sa divrem2(), newtonDivrem()
+CanonicalForm
+newtonDiv (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
+           const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
+                                   ///< which is monic in Variable (1)
+           const CanonicalForm& M  ///< [in] power of Variable (2)
+          );
+
+/// Karatsuba style modular multiplication for bivariate polynomials.
+///
+/// @return @a mulMod2 returns @a A * @a B mod @a M.
+CanonicalForm
+mulMod2 (const CanonicalForm& A, ///< [in] bivariate, compressed polynomial
+         const CanonicalForm& B, ///< [in] bivariate, compressed polynomial
+         const CanonicalForm& M  ///< [in] power of Variable (2)
+        );
+
+/// Karatsuba style modular multiplication for multivariate polynomials.
+///
+/// @return @a mulMod2 returns @a A * @a B mod @a MOD.
+CanonicalForm
+mulMod (const CanonicalForm& A, ///< [in] multivariate, compressed polynomial
+        const CanonicalForm& B, ///< [in] multivariate, compressed polynomial
+        const CFList& MOD       ///< [in] only contains powers of
+                                ///< Variables of level higher than 1
+       );
+
+/// reduce @a F modulo elements in @a M.
+///
+/// @return @a mod returns @a F modulo @a M
+CanonicalForm mod (const CanonicalForm& F, ///< [in] compressed polynomial
+                   const CFList& M         ///< [in] list containing only
+                                           ///< univariate polynomials
+                  );
+
+/// product of all elements in @a L modulo @a M via divide-and-conquer.
+///
+/// @return @a prodMod returns product of all elements in @a L modulo @a M.
+CanonicalForm
+prodMod (const CFList& L,       ///< [in] contains only bivariate, compressed
+                                ///< polynomials
+         const CanonicalForm& M ///< [in] power of Variable (2)
+        );
+
+/// product of all elements in @a L modulo @a M via divide-and-conquer.
+///
+/// @return @a prodMod returns product of all elements in @a L modulo @a M.
+CanonicalForm
+prodMod (const CFList& L, ///< [in] contains multivariate, compressed
+                          ///< polynomials
+         const CFList& M  ///< [in] contains only powers of Variables
+        );
+#endif
+/* FAC_MUL_H */
+

-- 
an open source computer algebra system



More information about the debian-science-commits mailing list