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

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


The following commit has been merged in the cleanedupstream branch:
commit 04290e40560003ed5ac61a60f046642627c79e35
Author: Martin Lee <martinlee84 at web.de>
Date:   Thu Feb 9 17:59:15 2012 +0100

    chg: more reorganization

diff --git a/factory/cf_gcd_smallp.cc b/factory/cf_gcd_smallp.cc
index d4e6bc9..e79b205 100755
--- a/factory/cf_gcd_smallp.cc
+++ b/factory/cf_gcd_smallp.cc
@@ -22,6 +22,7 @@
 #include "timing.h"
 
 #include "canonicalform.h"
+#include "algext.h"
 #include "cf_map.h"
 #include "cf_util.h"
 #include "templates/ftmpl_functions.h"
@@ -4026,8 +4027,8 @@ int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
   CFArray lcs= CFArray (2);
   lcs [0]= shiftedLCsEval1.getFirst();
   lcs [1]= shiftedLCsEval2.getFirst();
-  henselLift122 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
-                 lcs, false);
+  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
+                        lcs, false);
 
   for (CFListIterator i= factors; i.hasItem(); i++)
   {
@@ -4043,9 +4044,9 @@ int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
     liftBounds[0]= liftBound;
     for (int i= 1; i < U.level() - 1; i++)
       liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
-    factors= henselLift2 (UEval, factors, liftBounds, U.level() - 1, false,
-                          shiftedLCsEval1, shiftedLCsEval2, Pi, diophant,
-                          noOneToOne);
+    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
+                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
+                                  diophant, noOneToOne);
     delete [] liftBounds;
     if (noOneToOne)
       return 0;
diff --git a/factory/facFactorize.cc b/factory/facFactorize.cc
index e75e731..a7b94c0 100644
--- a/factory/facFactorize.cc
+++ b/factory/facFactorize.cc
@@ -477,8 +477,8 @@ precomputeLeadingCoeff (const CanonicalForm& LCF, const CFList& LCFFactors,
       CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
       CFArray Pi;
       CFList diophant;
-      henselLift122 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
-                     leadingCoeffs, false);
+      nonMonicHenselLift12 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
+                            leadingCoeffs, false);
 
       if (sqrfPartF.level() > 2)
       {
diff --git a/factory/facFqBivarUtil.cc b/factory/facFqBivarUtil.cc
index 9f48bd4..a596281 100644
--- a/factory/facFqBivarUtil.cc
+++ b/factory/facFqBivarUtil.cc
@@ -13,6 +13,7 @@
 #include <config.h>
 
 #include "cf_map.h"
+#include "algext.h"
 #include "cf_map_ext.h"
 #include "templates/ftmpl_functions.h"
 #include "ExtensionInfo.h"
diff --git a/factory/facFqFactorize.cc b/factory/facFqFactorize.cc
index e6351d0..dacc4c9 100644
--- a/factory/facFqFactorize.cc
+++ b/factory/facFqFactorize.cc
@@ -1655,8 +1655,8 @@ precomputeLeadingCoeff (const CanonicalForm& LCF, const CFList& LCFFactors,
       CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
       CFArray Pi;
       CFList diophant;
-      henselLift122 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
-                     leadingCoeffs, false);
+      nonMonicHenselLift12 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
+                            leadingCoeffs, false);
 
       if (sqrfPartF.level() > 2)
       {
diff --git a/factory/facHensel.cc b/factory/facHensel.cc
index bde2c1c..22fe683 100644
--- a/factory/facHensel.cc
+++ b/factory/facHensel.cc
@@ -20,9 +20,9 @@
 #include "debug.h"
 #include "timing.h"
 
+#include "algext.h"
 #include "facHensel.h"
 #include "facMul.h"
-#include "cf_util.h"
 #include "fac_util.h"
 #include "cf_algorithm.h"
 #include "cf_primes.h"
@@ -120,7 +120,7 @@ void tryDiophantine (CFList& result, const CanonicalForm& F,
   }
 }
 
-static inline
+static
 CFList mapinto (const CFList& L)
 {
   CFList result;
@@ -129,7 +129,7 @@ CFList mapinto (const CFList& L)
   return result;
 }
 
-static inline
+static
 int mod (const CFList& L, const CanonicalForm& p)
 {
   for (CFListIterator i= L; i.hasItem(); i++)
@@ -141,7 +141,7 @@ int mod (const CFList& L, const CanonicalForm& p)
 }
 
 
-static inline void
+static void
 chineseRemainder (const CFList & x1, const CanonicalForm & q1,
                   const CFList & x2, const CanonicalForm & q2,
                   CFList & xnew, CanonicalForm & qnew)
@@ -157,7 +157,7 @@ chineseRemainder (const CFList & x1, const CanonicalForm & q1,
   qnew= tmp2;
 }
 
-static inline
+static
 CFList Farey (const CFList& L, const CanonicalForm& q)
 {
   CFList result;
@@ -166,7 +166,7 @@ CFList Farey (const CFList& L, const CanonicalForm& q)
   return result;
 }
 
-static inline
+static
 CFList replacevar (const CFList& L, const Variable& a, const Variable& b)
 {
   CFList result;
@@ -354,7 +354,6 @@ void sortList (CFList& list, const Variable& x)
   }
 }
 
-static inline
 CFList diophantine (const CanonicalForm& F, const CFList& factors)
 {
   if (getCharacteristic() == 0)
@@ -664,9 +663,8 @@ henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
   return;
 }
 
-static inline
 CFList
-biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
+biDiophantine (const CanonicalForm& F, const CFList& factors, int d)
 {
   Variable y= F.mvar();
   CFList result;
@@ -765,10 +763,9 @@ biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
   }
 }
 
-static inline
 CFList
 multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
-                     const CFList& recResult, const CFList& M, const int d)
+                     const CFList& recResult, const CFList& M, int d)
 {
   Variable y= F.mvar();
   CFList result;
@@ -854,7 +851,7 @@ multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
   return result;
 }
 
-static inline
+static
 void
 henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
             const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
@@ -1096,7 +1093,7 @@ henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
 }
 
 CFList
-henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
+henselLift23 (const CFList& eval, const CFList& factors, int* l, CFList&
               diophant, CFArray& Pi, CFMatrix& M)
 {
   CFList buf= factors;
@@ -1163,8 +1160,7 @@ henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end,
 
 CFList
 henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
-            diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
-            lNew)
+            diophant, CFArray& Pi, CFMatrix& M, int lOld, int lNew)
 {
   diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
   int k= 0;
@@ -1204,8 +1200,8 @@ henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
 }
 
 CFList
