Changeset 08daea in git

Timestamp:

Nov 22, 2010, 11:12:46 AM (13 years ago)

Author:

Martin Lee <martinlee84@…>

Branches:

(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')

Children:

fb8ba2756f35787214f1a0b5702e31e4b853de97

Parents:

cde7084dc141670809afd29cafbd9c571bd9fc3a

Message:

new ezgcd, heuristic which gcd to choose and sparse modular gcd all for finite fields


git-svn-id: file:///usr/local/Singular/svn/trunk@13661 2c84dea3-7e68-4137-9b89-c4e89433aadc

Files:

• Singular/tesths.cc (modified) (1 diff)
• factory/cf_gcd.cc (modified) (2 diffs)
• factory/cf_gcd_smallp.cc (modified) (4 diffs)
• factory/cf_gcd_smallp.h (modified) (1 diff)
• factory/cf_map_ext.cc (modified) (1 diff)
• factory/cf_random.cc (modified) (1 diff)
• factory/facHensel.cc (modified) (2 diffs)
• factory/facHensel.h (modified) (1 diff)
• factory/ffreval.cc (modified) (1 diff)
• factory/ffreval.h (modified) (2 diffs)

Singular/tesths.cc

@@ -55 +55 @@
   //On(SW_USE_FF_MOD_GCD);
   //On(SW_USE_fieldGCD);
-  Off(SW_USE_EZGCD_P);
+  On(SW_USE_EZGCD_P);
   On(SW_USE_QGCD);
   Off(SW_USE_NTL_SORT); // may be changed by an command line option

factory/cf_gcd.cc

@@ -566 +566 @@
     else if (isOn( SW_USE_EZGCD_P ) && (!fc_and_gc_Univariate))
     {
-      if ( pe == 1 )
+      /*if ( pe == 1 )
         fc = fin_ezgcd( fc, gc );
       else if ( pe > 0 )// no variable at position 1
@@ -580 +580 @@
         gc = swapvar( gc, Variable(pe), Variable(1) );
         fc = swapvar( fin_ezgcd( fc, gc ), Variable(1), Variable(pe) );
-      }
+      }*/
+      fc= EZGCD_P (fc, gc);
     }
     else if (isOn(SW_USE_GCD_P))

factory/cf_gcd_smallp.cc

@@ -30 +30 @@
 #include "ftmpl_functions.h"
 #include "cf_random.h"
+#include "ffreval.h"
+#include "facHensel.h"
 
 // iinline helper functions:
@@ -46 +48 @@
 TIMING_DEFINE_PRINT(newton_interpolation);
 
+static const double log2exp= 1.442695041;
+
 /// compressing two polynomials F and G, M is used for compressing,
 /// N to reverse the compression
 static inline
 int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
-                CFMap & N, bool& topLevel)
+                CFMap & N, bool topLevel)
 {
   int n= tmax (F.level(), G.level());
@@ -205 +209 @@
   return 1;
 }
+
 
 int
@@ -1359 +1364 @@
 }
 
+CFArray
+solveVandermonde (const CFArray& M, const CFArray& A)
+{
+  int r= M.size();
+  ASSERT (A.size() == r, "vector does not have right size");
+
+  if (r == 1)
+  {
+    CFArray result= CFArray (1);
+    result [0]= A [0] / M [0];
+    return result;
+  }
+  // check solvability
+  bool notDistinct= false;
+  for (int i= 0; i < r - 1; i++)
+  {
+    for (int j= i + 1; j < r; j++)
+    {
+      if (M [i] == M [j])
+      {
+        notDistinct= true;
+        break;
+      }
+    }
+  }
+  if (notDistinct)
+    return CFArray();
+
+  CanonicalForm master= 1;
+  Variable x= Variable (1);
+  for (int i= 0; i < r; i++)
+    master *= x - M [i];
+  CFList Pj;
+  CanonicalForm tmp;
+  for (int i= 0; i < r; i++)
+  {
+    tmp= master/(x - M [i]);
+    tmp /= tmp (M [i], 1);
+    Pj.append (tmp);
+  }
+  CFArray result= CFArray (r);
+
+  CFListIterator j= Pj;
+  for (int i= 1; i <= r; i++, j++)
+  {
+    tmp= 0;
+    for (int l= 0; l < A.size(); l++)
+      tmp += A[l]*j.getItem()[l];
+    result[i - 1]= tmp;
+  }
+  return result;
+}
+
+CFArray
+solveGeneralVandermonde (const CFArray& M, const CFArray& A)
+{
+  int r= M.size();
+  ASSERT (c == r, "number of columns and rows not equal");
+  ASSERT (A.size() == r, "vector does not have right size");
+  if (r == 1)
+  {
+    CFArray result= CFArray (1);
+    result [0]= A[0] / M [0];
+    return result;
+  }
+  // check solvability
+  bool notDistinct= false;
+  for (int i= 0; i < r - 1; i++)
+  {
+    for (int j= i + 1; j < r; j++)
+    {
+      if (M [i] == M [j])
+      {
+        notDistinct= true;
+        break;
+      }
+    }
+  }
+  if (notDistinct)
+    return CFArray();
+
+  CanonicalForm master= 1;
+  Variable x= Variable (1);
+  for (int i= 0; i < r; i++)
+    master *= x - M [i];
+  master *= x;
+  CFList Pj;
+  CanonicalForm tmp;
+  for (int i= 0; i < r; i++)
+  {
+    tmp= master/(x - M [i]);
+    tmp /= tmp (M [i], 1);
+    Pj.append (tmp);
+  }
+
+  CFArray result= CFArray (r);
+
+  CFListIterator j= Pj;
+  for (int i= 1; i <= r; i++, j++)
+  {
+    tmp= 0;
+
+    for (int l= 1; l <= A.size(); l++)
+      tmp += A[l - 1]*j.getItem()[l];
+    result[i - 1]= tmp;
+  }
+  return result;
+}
+
+/// M in row echolon form, rk rank of M
+CFArray
+readOffSolution (const CFMatrix& M, const long rk)
+{
+  CFArray result= CFArray (rk);
+  CanonicalForm tmp1, tmp2, tmp3;
+  for (int i= rk; i >= 1; i--)
+  {
+    tmp3= 0;
+    tmp1= M (i, M.columns());
+    for (int j= M.columns() - 1; j >= 1; j--)
+    {
+      tmp2= M (i, j);
+      if (j == i)
+        break;
+      else
+        tmp3 += tmp2*result[j - 1];
+    }
+    result[i - 1]= (tmp1 - tmp3)/tmp2;
+  }
+  return result;
+}
+
+CFArray
+readOffSolution (const CFMatrix& M, const CFArray& L, const CFArray& partialSol)
+{
+  CFArray result= CFArray (M.rows());
+  CanonicalForm tmp1, tmp2, tmp3;
+  int k;
+  for (int i= M.rows(); i >= 1; i--)
+  {
+    tmp3= 0;
+    tmp1= L[i - 1];
+    k= 0;
+    for (int j= M.columns(); j >= 1; j--, k++)
+    {
+      tmp2= M (i, j);
+      if (j == i)
+        break;
+      else
+      {
+        if (k > partialSol.size() - 1)
+          tmp3 += tmp2*result[j - 1];
+        else
+          tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
+      }
+    }
+    result [i - 1]= (tmp1 - tmp3)/tmp2;
+  }
+  return result;
+}
+
+long
+gaussianElimFp (CFMatrix& M, CFArray& L)
+{
+  ASSERT (L.size() <= M.rows(), "dimension exceeded");
+  CFMatrix *N;
+  N= new CFMatrix (M.rows(), M.columns() + 1);
+
+  for (int i= 1; i <= M.rows(); i++)
+    for (int j= 1; j <= M.columns(); j++)
+      (*N) (i, j)= M (i, j);
+
+  int j= 1;
+  for (int i= 0; i < L.size(); i++, j++)
+    (*N) (j, M.columns() + 1)= L[i];
+  int p= getCharacteristic ();
+  zz_p::init (p);
+  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
+  long rk= gauss (*NTLN);
+
+  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
+
+  L= CFArray (M.rows());
+  for (int i= 0; i < M.rows(); i++)
+    L[i]= (*N) (i + 1, M.columns() + 1);
+  M= (*N) (1, M.rows(), 1, M.columns());
+  return rk;
+}
+
+long
+gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)
+{
+  ASSERT (L.size() <= M.rows(), "dimension exceeded");
+  CFMatrix *N;
+  N= new CFMatrix (M.rows(), M.columns() + 1);
+
+  for (int i= 1; i <= M.rows(); i++)
+    for (int j= 1; j <= M.columns(); j++)
+      (*N) (i, j)= M (i, j);
+
+  int j= 1;
+  for (int i= 0; i < L.size(); i++, j++)
+    (*N) (j, M.columns() + 1)= L[i];
+  int p= getCharacteristic ();
+  zz_p::init (p);
+  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+  zz_pE::init (NTLMipo);
+  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
+  long rk= gauss (*NTLN);
+
+  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
+
+  M= (*N) (1, M.rows(), 1, M.columns());
+  L= CFArray (M.rows());
+  for (int i= 0; i < M.rows(); i++)
+    L[i]= (*N) (i + 1, M.columns() + 1);
+  return rk;
+}
+
+CFArray
+solveSystemFp (const CFMatrix& M, const CFArray& L)
+{
+  ASSERT (L.size() <= M.rows(), "dimension exceeded");
+  CFMatrix *N;
+  N= new CFMatrix (M.rows(), M.columns() + 1);
+
+  for (int i= 1; i <= M.rows(); i++)
+    for (int j= 1; j <= M.columns(); j++)
+      (*N) (i, j)= M (i, j);
+
+  int j= 1;
+  for (int i= 0; i < L.size(); i++, j++)
+    (*N) (j, M.columns() + 1)= L[i];
+  int p= getCharacteristic ();
+  zz_p::init (p);
+  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
+  long rk= gauss (*NTLN);
+  if (rk != M.columns())
+    return CFArray();
+
+  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
+
+  CFArray A= readOffSolution (*N, rk);
+
+  return A;
+}
+
+CFArray
+solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
+{
+  ASSERT (L.size() <= M.rows(), "dimension exceeded");
+  CFMatrix *N;
+  N= new CFMatrix (M.rows(), M.columns() + 1);
+
+  for (int i= 1; i <= M.rows(); i++)
+    for (int j= 1; j <= M.columns(); j++)
+      (*N) (i, j)= M (i, j);
+  int j= 1;
+  for (int i= 0; i < L.size(); i++, j++)
+    (*N) (j, M.columns() + 1)= L[i];
+  int p= getCharacteristic ();
+  zz_p::init (p);
+  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+  zz_pE::init (NTLMipo);
+  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
+  long rk= gauss (*NTLN);
+  if (rk != M.columns())
+    return CFArray();
+
+  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
+
+  CFArray A= readOffSolution (*N, rk);
+
+  return A;
+}
+
+CFArray
+getMonoms (const CanonicalForm& F)
+{
+  if (F.inCoeffDomain())
+  {
+    CFArray result= CFArray (1);
+    result [0]= 1;
+    return result;
+  }
+  if (F.isUnivariate())
+  {
+    CFArray result= CFArray (size(F));
+    int j= 0;
+    for (CFIterator i= F; i.hasTerms(); i++, j++)
+      result[j]= power (F.mvar(), i.exp());
+    return result;
+  }
+  int numMon= size (F);
+  CFArray result= CFArray (numMon);
+  int j= 0;
+  CFArray recResult;
+  Variable x= F.mvar();
+  CanonicalForm powX;
+  for (CFIterator i= F; i.hasTerms(); i++)
+  {
+    powX= power (x, i.exp());
+    recResult= getMonoms (i.coeff());
+    for (int k= 0; k < recResult.size(); k++)
+      result[j+k]= powX*recResult[k];
+    j += recResult.size();
+  }
+  return result;
+}
+
+CFArray
+evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)
+{
+  if (F.inCoeffDomain())
+  {
+    CFArray result= CFArray (1);
+    result [0]= F;
+    return result;
+  }
+  if (F.isUnivariate())
+  {
+    ASSERT (evalPoints.length() == 1,
+            "expected an eval point with only one component");
+    CFArray result= CFArray (size(F));
+    int j= 0;
+    CanonicalForm evalPoint= evalPoints.getLast();
+    for (CFIterator i= F; i.hasTerms(); i++, j++)
+      result[j]= power (evalPoint, i.exp());
+    return result;
+  }
+  int numMon= size (F);
+  CFArray result= CFArray (numMon);
+  int j= 0;
+  CanonicalForm evalPoint= evalPoints.getLast();
+  CFList buf= evalPoints;
+  buf.removeLast();
+  CFArray recResult;
+  CanonicalForm powEvalPoint;
+  for (CFIterator i= F; i.hasTerms(); i++)
+  {
+    powEvalPoint= power (evalPoint, i.exp());
+    recResult= evaluateMonom (i.coeff(), buf);
+    for (int k= 0; k < recResult.size(); k++)
+      result[j+k]= powEvalPoint*recResult[k];
+    j += recResult.size();
+  }
+  return result;
+}
+
+CFArray
+evaluate (const CFArray& A, const CFList& evalPoints)
+{
+  CFArray result= A.size();
+  CanonicalForm tmp;
+  int k;
+  for (int i= 0; i < A.size(); i++)
+  {
+    tmp= A[i];
+    k= 1;
+    for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
+      tmp= tmp (j.getItem(), k);
+    result[i]= tmp;
+  }
+  return result;
+}
+
+CFList
+evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,
+                  CanonicalForm& Feval, CanonicalForm& Geval,
+                  const CanonicalForm& LCF, const bool& GF,
+                  const Variable& alpha, bool& fail, CFList& list
+                 )
+{
+  int k= tmax (F.level(), G.level()) - 1;
+  Variable x= Variable (1);
+  CFList result;
+  FFRandom genFF;
+  GFRandom genGF;
+  int p= getCharacteristic ();
+  int bound;
+  if (alpha != Variable (1))
+  {
+    bound= ipower (p, degree (getMipo(alpha)));
+    bound= ipower (bound, k);
+  }
+  else if (GF)
+  {
+    bound= ipower (p, getGFDegree());
+    bound= ipower (bound, k);
+  }
+  else
+    bound= ipower (p, k);
+
+  CanonicalForm random;
+  int j;
+  bool zeroOneOccured= false;
+  bool allEqual= false;
+  CanonicalForm buf;
+  do
+  {
+    random= 0;
+    // possible overflow if list.length() does not fit into a int
+    if (list.length() >= bound)
+    {
+      fail= true;
+      break;
+    }
+    for (int i= 0; i < k; i++)
+    {
+      if (GF)
+      {
+        result.append (genGF.generate());
+        random += result.getLast()*power (x, i);
+      }
+      else if (alpha != Variable(1))
+      {
+        AlgExtRandomF genAlgExt (alpha);
+        result.append (genAlgExt.generate());
+        random += result.getLast()*power (x, i);
+      }
+      else
+      {
+        result.append (genFF.generate());
+        random += result.getLast()*power (x, i);
+      }
+      if (result.getLast().isOne() || result.getLast().isZero())
+        zeroOneOccured= true;
+    }
+    if (find (list, random))
+    {
+      zeroOneOccured= false;
+      allEqual= false;
+      result= CFList();
+      continue;
+    }
+    if (zeroOneOccured)
+    {
+      list.append (random);
+      zeroOneOccured= false;
+      allEqual= false;
+      result= CFList();
+      continue;
+    }
+    // no zero at this point
+    if (k > 1)
+    {
+      allEqual= true;
+      CFIterator iter= random;
+      buf= iter.coeff();
+      iter++;
+      for (; iter.hasTerms(); iter++)
+        if (buf != iter.coeff())
+          allEqual= false;
+    }
+    if (allEqual)
+    {
+      list.append (random);
+      allEqual= false;
+      zeroOneOccured= false;
+      result= CFList();
+      continue;
+    }
+
+    Feval= F;
+    Geval= G;
+    CanonicalForm LCeval= LCF;
+    j= 1;
+    for (CFListIterator i= result; i.hasItem(); i++, j++)
+    {
+      Feval= Feval (i.getItem(), j);
+      Geval= Geval (i.getItem(), j);
+      LCeval= LCeval (i.getItem(), j);
+    }
+
+    if (LCeval.isZero())
+    {
+      if (!find (list, random))
+        list.append (random);
+      zeroOneOccured= false;
+      allEqual= false;
+      result= CFList();
+      continue;
+    }
+
+    if (list.length() >= bound)
+    {
+      fail= true;
+      break;
+    }
+  } while (find (list, random));
+
+  return result;
+}
+
+/// multiply two lists componentwise
+void mult (CFList& L1, const CFList& L2)
+{
+  ASSERT (L1.length() == L2.length(), "lists of the same size expected");
+
+  CFListIterator j= L2;
+  for (CFListIterator i= L1; i.hasItem(); i++, j++)
+    i.getItem() *= j.getItem();
+}
+
+void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,
+           CanonicalForm& Beval, const CFList& L)
+{
+  Aeval= A;
+  Beval= B;
+  int j= 1;
+  for (CFListIterator i= L; i.hasItem(); i++, j++)
+  {
+    Aeval= Aeval (i.getItem(), j);
+    Beval= Beval (i.getItem(), j);
+  }
+}
+
+CanonicalForm 
+monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
+                     const CanonicalForm& skeleton, const Variable& alpha,
+                     bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
+                    )
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+  if (F.isZero() && degree(G) > 0) return G/Lc(G);
+  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
+  else if (F.isZero() && G.isZero()) return F.genOne();
+  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
+  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
+  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
+  if (F == G) return F/Lc(F);
+
+  CFMap M,N;
+  int best_level= myCompress (A, B, M, N, false);
+
+  if (best_level == 0)
+    return B.genOne();
+
+  A= M(A);
+  B= M(B);
+
+  Variable x= Variable (1);
+  ASSERT (degree (A, y) == 0, "expected degree (F, 1) == 0");
+  ASSERT (degree (B, y) == 0, "expected degree (G, 1) == 0");
+
+  //univariate case 
+  if (A.isUnivariate() && B.isUnivariate())
+    return N (gcd (A, B));
+
+  CanonicalForm skel= M(skeleton);
+  CanonicalForm cA, cB;    // content of A and B
+  CanonicalForm ppA, ppB;    // primitive part of A and B
+  CanonicalForm gcdcAcB;
+  cA = uni_content (A);
+  cB = uni_content (B);
+  gcdcAcB= gcd (cA, cB);
+  ppA= A/cA;
+  ppB= B/cB;
+
+  CanonicalForm lcA, lcB;  // leading coefficients of A and B
+  CanonicalForm gcdlcAlcB;
+  lcA= uni_lcoeff (ppA);
+  lcB= uni_lcoeff (ppB);
+
+  if (fdivides (lcA, lcB))
+  {
+    if (fdivides (A, B))
+      return F/Lc(F);
+  }
+  if (fdivides (lcB, lcA))
+  {
+    if (fdivides (B, A))
+      return G/Lc(G);
+  }
+
+  gcdlcAlcB= gcd (lcA, lcB);
+  int skelSize= size (skel, skel.mvar());
+
+  int j= 0;
+  int biggestSize= 0;
+
+  for (CFIterator i= skel; i.hasTerms(); i++, j++)
+    biggestSize= tmax (biggestSize, size (i.coeff()));
+
+  CanonicalForm g, Aeval, Beval;
+
+  CFList evalPoints;
+  bool evalFail= false;
+  CFList list;
+  bool GF= false;
+  CanonicalForm LCA= LC (A);
+  CanonicalForm tmp;
+  CFArray gcds= CFArray (biggestSize);
+  CFList * pEvalPoints= new CFList [biggestSize];
+  Variable V_buf= alpha;
+  CFList source, dest;
+  CanonicalForm prim_elem, im_prim_elem;
+  for (int i= 0; i < biggestSize; i++)
+  {
+    if (i == 0)
+      evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
+                                    list);
+    else
+    {
+      mult (evalPoints, pEvalPoints [0]);
+      eval (A, B, Aeval, Beval, evalPoints);
+    }
+
+    if (evalFail)
+    {
+      if (V_buf != Variable (1))
+      {
+        do
+        {
+          int num_vars= tmin (getNumVars(A), getNumVars(B));
+          int d= totaldegree (A, Variable(2), Variable (A.level()));
+          d= tmin (d, totaldegree (B, Variable(2), Variable (B.level())));
+          Variable V_buf2= V_buf;
+          choose_extension (d, num_vars, V_buf2);
+          source= CFList();
+          dest= CFList();
+
+          bool prim_fail= false;
+          Variable V_buf3;
+          prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
+
+          ASSERT (!prim_fail, "failure in integer factorizer");
+          if (prim_fail)
+            ; //ERROR
+          else
+            im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
+
+          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
+          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
+
+          for (CFListIterator j= list; j.hasItem(); j++)
+            j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,
+                                im_prim_elem, source, dest);
+          for (int k= 0; k < i; k++)
+          {
+            for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
+              j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,
+                                  im_prim_elem, source, dest);
+            gcds[k]= mapUp (gcds[k], V_buf, V_buf2, prim_elem, im_prim_elem,
+                            source, dest);
+          }
+
+          if (alpha != Variable (1))
+          {
+            A= mapUp (A, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
+            B= mapUp (B, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
+          }
+          evalFail= false;
+          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
+                                        evalFail, list);
+        } while (evalFail);
+      }
+      else
+      {
+        CanonicalForm mipo;
+        int deg= 2;
+        do {
+          mipo= randomIrredpoly (deg, x);
+          V_buf= rootOf (mipo);
+          evalFail= false;
+          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
+                                        evalFail, list); 
+          deg++;
+        } while (evalFail);
+      }
+    }
+
+    g= gcd (Aeval, Beval);
+    g /= Lc (g);
+
+    if (degree (g) != degree (skel) || (skelSize != size (g)))
+    {
+      delete[] pEvalPoints;
+      fail= true;
+      return 0;
+    }
+    CFIterator l= skel;
+    for (CFIterator k= g; k.hasTerms(); k++, l++)
+    {
+      if (k.exp() != l.exp())
+      {
+        delete[] pEvalPoints;
+        fail= true;
+        return 0;
+      }
+    }
+    pEvalPoints[i]= evalPoints;
+    gcds[i]= g;
+
+    tmp= 0;
+    int j= 0;
+    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
+      tmp += k.getItem()*power (x, j);
+    list.append (tmp);
+  }
+
+  if (Monoms.size() == 0)
+    Monoms= getMonoms (skel);
+  if (coeffMonoms == NULL)
+    coeffMonoms= new CFArray [skelSize];
+  j= 0;
+  for (CFIterator i= skel; i.hasTerms(); i++, j++)
+    coeffMonoms[j]= getMonoms (i.coeff());
+
+  CFArray* pL= new CFArray [skelSize];
+  CFArray* pM= new CFArray [skelSize];
+  for (int i= 0; i < biggestSize; i++)
+  {
+    CFIterator l= gcds [i];
+    evalPoints= pEvalPoints [i];
+    for (int k= 0; k < skelSize; k++, l++)
+    {
+      if (i == 0)
+        pL[k]= CFArray (biggestSize);
+      pL[k] [i]= l.coeff();
+
+      if (i == 0)
+        pM[k]= evaluate (coeffMonoms [k], evalPoints);
+    }
+  }
+
+  CFArray solution;
+  CanonicalForm result= 0;
+  int ind= 0;
+  CFArray bufArray;
+  CFMatrix Mat;
+  for (int k= 0; k < skelSize; k++)
+  {
+    if (biggestSize != coeffMonoms[k].size())
+    {
+      bufArray= CFArray (coeffMonoms[k].size());
+      for (int i= 0; i < coeffMonoms[k].size(); i++)
+        bufArray [i]= pL[k] [i];
+      solution= solveGeneralVandermonde (pM [k], bufArray);
+    }
+    else
+      solution= solveGeneralVandermonde (pM [k], pL[k]);
+
+    if (solution.size() == 0)
+    {
+      delete[] pEvalPoints;
+      delete[] pM;
+      delete[] pL;
+      delete[] coeffMonoms;
+      fail= true;
+      return 0;
+    }
+    for (int l= 0; l < solution.size(); l++)
+      result += solution[l]*Monoms [ind + l];
+    ind += solution.size();
+  }
+
+  delete[] pEvalPoints;
+  delete[] pM;
+  delete[] pL;
+
+  if (alpha != Variable (1) && V_buf != alpha)
+  {
+    CFList u, v;
+    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
+  }
+
+  result= N(result);
+  if (fdivides (result, F) && fdivides (result, G))
+    return result;
+  else
+  {
+    delete[] coeffMonoms;
+    fail= true;
+    return 0;
+  }
+}
+
+CanonicalForm
+nonMonicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
+                        const CanonicalForm& skeleton, const Variable& alpha,
+                        bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
+                       )
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+  if (F.isZero() && degree(G) > 0) return G/Lc(G);
+  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
+  else if (F.isZero() && G.isZero()) return F.genOne();
+  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
+  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
+  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
+  if (F == G) return F/Lc(F);
+
+  CFMap M,N;
+  int best_level= myCompress (A, B, M, N, false);
+
+  if (best_level == 0)
+    return B.genOne();
+
+  A= M(A);
+  B= M(B);
+
+  Variable x= Variable (1);
+  ASSERT (degree (A, y) == 0, "expected degree (F, 1) == 0");
+  ASSERT (degree (B, y) == 0, "expected degree (G, 1) == 0");
+
+  //univariate case 
+  if (A.isUnivariate() && B.isUnivariate())
+    return N (gcd (A, B));
+
+  CanonicalForm skel= M(skeleton);
+
+  CanonicalForm cA, cB;    // content of A and B
+  CanonicalForm ppA, ppB;    // primitive part of A and B
+  CanonicalForm gcdcAcB;
+  cA = uni_content (A);
+  cB = uni_content (B);
+  gcdcAcB= gcd (cA, cB);
+  ppA= A/cA;
+  ppB= B/cB;
+
+  CanonicalForm lcA, lcB;  // leading coefficients of A and B
+  CanonicalForm gcdlcAlcB;
+  lcA= uni_lcoeff (ppA);
+  lcB= uni_lcoeff (ppB);
+
+  if (fdivides (lcA, lcB))
+  {
+    if (fdivides (A, B))
+      return F/Lc(F);
+  }
+  if (fdivides (lcB, lcA))
+  {
+    if (fdivides (B, A))
+      return G/Lc(G);
+  }
+
+  gcdlcAlcB= gcd (lcA, lcB);
+  int skelSize= size (skel, skel.mvar());
+
+  int j= 0;
+  int biggestSize= 0;
+  int bufSize;
+  int numberUni= 0;
+  for (CFIterator i= skel; i.hasTerms(); i++, j++)
+  {
+    bufSize= size (i.coeff());
+    biggestSize= tmax (biggestSize, bufSize);
+    numberUni += bufSize;
+  }
+  numberUni--;
+  numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
+  biggestSize= tmax (biggestSize , numberUni);
+
+  numberUni= biggestSize + size (LC(skel)) - 1;
+  int biggestSize2= tmax (numberUni, biggestSize);
+
+  CanonicalForm g, Aeval, Beval;
+
+  CFList evalPoints;
+  CFArray coeffEval;
+  bool evalFail= false;
+  CFList list;
+  bool GF= false;
+  CanonicalForm LCA= LC (A);
+  CanonicalForm tmp;
+  CFArray gcds= CFArray (biggestSize);
+  CFList * pEvalPoints= new CFList [biggestSize];
+  Variable V_buf= alpha;
+  CFList source, dest;
+  CanonicalForm prim_elem, im_prim_elem;
+  for (int i= 0; i < biggestSize; i++)
+  {
+    if (i == 0)
+    {
+      if (getCharacteristic() > 3)
+        evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
+                                    evalFail, list);
+      else
+        evalFail= true;
+
+      if (evalFail)
+      {
+        if (V_buf != Variable (1))
+        {
+          do
+          {
+            int num_vars= tmin (getNumVars(A), getNumVars(B));
+            int d= totaldegree (A, Variable(2), Variable (A.level()));
+            d= tmin (d, totaldegree (B, Variable(2), Variable (B.level())));
+            Variable V_buf2= V_buf;
+            choose_extension (d, num_vars, V_buf2);
+            source= CFList();
+            dest= CFList();
+
+            bool prim_fail= false;
+            Variable V_buf3;
+            prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
+
+            ASSERT (!prim_fail, "failure in integer factorizer");
+            if (prim_fail)
+              ; //ERROR
+            else
+              im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
+
+            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
+            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
+
+            for (CFListIterator i= list; i.hasItem(); i++) 
+              i.getItem()= mapUp (i.getItem(), V_buf, V_buf2, prim_elem,
+                                im_prim_elem, source, dest);
+            evalFail= false;
+            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
+                                          evalFail, list);
+          } while (evalFail);
+        }
+        else
+        {
+          CanonicalForm mipo;
+          int deg= 2;
+          do {
+            mipo= randomIrredpoly (deg, x);
+            V_buf= rootOf (mipo);
+            evalFail= false;
+            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
+                                          evalFail, list);
+            deg++;
+          } while (evalFail);
+        }
+      }
+    }
+    else
+    {
+      mult (evalPoints, pEvalPoints[0]);
+      eval (A, B, Aeval, Beval, evalPoints);
+    }
+
+    g= gcd (Aeval, Beval);
+    g /= Lc (g);
+
+    if (degree (g) != degree (skel) || (skelSize != size (g)))
+    {
+      delete[] pEvalPoints;
+      fail= true;
+      return 0;
+    }
+    CFIterator l= skel;
+    for (CFIterator k= g; k.hasTerms(); k++, l++)
+    {
+      if (k.exp() != l.exp())
+      {
+        delete[] pEvalPoints;
+        fail= true;
+        return 0;
+      }
+    }
+    pEvalPoints[i]= evalPoints;
+    gcds[i]= g;
+
+    tmp= 0;
+    int j= 0;
+    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
+      tmp += k.getItem()*power (x, j);
+    list.append (tmp);
+  }
+
+  if (Monoms.size() == 0)
+    Monoms= getMonoms (skel);
+
+  if (coeffMonoms == NULL)
+    coeffMonoms= new CFArray [skelSize];
+
+  j= 0;
+  for (CFIterator i= skel; i.hasTerms(); i++, j++)
+    coeffMonoms[j]= getMonoms (i.coeff());
+
+  int minimalColumnsIndex;
+  if (skelSize > 1)
+    minimalColumnsIndex= 1;
+  else
+    minimalColumnsIndex= 0;
+  int minimalColumns;
+
+  CFArray* pM= new CFArray [skelSize];
+  CFMatrix Mat;
+  // find the Matrix with minimal number of columns 
+  for (int i= 0; i < skelSize; i++)
+  {
+    pM[i]= CFArray (coeffMonoms[i].size());
+    if (i == 1)
+      minimalColumns= coeffMonoms[i].size();
+    if (i > 1)
+    {
+      if (minimalColumns > coeffMonoms[i].size())
+      {
+        minimalColumns= coeffMonoms[i].size();
+        minimalColumnsIndex= i;
+      }
+    }
+  }
+  CFMatrix* pMat= new CFMatrix [2];
+  pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
+  pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
+  CFArray* pL= new CFArray [skelSize];
+  for (int i= 0; i < biggestSize; i++)
+  {
+    CFIterator l= gcds [i];
+    evalPoints= pEvalPoints [i];
+    for (int k= 0; k < skelSize; k++, l++)
+    {
+      if (i == 0)
+        pL[k]= CFArray (biggestSize);
+      pL[k] [i]= l.coeff(); 
+
+      if (i == 0 && k != 0 && k != minimalColumnsIndex)
+        pM[k]= evaluate (coeffMonoms[k], evalPoints);
+      else if (i == 0 && (k == 0 || k == minimalColumnsIndex))
+        coeffEval= evaluate (coeffMonoms[k], evalPoints);
+      if (i > 0 && (k == 0 || k == minimalColumnsIndex))
+        coeffEval= evaluate (coeffMonoms[k], evalPoints);
+
+      if (k == 0)
+      {
+        if (pMat[k].rows() >= i + 1)
+        {
+          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
+            pMat[k] (i + 1, ind)= coeffEval[ind - 1];
+        }
+      }
+      if (k == minimalColumnsIndex)
+      {
+        if (pMat[1].rows() >= i + 1)
+        {
+          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
+            pMat[1] (i + 1, ind)= coeffEval[ind - 1];
+        }
+      }
+    }
+  }
+
+  CFArray solution;
+  CanonicalForm result= 0;
+  int ind= 1;
+  int matRows, matColumns;
+  matRows= pMat[1].rows();
+  matColumns= pMat[0].columns() - 1; 
+  matColumns += pMat[1].columns();
+
+  Mat= CFMatrix (matRows, matColumns);
+  for (int i= 1; i <= matRows; i++)
+    for (int j= 1; j <= pMat[1].columns(); j++)
+      Mat (i, j)= pMat[1] (i, j);
+
+  for (int j= pMat[1].columns() + 1; j <= matColumns; 
+       j++, ind++)
+  {
+    for (int i= 1; i <= matRows; i++) 
+      Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
+  }
+
+  CFArray firstColumn= CFArray (pMat[0].rows());
+  for (int i= 0; i < pMat[0].rows(); i++)
+    firstColumn [i]= pMat[0] (i + 1, 1); 
+
+  CFArray bufArray= pL[minimalColumnsIndex];
+
+  for (int i= 0; i < biggestSize; i++)
+    pL[minimalColumnsIndex] [i] *= firstColumn[i];
+
+  CFMatrix bufMat= pMat[1];
+  pMat[1]= Mat;
+
+  if (V_buf != x)
+    solution= solveSystemFq (pMat[1], 
+                             pL[minimalColumnsIndex], V_buf);
+  else
+    solution= solveSystemFp (pMat[1], 
+                             pL[minimalColumnsIndex]);
+
+  if (solution.size() == 0)
+  { //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
+    CFMatrix bufMat0= pMat[0];
+    delete [] pMat;
+    pMat= new CFMatrix [skelSize];
+    pL[minimalColumnsIndex]= bufArray; 
+    CFList* bufpEvalPoints;
+    CFArray bufGcds;
+    if (biggestSize != biggestSize2)
+    {
+      bufpEvalPoints= pEvalPoints;
+      pEvalPoints= new CFList [biggestSize2];
+      bufGcds= gcds;
+      gcds= CFArray (biggestSize2);
+      for (int i= 0; i < biggestSize; i++)
+      {
+        pEvalPoints[i]= bufpEvalPoints [i];
+        gcds[i]= bufGcds[i];
+      }
+      for (int i= 0; i < biggestSize2 - biggestSize; i++)
+      {
+        mult (evalPoints, pEvalPoints[0]);
+        eval (A, B, Aeval, Beval, evalPoints);
+        g= gcd (Aeval, Beval);
+        g /= Lc (g);
+        if (degree (g) != degree (skel) || (skelSize != size (g)))
+        {
+          delete[] pEvalPoints;
+          delete[] pMat;
+          delete[] pL;
+          delete[] coeffMonoms;
+          delete[] pM;
+          fail= true;
+          return 0;
+        }
+        CFIterator l= skel;
+        for (CFIterator k= g; k.hasTerms(); k++, l++)
+        {
+          if (k.exp() != l.exp())
+          {
+            delete[] pEvalPoints;
+            delete[] pMat;
+            delete[] pL;
+            delete[] coeffMonoms;
+            delete[] pM;
+            fail= true;
+            return 0;
+          }
+        }
+        pEvalPoints[i + biggestSize]= evalPoints;
+        gcds[i + biggestSize]= g;
+      }
+    }
+    for (int i= 0; i < biggestSize; i++)
+    {
+      CFIterator l= gcds [i];
+      evalPoints= pEvalPoints [i];
+      for (int k= 1; k < skelSize; k++, l++)
+      {
+        if (i == 0)
+          pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
+        if (k == minimalColumnsIndex)
+          continue;
+        coeffEval= evaluate (coeffMonoms[k], evalPoints);
+        if (pMat[k].rows() >= i + 1) 
+        {
+          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
+            pMat[k] (i + 1, ind)= coeffEval[ind - 1];
+        }
+      }
+    }
+    Mat= bufMat0;
+    pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
+    for (int j= 1; j <= Mat.rows(); j++)
+      for (int k= 1; k <= Mat.columns(); k++)
+        pMat [0] (j,k)= Mat (j,k);
+    Mat= bufMat;
+    for (int j= 1; j <= Mat.rows(); j++)
+      for (int k= 1; k <= Mat.columns(); k++)
+        pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
+    // write old matrix entries into new matrices
+    for (int i= 0; i < skelSize; i++)
+    {
+      bufArray= pL[i];
+      pL[i]= CFArray (biggestSize2);
+      for (int j= 0; j < bufArray.size(); j++)
+        pL[i] [j]= bufArray [j];
+    }
+    //write old vector entries into new and add new entries to old matrices
+    for (int i= 0; i < biggestSize2 - biggestSize; i++)
+    {
+      CFIterator l= gcds [i + biggestSize];
+      evalPoints= pEvalPoints [i + biggestSize];
+      for (int k= 0; k < skelSize; k++, l++)
+      {
+        pL[k] [i + biggestSize]= l.coeff(); 
+        coeffEval= evaluate (coeffMonoms[k], evalPoints);
+        if (pMat[k].rows() >= i + biggestSize + 1) 
+        { 
+          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
+            pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
+        }
+      }
+    }
+    // begin new
+    for (int i= 0; i < skelSize; i++)
+    {
+      if (pL[i].size() > 1)
+      {
+        for (int j= 2; j <= pMat[i].rows(); j++) 
+          pMat[i] (j, coeffMonoms[i].size() + j - 1)= 
+              -pL[i] [j - 1];
+      }
+    }
+
+    long rk;
+    matColumns= biggestSize2 - 1;
+    matRows= 0;
+    for (int i= 0; i < skelSize; i++)
+    {
+      if (V_buf == x)
+        rk= gaussianElimFp (pMat[i], pL[i]);
+      else
+        rk= gaussianElimFq (pMat[i], pL[i], V_buf);
+
+      if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
+      {
+        delete[] pEvalPoints;
+        delete[] pMat;
+        delete[] pL;
+        delete[] coeffMonoms;
+        delete[] pM;
+        fail= true;
+        return 0;
+      }
+      matRows += pMat[i].rows() - coeffMonoms[i].size();
+    }
+
+    CFMatrix bufMat;
+    Mat= CFMatrix (matRows, matColumns);
+    ind= 0;
+    bufArray= CFArray (matRows);
+    CFArray bufArray2;
+    for (int i= 0; i < skelSize; i++)
+    {
+      bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
+                       coeffMonoms[i].size() + 1, pMat[i].columns());
+
+      for (int j= 1; j <= bufMat.rows(); j++)
+        for (int k= 1; k <= bufMat.columns(); k++)
+          Mat (j + ind, k)= bufMat(j, k);
+      bufArray2= coeffMonoms[i].size();
+      for (int j= 1; j <= pMat[i].rows(); j++)
+      {
+        if (j > coeffMonoms[i].size())
+          bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
+        else 
+          bufArray2 [j - 1]= pL[i] [j - 1];
+      }
+      pL[i]= bufArray2;
+      ind += bufMat.rows();
+      pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
+    }
+
+    if (V_buf != x)
+      solution= solveSystemFq (Mat, bufArray, V_buf);
+    else
+      solution= solveSystemFp (Mat, bufArray);
+
+    if (solution.size() == 0)
+    {
+      delete[] pEvalPoints;
+      delete[] pMat;
+      delete[] pL;
+      delete[] coeffMonoms;
+      delete[] pM;
+      fail= true;
+      return 0;
+    }
+
+    ind= 0;
+    result= 0;
+    CFArray bufSolution;
+    for (int i= 0; i < skelSize; i++)
+    {
+      bufSolution= readOffSolution (pMat[i], pL[i], solution);
+      for (int i= 0; i < bufSolution.size(); i++, ind++)
+        result += Monoms [ind]*bufSolution[i];
+    }
+    if (alpha != Variable (1) && V_buf != alpha)
+    {
+      CFList u, v;
+      result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
+    }
+    result= N(result);
+    if (fdivides (result, F) && fdivides (result, G))
+      return result;
+    else
+    {
+      fail= true;
+      return 0;
+    }
+  } // end of deKleine, Monagan & Wittkopf
+
+  result += Monoms[0];
+  int ind2= 0, ind3= 2;
+  ind= 0;
+  for (int l= 0; l < minimalColumnsIndex; l++)
+    ind += coeffMonoms[l].size();
+  for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
+       l++, ind2++, ind3++)
+  {
+    result += solution[l]*Monoms [1 + ind2];
+    for (int i= 0; i < pMat[0].rows(); i++) 
+      firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
+  }
+  for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
+    result += solution[l]*Monoms [ind + l];
+  ind= coeffMonoms[0].size();
+  for (int k= 1; k < skelSize; k++)
+  {
+    if (k == minimalColumnsIndex)
+    {
+      ind += coeffMonoms[k].size();
+      continue;
+    }
+    if (k != minimalColumnsIndex)
+    {
+      for (int i= 0; i < biggestSize; i++)
+        pL[k] [i] *= firstColumn [i];
+    }
+
+    if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
+    {
+      bufArray= CFArray (coeffMonoms[k].size());
+      for (int i= 0; i < bufArray.size(); i++)
+        bufArray [i]= pL[k] [i];
+      solution= solveGeneralVandermonde (pM[k], bufArray);
+    }
+    else
+      solution= solveGeneralVandermonde (pM[k], pL[k]);
+
+    if (solution.size() == 0)
+    {
+      delete[] pEvalPoints;
+      delete[] pMat;
+      delete[] pL;
+      delete[] coeffMonoms;
+      delete[] pM;
+      fail= true;
+      return 0;
+    }
+    if (k != minimalColumnsIndex)
+    {
+      for (int l= 0; l < solution.size(); l++)
+        result += solution[l]*Monoms [ind + l];
+      ind += solution.size();
+    }
+  }
+
+  delete[] pEvalPoints;
+  delete[] pMat;
+  delete[] pL;
+  delete[] pM;
+
+  if (alpha != Variable (1) && V_buf != alpha)
+  {
+    CFList u, v;
+    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
+  }
+  result= N(result);
+
+  if (fdivides (result, F) && fdivides (result, G))
+    return result;
+  else
+  {
+    delete[] coeffMonoms;
+    fail= true;
+    return 0;
+  }
+}
+
+CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
+                           const Variable & alpha, CFList& l, bool& topLevel)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+  if (F.isZero() && degree(G) > 0) return G/Lc(G);
+  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
+  else if (F.isZero() && G.isZero()) return F.genOne();
+  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
+  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
+  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
+  if (F == G) return F/Lc(F);
+
+  CFMap M,N;
+  int best_level= myCompress (A, B, M, N, topLevel);
+
+  if (best_level == 0) return B.genOne();
+
+  A= M(A);
+  B= M(B);
+
+  Variable x= Variable (1);
+
+  //univariate case 
+  if (A.isUnivariate() && B.isUnivariate())
+    return N (gcd (A, B));
+
+  int substitute= substituteCheck (A, B);
+
+  if (substitute > 1)
+    subst (A, B, A, B, substitute);
+
+  CanonicalForm cA, cB;    // content of A and B
+  CanonicalForm ppA, ppB;    // primitive part of A and B
+  CanonicalForm gcdcAcB;
+  if (topLevel)
+  {
+    if (best_level <= 2)
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
+    else 
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 
+  }
+  else
+  {
+    cA = uni_content (A);
+    cB = uni_content (B);
+    gcdcAcB= gcd (cA, cB);
+    ppA= A/cA;
+    ppB= B/cB;
+  }
+
+  CanonicalForm lcA, lcB;  // leading coefficients of A and B
+  CanonicalForm gcdlcAlcB;
+  lcA= uni_lcoeff (ppA);
+  lcB= uni_lcoeff (ppB);
+
+  if (fdivides (lcA, lcB))
+  {
+    if (fdivides (A, B))
+      return F/Lc(F);
+  }
+  if (fdivides (lcB, lcA))
+  {
+    if (fdivides (B, A))
+      return G/Lc(G);
+  }
+
+  gcdlcAlcB= gcd (lcA, lcB);
+
+  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
+  int d0;
+
+  if (d == 0) 
+  {
+    if (substitute > 1)
+      return N(reverseSubst (gcdcAcB, substitute));
+    else
+      return N(gcdcAcB);
+  }
+  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
+
+  if (d0 < d)
+    d= d0;
+
+  if (d == 0)
+  {
+    if (substitute > 1)
+      return N(reverseSubst (gcdcAcB, substitute));
+    else
+      return N(gcdcAcB);
+  }
+
+  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
+  CanonicalForm newtonPoly= 1;
+  m= gcdlcAlcB;
+  G_m= 0;
+  H= 0;
+  bool fail= false;
+  topLevel= false;
+  bool inextension= false;
+  Variable V_buf= alpha;
+  CanonicalForm prim_elem, im_prim_elem;
+  CFList source, dest;
+  do // first do
+  {
+    random_element= randomElement (m, V_buf, l, fail);
+    if (random_element == 0 && !fail)
+    {
+      if (!find (l, random_element))
+        l.append (random_element);
+      continue;
+    }
+    if (fail)
+    {
+      source= CFList();
+      dest= CFList();
+      int num_vars= tmin (getNumVars(A), getNumVars(B));
+      num_vars--;
+
+      choose_extension (d, num_vars, V_buf);
+      bool prim_fail= false;
+      Variable V_buf2;
+      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
+
+      ASSERT (!prim_fail, "failure in integer factorizer");
+      if (prim_fail)
+        ; //ERROR
+      else
+        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
+
+      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
+      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
+      inextension= true;
+      for (CFListIterator i= l; i.hasItem(); i++) 
+        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
+                             im_prim_elem, source, dest);
+      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
+                          source, dest);
+      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
+                         source, dest);
+
+      fail= false;
+      random_element= randomElement (m, V_buf, l, fail );
+      DEBOUTLN (cerr, "fail= " << fail);
+      CFList list;
+      TIMING_START (gcd_recursion);
+      G_random_element=
+      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
+                        list, topLevel);
+      TIMING_END_AND_PRINT (gcd_recursion,
+                            "time for recursive call: ");
+      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+    }
+    else
+    {
+      CFList list;
+      TIMING_START (gcd_recursion);
+      G_random_element=
+      sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
+                        list, topLevel);
+      TIMING_END_AND_PRINT (gcd_recursion,
+                            "time for recursive call: ");
+      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+    }
+
+    d0= totaldegree (G_random_element, Variable(2),
+                     Variable (G_random_element.level()));
+    if (d0 == 0)
+    {
+      if (substitute > 1)
+        return N(reverseSubst (gcdcAcB, substitute));
+      else
+        return N(gcdcAcB);
+    }
+    if (d0 >  d)
+    {
+      if (!find (l, random_element))
+        l.append (random_element);
+      continue;
+    }
+
+    G_random_element=
+    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
+    * G_random_element;
+
+    skeleton= G_random_element;
+    d0= totaldegree (G_random_element, Variable(2), 
+                     Variable(G_random_element.level()));
+    if (d0 <  d)
+    {
+      m= gcdlcAlcB;
+      newtonPoly= 1;
+      G_m= 0;
+      d= d0;
+    }
+
+    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
+    if (uni_lcoeff (H) == gcdlcAlcB)
+    {
+      cH= uni_content (H);
+      ppH= H/cH;
+      if (inextension)
+      {
+        CFList u, v;
+        //maybe it's better to test if ppH is an element of F(\alpha) before
+        //mapping down
+        DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
+        ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
+        ppH /= Lc(ppH);
+        DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
+        if (fdivides (ppH, A) && fdivides (ppH, B))
+        {
+          if (substitute > 1)
+          {
+            ppH= reverseSubst (ppH, substitute);
+            gcdcAcB= reverseSubst (gcdcAcB, substitute);
+          }
+          return N(gcdcAcB*ppH);
+        }
+      }
+      else if (fdivides (ppH, A) && fdivides (ppH, B))
+      {
+        if (substitute > 1)
+        {
+          ppH= reverseSubst (ppH, substitute);
+          gcdcAcB= reverseSubst (gcdcAcB, substitute);
+        }
+        return N(gcdcAcB*ppH);
+      }
+    }
+    G_m= H;
+    newtonPoly= newtonPoly*(x - random_element);
+    m= m*(x - random_element);
+    if (!find (l, random_element))
+      l.append (random_element);
+
+    if (getCharacteristic () > 3 && size (skeleton) < 100
+        || size (skeleton) < 10)
+    {
+      CFArray Monoms;
+      CFArray *coeffMonoms= NULL;
+      do //second do
+      {
+        random_element= randomElement (m, V_buf, l, fail);
+        if (random_element == 0 && !fail)
+        {
+          if (!find (l, random_element))
+            l.append (random_element);
+          continue;
+        }
+        if (fail)
+        {
+          source= CFList();
+          dest= CFList();
+          int num_vars= tmin (getNumVars(A), getNumVars(B));
+          num_vars--;
+
+          choose_extension (d, num_vars, V_buf);
+          bool prim_fail= false;
+          Variable V_buf2;
+          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
+
+          ASSERT (!prim_fail, "failure in integer factorizer");
+          if (prim_fail)
+            ; //ERROR
+          else
+            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
+
+          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
+          DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
+          inextension= true;
+          for (CFListIterator i= l; i.hasItem(); i++) 
+            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
+                                im_prim_elem, source, dest);
+          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
+                              source, dest);
+          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+
+          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
+                            source, dest);
+
+          fail= false;
+          random_element= randomElement (m, V_buf, l, fail);
+          DEBOUTLN (cerr, "fail= " << fail);
+          CFList list;
+          TIMING_START (gcd_recursion);
+
+          //sparseInterpolation
+          bool sparseFail= false;
+          if (LC (skeleton).inCoeffDomain())
+            G_random_element= 
+            monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
+                                skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
+          else
+            G_random_element=
+            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
+                                    skeleton, V_buf, sparseFail, coeffMonoms,
+                                    Monoms);
+          if (sparseFail)
+            break;
+
+          TIMING_END_AND_PRINT (gcd_recursion,
+                                "time for recursive call: ");
+          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+        }
+        else
+        {
+          CFList list;
+          TIMING_START (gcd_recursion);
+          bool sparseFail= false;
+          if (LC (skeleton).inCoeffDomain())
+            G_random_element= 
+            monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
+                                skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
+          else
+            G_random_element=
+            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
+                                    skeleton, V_buf, sparseFail, coeffMonoms, 
+                                    Monoms);
+          if (sparseFail)
+            break;
+
+          TIMING_END_AND_PRINT (gcd_recursion,
+                                "time for recursive call: ");
+          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+        }
+
+        d0= totaldegree (G_random_element, Variable(2),
+                        Variable (G_random_element.level()));
+        if (d0 == 0)
+        {
+          if (substitute > 1)
+            return N(reverseSubst (gcdcAcB, substitute));
+          else
+            return N(gcdcAcB);
+        }
+        if (d0 >  d)
+        {
+          if (!find (l, random_element))
+            l.append (random_element);
+          continue;
+        }
+
+        G_random_element=
+        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
+        * G_random_element;
+
+        d0= totaldegree (G_random_element, Variable(2),
+                        Variable(G_random_element.level()));
+        if (d0 <  d)
+        {
+          m= gcdlcAlcB;
+          newtonPoly= 1;
+          G_m= 0;
+          d= d0;
+        }
+
+        TIMING_START (newton_interpolation);
+        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
+        TIMING_END_AND_PRINT (newton_interpolation,
+                              "time for newton interpolation: ");
+
+        //termination test
+        if (uni_lcoeff (H) == gcdlcAlcB)
+        {
+          cH= uni_content (H);
+          ppH= H/cH;
+          if (inextension)
+          {
+            CFList u, v;
+            //maybe it's better to test if ppH is an element of F(\alpha) before
+            //mapping down
+            DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
+            ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
+            ppH /= Lc(ppH);
+            DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
+            if (fdivides (ppH, A) && fdivides (ppH, B))
+            {
+              if (substitute > 1)
+              {
+                ppH= reverseSubst (ppH, substitute);
+                gcdcAcB= reverseSubst (gcdcAcB, substitute);
+              }
+              return N(gcdcAcB*ppH);
+            }
+          }
+          else if (fdivides (ppH, A) && fdivides (ppH, B))
+          {
+            if (substitute > 1)
+            {
+              ppH= reverseSubst (ppH, substitute);
+              gcdcAcB= reverseSubst (gcdcAcB, substitute);
+            }
+            return N(gcdcAcB*ppH);
+          }
+        }
+
+        G_m= H;
+        newtonPoly= newtonPoly*(x - random_element);
+        m= m*(x - random_element);
+        if (!find (l, random_element))
+          l.append (random_element);
+
+      } while (1);
+    }
+  } while (1);
+}
+
+CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
+                           bool& topLevel, CFList& l)
+{
+  CanonicalForm A= F;
+  CanonicalForm B= G;
+  if (F.isZero() && degree(G) > 0) return G/Lc(G);
+  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
+  else if (F.isZero() && G.isZero()) return F.genOne();
+  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
+  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
+  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
+  if (F == G) return F/Lc(F);
+
+  CFMap M,N;
+  int best_level= myCompress (A, B, M, N, topLevel);
+
+  if (best_level == 0) return B.genOne();
+
+  A= M(A);
+  B= M(B);
+
+  Variable x= Variable (1);
+
+  //univariate case 
+  if (A.isUnivariate() && B.isUnivariate())
+    return N (gcd (A, B));
+
+  int substitute= substituteCheck (A, B);
+
+  if (substitute > 1)
+    subst (A, B, A, B, substitute);
+
+  CanonicalForm cA, cB;    // content of A and B
+  CanonicalForm ppA, ppB;    // primitive part of A and B
+  CanonicalForm gcdcAcB;
+  if (topLevel)
+  {
+    if (best_level <= 2)
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
+    else
+      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);
+  }
+  else
+  {
+    cA = uni_content (A);
+    cB = uni_content (B);
+    gcdcAcB= gcd (cA, cB);
+    ppA= A/cA;
+    ppB= B/cB;
+  }
+
+  CanonicalForm lcA, lcB;  // leading coefficients of A and B
+  CanonicalForm gcdlcAlcB;
+  lcA= uni_lcoeff (ppA);
+  lcB= uni_lcoeff (ppB);
+
+  if (fdivides (lcA, lcB))
+  {
+    if (fdivides (A, B))
+      return F/Lc(F);
+  }
+  if (fdivides (lcB, lcA))
+  {
+    if (fdivides (B, A))
+      return G/Lc(G);
+  }
+
+  gcdlcAlcB= gcd (lcA, lcB);
+
+  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
+  int d0;
+
+  if (d == 0)
+  {
+    if (substitute > 1)
+      return N(reverseSubst (gcdcAcB, substitute));
+    else
+      return N(gcdcAcB);
+  }
+  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
+
+  if (d0 < d)
+    d= d0;
+
+  if (d == 0)
+  {
+    if (substitute > 1)
+      return N(reverseSubst (gcdcAcB, substitute));
+    else
+      return N(gcdcAcB);
+  }
+
+  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
+  CanonicalForm newtonPoly= 1;
+  m= gcdlcAlcB;
+  G_m= 0;
+  H= 0;
+  bool fail= false;
+  bool topLevel2= topLevel;
+  int loops= 0;
+  topLevel= false;
+  bool inextension= false;
+  bool inextensionextension= false;
+  Variable V_buf, alpha;
+  CanonicalForm prim_elem, im_prim_elem;
+  CFList source, dest;
+  do //first do
+  {
+    if (inextension)
+      random_element= randomElement (m, alpha, l, fail);
+    else
+      random_element= FpRandomElement (m, l, fail);
+    if (random_element == 0 && !fail)
+    {
+      if (!find (l, random_element))
+        l.append (random_element);
+      continue;
+    }
+
+    if (!fail && !inextension)
+    {
+      CFList list;
+      TIMING_START (gcd_recursion);
+      G_random_element=
+      sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
+                   list);
+      TIMING_END_AND_PRINT (gcd_recursion,
+                            "time for recursive call: ");
+      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+    }
+    else if (!fail && inextension)
+    {
+      CFList list;
+      TIMING_START (gcd_recursion);
+      G_random_element=
+      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
+                        list, topLevel);
+      TIMING_END_AND_PRINT (gcd_recursion,
+                            "time for recursive call: ");
+      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+    }
+    else if (fail && !inextension)
+    {
+      source= CFList();
+      dest= CFList();
+      CFList list;
+      CanonicalForm mipo;
+      int deg= 2;
+      do
+      {
+        mipo= randomIrredpoly (deg, x);
+        alpha= rootOf (mipo);
+        inextension= true;
+        fail= false;
+        random_element= randomElement (m, alpha, l, fail); 
+        deg++;
+      } while (fail);
+      list= CFList();
+      TIMING_START (gcd_recursion);
+      G_random_element=
+      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
+                        list, topLevel);
+      TIMING_END_AND_PRINT (gcd_recursion,
+                            "time for recursive call: ");
+      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+    }
+    else if (fail && inextension)
+    {
+      source= CFList();
+      dest= CFList();
+      int num_vars= tmin (getNumVars(A), getNumVars(B));
+      num_vars--;
+      V_buf= alpha;
+      choose_extension (d, num_vars, V_buf);
+      bool prim_fail= false;
+      Variable V_buf2;
+      CanonicalForm prim_elem, im_prim_elem;
+      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
+
+      ASSERT (!prim_fail, "failure in integer factorizer");
+      if (prim_fail)
+        ; //ERROR
+      else
+        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
+
+      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
+      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
+
+      inextensionextension= true;
+      for (CFListIterator i= l; i.hasItem(); i++)
+        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
+                             im_prim_elem, source, dest);
+      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
+                          source, dest);
+      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
+                         source, dest);
+      fail= false;
+      random_element= randomElement (m, V_buf, l, fail );
+      DEBOUTLN (cerr, "fail= " << fail);
+      CFList list;
+      TIMING_START (gcd_recursion);
+      G_random_element=
+      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
+                        list, topLevel);
+      TIMING_END_AND_PRINT (gcd_recursion,
+                            "time for recursive call: ");
+      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+      alpha= V_buf;
+    }
+
+    d0= totaldegree (G_random_element, Variable(2),
+                     Variable (G_random_element.level()));
+    if (d0 == 0)
+    {
+      if (substitute > 1)
+        return N(reverseSubst (gcdcAcB, substitute));
+      else
+        return N(gcdcAcB);
+    } 
+    if (d0 >  d)
+    {
+      if (!find (l, random_element))
+        l.append (random_element);
+      continue;
+    }
+
+    G_random_element=
+    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
+    * G_random_element;
+
+    skeleton= G_random_element;
+
+    d0= totaldegree (G_random_element, Variable(2),
+                     Variable(G_random_element.level()));
+    if (d0 <  d)
+    {
+      m= gcdlcAlcB;
+      newtonPoly= 1;
+      G_m= 0;
+      d= d0;
+    }
+
+    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
+
+    if (uni_lcoeff (H) == gcdlcAlcB)
+    {
+      cH= uni_content (H);
+      ppH= H/cH;
+      ppH /= Lc (ppH);
+      DEBOUTLN (cerr, "ppH= " << ppH);
+
+      if (fdivides (ppH, A) && fdivides (ppH, B))
+      {
+        if (substitute > 1)
+        {
+          ppH= reverseSubst (ppH, substitute);
+          gcdcAcB= reverseSubst (gcdcAcB, substitute);
+        }
+        return N(gcdcAcB*ppH);
+      }
+    }
+    G_m= H;
+    newtonPoly= newtonPoly*(x - random_element);
+    m= m*(x - random_element);
+    if (!find (l, random_element))
+      l.append (random_element);
+
+    if ((getCharacteristic() > 3 && size (skeleton) < 100) ||
+         size (skeleton) < 10)
+    {
+      CFArray Monoms;
+      CFArray* coeffMonoms= NULL;
+
+      do //second do
+      {
+        if (inextension)
+          random_element= randomElement (m, alpha, l, fail);
+        else
+          random_element= FpRandomElement (m, l, fail);
+        if (random_element == 0 && !fail)
+        {
+          if (!find (l, random_element))
+            l.append (random_element);
+          continue;
+        }
+
+        bool sparseFail= false;
+        if (!fail && !inextension)
+        {
+          CFList list;
+          TIMING_START (gcd_recursion);
+          if (LC (skeleton).inCoeffDomain())
+            G_random_element=
+            monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
+                                skeleton, Variable(1), sparseFail, coeffMonoms,
+                                Monoms);
+          else
+            G_random_element=
+            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
+                                    skeleton, Variable (1), sparseFail,
+                                    coeffMonoms, Monoms);
+          TIMING_END_AND_PRINT (gcd_recursion,
+                                "time for recursive call: ");
+          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+        }
+        else if (!fail && inextension)
+        {
+          CFList list;
+          TIMING_START (gcd_recursion);
+          if (LC (skeleton).inCoeffDomain())
+            G_random_element=
+            monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
+                                skeleton, alpha, sparseFail, coeffMonoms,
+                                Monoms);
+          else
+            G_random_element=
+            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
+                                   skeleton, alpha, sparseFail, coeffMonoms,
+                                   Monoms);
+          TIMING_END_AND_PRINT (gcd_recursion,
+                                "time for recursive call: ");
+          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+        }
+        else if (fail && !inextension)
+        {
+          source= CFList();
+          dest= CFList();
+          CFList list;
+          CanonicalForm mipo;
+          int deg= 2;
+          do 
+          {
+            mipo= randomIrredpoly (deg, x);
+            alpha= rootOf (mipo);
+            inextension= true;
+            fail= false;
+            random_element= randomElement (m, alpha, l, fail);
+            deg++;
+          } while (fail);
+          list= CFList();
+          TIMING_START (gcd_recursion);
+          if (LC (skeleton).inCoeffDomain())
+            G_random_element=
+            monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
+                                skeleton, alpha, sparseFail, coeffMonoms,
+                                Monoms);
+          else
+            G_random_element=
+            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
+                                   skeleton, alpha, sparseFail, coeffMonoms,
+                                   Monoms);
+          TIMING_END_AND_PRINT (gcd_recursion,
+                                "time for recursive call: ");
+          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+        }
+        else if (fail && inextension)
+        {
+          source= CFList();
+          dest= CFList();
+          int num_vars= tmin (getNumVars(A), getNumVars(B));
+          num_vars--;
+          V_buf= alpha;
+          choose_extension (d, num_vars, V_buf);
+          bool prim_fail= false;
+          Variable V_buf2;
+          CanonicalForm prim_elem, im_prim_elem;
+          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
+
+          ASSERT (!prim_fail, "failure in integer factorizer");
+          if (prim_fail)
+            ; //ERROR
+          else
+            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
+
+          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
+          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
+
+          inextensionextension= true;
+          for (CFListIterator i= l; i.hasItem(); i++)
+            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
+                                im_prim_elem, source, dest);
+          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
+                              source, dest);
+          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
+          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
+                            source, dest);
+          fail= false;
+          random_element= randomElement (m, V_buf, l, fail );
+          DEBOUTLN (cerr, "fail= " << fail);
+          CFList list;
+          TIMING_START (gcd_recursion);
+          if (LC (skeleton).inCoeffDomain())
+            G_random_element=
+            monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
+                                skeleton, V_buf, sparseFail, coeffMonoms,
+                                Monoms);
+          else
+            G_random_element=
+            nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
+                                    skeleton, V_buf, sparseFail,
+                                    coeffMonoms, Monoms);
+          TIMING_END_AND_PRINT (gcd_recursion,
+                                "time for recursive call: ");
+          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
+          alpha= V_buf;
+        }
+
+        if (sparseFail)
+          break;
+
+        d0= totaldegree (G_random_element, Variable(2),
+                        Variable (G_random_element.level()));
+        if (d0 == 0)
+        {
+          if (substitute > 1)
+            return N(reverseSubst (gcdcAcB, substitute));
+          else
+            return N(gcdcAcB);
+        }
+        if (d0 >  d)
+        {
+          if (!find (l, random_element))
+            l.append (random_element);
+          continue;
+        }
+
+        G_random_element=
+        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
+        * G_random_element;
+
+        d0= totaldegree (G_random_element, Variable(2),
+                        Variable(G_random_element.level()));
+        if (d0 <  d)
+        {
+          m= gcdlcAlcB;
+          newtonPoly= 1;
+          G_m= 0;
+          d= d0;
+        }
+
+        TIMING_START (newton_interpolation);
+        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
+        TIMING_END_AND_PRINT (newton_interpolation,
+                              "time for newton interpolation: ");
+
+        //termination test
+        if (uni_lcoeff (H) == gcdlcAlcB)
+        {
+          cH= uni_content (H);
+          ppH= H/cH;
+          ppH /= Lc (ppH);
+          DEBOUTLN (cerr, "ppH= " << ppH);
+          if (fdivides (ppH, A) && fdivides (ppH, B))
+          {
+            if (substitute > 1)
+            {
+              ppH= reverseSubst (ppH, substitute);
+              gcdcAcB= reverseSubst (gcdcAcB, substitute);
+            }
+            return N(gcdcAcB*ppH);
+          }
+        }
+
+        G_m= H;
+        newtonPoly= newtonPoly*(x - random_element);
+        m= m*(x - random_element);
+        if (!find (l, random_element))
+          l.append (random_element);
+
+      } while (1); //end of second do
+    }
+  } while (1); //end of first do
+}
+
+static inline
+int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
+                    CFMap & N, int& both_non_zero)
+{
+  int n= tmax (F.level(), G.level());
+  int * degsf= new int [n + 1];
+  int * degsg= new int [n + 1];
+
+  for (int i = 0; i <= n; i++)
+    degsf[i]= degsg[i]= 0;
+
+  degsf= degrees (F, degsf);
+  degsg= degrees (G, degsg);
+
+  both_non_zero= 0;
+  int f_zero= 0;
+  int g_zero= 0;
+  int both_zero= 0;
+
+  for (int i= 1; i <= n; i++)
+  {
+    if (degsf[i] != 0 && degsg[i] != 0)
+    {
+      both_non_zero++;
+      continue;
+    }
+    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
+    {
+      f_zero++;
+      continue;
+    }
+    if (degsg[i] == 0 && degsf[i] && i <= F.level())
+    {
+      g_zero++;
+      continue;
+    }
+  }
+
+  if (both_non_zero == 0) return 0;
+
+  // map Variables which do not occur in both polynomials to higher levels
+  int k= 1;
+  int l= 1;
+  for (int i= 1; i <= n; i++)
+  {
+    if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
+    {
+      if (k + both_non_zero != i)
+      {
+        M.newpair (Variable (i), Variable (k + both_non_zero));
+        N.newpair (Variable (k + both_non_zero), Variable (i));
+      }
+      k++;
+    }
+    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
+    {
+      if (l + g_zero + both_non_zero != i)
+      {
+        M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
+        N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
+      }
+      l++;
+    }
+  }
+
+  // sort Variables x_{i} in decreasing order of
+  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
+  int m= tmin (F.level(), G.level());
+  int max_min_deg;
+  k= both_non_zero;
+  l= 0;
+  int i= 1;
+  while (k > 0)
+  {
+    max_min_deg= tmin (degsf[i], degsg[i]);
+    while (max_min_deg == 0)
+    {
+      i++;
+      max_min_deg= tmin (degsf[i], degsg[i]);
+    }
+    for (int j= i + 1; j <= m; j++)
+    {
+      if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
+          (tmin (degsf[j], degsg[j]) != 0))
+      {
+        max_min_deg= tmin (degsf[j], degsg[j]);
+        l= j;
+      }
+    }
+
+    if (l != 0)
+    {
+      if (l != k)
+      {
+        M.newpair (Variable (l), Variable(k));
+        N.newpair (Variable (k), Variable(l));
+        degsf[l]= 0;
+        degsg[l]= 0;
+        l= 0;
+      }
+      else
+      {
+        degsf[l]= 0;
+        degsg[l]= 0;
+        l= 0;
+      }
+    }
+    else if (l == 0)
+    {
+      if (i != k)
+      {
+        M.newpair (Variable (i), Variable (k));
+        N.newpair (Variable (k), Variable (i));
+        degsf[i]= 0;
+        degsg[i]= 0;
+      }
+      else
+      {
+        degsf[i]= 0;
+        degsg[i]= 0;
+      }
+      i++;
+    }
+    k--;
+  }
+
+  delete [] degsf;
+  delete [] degsg;
+
+  return both_non_zero;
+}
+
+static inline
+CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
+                            const CFList& evaluation)
+{
+  CanonicalForm A= F;
+  int k= 2;
+  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
+    A= A (Variable (k) + i.getItem(), k);
+
+  CanonicalForm buf= A;
+  Feval= CFList();
+  Feval.append (buf);
+  for (k= evaluation.length() + 1; k > 2; k--)
+  {
+    buf= mod (buf, Variable (k));
+    Feval.insert (buf);
+  }
+  return A;
+}
+
+static inline
+CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation)
+{
+  int l= evaluation.length() + 1;
+  CanonicalForm result= F;
+  CFListIterator j= evaluation;
+  for (int i= 2; i < l + 1; i++, j++)
+  {
+    if (F.level() < i)
+      continue;
+    result= result (Variable (i) - j.getItem(), i);
+  }
+  return result;
+}
+
+static inline
+bool Hensel_P (const CanonicalForm & U, CFArray & G, const Evaluation & A,
+                const Variable & x, const CFArray& LCs )
+{
+  CFList factors;
+  factors.append (G[1]);
+  factors.append (G[2]);
+  CFList evaluation;
+  for (int i= A.min(); i <= A.max(); i++)
+    evaluation.append (A [i]);
+  CFList UEval;
+  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
+  CFArray shiftedLCs= CFArray (2);
+  CFList shiftedLCsEval1, shiftedLCsEval2;
+  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
+  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
+  CanonicalForm LcUEval= LC (UEval.getFirst(), x);
+  factors.insert (1);
+  int liftBound= degree (UEval.getLast(), 2) + 1;
+  CFArray Pi;
+  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
+  CFList diophant;
+  CFArray lcs= CFArray (2);
+  lcs [0]= shiftedLCsEval1.getFirst();
+  lcs [1]= shiftedLCsEval2.getFirst();
+  henselLift122 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
+                 lcs, false);
+
+  bool noChange= true;
+  for (CFListIterator i= factors; i.hasItem(); i++)
+  {
+    if (degree (i.getItem(), 2) != 0)
+      noChange= false;
+  }
+  if (noChange)
+    return !noChange;
+  int * liftBounds;
+  if (U.level() > 2)
+  {
+    liftBounds= new int [U.level() - 1];
+    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);
+  }
+  G[1]= factors.getFirst();
+  G[2]= factors.getLast();
+  G[1]= myReverseShift (G[1], evaluation);
+  G[2]= myReverseShift (G[2], evaluation);
+  return true;
+}
+
+static inline
+bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
+                 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
+                 FFREvaluation & b, int delta, int degF, int degG, int maxeval,
+                 int & count)
+{
+  if( count == 0 && delta != 0)
+  {
+    if( count++ > maxeval )
+      return false;
+  }
+  if (count > 0)
+  {
+    b.nextpoint();
+    if (count++ > maxeval)
+      return false;
+  }
+  while( true )
+  {
+    Fb = b( F );
+    if( degree( Fb, 1 ) == degF )
+    {
+      Gb = b( G );
+      if( degree( Gb, 1 ) == degG )
+      {
+        Db = gcd( Fb, Gb );
+        if( delta > 0 )
+        {
+          if( degree( Db, 1 ) <= delta )
+            return true;
+        }
+        else
+          return true;
+      }
+    }
+    b.nextpoint();
+    if( count++ > maxeval )
+      return false;
+  }
+}
+
+// parameters for heuristic
+static int maxNumEval= 200;
+static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
+static int sizePerVars2= 290; //try psr gcd if size/#variables is less
+
+/// Extended Zassenhaus GCD for finite fields
+CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
+{
+  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
+  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
+  else if (FF.isZero() && GG.isZero()) return FF.genOne();
+  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
+  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
+  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
+  if (FF == GG) return FF/Lc(FF);
+
+  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
+                lcD;
+  CFArray DD( 1, 2 ), lcDD( 1, 2 );
+  int degF, degG, delta, count;
+  int maxeval;
+  maxeval= tmin((getCharacteristic()/
+                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
+  count= 0; // number of eval. used
+  FFREvaluation b, bt;
+  bool gcdfound = false;
+  Variable x = Variable(1);
+
+  F= FF;
+  G= GG;
+
+  CFMap M,N;
+  int smallestDegLev;
+  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
+
+  if (best_level == 0) return G.genOne();
+
+  F= M (F);
+  G= M (G);
+
+  f = content( F, x ); g = content( G, x );
+  d = gcd( f, g );
+  F /= f; G /= g;
+
+  if( F.isUnivariate() && G.isUnivariate() )
+  {
+    if( F.mvar() == G.mvar() )
+      d *= gcd( F, G );
+    return N (d);
+  }
+
+  int maxNumVars= tmax (getNumVars (F), getNumVars (G));
+  Variable a,bv;
+  int sizeF= size (F);
+  int sizeG= size (G);
+
+  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
+  {
+    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
+      return N (d*GCD_Fp_extension (F, G, a));
+    else if (CFFactory::gettype() == GaloisFieldDomain)
+      return N (d*GCD_GF (F, G));
+    else
+      return N (d*GCD_small_p (F, G));
+  }
+
+  if( gcd_test_one( F, G, false ) )
+  {
+    return N (d);
+  }
+
+  lcF = LC( F, x ); lcG = LC( G, x );
+  lcD = gcd( lcF, lcG );
+
+  delta = 0;
+  degF = degree( F, x ); degG = degree( G, x );
+  if(hasFirstAlgVar(F,a))
+  {
+    if(hasFirstAlgVar(G,bv))
+    {
+      if(bv.level() > a.level())
+        a = bv;
+    }
+    b = FFREvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
+  }
+  else // F not in extension
+  {
+    if(hasFirstAlgVar(G,a))
+      b = FFREvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
+    else
+    { // both not in extension given by algebraic variable
+      if (CFFactory::gettype() != GaloisFieldDomain)
+        b = FFREvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
+      else
+        b = FFREvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
+    }
+  }
+  CanonicalForm cand;
+  while( !gcdfound )
+  {
+    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count ))
+    { // too many eval. used --> try another method
+      if (sizeF/maxNumVars < sizePerVars2 && sizeG/maxNumVars < sizePerVars2)
+      {
+        CanonicalForm dummy1, dummy2;
+        Variable y= Variable (tmax (F.level(), G.level()));
+        Variable z= Variable (smallestDegLev);
+        dummy1= swapvar (F, x, y);
+        dummy2= swapvar (G, x, y);
+        if (size (psr (dummy1, dummy2, y)) <= tmin (size (F), size (G)))
+        {
+          Off (SW_USE_EZGCD_P);
+          CanonicalForm result= swapvar (gcd (dummy1,dummy2), x, y);
+          On (SW_USE_EZGCD_P);
+          return N (d*result);
+        }
+      }
+      if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
+        return N (d*sparseGCDFq (F, G, a));
+      if (CFFactory::gettype() == GaloisFieldDomain)
+        return N (d*GCD_GF (F, G));
+      else
+        return N (d*sparseGCDFp (F,G));
+    }
+    delta = degree( Db );
+    if( delta == 0 )
+      return N (d);
+    while( true )
+    {
+      bt = b;
+      if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count ))
+      { // too many eval. used --> try another method
+        if (sizeF/maxNumVars < sizePerVars2 && sizeG/maxNumVars < sizePerVars2)
+        {
+          CanonicalForm dummy1, dummy2;
+          Variable y= Variable (tmax (F.level(), G.level()));
+          Variable z= Variable (smallestDegLev);
+          dummy1= swapvar (F, x, y);
+          dummy2= swapvar (G, x, y);
+          if (size (psr (dummy1, dummy2, y)) <= tmin (size (F), size (G)))
+          {
+            Off (SW_USE_EZGCD_P);
+            CanonicalForm result= swapvar (gcd (dummy1,dummy2), x, y);
+            On (SW_USE_EZGCD_P);
+            return N (d*result);
+          }
+        }
+        if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
+          return N (d*sparseGCDFq (F, G, a));
+        if (CFFactory::gettype() == GaloisFieldDomain)
+          return N (d*GCD_GF (F, G));
+        else
+          return N (d*sparseGCDFp (F,G));
+      }
+      int dd = degree( Dbt );
+      if( dd == 0 )
+        return N (d);
+      if( dd == delta )
+        break;
+      if( dd < delta )
+      {
+        delta = dd;
+        b = bt;
+        Db = Dbt; Fb = Fbt; Gb = Gbt;
+      }
+    }
+    if( degF <= degG && delta == degF && fdivides( F, G ) )
+      return N (d*F);
+    if( degG < degF && delta == degG && fdivides( G, F ) )
+      return N (d*G);
+    if( delta != degF && delta != degG )
+    {
+      bool B_is_F;
+      CanonicalForm xxx1, xxx2;
+      CanonicalForm buf;
+      DD[1] = Fb / Db;
+      buf= Gb/Db;
+      xxx1 = gcd( DD[1], Db );
+      xxx2 = gcd( buf, Db );
+      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
+          (size (F) <= size (G))) 
+          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
+      {
+        B = F;
+        DD[2] = Db;
+        lcDD[1] = lcF;
+        lcDD[2] = lcD;
+        B_is_F = true;
+      }
+      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
+               (size (G) < size (F))) 
+               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
+      {
+        DD[1] = buf;
+        B = G;
+        DD[2] = Db;
+        lcDD[1] = lcG;
+        lcDD[2] = lcD;
+        B_is_F = false;
+      }
+      else // special case handling
+      {
+        if (sizeF/maxNumVars < sizePerVars2 && sizeG/maxNumVars < sizePerVars2)
+        {
+          CanonicalForm dummy1, dummy2;
+          Variable y= Variable (tmax (F.level(), G.level()));
+          Variable z= Variable (smallestDegLev);
+          dummy1= swapvar (F, x, y);
+          dummy2= swapvar (G, x, y);
+          if (size (psr (dummy1, dummy2, y)) <= tmin (size (F), size (G)))
+          {
+            Off (SW_USE_EZGCD_P);
+            CanonicalForm result= swapvar (gcd (dummy1,dummy2), x, y);
+            On (SW_USE_EZGCD_P);
+            return N (d*result);
+          }
+        }
+        if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
+          return N (d*sparseGCDFq (F, G, a));
+        if (CFFactory::gettype() == GaloisFieldDomain)
+          return N (d*GCD_GF (F, G));
+        else
+          return N (d*sparseGCDFp (F,G));
+      }
+      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
+      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
+
+      gcdfound= Hensel_P (B*lcD, DD, b, x, lcDD);
+
+      if( gcdfound )
+      {
+        cand = DD[2] / content( DD[2], Variable(1) );
+        if( B_is_F )
+          gcdfound = fdivides( cand, G );
+        else
+          gcdfound = fdivides( cand, F );
+      }
+    }
+    delta= 0;
+  }
+  return N (d*cand);
+}
+
+
 #endif