-henselLift (const CFList& eval, const CFList& factors, const int* l, const int
-            lLength, bool sort)
+henselLift (const CFList& eval, const CFList& factors, int* l, int lLength,
+            bool sort)
 {
   CFList diophant;
   CFList buf= factors;
@@ -1238,10 +1234,12 @@ henselLift (const CFList& eval, const CFList& factors, const int* l, const int
   return result;
 }
 
+// nonmonic
+
 void
-henselStep122 (const CanonicalForm& F, const CFList& factors,
-              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
-              CFArray& Pi, int j, const CFArray& LCs)
+nonMonicHenselStep12 (const CanonicalForm& F, const CFList& factors,
+                      CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
+                      CFArray& Pi, int j, const CFArray& LCs)
 {
   Variable x= F.mvar();
   CanonicalForm xToJ= power (x, j);
@@ -1450,8 +1448,9 @@ henselStep122 (const CanonicalForm& F, const CFList& factors,
 }
 
 void
-henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
-              CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
+nonMonicHenselLift12 (const CanonicalForm& F, CFList& factors, int l,
+                      CFArray& Pi, CFList& diophant, CFMatrix& M,
+                      const CFArray& LCs, bool sort)
 {
   if (sort)
     sortList (factors, Variable (1));
@@ -1517,7 +1516,7 @@ henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
   }
 
   for (i= 1; i < l; i++)
-    henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
+    nonMonicHenselStep12 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
 
   factors= CFList();
   for (i= 0; i < bufFactors.size(); i++)
@@ -1528,7 +1527,6 @@ henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
 
 /// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
 /// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
-static inline
 CFList
 diophantine (const CFList& recResult, const CFList& factors,
              const CFList& products, const CFList& M, const CanonicalForm& E,
@@ -1607,11 +1605,11 @@ diophantine (const CFList& recResult, const CFList& factors,
   return result;
 }
 
-static inline
 void
-henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
-            const CFList& diophant, CFMatrix& M, CFArray& Pi,
-            const CFList& products, int j, const CFList& MOD, bool& noOneToOne)
+nonMonicHenselStep (const CanonicalForm& F, const CFList& factors,
+                    CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
+                    CFArray& Pi, const CFList& products, int j,
+                    const CFList& MOD, bool& noOneToOne)
 {
   CanonicalForm E;
   CanonicalForm xToJ= power (F.mvar(), j);
@@ -1847,9 +1845,9 @@ CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
 }
 
 CFList
-henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
-              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
-              const CFList& LCs2, bool& bad)
+nonMonicHenselLift232(const CFList& eval, const CFList& factors, int* l, CFList&
+                      diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
+                      const CFList& LCs2, bool& bad)
 {
   CFList buf= factors;
   int k= 0;
@@ -1897,7 +1895,8 @@ henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
 
   for (int d= 1; d < l[1]; d++)
   {
-    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
+    nonMonicHenselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, products,
+                        d, MOD, bad);
     if (bad)
       return CFList();
   }
@@ -1909,9 +1908,10 @@ henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
 
 
 CFList
-henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
-            diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
-            lNew, const CFList& LCs1, const CFList& LCs2, bool& bad)
+nonMonicHenselLift2 (const CFList& F, const CFList& factors, const CFList& MOD,
+                    CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
+                    int& lNew, const CFList& LCs1, const CFList& LCs2, bool& bad
+                   )
 {
   int k= 0;
   CFArray bufFactors= CFArray (factors.length());
@@ -1958,7 +1958,8 @@ henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
 
   for (int d= 1; d < lNew; d++)
   {
-    henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
+    nonMonicHenselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, products,
+                        d, MOD, bad);
     if (bad)
       return CFList();
   }
@@ -1970,9 +1971,9 @@ henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
 }
 
 CFList
-henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
-             lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
-             const CFArray& Pi, const CFList& diophant, bool& bad)
+nonMonicHenselLift2 (const CFList& eval, const CFList& factors, int* l, int
+                    lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
+                    const CFArray& Pi, const CFList& diophant, bool& bad)
 {
   CFList bufDiophant= diophant;
   CFList buf= factors;
@@ -1980,8 +1981,8 @@ henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
     sortList (buf, Variable (1));
   CFArray bufPi= Pi;
   CFMatrix M= CFMatrix (l[1], factors.length());
-  CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2,
-                               bad);
+  CFList result= 
+    nonMonicHenselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2, bad);
   if (bad)
     return CFList();
 
@@ -2009,8 +2010,8 @@ henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
     bufLCs1.append (jj.getItem());
     bufLCs2.append (jjj.getItem());
     M= CFMatrix (l[i], factors.length());
-    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
-                         l[i], bufLCs1, bufLCs2, bad);
+    result= nonMonicHenselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M,
+                                 l[i - 1], l[i], bufLCs1, bufLCs2, bad);
     if (bad)
       return CFList();
     MOD.append (power (Variable (i + 2), l[i]));
@@ -2107,7 +2108,8 @@ nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
   MOD.insert (yToL);
   for (int d= 1; d < liftBound; d++)
   {
-    henselStep2 (F, factors, bufFactors, diophant, M, Pi, products, d, MOD, bad);
+    nonMonicHenselStep (F, factors, bufFactors, diophant, M, Pi, products, d,
+                        MOD, bad);
     if (bad)
       return CFList();
   }
@@ -2120,7 +2122,7 @@ nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
 
 CFList
 nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