factory/cf_gcd_smallp.h

@@ -63 +63 @@
 }
 
+CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
+                           bool& topLevel, CFList& l);
+
+/// Zippel's sparse GCD over Fp
+static inline
+CanonicalForm sparseGCDFp (const CanonicalForm& A, const CanonicalForm& B)
+{
+  ASSERT (CFFactory::gettype() == GaloisFieldDomain, 
+          "GF as base field expected");
+  CFList list;
+  bool topLevel= true;
+  return sparseGCDFp (A, B, topLevel, list);
+}
+
+/// Zippel's sparse GCD over Fq
+CanonicalForm
+sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
+             const Variable& alpha, CFList& l, bool& topLevel);
+
+static inline
+CanonicalForm sparseGCDFq (const CanonicalForm& A, const CanonicalForm& B,
+                           const Variable& alpha)
+{
+  CFList list;
+  bool topLevel= true;
+  return sparseGCDFq (A, B, alpha, list, topLevel);
+}
+
+/// extended Zassenhaus GCD
+CanonicalForm
+EZGCD_P (const CanonicalForm& A, const CanonicalForm& B);
+
 CanonicalForm
 randomIrredpoly (int i, const Variable & x) ;

factory/cf_map_ext.cc

@@ -340 +340 @@
       return 0;
   } while (1);
-  zz_pE::init (NTL_mipo);
-  zz_pEX NTL_alpha_mipo= convertFacCF2NTLzz_pEX (mipo, NTL_mipo);
-  zz_pE root= FindRoot (NTL_alpha_mipo);
+  zz_pX alpha_mipo= convertFacCF2NTLzzpX (mipo);
+  zz_pE::init (alpha_mipo);
+  zz_pEX NTL_beta_mipo= to_zz_pEX (NTL_mipo);
+  zz_pE root= FindRoot (NTL_beta_mipo);
   return convertNTLzzpE2CF (root, alpha);
 }

factory/cf_random.cc

@@ -65 +65 @@
 CanonicalForm GFRandom::generate () const
 {
-    return CanonicalForm( int2imm_gf( factoryrandom( gf_q ) ) );
+    int i= factoryrandom( gf_q );
+    if ( i == gf_q1 ) i++;
+    return CanonicalForm( int2imm_gf( i ) );
 }
 

factory/facHensel.cc

@@ -26 +26 @@
 #include "facHensel.h"
 #include "cf_util.h"
+#include "fac_util.h"
+#include "cf_algorithm.h"
 
 #ifdef HAVE_NTL
@@ -1553 +1555 @@
 }
 
+void
+henselStep122 (const CanonicalForm& F, const CFList& factors,
+              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
+              CFArray& Pi, int j, const CFArray& LCs)
+{
+  Variable x= F.mvar();
+  CanonicalForm xToJ= power (x, j);
+
+  CanonicalForm E;
+  // compute the error
+  if (degree (Pi [factors.length() - 2], x) > 0)
+    E= F[j] - Pi [factors.length() - 2] [j];
+  else
+    E= F[j];
+
+  CFArray buf= CFArray (diophant.length());
+
+  int k= 0;
+  CanonicalForm remainder;
+  // actual lifting
+  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
+  {
+    if (degree (bufFactors[k], x) > 0)
+      remainder= modNTL (E, bufFactors[k] [0]);
+    else
+      remainder= modNTL (E, bufFactors[k]);
+    buf[k]= mulNTL (i.getItem(), remainder);
+    if (degree (bufFactors[k], x) > 0)
+      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
+    else
+      buf[k]= modNTL (buf[k], bufFactors[k]); 
+  }
+
+  for (k= 0; k < factors.length(); k++)
+    bufFactors[k] += xToJ*buf[k];
+
+  // update Pi [0]
+  int degBuf0= degree (bufFactors[0], x);
+  int degBuf1= degree (bufFactors[1], x);
+  if (degBuf0 > 0 && degBuf1 > 0)
+  {
+    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
+    if (j + 2 <= M.rows())
+      M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
+  }
+  CanonicalForm uIZeroJ;
+  if (degBuf0 > 0 && degBuf1 > 0)
+    uIZeroJ= mulNTL(bufFactors[0][0],buf[1])+mulNTL (bufFactors[1][0], buf[0]);
+  else if (degBuf0 > 0)
+    uIZeroJ= mulNTL (buf[0], bufFactors[1]);
+  else if (degBuf1 > 0)
+    uIZeroJ= mulNTL (bufFactors[0], buf [1]);
+  else
+    uIZeroJ= 0;
+  Pi [0] += xToJ*uIZeroJ;
+
+  CFArray tmp= CFArray (factors.length() - 1);
+  for (k= 0; k < factors.length() - 1; k++)
+    tmp[k]= 0;
+  CFIterator one, two;
+  one= bufFactors [0];
+  two= bufFactors [1];
+  if (degBuf0 > 0 && degBuf1 > 0)
+  {
+    while (one.exp() > j) one++;
+    while (two.exp() > j) two++;
+    for (k= 1; k <= (int) ceil (j/2.0); k++)
+    {
+      if (k != j - k + 1)
+      {
+        if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
+        {
+          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
+                      two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
+          one++;
+          two++;
+        }
+        else if (one.exp() == j - k + 1)
+        {
+          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
+                      M (k + 1, 1);
+          one++;
+        }
+        else if (two.exp() == j - k + 1)
+        {
+          tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
+                    M (k + 1, 1);
+          two++;
+        }
+      }
+      else
+        tmp[0] += M (k + 1, 1);
+    }
+
+    if (j + 2 <= M.rows())
+    {
+      if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
+        tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
+                           (bufFactors [1] [j + 1] + bufFactors [1] [0]))
+                           - M(1,1) - M (j + 2,1);
+      else if (degBuf0 >= j + 1)
+        tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
+      else if (degBuf1 >= j + 1)
+        tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
+    }
+  }
+  Pi [0] += tmp[0]*xToJ*F.mvar();
+
+  /*// update Pi [l]
+  int degPi, degBuf;
+  for (int l= 1; l < factors.length() - 1; l++)
+  {
+    degPi= degree (Pi [l - 1], x);
+    degBuf= degree (bufFactors[l + 1], x);
+    if (degPi > 0 && degBuf > 0)
+      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
+    if (j == 1)
+    {
+      if (degPi > 0 && degBuf > 0)
+        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
+                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
+                  M (1, l + 1));
+      else if (degPi > 0)
+        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
+      else if (degBuf > 0)
+        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
+    }
+    else
+    {
+      if (degPi > 0 && degBuf > 0)
+      {
+        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
+        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]); 
+      }
+      else if (degPi > 0)
+        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
+      else if (degBuf > 0)
+      {
+        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
+        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
+      }
+      Pi[l] += xToJ*uIZeroJ;
+    }
+    one= bufFactors [l + 1];
+    two= Pi [l - 1];
+    if (two.exp() == j + 1)
+    {
+      if (degBuf > 0 && degPi > 0)
+      {
+          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
+          two++;
+      }
+      else if (degPi > 0)
+      {
+          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
+          two++;
+      }
+    }
+    if (degBuf > 0 && degPi > 0)
+    {
+      for (k= 1; k <= (int) ceil (j/2.0); k++)
+      {
+        if (k != j - k + 1)
+        {
+          if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
+          {
+            tmp[l] += mulNTL ((bufFactors[l + 1][k]+one.coeff()),(Pi[l-1] [k] +
+                        two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
+            one++;
+            two++;
+          }
+          else if (one.exp() == j - k + 1)
+          {
+            tmp[l] += mulNTL ((bufFactors[l+1][k]+one.coeff()),Pi[l - 1] [k]) -
+                        M (k + 1, l + 1);
+            one++;
+          }
+          else if (two.exp() == j - k + 1)
+          {
+            tmp[l] += mulNTL (bufFactors[l+1][k],(Pi[l-1] [k] + two.coeff())) -
+                      M (k + 1, l + 1);
+            two++;
+          }
+        }
+        else
+          tmp[l] += M (k + 1, l + 1);
+      }
+    }
+    Pi[l] += tmp[l]*xToJ*F.mvar();
+  }*/
+  return;
+}
+
+void
+henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
+              CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
+{
+  if (sort)
+    sortList (factors, Variable (1));
+  Pi= CFArray (factors.length() - 2);
+  CFList bufFactors2= factors;
+  bufFactors2.removeFirst();
+  diophant.removeFirst();
+  CFListIterator iter= diophant;
+  CanonicalForm s,t;
+  extgcd (bufFactors2.getFirst(), bufFactors2.getLast(), s, t);
+  diophant= CFList();
+  diophant.append (t);
+  diophant.append (s);
+  DEBOUTLN (cerr, "diophant= " << diophant);
+
+  CFArray bufFactors= CFArray (bufFactors2.length());
+  int i= 0;
+  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
+    bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
+
+  Variable x= F.mvar();
+  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
+  {
+    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
+    Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
+                        mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
+  }
+  else if (degree (bufFactors[0], x) > 0)
+  {
+    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
+    Pi [0]= M (1, 1) +
+            mulNTL (bufFactors [0] [1], bufFactors[1])*x;
+  }
+  else if (degree (bufFactors[1], x) > 0)
+  {
+    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
+    Pi [0]= M (1, 1) +
+            mulNTL (bufFactors [0], bufFactors[1] [1])*x;
+  }
+  else
+  {
+    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
+    Pi [0]= M (1, 1);
+  }
+
+  for (i= 1; i < l; i++)
+    henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
+
+  factors= CFList();
+  for (i= 0; i < bufFactors.size(); i++)
+    factors.append (bufFactors[i]);
+  return;
+}
+
+/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$ 
+/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
+static inline
+CFList
+diophantine (const CFList& recResult, const CFList& factors, 
+             const CFList& products, const CFList& M, const CanonicalForm& E)
+{
+  if (M.isEmpty())
+  {
+    CFList result;
+    CFListIterator j= factors;
+    CanonicalForm buf;
+    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
+    {
+      ASSERT (E.isUnivariate() || E.inCoeffDomain(),
+              "constant or univariate poly expected");
+      ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
+              "constant or univariate poly expected");
+      ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
+              "constant or univariate poly expected");
+      buf= mulNTL (E, i.getItem());
+      result.append (modNTL (buf, j.getItem()));
+    }
+    return result;
+  }
+  Variable y= M.getLast().mvar();
+  CFList bufFactors= factors;
+  for (CFListIterator i= bufFactors; i.hasItem(); i++)
+    i.getItem()= mod (i.getItem(), y);
+  CFList bufProducts= products;
+  for (CFListIterator i= bufProducts; i.hasItem(); i++)
+    i.getItem()= mod (i.getItem(), y);
+  CFList buf= M;
+  buf.removeLast();
+  CanonicalForm bufE= mod (E, y);
+  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
+                                      bufE);
+
+  CanonicalForm e= E;
+  CFListIterator j= products;
+  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
+    e -= mulMod (i.getItem(), j.getItem(), M);
+
+  CFList result= recDiophantine;
+  int d= degree (M.getLast());
+  CanonicalForm coeffE;
+  for (int i= 1; i < d; i++)
+  {
+    if (degree (e, y) > 0)
+      coeffE= e[i];
+    else
+      coeffE= 0;
+    if (!coeffE.isZero())
+    {
+      CFListIterator k= result;
+      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
+                                   coeffE);
+      CFListIterator l= products;
+      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
+      {
+        k.getItem() += j.getItem()*power (y, i);
+        e -= mulMod (j.getItem()*power (y, i), l.getItem(), M);
+      }
+    }
+    if (e.isZero())
+      break;
+  }
+  return result;
+}
+
+// non monic case
+static inline
+CFList
+multiRecDiophantine2 (const CanonicalForm& F, const CFList& factors,
+                     const CFList& recResult, const CFList& M, const int d)
+{
+  Variable y= F.mvar();
+  CFList result;
+  CFListIterator i;
+  CanonicalForm e= 1;
+  CFListIterator j= factors;
+  CFList p;
+  CFArray bufFactors= CFArray (factors.length());
+  CanonicalForm yToD= power (y, d);
+  int k= 0;
+  for (CFListIterator i= factors; i.hasItem(); i++, k++)
+    bufFactors [k]= i.getItem();
+  CanonicalForm b;
+  CFList buf= M;
+  buf.removeLast();
+  buf.append (yToD);
+  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
+  {
+    b= 1;
+    if (fdivides (bufFactors[k], F))
+      b= F/bufFactors[k];
+    else 
+    {
+      for (int l= 0; l < factors.length(); l++)
+      {
+        if (l == k)
+          continue;
+        else
+        {
+          b= mulMod (b, bufFactors[l], buf);
+        }
+      }
+    }
+    p.append (b);
+  }
+  j= p;
+  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
+    e -= mulMod (i.getItem(), j.getItem(), M);
+  if (e.isZero())
+    return recResult;
+  CanonicalForm coeffE;
+  result= recResult;
+  CanonicalForm g;
+  buf.removeLast();
+  for (int i= 1; i < d; i++)
+  {
+    if (degree (e, y) > 0)
+      coeffE= e[i];
+    else
+      coeffE= 0;
+    if (!coeffE.isZero())
+    {
+      CFListIterator k= result;
+      CFListIterator l= p;
+      j= recResult;
+      int ii= 0;
+      CanonicalForm dummy;
+      CFList recDiophantine;
+      CFList buf2, buf3;
+      buf2= factors;
+      buf3= p;
+      for (CFListIterator iii= buf2; iii.hasItem(); iii++)
+        iii.getItem()= mod (iii.getItem(), y);
+      for (CFListIterator iii= buf3; iii.hasItem(); iii++)
+        iii.getItem()= mod (iii.getItem(), y);
+      recDiophantine= diophantine (recResult, buf2, buf3, buf, coeffE);
+      CFListIterator iter= recDiophantine;
+      for (; j.hasItem(); j++, k++, l++, ii++, iter++)
+      {
+        k.getItem() += iter.getItem()*power (y, i);
+        e -= mulMod (iter.getItem()*power (y, i), l.getItem(), M);
+      }
+    }
+    if (e.isZero())
+      break;
+  }
+
+#ifdef DEBUGOUTPUT
+    CanonicalForm test= 0;
+    j= p; 
+    for (CFListIterator i= result; i.hasItem(); i++, j++)
+      test += mod (i.getItem()*j.getItem(), power (y, d));
+    DEBOUTLN (cerr, "test= " << test);
+#endif
+  return result;
+}
+
+// so far only tested for two! factor Hensel lifting
+static inline
+void
+henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
+            const CFList& diophant, CFMatrix& M, CFArray& Pi,
+            const CFList& products, int j, const CFList& MOD)
+{
+  CanonicalForm E;
+  CanonicalForm xToJ= power (F.mvar(), j);
+  Variable x= F.mvar();
+
+  // compute the error
+#ifdef DEBUGOUTPUT
+    CanonicalForm test= 1;
+    for (int i= 0; i < factors.length(); i++)
+    {
+      if (i == 0)
+        test *= mod (bufFactors [i], power (x, j));
+      else
+        test *= bufFactors[i];
+    }
+    test= mod (test, power (x, j));
+    test= mod (test, MOD);
+    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
+    DEBOUTLN (cerr, "test= " << test2);
+#endif
+
+  if (degree (Pi [factors.length() - 2], x) > 0)
+    E= F[j] - Pi [factors.length() - 2] [j];
+  else
+    E= F[j];
+
+  CFArray buf= CFArray (diophant.length());
+
+  // actual lifting
+  CFList diophantine2= diophantine (diophant, factors, products, MOD, E);
+
+  int k= 0;
+  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
+  {
+    buf[k]= i.getItem();
+    bufFactors[k] += xToJ*i.getItem();
+  }
+
+  // update Pi [0]
+  int degBuf0= degree (bufFactors[0], x);
+  int degBuf1= degree (bufFactors[1], x);
+  if (degBuf0 > 0 && degBuf1 > 0)
+  {
+    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
+    if (j + 2 <= M.rows())
+      M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
+  }
+  CanonicalForm uIZeroJ;
+
+    if (degBuf0 > 0 && degBuf1 > 0)
+      uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
+               mulMod (bufFactors[1] [0], buf[0], MOD);
+    else if (degBuf0 > 0)
+      uIZeroJ= mulMod (buf[0], bufFactors[1], MOD);
+    else if (degBuf1 > 0)
+      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
+    else
+      uIZeroJ= 0;
+    Pi [0] += xToJ*uIZeroJ; 
+
+  CFArray tmp= CFArray (factors.length() - 1);
+  for (k= 0; k < factors.length() - 1; k++)
+    tmp[k]= 0;
+  CFIterator one, two;
+  one= bufFactors [0];
+  two= bufFactors [1];
+  if (degBuf0 > 0 && degBuf1 > 0)
+  {
+    while (one.exp() > j) one++;
+    while (two.exp() > j) two++;
+    for (k= 1; k <= (int) ceil (j/2.0); k++)
+    {
+      if (k != j - k + 1)
+      {
+        if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
+        {
+          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
+                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
+                    M (j - k + 2, 1);
+          one++;
+          two++;
+        }
+        else if (one.exp() == j - k + 1)
+        {
+          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
+                    bufFactors[1] [k], MOD) - M (k + 1, 1);
+          one++;
+        }
+        else if (two.exp() == j - k + 1)
+        {
+          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
+                    two.coeff()), MOD) - M (k + 1, 1);
+          two++;
+        }
+      }
+      else
+      {
+        tmp[0] += M (k + 1, 1);
+      }
+    }
+
+    if (j + 2 <= M.rows())
+    {
+      if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
+        tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
+                           (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
+                   - M(1,1) - M (j + 2,1);
+      else if (degBuf0 >= j + 1)
+        tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
+      else if (degBuf1 >= j + 1)
+        tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
+    }
+  }
+  Pi [0] += tmp[0]*xToJ*F.mvar();
+
+  // update Pi [l]
+  int degPi, degBuf;
+  for (int l= 1; l < factors.length() - 1; l++)
+  {
+    degPi= degree (Pi [l - 1], x);
+    degBuf= degree (bufFactors[l + 1], x);
+    if (degPi > 0 && degBuf > 0)
+      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
+    if (j == 1)
+    {
+      if (degPi > 0 && degBuf > 0)
+        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
+                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
+                  M (1, l + 1));
+      else if (degPi > 0)
+        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
+      else if (degBuf > 0)
+        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
+    }
+    else
+    {
+      if (degPi > 0 && degBuf > 0)
+      {
+        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
+        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
+      }
+      else if (degPi > 0)
+        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
+      else if (degBuf > 0)
+      {
+        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
+        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
+      }
+      Pi[l] += xToJ*uIZeroJ;
+    }
+    one= bufFactors [l + 1];
+    two= Pi [l - 1];
+    if (two.exp() == j + 1)
+    {
+      if (degBuf > 0 && degPi > 0)
+      {
+          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
+          two++;
+      }
+      else if (degPi > 0)
+      {
+          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
+          two++;
+      }
+    }
+    if (degBuf > 0 && degPi > 0)
+    {
+      for (k= 1; k <= (int) ceil (j/2.0); k++)
+      {
+        if (k != j - k + 1)
+        {
+          if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
+          {
+            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
+                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
+                      M (j - k + 2, l + 1);
+            one++;
+            two++;
+          }
+          else if (one.exp() == j - k + 1)
+          {
+            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
+                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
+            one++;
+          }
+          else if (two.exp() == j - k + 1)
+          {
+            tmp[l] += mulMod (bufFactors[l + 1] [k],
+                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
+            two++;
+          }
+        }
+        else
+          tmp[l] += M (k + 1, l + 1);
+      }
+    }
+    Pi[l] += tmp[l]*xToJ*F.mvar();
+  }
+
+  return;
+}
+
+// wrt. Variable (1)
+CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
+{
+  if (degree (F, 1) <= 0)
+   return c;
+  else
+  {
+    CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
+    result += (swapvar (c, Variable (F.level() + 1), Variable (1))
+              - LC (result))*power (result.mvar(), degree (result));
+    return swapvar (result, Variable (F.level() + 1), Variable (1));
+  }
+}
+
+// so far only for two factor Hensel lifting
+CFList
+henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
+              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
+              const CFList& LCs2)
+{
+  CFList buf= factors;
+  int k= 0;
+  int liftBoundBivar= l[k];
+  CFList bufbuf= factors;
+  Variable v= Variable (2);
+
+  CFList MOD;
+  MOD.append (power (Variable (2), liftBoundBivar));
+  CFArray bufFactors= CFArray (factors.length());
+  k= 0;
+  CFListIterator j= eval;
+  j++;
+  CFListIterator iter1= LCs1;
+  CFListIterator iter2= LCs2;
+  iter1++;
+  iter2++;
+  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
+  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
+
+  CFListIterator i= buf;
+  i++;
+  Variable y= j.getItem().mvar();
+  if (y.level() != 3)
+    y= Variable (3);
+
+  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
+  M (1, 1)= Pi[0];
+  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
+    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
+                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
+  else if (degree (bufFactors[0], y) > 0)
+    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
+  else if (degree (bufFactors[1], y) > 0)
+    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
+
+  CFList products;
+  for (int i= 0; i < bufFactors.size(); i++)
+  {
+    if (degree (bufFactors[i], y) > 0)
+      products.append (M (1, 1)/bufFactors[i] [0]);
+    else
+      products.append (M (1, 1)/bufFactors[i]);
+  }
+
+  for (int d= 1; d < l[1]; d++)
+    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD);
+  CFList result;
+  for (k= 0; k < factors.length(); k++)
+    result.append (bufFactors[k]);
+  return result;
+}
+
+
+CFList
+henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
+            diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
+            lNew, const CFList& LCs1, const CFList& LCs2)
+{
+  int k= 0;
+  CFArray bufFactors= CFArray (factors.length());
+  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
+  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
+  CFList buf= factors;
+  Variable y= F.getLast().mvar();
+  Variable x= F.getFirst().mvar();
+  CanonicalForm xToLOld= power (x, lOld);
+  Pi [0]= mod (Pi[0], xToLOld);
+  M (1, 1)= Pi [0];
+
+  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
+    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) + 
+                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
+  else if (degree (bufFactors[0], y) > 0)
+    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
+  else if (degree (bufFactors[1], y) > 0)
+    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
+
+  CFList products;
+  for (int i= 0; i < bufFactors.size(); i++)
+  {
+    if (degree (bufFactors[i], y) > 0)
+    {
+      ASSERT (fdivides (bufFactors[i][0], M (1, 1)), "expected exact division");
+      products.append (M (1, 1)/bufFactors[i] [0]);
+    }
+    else
+    {
+      ASSERT (fdivides (bufFactors[i], M (1, 1)), "expected exact division");
+      products.append (M (1, 1)/bufFactors[i]);
+    }
+  }
+
+  for (int d= 1; d < lNew; d++)
+    henselStep2 (F.getLast(),buf,bufFactors, diophant, M, Pi, products, d, MOD);
+
+  CFList result;
+  for (k= 0; k < factors.length(); k++)
+    result.append (bufFactors[k]);
+  return result;
+}
+
+// so far only for two! factor Hensel lifting
+CFList
+henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
+             lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
+             const CFArray& Pi, const CFList& diophant)
+{
+  CFList bufDiophant= diophant;
+  CFList buf= factors;
+  if (sort)
+    sortList (buf, Variable (1));
+  CFArray bufPi= Pi;
+  CFMatrix M= CFMatrix (l[1], factors.length());
+  CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2);
+  if (eval.length() == 2)
+    return result;
+  CFList MOD;
+  for (int i= 0; i < 2; i++)
+    MOD.append (power (Variable (i + 2), l[i]));
+  CFListIterator j= eval;
+  j++;
+  CFList bufEval;
+  bufEval.append (j.getItem());
+  j++;
+  CFListIterator jj= LCs1;
+  CFListIterator jjj= LCs2;
+  CFList bufLCs1, bufLCs2;
+  jj++, jjj++;
+  bufLCs1.append (jj.getItem());
+  bufLCs2.append (jjj.getItem());
+  jj++, jjj++;
+
+  for (int i= 2; i <= lLength && j.hasItem(); i++, j++, jj++, jjj++)
+  {
+    bufEval.append (j.getItem());
+    bufLCs1.append (jj.getItem());
+    bufLCs2.append (jjj.getItem());
+    M= CFMatrix (l[i], factors.length());
+    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
+                         l[i], bufLCs1, bufLCs2);
+    MOD.append (power (Variable (i + 2), l[i]));
+    bufEval.removeFirst();
+    bufLCs1.removeFirst();
+    bufLCs2.removeFirst();
+  }
+  return result;
+}
+
 #endif
 /* HAVE_NTL */