-                    CFList& diophant, CFArray& Pi, CFMatrix& M, const int lOld,
+                    CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
                     int& lNew, const CFList& MOD, bool& noOneToOne
                    )
 {
@@ -2187,8 +2189,8 @@ nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
 
   for (int d= 1; d < lNew; d++)
   {
-    henselStep2 (F.getLast(), factors, bufFactors, diophant, M, Pi, products, d,
-                 MOD, noOneToOne);
+    nonMonicHenselStep (F.getLast(), factors, bufFactors, diophant, M, Pi,
+                        products, d, MOD, noOneToOne);
     if (noOneToOne)
       return CFList();
   }
diff --git a/factory/facHensel.h b/factory/facHensel.h
index e3bf585..7cf857c 100644
--- a/factory/facHensel.h
+++ b/factory/facHensel.h
@@ -22,9 +22,6 @@
 #include "timing.h"
 
 #include "canonicalform.h"
-#include "cf_iter.h"
-#include "templates/ftmpl_functions.h"
-#include "algext.h"
 
 /// sort a list of polynomials by their degree in @a x.
 ///
@@ -83,7 +80,7 @@ henselLift23 (const CFList& eval,    ///< [in] contains compressed, bivariate
               const CFList& factors, ///< [in] monic bivariate factors,
                                      ///< including leading coefficient
                                      ///< as first element.
-              const int* l,          ///< [in] l[0]: precision of bivariate
+              int* l,          ///< [in] l[0]: precision of bivariate
                                      ///< lifting, l[1] as above
               CFList& diophant,      ///< [in,out] returns the result of
                                      ///< biDiophantine()
@@ -127,8 +124,8 @@ henselLift (const CFList& F,       ///< [in] two compressed, multivariate
             CFList& diophant,      ///< [in,out] result of multiRecDiophantine()
             CFArray& Pi,           ///< [in,out] stores intermediate results
             CFMatrix& M,           ///< [in,out] stores intermediate results
-            const int lOld,       ///< [in] lifting precision of F.mvar()
-            const int lNew        ///< [in] lifting precision of G.mvar()
+            int lOld,       ///< [in] lifting precision of F.mvar()
+            int lNew        ///< [in] lifting precision of G.mvar()
            );
 
 /// Hensel lifting from bivariate to multivariate
@@ -144,16 +141,15 @@ henselLift (const CFList& eval,    ///< [in] a list of polynomials the last
                                    ///< compressed bivariate poly.
             const CFList& factors, ///< [in] bivariate factors, including
                                    ///< leading coefficient
-            const int* l,          ///< [in] lifting bounds
-            const int lLength,     ///< [in] length of l
+            int* l,          ///< [in] lifting bounds
+            int lLength,     ///< [in] length of l
             bool sort= true        ///< [in] sort factors by degree in
                                    ///< Variable(1)
            );
 
-/// two factor Hensel lifting from univariate to bivariate, factors need not to
-/// be monic
+/// Hensel lifting from univariate to bivariate, factors need not to be monic
 void
-henselLift122 (const CanonicalForm& F,///< [in] a bivariate poly
+nonMonicHenselLift12 (const CanonicalForm& F,///< [in] a bivariate poly
                CFList& factors,       ///< [in, out] a list of univariate polys
                                       ///< lifted factors
                int l,                 ///< [in] lift bound
@@ -170,7 +166,7 @@ henselLift122 (const CanonicalForm& F,///< [in] a bivariate poly
 ///
 /// @return @a henselLift122 returns a list of lifted factors
 CFList
-henselLift2 (const CFList& eval,   ///< [in] a list of polynomials the last
+nonMonicHenselLift2 (const CFList& eval,///< [in] a list of polynomials the last
                                    ///< element is a compressed multivariate
                                    ///< poly, last but one element equals the
                                    ///< last elements modulo its main variable
@@ -178,7 +174,7 @@ henselLift2 (const CFList& eval,   ///< [in] a list of polynomials the last
                                    ///< compressed bivariate poly.
              const CFList& factors,///< [in] bivariate factors
              int* l,               ///< [in] lift bounds
-             const int lLength,    ///< [in] length of l
+             int lLength,          ///< [in] length of l
              bool sort,            ///< [in] if true factors are sorted by
                                    ///< their degree in Variable(1)
              const CFList& LCs1,   ///< [in] a list of evaluated LC of first
diff --git a/factory/facMul.cc b/factory/facMul.cc
index f022ac5..27e3b2f 100644
--- a/factory/facMul.cc
+++ b/factory/facMul.cc
@@ -16,6 +16,7 @@
 #include "canonicalform.h"
 #include "facMul.h"
 #include "algext.h"
+#include "cf_util.h"
 #include "templates/ftmpl_functions.h"
 
 #ifdef HAVE_NTL
@@ -26,6 +27,8 @@
 #include "FLINTconvert.h"
 #endif
 
+// univariate polys
+
 #ifdef HAVE_FLINT
 void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
 {
@@ -47,8 +50,8 @@ void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
 
 
 CanonicalForm
-reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha,
-              const CanonicalForm& den)
+reverseSubstQa (const fmpz_poly_t F, int d, const Variable& alpha,
+                const CanonicalForm& den)
 {
   Variable x= Variable (1);
 
@@ -110,7 +113,7 @@ mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G,
   fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
 
   denA *= denB;
-  A= reverseSubst (FLINTA, d, alpha, denA);
+  A= reverseSubstQa (FLINTA, d, alpha, denA);
 
   fmpz_poly_clear (FLINTA);
   fmpz_poly_clear (FLINTB);
@@ -221,7 +224,7 @@ mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G,
   fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
 
   denA *= denB;
-  A= reverseSubst (FLINTA, d, alpha, denA);
+  A= reverseSubstQa (FLINTA, d, alpha, denA);
   fmpz_poly_clear (FLINTA);
   fmpz_poly_clear (FLINTB);
   return A;
@@ -569,6 +572,10 @@ divNTL (const CanonicalForm& F, const CanonicalForm& G)
   return result;
 }
 
+// end univariate polys
+//*************************
+// bivariate polys
+
 #ifdef HAVE_FLINT
 void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
 {
@@ -591,34 +598,6 @@ void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
   _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);
@@ -663,174 +642,149 @@ void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
   fmpq_clear (coeff);
   _fmpq_poly_normalise (result);
 }
-#endif
 
-zz_pX kronSubFp (const CanonicalForm& A, int d)
+void
+kronSubReciproFp (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
+                  int d)
 {
   int degAy= degree (A);
-  zz_pX result;
-  result.rep.SetLength (d*(degAy + 1));
+  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));
 
-  zz_p *resultp;
-  resultp= result.rep.elts();
-  zz_pX buf;
-  zz_p *bufp;
-  int j, k, bufRepLength;
+  nmod_poly_t buf;
 
+  int k, kk, j, bufRepLength;
   for (CFIterator i= A; i.hasTerms(); i++)
   {
-    if (i.coeff().inCoeffDomain())
-      buf= convertFacCF2NTLzzpX (i.coeff());
-    else
-      buf= convertFacCF2NTLzzpX (i.coeff());
+    convertFacCF2nmod_poly_t (buf, i.coeff());
 
     k= i.exp()*d;
-    bufp= buf.rep.elts();
-    bufRepLength= (int) buf.rep.length();
+    kk= (degAy - i.exp())*d;
+    bufRepLength= (int) nmod_poly_length (buf);
     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));
+      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()
+                                       )
+                             );
     }
-    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];
+    nmod_poly_clear (buf);
   }
-  result.normalize();
-
-  return result;
+  _nmod_poly_normalise (subA1);
+  _nmod_poly_normalise (subA2);
 }
 
 void
-kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
-                const Variable& alpha)
+kronSubReciproQ (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
+                 int d)
 {
   int degAy= degree (A);
-  subA1.rep.SetLength ((long) d*(degAy + 2));
-  subA2.rep.SetLength ((long) d*(degAy + 2));
+  fmpz_poly_init2 (subA1, d*(degAy + 2));
+  fmpz_poly_init2 (subA2, 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;
+  fmpz_poly_t buf;
+  fmpz_t coeff1, coeff2;
 
+  int k, kk, j, 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);
+    convertFacCF2Fmpz_poly_t (buf, i.coeff());
 
     k= i.exp()*d;
     kk= (degAy - i.exp())*d;
-    bufp= buf.rep.elts();
-    bufRepLength= (int) buf.rep.length();
+    bufRepLength= (int) fmpz_poly_length (buf);
     for (j= 0; j < bufRepLength; j++)
     {
-      subA1p [j + k] += bufp [j];
-      subA2p [j + kk] += bufp [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);
   }
-  subA1.normalize();
-  subA2.normalize();
+  fmpz_clear (coeff1);
+  fmpz_clear (coeff2);
+  _fmpz_poly_normalise (subA1);
+  _fmpz_poly_normalise (subA2);
 }
 
-void
-kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
+CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
 {
-  int degAy= degree (A);
-  subA1.rep.SetLength ((long) d*(degAy + 2));
-  subA2.rep.SetLength ((long) d*(degAy + 2));
+  Variable y= Variable (2);
+  Variable x= Variable (1);
 
-  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;
+  fmpz_poly_t f;
+  fmpz_poly_init (f);
+  fmpz_poly_set (f, F);
 
-  for (CFIterator i= A; i.hasTerms(); i++)
+  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)
   {
-    buf= convertFacCF2NTLzzpX (i.coeff());
+    degfSubK= degf - k;
+    if (degfSubK >= d)
+      repLength= d;
+    else
+      repLength= degfSubK + 1;
 
-    k= i.exp()*d;
-    kk= (degAy - i.exp())*d;
-    bufp= buf.rep.elts();
-    bufRepLength= (int) buf.rep.length();
-    for (j= 0; j < bufRepLength; j++)
+    fmpz_poly_init2 (buf, repLength);
+    fmpz_init (coeff);
+    for (j= 0; j < repLength; j++)
     {
-      subA1p [j + k] += bufp [j];
-      subA2p [j + kk] += bufp [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);
   }
-  subA1.normalize();
-  subA2.normalize();
+  fmpz_poly_clear (f);
+
+  return result;
 }
 
 CanonicalForm
-reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
-              const Variable& alpha)
+reverseSubstReciproFp (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
 {
   Variable y= Variable (2);
   Variable x= Variable (1);
 
-  zz_pEX f= F;
-  zz_pEX g= G;
-  int degf= deg(f);
-  int degg= deg(g);
+  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);
 
-  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));
+  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));
 
-  zz_pE *gp= g.rep.elts();
-  zz_pE *fp= f.rep.elts();
   CanonicalForm result= 0;
   int i= 0;
   int lf= 0;
@@ -838,8 +792,6 @@ reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int 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)
@@ -848,14 +800,13 @@ reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
       repLengthBuf1= 0;
     else
       repLengthBuf1= degfSubLf + 1;
-    buf1.rep.SetLength((long) repLengthBuf1);
+    nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
 
-    buf1p= buf1.rep.elts();
     for (ind= 0; ind < repLengthBuf1; ind++)
-      buf1p [ind]= fp [ind + lf];
-    buf1.normalize();
+      nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
+    _nmod_poly_normalise (buf1);
 
-    repLengthBuf1= buf1.rep.length();
+    repLengthBuf1= nmod_poly_length (buf1);
 
     if (deggSubLg >= d - 1)
       repLengthBuf2= d - 1;
@@ -864,27 +815,23 @@ reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
     else
       repLengthBuf2= deggSubLg + 1;
 
-    buf2.rep.SetLength ((long) repLengthBuf2);
-    buf2p= buf2.rep.elts();
+    nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
     for (ind= 0; ind < repLengthBuf2; ind++)
-      buf2p [ind]= gp [ind + lg];
-    buf2.normalize();
+      nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
 
-    repLengthBuf2= buf2.rep.length();
+    _nmod_poly_normalise (buf2);
+    repLengthBuf2= nmod_poly_length (buf2);
 
-    buf3.rep.SetLength((long) repLengthBuf2 + d);
-    buf3p= buf3.rep.elts();
-    buf2p= buf2.rep.elts();
-    buf1p= buf1.rep.elts();
+    nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
     for (ind= 0; ind < repLengthBuf1; ind++)
-      buf3p [ind]= buf1p [ind];
+      nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
     for (ind= repLengthBuf1; ind < d; ind++)
-      buf3p [ind]= zzpEZero;
+      nmod_poly_set_coeff_ui (buf3, ind, 0);
     for (ind= 0; ind < repLengthBuf2; ind++)
-      buf3p [ind + d]= buf2p [ind];
-    buf3.normalize();
+      nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
+    _nmod_poly_normalise (buf3);
 
-    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
+    result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
     i++;
 
 
@@ -894,55 +841,67 @@ reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
     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];
+        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;
-
-    buf2p= buf2.rep.elts();
+    }
     if (degfSubLf >= 0)
     {
       for (ind= 0; ind < repLengthBuf2; ind++)
-        fp [ind + lf] -= buf2p [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 zz_pX& F, const zz_pX& G, int d, int k)
+reverseSubstReciproQ (const fmpz_poly_t F, const fmpz_poly_t 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);
+  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);
 
-  zz_pX buf1;
-  zz_pX buf2;
-  zz_pX buf3;
 
-  zz_p *buf1p;
-  zz_p *buf2p;
-  zz_p *buf3p;
+  fmpz_poly_t buf1,buf2, buf3;
 