factory/facHensel.h

@@ -364 +364 @@
            );
 
+/// two factor Hensel lifting from univariate to bivariate, factors need not to
+/// be monic
+void
+henselLift122 (const CanonicalForm& F,///< [in] a bivariate poly
+               CFList& factors,       ///< [in, out] a list of univariate polys
+                                      ///< lifted factors
+               int l,                 ///< [in] lift bound
+               CFArray& Pi,           ///< [in, out] stores intermediate results
+               CFList& diophant,      ///< [in, out] result of diophantine
+               CFMatrix& M,           ///< [in, out] stores intermediate results
+               const CFArray& LCs,    ///< [in] leading coefficients
+               bool sort              ///< [in] if true factors are sorted by
+                                      ///< their degree
+              );
+
+/// two factor Hensel lifting from bivariate to multivariate, factors need not
+/// to be monic
+///
+/// @return @a henselLift122 returns a list of lifted factors
+CFList
+henselLift2 (const CFList& eval,   ///< [in] a list of polynomials the last
+                                   ///< element is a compressed multivariate
+                                   ///< poly, last but one element equals the
+                                   ///< last elements modulo its main variable
+                                   ///< and so on. The first element is a
+                                   ///< compressed bivariate poly.
+             const CFList& factors,///< [in] bivariate factors
+             int* l,               ///< [in] lift bounds
+             const int lLength,    ///< [in] length of l
+             bool sort,            ///< [in] if true factors are sorted by
+                                   ///< their degree in Variable(1)
+             const CFList& LCs1,   ///< [in] a list of evaluated LC of first
+                                   ///< factor
+             const CFList& LCs2,   ///< [in] a list of evaluated LC of second
+                                   ///< factor
+             const CFArray& Pi,    ///< [in] intermediate result
+             const CFList& diophant///< [in] result of diophantine
+            );
+
 #endif
 /* FAC_HENSEL_H */