-  if (f.rep.length() < (long) d*(k+1)) //zero padding
-    f.rep.SetLength ((long)d*(k+1));
+  if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
+    fmpz_poly_fit_length (f,(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;
@@ -950,7 +909,7 @@ reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
   int degfSubLf= degf;
   int deggSubLg= degg-lg;
   int repLengthBuf2, repLengthBuf1, ind, tmp;
-  zz_p zzpZero= zz_p();
+  fmpz_t tmp1, tmp2;
   while (degf >= lf || lg >= 0)
   {
     if (degfSubLf >= d)
@@ -959,14 +918,17 @@ reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
       repLengthBuf1= 0;
     else
       repLengthBuf1= degfSubLf + 1;
-    buf1.rep.SetLength((long) repLengthBuf1);
 
-    buf1p= buf1.rep.elts();
+    fmpz_poly_init2 (buf1, repLengthBuf1);
+
     for (ind= 0; ind < repLengthBuf1; ind++)
-      buf1p [ind]= fp [ind + lf];
-    buf1.normalize();
+    {
+      fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
+      fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
+    }
+    _fmpz_poly_normalise (buf1);
 
-    repLengthBuf1= buf1.rep.length();
+    repLengthBuf1= fmpz_poly_length (buf1);
 
     if (deggSubLg >= d - 1)
       repLengthBuf2= d - 1;
@@ -975,29 +937,33 @@ reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
     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();
+    fmpz_poly_init2 (buf2, repLengthBuf2);
 
-    repLengthBuf2= buf2.rep.length();
+    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);
 
-    buf3.rep.SetLength((long) repLengthBuf2 + d);
-    buf3p= buf3.rep.elts();
-    buf2p= buf2.rep.elts();
-    buf1p= buf1.rep.elts();
+    fmpz_poly_init2 (buf3, repLengthBuf2 + d);
     for (ind= 0; ind < repLengthBuf1; ind++)
-      buf3p [ind]= buf1p [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++)
-      buf3p [ind]= zzpZero;
+      fmpz_poly_set_coeff_ui (buf3, ind, 0);
     for (ind= 0; ind < repLengthBuf2; ind++)
-      buf3p [ind + d]= buf2p [ind];
-    buf3.normalize();
+    {
+      fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
+      fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
+    }
+    _fmpz_poly_normalise (buf3);
 
-    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
+    result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
     i++;
 
 
@@ -1007,80 +973,63 @@ reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
     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];
+      {
+        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;
-
-    buf2p= buf2.rep.elts();
+    }
     if (degfSubLf >= 0)
     {
       for (ind= 0; ind < repLengthBuf2; ind++)
-        fp [ind + lf] -= buf2p [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);
   }
 
-  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;
-  }
+  fmpz_poly_clear (f);
+  fmpz_poly_clear (g);
+  fmpz_clear (tmp1);
+  fmpz_clear (tmp2);
 
   return result;
 }
 
-CanonicalForm reverseSubstFp (const zz_pX& F, int d)
+CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
 {
   Variable y= Variable (2);
   Variable x= Variable (1);
 
-  zz_pX f= F;
-  zz_p *fp= f.rep.elts();
+  nmod_poly_t f;
+  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
+  nmod_poly_set (f, F);
 
-  zz_pX buf;
-  zz_p *bufp;
+  nmod_poly_t buf;
   CanonicalForm result= 0;
   int i= 0;
-  int degf= deg(f);
+  int degf= nmod_poly_degree(f);
   int k= 0;
   int degfSubK, repLength, j;
   while (degf >= k)
@@ -1091,56 +1040,67 @@ CanonicalForm reverseSubstFp (const zz_pX& F, int d)
     else
       repLength= degfSubK + 1;
 
-    buf.rep.SetLength ((long) repLength);
-    bufp= buf.rep.elts();
+    nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
     for (j= 0; j < repLength; j++)
-      bufp [j]= fp [j + k];
-    buf.normalize();
+      nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
+    _nmod_poly_normalise (buf);
 
-    result += convertNTLzzpX2CF (buf, x)*power (y, i);
+    result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
     i++;
     k= d*i;
+    nmod_poly_clear (buf);
   }
+  nmod_poly_clear (f);
 
   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)
+mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
+                    CanonicalForm& M)
 {
-  int d1= degree (F, 1) + degree (G, 1) + 1;
+  int d1= tmax (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);
+  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);
+  kronSubReciproFp (F1, F2, F, d1);
+
+  nmod_poly_t G1, G2;
+  nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
+  nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
+  kronSubReciproFp (G1, G2, G, d1);
 
   int k= d1*degree (M);
-  MulTrunc (F1, F1, G1, (long) k);
+  nmod_poly_mullow (F1, F1, G1, (long) k);
 
-  int degtailF= degree (tailcoeff (F), 1);
+  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 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);
 
-  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
-  return reverseSubst (F1, F2, d1, d2);
+
+  CanonicalForm result= reverseSubstReciproFp (F1, F2, d1, d2);
+
+  nmod_poly_clear (F1);
+  nmod_poly_clear (F2);
+  nmod_poly_clear (G1);
+  nmod_poly_clear (G2);
+  return result;
 }
 
-//Kronecker substitution
 CanonicalForm
-mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
-              CanonicalForm& M)
+mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
+                CanonicalForm& M)
 {
   CanonicalForm A= F;
   CanonicalForm B= G;
@@ -1153,244 +1113,264 @@ mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
   int d2= tmax (degAy, degBy);
 
   if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
-    return mulMod2NTLFpReci (A, B, M);
+    return mulMod2FLINTFpReci (A, B, M);
 
-  zz_pX NTLA= kronSubFp (A, d1);
-  zz_pX NTLB= kronSubFp (B, d1);
+  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);
-  MulTrunc (NTLA, NTLA, NTLB, (long) k);
+  nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
 
-  A= reverseSubstFp (NTLA, d1);
+  A= reverseSubstFp (FLINTA, d1);
 
+  nmod_poly_clear (FLINTA);
+  nmod_poly_clear (FLINTB);
   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;
+mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
+                    CanonicalForm& M)
+{
+  int d1= tmax (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);
+  fmpz_poly_t F1, F2;
+  fmpz_poly_init (F1);
+  fmpz_poly_init (F2);
+  kronSubReciproQ (F1, F2, F, d1);
+
+  fmpz_poly_t G1, G2;
+  fmpz_poly_init (G1);
+  fmpz_poly_init (G2);
+  kronSubReciproQ (G1, G2, G, d1);
 
   int k= d1*degree (M);
-  MulTrunc (F1, F1, G1, (long) k);
+  fmpz_poly_mullow (F1, F1, G1, (long) k);
 
-  int degtailF= degree (tailcoeff (F), 1);
+  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 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);
 
-  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
-  return reverseSubst (F1, F2, d1, d2, alpha);
-}
+  CanonicalForm result= reverseSubstReciproQ (F1, F2, d1, d2);
 
-#ifdef HAVE_FLINT
-CanonicalForm
-mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
-                CanonicalForm& M);
-#endif
+  fmpz_poly_clear (F1);
+  fmpz_poly_clear (F2);
+  fmpz_poly_clear (G1);
+  fmpz_poly_clear (G2);
+  return result;
+}
 
 CanonicalForm
-mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
-              CanonicalForm& M)
+mulMod2FLINTQ (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 degAx= degree (A, 1);
+  int degBx= degree (B, 1);
+  int d1= degAx + 1 + degBx;
 
-    int k= d1*degree (M);
+  CanonicalForm f= bCommonDen (F);
+  CanonicalForm g= bCommonDen (G);
+  A *= f;
+  B *= g;
 
-    MulTrunc (NTLA, NTLA, NTLB, (long) k);
+  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);
 
-    A= reverseSubst (NTLA, d1, alpha);
+  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
+  A= reverseSubstQ (FLINTA, d1);
+  fmpz_poly_clear (FLINTA);
+  fmpz_poly_clear (FLINTB);
+  return A/(f*g);
+}
 
-    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)
+zz_pX kronSubFp (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));
+  zz_pX result;
+  result.rep.SetLength (d*(degAy + 1));
 
-  nmod_poly_t buf;
+  zz_p *resultp;
+  resultp= result.rep.elts();
+  zz_pX buf;
+  zz_p *bufp;
+  int j, k, bufRepLength;
 
-  int k, kk, j, bufRepLength;
   for (CFIterator i= A; i.hasTerms(); i++)
   {
-    convertFacCF2nmod_poly_t (buf, i.coeff());
+    if (i.coeff().inCoeffDomain())
+      buf= convertFacCF2NTLzzpX (i.coeff());
+    else
+      buf= convertFacCF2NTLzzpX (i.coeff());
 
     k= i.exp()*d;
-    kk= (degAy - i.exp())*d;
-    bufRepLength= (int) nmod_poly_length (buf);
+    bufp= buf.rep.elts();
+    bufRepLength= (int) buf.rep.length();
     for (j= 0; j < bufRepLength; j++)
+      resultp [j + k]= bufp [j];
+  }
+  result.normalize();
+
+  return result;
+}
+
+zz_pEX kronSubFq (const CanonicalForm& A, int d, const Variable& alpha)
+{
+  int degAy= degree (A);
+  zz_pEX result;
+  result.rep.SetLength (d*(degAy + 1));
+
+  Variable v;
+  zz_pE *resultp;
+  resultp= result.rep.elts();
+  zz_pEX buf1;
+  zz_pE *buf1p;
+  zz_pX buf2;
+  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+  int j, k, buf1RepLength;
+
+  for (CFIterator i= A; i.hasTerms(); i++)
+  {
+    if (i.coeff().inCoeffDomain())
     {
-      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()
-                                       )
-                             );
+      buf2= convertFacCF2NTLzzpX (i.coeff());
+      buf1= to_zz_pEX (to_zz_pE (buf2));
     }
-    nmod_poly_clear (buf);
+    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];
   }
-  _nmod_poly_normalise (subA1);
-  _nmod_poly_normalise (subA2);
+  result.normalize();
+
+  return result;
 }
 
 void
-kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
-                int d)
+kronSubReciproFq (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
+                  const Variable& alpha)
 {
   int degAy= degree (A);
-  fmpz_poly_init2 (subA1, d*(degAy + 2));
-  fmpz_poly_init2 (subA2, d*(degAy + 2));
+  subA1.rep.SetLength ((long) d*(degAy + 2));
+  subA2.rep.SetLength ((long) d*(degAy + 2));
 
-  fmpz_poly_t buf;
-  fmpz_t coeff1, coeff2;
+  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;
 
-  int k, kk, j, bufRepLength;
   for (CFIterator i= A; i.hasTerms(); i++)
   {
-    convertFacCF2Fmpz_poly_t (buf, i.coeff());
+    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;
-    bufRepLength= (int) fmpz_poly_length (buf);
+    bufp= buf.rep.elts();
+    bufRepLength= (int) buf.rep.length();
     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);
+      subA1p [j + k] += bufp [j];
+      subA2p [j + kk] += bufp [j];
     }
-    fmpz_poly_clear (buf);
   }
-  fmpz_clear (coeff1);
-  fmpz_clear (coeff2);
-  _fmpz_poly_normalise (subA1);
-  _fmpz_poly_normalise (subA2);
+  subA1.normalize();
+  subA2.normalize();
 }
 
-CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
+void
+kronSubReciproFp (zz_pX& subA1, zz_pX& subA2, const CanonicalForm& A, int d)
 {
-  Variable y= Variable (2);
-  Variable x= Variable (1);
+  int degAy= degree (A);
+  subA1.rep.SetLength ((long) d*(degAy + 2));
+  subA2.rep.SetLength ((long) d*(degAy + 2));
 
-  fmpz_poly_t f;
-  fmpz_poly_init (f);
-  fmpz_poly_set (f, F);
+  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;
 
-  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)
+  for (CFIterator i= A; i.hasTerms(); i++)
   {
-    degfSubK= degf - k;
-    if (degfSubK >= d)
-      repLength= d;
-    else
-      repLength= degfSubK + 1;
+    buf= convertFacCF2NTLzzpX (i.coeff());
 
-    fmpz_poly_init2 (buf, repLength);
-    fmpz_init (coeff);
-    for (j= 0; j < repLength; j++)
+    k= i.exp()*d;
+    kk= (degAy - i.exp())*d;
+    bufp= buf.rep.elts();
+    bufRepLength= (int) buf.rep.length();
+    for (j= 0; j < bufRepLength; j++)
     {
-      fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
-      fmpz_poly_set_coeff_fmpz (buf, j, coeff);
+      subA1p [j + k] += bufp [j];
+      subA2p [j + kk] += bufp [j];
     }
-    _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;
+  subA1.normalize();
+  subA2.normalize();
 }
 
 CanonicalForm
-reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
+reverseSubstReciproFq (const zz_pEX& F, const zz_pEX& G, int d, int k,
+                       const Variable& alpha)
 {
   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);
-
+  zz_pEX f= F;
+  zz_pEX g= G;
+  int degf= deg(f);
+  int degg= deg(g);
 