factory/ffreval.cc

@@ -59 +59 @@
 }
 
-/* ------------------------------------------------------------------------ */
-/* forward declarations: fin_ezgcd stuff*/
-static bool findeval_P( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG, int bound, int & counter );
-static void solveF_P ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, int r, CFArray & a );
-static CanonicalForm derivAndEval_P ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a );
-//static CanonicalForm checkCombination_P ( const CanonicalForm & W, const Evaluation & a, const IteratedFor & e, int k );
-static CFArray findCorrCoeffs_P ( const CFArray & P, const CFArray & Q, const CFArray & P0, const CFArray & Q0, const CFArray & S, const CFArray & T, const CanonicalForm & C, const Evaluation & I, int r, int k, int h, int * n );
-static bool liftStep_P ( CFArray & P, int k, int r, int t, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h );
-static bool Hensel_P ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const Variable & x );
-static bool hasFirstAlgVar( const CanonicalForm & f, Variable & a );
-
-CanonicalForm fin_ezgcd( const CanonicalForm & FF, const CanonicalForm & GG )
-{
-  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD;
-  CFArray DD( 1, 2 ), lcDD( 1, 2 );
-  int degF, degG, delta, count, maxeval;
-  maxeval = getCharacteristic(); // bound on the number of eval. to use
-  count = 0; // number of eval. used
-  REvaluation b, bt;
-  bool gcdfound = false;
-  Variable x = Variable(1);
-  f = content( FF, x ); g = content( GG, x ); d = gcd( f, g );
-  F = FF / f; G = GG / g;
-  if( F.isUnivariate() && G.isUnivariate() )
-  {
-    if( F.mvar() == G.mvar() )
-      d *= gcd( F, G );
-    return d;
-  }
-  if( gcd_test_one( F, G, false ) )
-    return d;
-
-  lcF = LC( F, x ); lcG = LC( G, x );
-  lcD = gcd( lcF, lcG );
-  delta = 0;
-  degF = degree( F, x ); degG = degree( G, x );
-  Variable a,bv;
-  if(hasFirstAlgVar(F,a))
-  {
-    if(hasFirstAlgVar(G,bv))
-    {
-      if(bv.level() > a.level())
-        a = bv;
-    }
-    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
-  }
-  else // F not in extension
-  {
-    if(hasFirstAlgVar(G,a))
-      b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
-    else // both not in extension
-      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
-  }
-  while( !gcdfound )
-  {
-    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count ))
-    { // too many eval. used --> try another method
-      Off( SW_USE_EZGCD_P );
-      d *= gcd( F, G );
-      On( SW_USE_EZGCD_P );
-      return d;
-    }
-    delta = degree( Db );
-    if( delta == 0 )
-      return d;
-    while( true )
-    {
-      bt = b;
-      if( !findeval_P( F, G, Fbt, Gbt, Dbt, bt, delta+1, degF, degG, maxeval, count ))
-      { // too many eval. used --> try another method
-        Off( SW_USE_EZGCD_P );
-        d *= gcd( F, G );
-        On( SW_USE_EZGCD_P );
-        return d;
-      }
-      int dd = degree( Dbt );
-      if( dd == 0 )
-        return d;
-      if( dd == delta )
-        break;
-      if( dd < delta )
-      {
-        delta = dd;
-        b = bt;
-        Db = Dbt; Fb = Fbt; Gb = Gbt;
-      }
-    }
-    if( degF <= degG && delta == degF && fdivides( F, G ) )
-      return d*F;
-    if( degG < degF && delta == degG && fdivides( G, F ) )
-      return d*G;
-    if( delta != degF && delta != degG )
-    {
-      bool B_is_F;
-      CanonicalForm xxx;
-      DD[1] = Fb / Db;
-      xxx = gcd( DD[1], Db );
-      if( xxx.inCoeffDomain() )
-      {
-        B = F;
-        DD[2] = Db;
-        lcDD[1] = lcF;
-        lcDD[2] = lcF;
-        B *= lcF;
-        B_is_F = true;
-      }
-      else
-      {
-        DD[1] = Gb / Db;
-        xxx = gcd( DD[1], Db );
-        if( xxx.inCoeffDomain() )
-        {
-          B = G;
-          DD[2] = Db;
-          lcDD[1] = lcG;
-          lcDD[2] = lcG;
-          B *= lcG;
-          B_is_F = false;
-        }
-        else // special case handling
-        {
-          Off( SW_USE_EZGCD_P );
-          d *= gcd( F, G ); // try another method
-          On( SW_USE_EZGCD_P );
-          return d;
-        }
-      }
-      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
-      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
-
-      gcdfound = Hensel_P( B, DD, lcDD, b, x );
-
-      if( gcdfound )
-      {
-        CanonicalForm cand = DD[2] / content( DD[2], Variable(1) );
-        if( B_is_F )
-          gcdfound = fdivides( cand, G );
-        else
-          gcdfound = fdivides( cand, F );
-      }
-    }
-    delta++;
-  }
-  return d * DD[2] / content( DD[2],Variable(1) );
-}
-
-
-static bool hasFirstAlgVar( const CanonicalForm & f, Variable & a )
-{
-  if( f.inBaseDomain() ) // f has NO alg. variable
-    return false;
-  if( f.level()<0 ) // f has only alg. vars, so take the first one
-  {
-    a = f.mvar();
-    return true;
-  }
-  for(CFIterator i=f; i.hasTerms(); i++)
-    if( hasFirstAlgVar( i.coeff(), a ))
-      return true; // 'a' is already set
-  return false;
-}
-
-
-static bool findeval_P( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG, int maxeval, int & count )
-{
-  if( delta != 0 )
-  {
-    b.nextpoint();
-    if( count++ > maxeval )
-      return false; // too many eval. used
-  }
-  while( true )
-  {
-    Fb = b( F );
-    if( degree( Fb ) == degF )
-    {
-      Gb = b( G );
-      if( degree( Gb ) == degG )
-      {
-        Db = gcd( Fb, Gb );
-        if( delta > 0 )
-        {
-          if( degree( Db ) < delta )
-            return true;
-        }
-        else
-          return true;
-      }
-    }
-    b.nextpoint();
-    if( count++ > maxeval ) // too many eval. used
-      return false;
-  }
-}
-
-
-static void solveF_P ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, int r, CFArray & a )
-{
-  CanonicalForm bb, b = C;
-  for ( int j = 1; j < r; j++ )
-  {
-    divrem( S[j] * b, Q[j], a[j], bb );
-    a[j] = T[j] * b + a[j] * P[j];
-    b = bb;
-  }
-  a[r] = b;
-}
-
-
-static CanonicalForm derivAndEval_P ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a )
-{
-  if ( n == 0 )
-    return f( a, x );
-  if ( f.degree( x ) < n )
-    return 0;
-  CFIterator i;
-  CanonicalForm sum = 0, fact;
-  int min, j;
-  Variable v = Variable( f.level() + 1 );
-  for ( i = swapvar( f, x, v); i.hasTerms() && i.exp() >= n; i++ )
-  {
-    fact = 1;
-    min = i.exp() - n;
-    for ( j = i.exp(); j > min; j-- )
-      fact *= j;
-    sum += fact * i.coeff() * power( v, min );
-  }
-  return sum( a, v );
-}
-
-
-static CFArray findCorrCoeffs_P ( const CFArray & P, const CFArray & Q, const CFArray & P0, const CFArray & Q0, const CFArray & S, const CFArray & T, const CanonicalForm & C, const Evaluation & I, int r, int k, int h, int * n )
-{
-  bool what;
-  int i, j, m;
-  CFArray A(1,r), a(1,r);
-  CanonicalForm C0, Dm, g, prodcomb;
-  C0 = I( C, 2, k-1 );
-  solveF_P( P0, Q0, S, T, 1, r, a );
-  for ( i = 1; i <= r; i++ )
-    A[i] = (a[i] * C0) % P0[i];
-  for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ )
-  {
-    Dm = -C;
-    prodcomb = 1;
-    for ( i = 1; i <= r; i++ )
-    {
-      Dm += prodcomb * A[i] * Q[i];
-      prodcomb *= P[i];
-    }
-    if ( Dm != 0 )
-    {
-      if ( k == 2 )
-      {
-        solveF_P( P0, Q0, S, T, Dm, r, a );
-        for ( i = 1; i <= r; i++ )
-          A[i] -= a[i];
-      }
-      else
-      {
-        IteratedFor e( 2, k-1, m+1 );
-        while ( e.iterations_left() )
-        {
-          j = 0;
-          what = true;
-          for ( i = 2; i < k; i++ )
-          {
-            j += e[i];
-            if ( e[i] > n[i] )
-            {
-              what = false;
-              break;
-            }
-          }
-          if ( what && j == m+1 )
-          {
-            g = Dm;
-            for ( i = k-1; i >= 2 && ! g.isZero(); i-- )
-              g = derivAndEval_P( g, e[i], Variable( i ), I[i] );
-            if ( ! g.isZero() )
-            {
-              prodcomb = 1;
-              for ( i = 2; i < k; i++ )
-                for ( j = 2; j <= e[i]; j++ )
-                  prodcomb *= j;
-              g /= prodcomb;
-              if( ! (g.mvar() > Variable(1)) )
-              {
-                prodcomb = 1;
-                for ( i = k-1; i > 1; i-- )
-                  prodcomb *= binomialpower( Variable(i), -I[i], e[i] );
-                solveF_P( P0, Q0, S, T, g, r, a );
-                for ( i = 1; i <= r; i++ )
-                  A[i] -= a[i] * prodcomb;
-              }
-            }
-          }
-          e++;
-        }
-      }
-    }
-  }
-  return A;
-}
-
-
-static bool liftStep_P ( CFArray & P, int k, int r, int t, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h )
-{
-  CFArray K( 1, r ), Q( 1, r ), Q0( 1, r ), P0( 1, r ), S( 1, r ), T( 1, r ), alpha( 1, r );
-  CanonicalForm Rm, C, D, xa = Variable(k) - A[k];
-  CanonicalForm xa1 = xa, xa2 = xa*xa;
-  int i, m;
-  for ( i = 1; i <= r; i++ )
-  {
-    Variable vm = Variable( t + 1 );
-    Variable v1 = Variable(1);
-    K[i] = swapvar( replaceLc( swapvar( P[i], v1, vm ), swapvar( A( lcG[i], k+1, t ), v1, vm ) ), v1, vm );
-    P[i] = A( K[i], k, t );
-  }
-  Q[r] = 1;
-  for ( i = r; i > 1; i-- )
-  {
-    Q[i-1] = Q[i] * P[i];
-    P0[i] = A( P[i], 2, k-1 );
-    Q0[i] = A( Q[i], 2, k-1 );
-    (void) extgcd( P0[i], Q0[i], S[i], T[i] );
-  }
-  P0[1] = A( P[1], 2, k-1 );
-  Q0[1] = A( Q[1], 2, k-1 );
-  (void) extgcd( P0[1], Q0[1], S[1], T[1] );
-  for ( m = 1; m <= n[k]+1; m++ )
-  {
-    Rm = prod( K ) - Uk;
-    for ( i = 2; i < k; i++ )
-      Rm.mod( binomialpower( Variable(i), -A[i], n[i]+1 ) );
-    if ( mod( Rm, xa2 ) != 0 )
-    {
-      C = derivAndEval_P( Rm, m, Variable( k ), A[k] );
-      D = 1;
-      for ( i = 2; i <= m; i++ )
-        D *= i;
-      C /= D;
-      alpha = findCorrCoeffs_P( P, Q, P0, Q0, S, T, C, A, r, k, h, n );
-      for ( i = 1; i <= r; i++ )
-        K[i] -= alpha[i] * xa1;
-    }
-    xa1 = xa2;
-    xa2 *= xa;
-  }
-  for ( i = 1; i <= r; i++ )
-    P[i] = K[i];
-  return prod( K ) - Uk == 0;
-}
-
-
-static bool Hensel_P ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const Variable & x )
-{
-  int k, i, h, t = A.max();
-  bool goodeval = true;
-  CFArray Uk( A.min(), A.max() );
-  int * n = new int[t+1];
-  Uk[t] = U;
-  for ( k = t-1; k > 1; k-- )
-  {
-    Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) );
-    n[k] = degree( Uk[k], Variable( k ) );
-  }
-  for ( k = A.min(); goodeval && (k <= t); k++ )
-  {
-    h = totaldegree( Uk[k], Variable(A.min()), Variable(k-1) );
-    for ( i = A.min(); i <= k; i++ )
-      n[i] = degree( Uk[k], Variable(i) );
-    goodeval = liftStep_P( G, k, G.max(), t, A, lcG, Uk[k], n, h );
-  }
-  delete[] n;
-  return goodeval;
-}

factory/ffreval.h

@@ -5 +5 @@
 
 #include "canonicalform.h"
-#include "cf_eval.h"
-
+#include "cf_reval.h"
 
 class FFREvaluation : public REvaluation
@@ -17 +16 @@
     {
       for( int i=min; i<=max; i++ )
-        values[i] = start[i] = gen->generate();  //generate random point
+        values[i] = start[i] = 0; // start with 0
+
+      nextpoint();
+    }
+    FFREvaluation( int min, int max, const GFRandom & sample ) : REvaluation( min, max, sample ), start( min, max )
+    {
+      for( int i=min; i<=max; i++ )
+        values[i] = start[i] = 0; // start with 0
+
+      nextpoint();
+    }
+    FFREvaluation( int min, int max, const AlgExtRandomF & sample ) : REvaluation( min, max, sample ), start( min, max )
+    {
+      for( int i=min; i<=max; i++ )
+        values[i] = start[i] = 0; // start with 0
 
       nextpoint();
     }
     FFREvaluation& operator= ( const FFREvaluation & e );
+
     void nextpoint();
     bool hasNext();
 };
 
-CanonicalForm fin_ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG );
-
 #endif