-  nmod_poly_t buf1,buf2, buf3;
+  zz_pEX buf1;
+  zz_pEX buf2;
+  zz_pEX buf3;
 
-  if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
-    nmod_poly_fit_length (f,(long)d*(k+1));
+  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;
@@ -1398,6 +1378,8 @@ reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int 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)
@@ -1406,13 +1388,14 @@ reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
       repLengthBuf1= 0;
     else
       repLengthBuf1= degfSubLf + 1;
-    nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
+    buf1.rep.SetLength((long) repLengthBuf1);
 
+    buf1p= buf1.rep.elts();
     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);
+      buf1p [ind]= fp [ind + lf];
+    buf1.normalize();
 
-    repLengthBuf1= nmod_poly_length (buf1);
+    repLengthBuf1= buf1.rep.length();
 
     if (deggSubLg >= d - 1)
       repLengthBuf2= d - 1;
@@ -1421,23 +1404,27 @@ reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
     else
       repLengthBuf2= deggSubLg + 1;
 
-    nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
+    buf2.rep.SetLength ((long) repLengthBuf2);
+    buf2p= buf2.rep.elts();
     for (ind= 0; ind < repLengthBuf2; ind++)
-      nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
+      buf2p [ind]= gp [ind + lg];
+    buf2.normalize();
 
-    _nmod_poly_normalise (buf2);
-    repLengthBuf2= nmod_poly_length (buf2);
+    repLengthBuf2= buf2.rep.length();
 
-    nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
+    buf3.rep.SetLength((long) repLengthBuf2 + d);
+    buf3p= buf3.rep.elts();
+    buf2p= buf2.rep.elts();
+    buf1p= buf1.rep.elts();
     for (ind= 0; ind < repLengthBuf1; ind++)
-      nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
+      buf3p [ind]= buf1p [ind];
     for (ind= repLengthBuf1; ind < d; ind++)
-      nmod_poly_set_coeff_ui (buf3, ind, 0);
+      buf3p [ind]= zzpEZero;
     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);
+      buf3p [ind + d]= buf2p [ind];
+    buf3.normalize();
 
-    result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
+    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
     i++;
 
 
@@ -1447,67 +1434,55 @@ reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
     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++)
-        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()
-                                         )
-                               );
+        gp [ind + lg] -= buf1p [ind];
     }
+
     if (lg < 0)
-    {
-      nmod_poly_clear (buf1);
-      nmod_poly_clear (buf2);
-      nmod_poly_clear (buf3);
       break;
-    }
+
+    buf2p= buf2.rep.elts();
     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()
-                                         )
-                               );
+        fp [ind + lf] -= buf2p [ind];
     }
-    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)
+reverseSubstReciproFp (const zz_pX& F, const zz_pX& 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);
+  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;
 
-  fmpz_poly_t buf1,buf2, buf3;
+  zz_p *buf1p;
+  zz_p *buf2p;
+  zz_p *buf3p;
 
-  if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
-    fmpz_poly_fit_length (f,(long)d*(k+1));
+  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;
@@ -1515,7 +1490,7 @@ reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
   int degfSubLf= degf;
   int deggSubLg= degg-lg;
   int repLengthBuf2, repLengthBuf1, ind, tmp;
-  fmpz_t tmp1, tmp2;
+  zz_p zzpZero= zz_p();
   while (degf >= lf || lg >= 0)
   {
     if (degfSubLf >= d)
@@ -1524,17 +1499,14 @@ reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
       repLengthBuf1= 0;
     else
       repLengthBuf1= degfSubLf + 1;
+    buf1.rep.SetLength((long) repLengthBuf1);
 
-    fmpz_poly_init2 (buf1, repLengthBuf1);
-
+    buf1p= buf1.rep.elts();
     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);
+      buf1p [ind]= fp [ind + lf];
+    buf1.normalize();
 
-    repLengthBuf1= fmpz_poly_length (buf1);
+    repLengthBuf1= buf1.rep.length();
 
     if (deggSubLg >= d - 1)
       repLengthBuf2= d - 1;
@@ -1543,33 +1515,29 @@ reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
     else
       repLengthBuf2= deggSubLg + 1;
 
-    fmpz_poly_init2 (buf2, repLengthBuf2);
-
+    buf2.rep.SetLength ((long) repLengthBuf2);
+    buf2p= buf2.rep.elts();
     for (ind= 0; ind < repLengthBuf2; ind++)
-    {
-      fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
-      fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
-    }
+      buf2p [ind]= gp [ind + lg];
 
-    _fmpz_poly_normalise (buf2);
-    repLengthBuf2= fmpz_poly_length (buf2);
+    buf2.normalize();
 
-    fmpz_poly_init2 (buf3, repLengthBuf2 + d);
+    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++)
-    {
-      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);
-    }
+      buf3p [ind]= buf1p [ind];
     for (ind= repLengthBuf1; ind < d; ind++)
-      fmpz_poly_set_coeff_ui (buf3, ind, 0);
+      buf3p [ind]= zzpZero;
     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);
+      buf3p [ind + d]= buf2p [ind];
+    buf3.normalize();
 
-    result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
+    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
     i++;
 
 
@@ -1579,63 +1547,80 @@ reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
     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++)
-      {
-        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);
-      }
+        gp [ind + lg] -= buf1p [ind];
     }
     if (lg < 0)
-    {
-      fmpz_poly_clear (buf1);
-      fmpz_poly_clear (buf2);
-      fmpz_poly_clear (buf3);
       break;
-    }
+
+    buf2p= buf2.rep.elts();
     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);
-      }
+        fp [ind + lf] -= buf2p [ind];
     }
-    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 reverseSubstFq (const zz_pEX& F, int d, const Variable& alpha)
+{
+  Variable y= Variable (2);
+  Variable x= Variable (1);
+
+  zz_pEX f= F;
+  zz_pE *fp= f.rep.elts();
+
+  zz_pEX buf;
+  zz_pE *bufp;
+  CanonicalForm result= 0;
+  int i= 0;
+  int degf= deg(f);
+  int k= 0;
+  int degfSubK, repLength, j;
+  while (degf >= k)
+  {
+    degfSubK= degf - k;
+    if (degfSubK >= d)
+      repLength= d;
+    else
+      repLength= degfSubK + 1;
+
+    buf.rep.SetLength ((long) repLength);
+    bufp= buf.rep.elts();
+    for (j= 0; j < repLength; j++)
+      bufp [j]= fp [j + k];
+    buf.normalize();
+
+    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
+    i++;
+    k= d*i;
+  }
 
   return result;
 }
 
-CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
+CanonicalForm reverseSubstFp (const zz_pX& 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);
+  zz_pX f= F;
+  zz_p *fp= f.rep.elts();
 
-  nmod_poly_t buf;
+  zz_pX buf;
+  zz_p *bufp;
   CanonicalForm result= 0;
   int i= 0;
-  int degf= nmod_poly_degree(f);
+  int degf= deg(f);
   int k= 0;
   int degfSubK, repLength, j;
   while (degf >= k)
@@ -1646,67 +1631,56 @@ CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
     else
       repLength= degfSubK + 1;
 
-    nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
+    buf.rep.SetLength ((long) repLength);
+    bufp= buf.rep.elts();
     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);
+      bufp [j]= fp [j + k];
+    buf.normalize();
 
-    result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
+    result += convertNTLzzpX2CF (buf, x)*power (y, i);
     i++;
     k= d*i;
-    nmod_poly_clear (buf);
   }
-  nmod_poly_clear (f);
 
   return result;
 }
 
+// assumes input to be reduced mod M and to be an element of Fq not Fp
 CanonicalForm
-mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
-                    CanonicalForm& M)
+mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
+                  CanonicalForm& M)
 {
-  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
+  int d1= 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);
+  zz_pX F1, F2;
+  kronSubReciproFp (F1, F2, F, d1);
+  zz_pX G1, G2;
+  kronSubReciproFp (G1, G2, G, d1);
 
   int k= d1*degree (M);
-  nmod_poly_mullow (F1, F1, G1, (long) k);
+  MulTrunc (F1, F1, G1, (long) k);
 
-  int degtailF= degree (tailcoeff (F), 1);;
+  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);
 
-  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);
+  reverse (F2, F2);
+  reverse (G2, G2);
+  MulTrunc (F2, F2, G2, b + 1);
+  reverse (F2, F2, b);
 
-  nmod_poly_clear (F1);
-  nmod_poly_clear (F2);
-  nmod_poly_clear (G1);
-  nmod_poly_clear (G2);
-  return result;
+  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
+  return reverseSubstReciproFp (F1, F2, d1, d2);
 }
 
+//Kronecker substitution
 CanonicalForm
-mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
-                CanonicalForm& M)
+mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
+              CanonicalForm& M)
 {
   CanonicalForm A= F;
   CanonicalForm B= G;
@@ -1719,97 +1693,100 @@ mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
   int d2= tmax (degAy, degBy);
 
   if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
-    return mulMod2FLINTFpReci (A, B, M);
+    return mulMod2NTLFpReci (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);
+  zz_pX NTLA= kronSubFp (A, d1);
+  zz_pX NTLB= kronSubFp (B, d1);
 
   int k= d1*degree (M);
-  nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
+  MulTrunc (NTLA, NTLA, NTLB, (long) k);
 
-  A= reverseSubstFp (FLINTA, d1);
+  A= reverseSubstFp (NTLA, d1);
 
-  nmod_poly_clear (FLINTA);
-  nmod_poly_clear (FLINTB);
   return A;
 }
 
+// assumes input to be reduced mod M and to be an element of Fq not Fp
 CanonicalForm
-mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
-                    CanonicalForm& M)
+mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
+                  CanonicalForm& M, const Variable& alpha)
 {
-  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
+  int d1= 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);
+  zz_pEX F1, F2;
+  kronSubReciproFq (F1, F2, F, d1, alpha);
+  zz_pEX G1, G2;
+  kronSubReciproFq (G1, G2, G, d1, alpha);
 
   int k= d1*degree (M);
-  fmpz_poly_mullow (F1, F1, G1, (long) k);
+  MulTrunc (F1, F1, G1, (long) k);
 
-  int degtailF= degree (tailcoeff (F), 1);;
+  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);
 
-  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);
+  reverse (F2, F2);
+  reverse (G2, G2);
+  MulTrunc (F2, F2, G2, b + 1);
+  reverse (F2, F2, b);
 
-  fmpz_poly_clear (F1);
-  fmpz_poly_clear (F2);
-  fmpz_poly_clear (G1);
-  fmpz_poly_clear (G2);
-  return result;
+  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
+  return reverseSubstReciproFq (F1, F2, d1, d2, alpha);
 }
 
+#ifdef HAVE_FLINT
 CanonicalForm
-mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
-                CanonicalForm& M)
+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;
 
-  int degAx= degree (A, 1);
-  int degBx= degree (B, 1);
-  int d1= degAx + 1 + degBx;
+  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);
 
-  CanonicalForm f= bCommonDen (F);
-  CanonicalForm g= bCommonDen (G);
-  A *= f;
-  B *= g;
+    int degMipo= degree (getMipo (alpha));
+    if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
+        (2*degAy > degree (M)))
+      return mulMod2NTLFqReci (A, B, M, alpha);
 
-  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);
+    zz_pEX NTLA= kronSubFq (A, d1, alpha);
+    zz_pEX NTLB= kronSubFq (B, d1, alpha);
 
-  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
-  A= reverseSubstQ (FLINTA, d1);
-  fmpz_poly_clear (FLINTA);
-  fmpz_poly_clear (FLINTB);
-  return A/(f*g);
-}
+    int k= d1*degree (M);
+
+    MulTrunc (NTLA, NTLA, NTLB, (long) k);
 
+    A= reverseSubstFq (NTLA, d1, alpha);
+
+    return A;
+  }
+  else
+#ifdef HAVE_FLINT
+    return mulMod2FLINTFp (A, B, M);
+#else
+    return mulMod2NTLFp (A, B, M);
 #endif
+}
 
 CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
                        const CanonicalForm& M)
@@ -1886,6 +1863,10 @@ CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
   DEBOUTLN (cerr, "fatal end in mulMod2");
 }
 
+// end bivariate polys
+//**********************
+// multivariate polys
+
 CanonicalForm mod (const CanonicalForm& F, const CFList& M)
 {
   CanonicalForm A= F;
@@ -2030,6 +2011,10 @@ CanonicalForm prodMod (const CFList& L, const CFList& M)
   }
 }
 
+// end multivariate polys
+//***************************
+// division
+
 CanonicalForm reverse (const CanonicalForm& F, int d)
 {
   if (d == 0)
@@ -2509,4 +2494,6 @@ void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
   return;
 }
 
+// end division
+
 #endif

-- 
an open source computer algebra system



More information about the debian-science-commits mailing list